/*---------------------------------------------
                    STIMULATION
  ---------------------------------------------*/


gna =  .025
AXONM = 5
gkdr = 0.01
//celsius = 35.0  
KMULT =  0.03
KMULTP = 0.03
ghd=0.00005

proc canalAllMigliore05(){ 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)
            }
		}
    }
    
    access eqSection[0]
    distance()    

    forsec "eqSection[1]" {   
        insert nax gbar_nax=gna * AXONM*surfFact.x(1)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(1)
        insert kap gkabar_kap = KMULTP*surfFact.x(1)
    }    
    
    
    forsec "eqSection[0]" {   
		insert hd ghdbar_hd=ghd*surfFact.x(0)	vhalfl_hd=-73
        insert na3 gbar_na3=gna*surfFact.x(0)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(0)
        insert kap gkabar_kap = KMULTP*surfFact.x(0)
    } 
    
  
        
    forall if (issection("eqSection[2]")||issection("eqSection[4]")) {
        if (issection("eqSection[2]")) {ii=2}else{ii=4}
    	insert ds
		insert hd ghdbar_hd=ghd*surfFact.x(ii)
        insert na3 gbar_na3=gna*surfFact.x(ii)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(ii)
		insert kap gkabar_kap=0*surfFact.x(ii)
		insert kad gkabar_kad=0*surfFact.x(ii)
        maxDist = distance(1)
//            print "maxDist: ",maxDist
        maxDist = 100
        
		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/maxDist)*surfFact.x(ii)
            {	vhalfl_hd=-73
                gkabar_kap(x) = KMULTP*(1+xdist/maxDist)*surfFact.x(ii)}
		}    
    }
    
    forall if (issection("eqSection[3]")||issection("eqSection[5]")||issection("eqSection[6]")||issection("eqSection[9]")||issection("eqSection[10]")) {
        if (issection("eqSection[3]")) {ii=3}else if (issection("eqSection[5]")){ii=5}else if (issection("eqSection[6]")){ii=6}else if (issection("eqSection[9]")){ii=9}else{ii=10}
    	insert ds
		insert hd ghdbar_hd=ghd*surfFact.x(ii)
        insert na3 gbar_na3=gna*surfFact.x(ii)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(ii)
		insert kap gkabar_kap=0*surfFact.x(ii)
		insert kad gkabar_kad=0*surfFact.x(ii)
        minDist = distance(0)
//            print "minDist: ",minDist
        minDist = 100
        
		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/minDist)*surfFact.x(ii)
            
            {   vhalfl_hd=-81
            	gkabar_kad(x) = KMULT*(1+xdist/minDist)*surfFact.x(ii)}
		}    
    }    


    forsec "eqSection[7]" {
		insert hd ghdbar_hd=ghd*surfFact.x(7) vhalfl_hd=-73
        insert na3 gbar_na3=gna*surfFact.x(7)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(7)
		insert kap gkabar_kap=0*surfFact.x(7)
		insert kad gkabar_kad=0*surfFact.x(7)
        maxDist = distance(1)
//            print "maxDist: ",maxDist
        maxDist = 100
        
		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/maxDist)*surfFact.x(7)
            {
				vhalfl_hd=-73
            	gkabar_kap(x) = KMULTP*(1+xdist/maxDist)*surfFact.x(7)
            }
		}
    }
 
    forsec "eqSection[8]" {
		insert hd ghdbar_hd=ghd*surfFact.x(8) vhalfl_hd=-73
        insert na3 gbar_na3=gna*surfFact.x(8)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(8)
		insert kap gkabar_kap=0*surfFact.x(8)
		insert kad gkabar_kad=0*surfFact.x(8)
        minDist = distance(0)
//            print "minDist: ",minDist
        minDist = 100
        
		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            //print "xdist: ",xdist
            ghdbar_hd(x) = ghd*(1+3*xdist/minDist)*surfFact.x(8)
            {
				vhalfl_hd=-81
            	gkabar_kad(x) = KMULT*(1+xdist/minDist)*surfFact.x(8)
            	//print "gkabar_kad(x): ",gkabar_kad(x)
            }
		}
    } 
    
    forsec "eqSection[11]" {
		insert hd ghdbar_hd=ghd*surfFact.x(11) vhalfl_hd=-73
        insert na3 gbar_na3=gna*surfFact.x(11)
        insert kdr gkdrbar_kdr=gkdr*surfFact.x(11)
		insert kap gkabar_kap=0*surfFact.x(11)
		insert kad gkabar_kad=0*surfFact.x(11)
        minDist = distance(0)
//            print "minDist: ",minDist
        minDist = 100
        
		for (x) if (x>0 && x<1) { 
            xdist = distance(x)
            ghdbar_hd(x) = ghd*(1+3*xdist/minDist)*surfFact.x(11)
            {
				vhalfl_hd=-81
            	gkabar_kad(x) = KMULT*(1+xdist/minDist)*surfFact.x(11)
            }
		}
    }         

}


