/*---------------------------------------------------------------------

   Written by Alessandro Limongiello
   University "Federico II" of Naples, 2012

-----------------------------------------------------------------------*/

load_file("nrngui.hoc")

/*---------------------------------------------
                   USEFUL DATA
  ---------------------------------------------*/
  
create soma[1] 
create eqSection[1]  
                                                       	
xopen("fixnseg.hoc")            

Vrest = -65
Ra0 = 150                       //Ohm*cm
RaAx = 50
g_pas0 = 1/28000                //Siemens/cm^2   1/30000
cm0 = 1
freq=50
factRm = 1
factCm = 1
factRa = 1
celsius = 34

xopen("useful&InitProc.hoc") 
xopen("stimulation1.hoc") 

//////   Spike counter//////////

objref apcOr

proc insert_APCor() {
    apcOr = new APCount(0.5)
    apcOr.thresh = -10
}


/*---------------------------------------------
                      MAIN
  ---------------------------------------------*/

NSeed = 1
maxStim = 0

objref vbox3

vbox3 = new VBox()    
vbox3.intercept(1)
xpanel("Input Data:")
xvalue("Number of Seed","NSeed",0)
xbutton("Start Simulation","go()")
xpanel()
vbox3.intercept(0)
vbox3.map("Input Data", 0, 120, 160, 60)

synType = 0
objref vecNSpikeOr,vecMuNSpike

vecNSpikeOr = new Vector()
 

/////////////////////////
/* simulation control */
/////////////////////////
strdef morphName
objref f, fLog

proc fileTransfer(){ 

    
    
    strdef strAPPP
    sprint(strAPPP,"mkdir ./MaxStim_Output")

    system(strAPPP)
    
    strdef strFileRisNSpike
    sprint(strFileRisNSpike,"./MaxStim_Output/NSpike_%s.xls",morphName)
    strdef strFileRisNSpikeLog
    sprint(strFileRisNSpikeLog,"./MaxStim_Output/filelog_%s.txt",morphName)   
    
    
    f = new File(strFileRisNSpike) 
    f.wopen(strFileRisNSpike)
    f.printf("synNum\tmuSpikesOr\tsdSpikesOr\n")
    f.close(strFileRisNSpike)
      
    
    fLog = new File(strFileRisNSpikeLog) 
    fLog.wopen(strFileRisNSpikeLog)
    fLog.printf("Log file:\n")
    fLog.close(strFileRisNSpikeLog) 
}



optPrint = 0
tstop = 500 //500

objref synNumVect
synNumVect = new Vector(0,0)

for (j = 2; j <= 25; j = j + 1){
    synNumVect.append(20*j)
}


objref strobj
strobj = new StringFunctions()


proc initialize(){

    t = 0
     
    AleFunctionalSetting()
    cellAnalysis()
    clusteriseAleFunz1()  
    noise0 = 1
    
    insertSyn1(synNum,positionSeed,synType,clusterList)

    canalAllMigliore05Morph()    
    initMigliore05()    
    
}


proc integrateMorph() {

    fadvance()
    objref apcOr
    access soma[0]
    insert_APCor() 
        
    continuerun(tstop)    

}