proc initMigliore05(){
     
	t=0
    forall {
        v=Vrest
        if (ismembrane("nax") || ismembrane("na3")) {ena=55}
        if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {ek=-90}
        if (ismembrane("hd") ) {ehd_hd=-30}
	}
	finitialize(Vrest)
    fcurrent()

    forall {
    	for (x) {
        	if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)}
        	if (ismembrane("hd")) {e_pas(x)=e_pas(x)+i_hd(x)/g_pas(x)}
    		}
	}  
	
}	
	

objref stimApp,synApp,netConApp,stimEqApp,synEqApp,netConEqApp

tau10 = 0.4
tau20 = 0.6
intervalSynSeed = 0
interval_mean = 18.75
interval_sd = 2
syn_asyn = 0
gmaxSynSeed = 0
poissonSeed = 987651119
noise0=0

objref randNum,randInterval,randGmax0,randPoss
objectvar StimList,Exp2SynList,NetConList,StimEqList,Exp2SynEqList,NetConEqList

objref distanceVect
distanceVect = new Vector(1000)

objref forbiddenVect


gmax0 = 2e-3
gmax0_sd = 1e-3

objref rslist,rs
rslist = new List()  //random list


proc insertSyn1(){ local sec, secApp, possApp, i, intervalApp
 
    clusterListApp = $o4
    StimList = new List()
    rslist = new List()  //random list
    Exp2SynList = new List()
    NetConList = new List()
    randNum = new Random($2)
    secApp = randNum.discunif(0,NSEC-1)
    randPoss = new Random($2)
    possApp = randPoss.uniform(0, 1)
    
    randInterval = new Random(intervalSynSeed)
    intervalApp = randInterval.normal(interval_mean,interval_sd^2)
    randGmax0 = new Random(gmaxSynSeed)
    gmax0App = randGmax0.normal(gmax0,(gmax0_sd)^2)  
    
    access soma 
    distance()
    sec = 0
    forbiddenVect = new Vector()
    forall { 
       if (issection("eqSection.*")||issection("soma[0]")||issection("axon[0]")){
           forbiddenVect.append(sec)
       }    
       sec=sec+1
    }
    
    synSecMorphVect = new Vector($1)
    synPosMorphVect = new Vector($1)
    synGmaxMorphVect = new Vector($1)    
        
    for (i = 0; i < $1; i = i + 1){
        secApp = randNum.repick()
        while (forbiddenVect.contains(secApp)){secApp = randNum.repick()}
        if(optPrint) print secApp, i
        sec = 0 
        forall {
          if (sec==secApp) {
            if (syn_asyn==0) intervalApp = randInterval.repick()
            possApp = randPoss.repick()
            gmax0App = randGmax0.repick()
            if (gmax0App<0) gmax0App=0
            synSecMorphVect.x(i) = secApp
            synPosMorphVect.x(i) = possApp 
            synGmaxMorphVect.x(i) = gmax0App   
            
            if(optPrint) print "ALE0: ",gmax0App,secApp                           
            distanceVect.set(i,distance(0.5))
            stimApp = new NetStim(.5)
            stimApp.number=50
            stimApp.interval=intervalApp
            stimApp.start=stimApp.interval/2-stimApp.interval/4
            stimApp.noise=noise0
            
            stimApp.seed(poissonSeed)
            StimList.append(stimApp)       
            
            
            synApp = new Exp2Syn(possApp)
            Exp2SynList.append(synApp)
            synApp.tau1=tau10
            synApp.tau2=tau20
            synApp.e=0
            if (syn_asyn==0) {
               netConApp = new NetCon(StimList.o[i], Exp2SynList.o[i], 0, 0, gmax0App)   
            }else{
               netConApp = new NetCon(StimList.o[0], Exp2SynList.o[i], 0, 0, gmax0App)
            }
            NetConList.append(netConApp)
          }
          sec=sec+1
        }   
    }   
}


possEq0 = 0.999
factTau1Syn = 1
factTau2Syn = 1
factGmaxSyn = 1


factGmaxSyn0 = 1
factGmaxSyn1 = 0


proc createEqSyn(){local j
    StimEqList = new List()
    Exp2SynEqList = new List()
    NetConEqList = new List()  

    for (j = 0; j < $1; j = j + 1){     
        stimEqApp = new NetStim(.5)
        stimEqApp.number = StimList.o[j].number
        stimEqApp.interval=StimList.o[j].interval
        stimEqApp.start=StimList.o[j].start
        stimEqApp.noise=StimList.o[j].noise
        stimEqApp.seed(987651119)                    
        StimEqList.append(stimEqApp) 

        synEqApp = new Exp2Syn(.5)
        synEqApp.e = 0  
        Exp2SynEqList.append(synEqApp)  
        netConEqApp = new NetCon(StimEqList.o(j),Exp2SynEqList.o(j),0,0,gmax0)   
        NetConEqList.append(netConEqApp)      
    }     
  
}


objectvar alreadyVisited

proc changeSurf(){local sec,j,i,index,possApp
    clusterListApp = $o1
    nCluster = clusterListApp.count
    strdef sname
    alreadyVisited = new Vector(0,0)    
    nsyn = $2  
    
    for (j = 0; j < nsyn; j = j + 1){
        possApp=Exp2SynList.o(j).get_loc()

        sname = secname()
    	sec = 0
    	forall if (!issection("eqSection.*")){
            if (issection(sname)) {break}
            sec=sec+1
        }
        
        for (i = 0; i < nCluster; i = i + 1){
            if (clusterListApp.o(i).contains(sec)) {break}
        }        
        
        if (!alreadyVisited.contains(sec)) {
           alreadyVisited.append(sec)
           orSurf2.set(i,orSurf2.x(i)+secSurf.x(sec))   
        }        
        pop_section()
    }
    
    if(optPrint) print "__________surfFact after synapse stimulus__________"
    
    for (i = 0; i < nCluster; i = i + 1){
        if (SURFFACT_PARTIAL_METHOD == 2 ){
           if (orSurf2.x(i)!=0){
              surfFact.set(i,surfFactVect.x(i)*orSurf2.x(i)/eqSurf.x(i))
              surfFactRmCm.set(i,surfFactVect.x(i)*orSurf2.x(i)/eqSurf.x(i))
           }else{
              surfFact.set(i,surfFactVect.x(i)*orSurf1.x(i)/eqSurf.x(i))
              surfFactRmCm.set(i,surfFactVect.x(i)*orSurf1.x(i)/eqSurf.x(i))
           }
        }else if (SURFFACT_PARTIAL_METHOD > 2 ){
           if (orSurf2.x(i)!=0){
              surfFact.set(i,surfFactVect.x(i)*orSurf2.x(i)/eqSurf.x(i))
           }else{
              surfFact.set(i,surfFactVect.x(i)*orSurf1.x(i)/eqSurf.x(i))
           }              
        }
        if(optPrint) print surfFact.x(i)
    }    

}



objref strobj
strobj = new StringFunctions()

maxStimExp1 = 0.75
maxStimExp2 = 0.6