proc go(){local ii, j, idum, idum1, muOr, sdOr, muEq, sdEq, TmuOr, TsdOr, TmuEq, TsdEq localobj q,q1
    q = new Random()
    q1 = new Random()
    positionSeed = q.discunif(0,9943799)
    poissonSeed = q1.discunif(0,6865119)
    mcell_ran4_init(1)
    q1.MCellRan4(1) 
    q.MCellRan4(1)
    idum = q.seq()
    idum1 = q1.seq()
             
    init()
    readingFile()   
    sprint(morphName,"%s",strFile)
    strobj.right(morphName, 15)    
    fileTransfer()
    vecMuNSpike = new Vector()
    
    for (j = 0; j < synNumVect.size(); j = j + 1){
        positionSeed = q.seq(idum)
        poissonSeed = q1.seq(idum1)    
        {muOr=0 sdOr=0 muEq=0 sdEq=0 TmuOr=0 TsdOr=0 TmuEq=0 TsdEq=0}
        synNum = synNumVect.x(j)
                 
        
        printf ("\n--------------------------------------------------\n")  
        printf("Morphology: %s\n", morphName)
        printf("synNum: %d\n", synNumVect.x(j))
        fLog.aopen(strFileRisNSpikeLog)
        fLog.printf ("\n--------------------------------------------------\n")  
        fLog.printf("Morphology: %s\n", morphName)
        fLog.printf("synNum: %d\n", synNumVect.x(j))                
        fLog.close(strFileRisNSpikeLog)
                                            
        for (ii = 0; ii < NSeed; ii = ii + 1){
            positionSeed = q.repick()
            poissonSeed = q1.repick()
            
            initialize() 
            printf("positionSeed: %d\n", positionSeed)
            printf("poissonSeed: %d\n", poissonSeed)  
            fLog.aopen(strFileRisNSpikeLog)
            fLog.printf("positionSeed: %d\n", positionSeed)
            fLog.printf("poissonSeed: %d\n", poissonSeed)                                     
            fLog.close(strFileRisNSpikeLog)                    
                
            goMorph()                       
            vecNSpikeOr.append(apcOr.n)   
                                                                
            printf("NSpike Morfologico: %d\n", vecNSpikeOr.x[ii])                    
            fLog.aopen(strFileRisNSpikeLog)
            fLog.printf("NSpike Morfologico: %d\n", vecNSpikeOr.x[ii])                                   
            fLog.close(strFileRisNSpikeLog)                     
                                                          
         }  
         
         if (vecNSpikeOr.size()>0){
            muOr = vecNSpikeOr.mean() 
            if (vecNSpikeOr.size()>1) sdOr = vecNSpikeOr.stdev()
            vecNSpikeOr.remove(0, vecNSpikeOr.size()-1)
         }
         vecMuNSpike.append(muOr)

         print "\n\n%%%%%%%%%%%%%%%%%%%%%%%"
         print "synNum\tmuSpikesOr\tsdSpikesOr"
         print synNum,"\t",muOr,"\t",sdOr
         print "%%%%%%%%%%%%%%%%%%%%%%%"
         fLog.aopen(strFileRisNSpikeLog)
         fLog.printf("\n\n%%%%%%%%%%%%%%%%%%%%%%%\n")
         fLog.printf("synNum\tmuSpikesOr\tsdSpikesOr\n")
         fLog.printf("%d\t%2.2f\t%2.2f\n",synNum,muOr,sdOr)
         fLog.printf("%%%%%%%%%%%%%%%%%%%%%%%\n")
         fLog.close(strFileRisNSpikeLog)          
         f.aopen(strFileRisNSpike)
         f.printf("%d\t%2.2f\t%2.2f\t\n",synNum,muOr,sdOr)
         f.close(strFileRisNSpike)        

     }
     
     maxStim = synNumVect.x(vecMuNSpike.max_ind())
     print "MaxStim Value: ",maxStim
        
}

proc goMorph(){ 
    objectvar StimEqList,Exp2SynEqList,NetConEqList
    objref stimEqApp,synEqApp,netConEqApp
    
    integrateMorph()
}



proc canalAllMigliore05Morph(){ local ii,i

    forall if (issection("axon.*")) {   
        insert nax gbar_nax=gna * AXONM
        insert kdr gkdrbar_kdr=gkdr
        insert kap gkabar_kap = KMULTP
    }
    
    forall if (issection("soma.*")) {  
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
        insert na3 gbar_na3=gna
        insert kdr gkdrbar_kdr=gkdr
        insert kap gkabar_kap = KMULTP
    }
    
    access soma
    distance()

    forall if (issection("dendrite.*")) {
		insert hd ghdbar_hd=ghd vhalfl_hd=-73
        insert na3 gbar_na3=gna
        insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/100)
            if (xdist > 100){
				vhalfl_hd=-81
            	gkabar_kad(x) = KMULT*(1+xdist/100)
            }else {
				vhalfl_hd=-73
            	gkabar_kap(x) = KMULTP*(1+xdist/100)
            }
		}
    }
                    
    forall if (issection("apical_dendrite.*")) {
    	insert ds
		insert hd ghdbar_hd=ghd
        insert na3 gbar_na3=gna
        insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/100)
            if (xdist > 100){
			    vhalfl_hd=-81
                gkabar_kad(x) = KMULT*(1+xdist/100)
            }else{
				vhalfl_hd=-73
                gkabar_kap(x) = KMULTP*(1+xdist/100)
            }
		}
    }
    
    
    forall if (issection("user5.*")) {
    	insert ds
		insert hd ghdbar_hd=ghd
        insert na3 gbar_na3=gna
        insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/100)
            if (xdist > 100){
				vhalfl_hd=-81
            	gkabar_kad(x) = KMULT*(1+xdist/100)
            }else{
				vhalfl_hd=-73
                gkabar_kap(x) = KMULTP*(1+xdist/100)
            }
		}
    }

    
}