proc findSyn(){ local sec,j,i,index,possApp
    clusterListApp = $o1
    nCluster = clusterListApp.count
    
    strdef sname
    StimEqList = new List()
    Exp2SynEqList = new List()
    NetConEqList = new List() 

    if(optPrint) {print "Synapsi\tSection number\tSection name\tCluster\tpathRi\tpathRi0\tRi(pathRi- pathRi0)\
          \tposs\tpathRiEq\tpathRiEq0\tRiEq (pathRiEq - pathRiEq0)\tpossEq\tfactGmaxSyn"}
    
    nsyn = $3
    if ($4==1){ nsyn= NSEC-2 }
    
    synSecVect = new Vector(nsyn)
    synPosVect = new Vector(nsyn)
    synGmaxVect = new Vector(nsyn)
    
    for (j = 0; j < nsyn; j = j + 1){
             
        possApp=Exp2SynList.o(j).get_loc()
        sname = secname()
       	sec = 0
    	
    	forall if (!issection("eqSection.*")){
            if (issection(sname)) {/*print "section number: ",sec*/ break}
            sec=sec+1
        }

        for (i = 0; i < nCluster; i = i + 1){
            if (clusterListApp.o(i).contains(sec)) {/*print "EqSection number: ",i*/ break}
        }
        synSecVect.x(j)=i
            

        if (withLocation==1){    
           locationSearch2(sec,possApp,i)
           possEq = eqPoss
        }else {possEq = possEq0}
        
        synPosVect.x(j)=possEq
        
        maxmax = maxStim
        if (nsyn<=maxmax){
           factGmaxSyn0 = (nsyn/maxmax)^maxStimExp1
        }else{
           factGmaxSyn0 = 1 + ((nsyn - maxmax)/maxmax)^maxStimExp2
        }

        if (maxStimOnOff==0) factGmaxSyn0 = 1        
        
        if ($2==0){
           if (withLocation==0){
              factGmaxSyn = factGmaxSyn0*eqsecpathri.x(i)/secpathri.x(sec)
           }else{
              factGmaxSyn = factGmaxSyn0*eqLocalpathri/Localpathri
              if(optPrint) print "pappa: ",eqLocalpathri,Localpathri,eqLocalpathri/Localpathri,factGmaxSyn
           }      
           factTauSyn = 1
        } else if ($2==1){
           if (withLocation==0){
              factGmaxSyn = factGmaxSyn0*secpathri.x(sec)^(factGmaxSyn1)
           }else{
              factGmaxSyn = factGmaxSyn0*Localpathri^(factGmaxSyn1)
           }
           factTauSyn = 1 
        } else if ($2==2){
           factGmaxSyn = 1
           factTauSyn = 1 
        }       
        
    
        if(optPrint) {print j+1,sec,secname(),"\t",i,secpathri.x(sec),secpathri0.x(sec),\
              secpathri.x(sec)-secpathri0.x(sec),possApp,eqsecpathri.x(i),eqsecpathri0.x(i),\
              eqsecpathri.x(i)-eqsecpathri0.x(i),possEq,factGmaxSyn}
   

        StimEqList.append(StimList.o(j)) 
         
        
        synEqApp = new Exp2Syn(possEq)
        synEqApp.tau1 = tau10*factTau1Syn
        synEqApp.tau2 = tau20*factTau2Syn
        synEqApp.e = 0  
        Exp2SynEqList.append(synEqApp)  

        ww =   NetConList.o(j).weight   
        
        
        if (syn_asyn==0) {
           netConEqApp = new NetCon(StimEqList.o(j),Exp2SynEqList.o(j),0,0,factGmaxSyn*ww)   
        }else{
           netConEqApp = new NetCon(StimEqList.o(0),Exp2SynEqList.o(j),0,0,factGmaxSyn*ww)
        }
        synGmaxVect.x(j)=factGmaxSyn*ww
        NetConEqList.append(netConEqApp)  
        
        eqSection[i] synEqApp.loc(possEq)


        pop_section()


    }


}


proc locationSearch2(){local nCluster,secNum,orPoss,j
    {secNum=$1 orPoss=$2 j=$3}
    
   
    if(optPrint) print "______locationSearch________"
    
    
    Localpathri=orPoss*secpathri.x(secNum)+(1-orPoss)*secpathri0.x(secNum)
    eqLocalpathri=(Localpathri-orMinpathri.x(j))*(eqsecpathri.x(j)-eqsecpathri0.x(j))\
                   /(orMaxpathri.x(j)-orMinpathri.x(j)) + eqsecpathri0.x(j)
    eqPoss = (eqLocalpathri-eqsecpathri0.x(j))/(eqsecpathri.x(j)-eqsecpathri0.x(j))


    if (eqPoss<=0) eqPoss=0.001
    if (eqPoss>=1) eqPoss=0.999
    if(optPrint) print "EqSection number: ",j,"eqPoss: ",eqPoss

}