load_file("kinetics.hoc")

proc choose_stim_site(){local i,avg_diam,totalL,curL,arcl,l1,l2,arcd localobj cell,rvp,sec_list,sn
    // choose the "stimulation site" of the crux of the axon for the cell
    // which is always 1/4 of the electrotonic length away from the crux
    cell = $o1
    kinetics = $2

    rvp = new RangeVarPlot("v")
    cell.axon[cell.junction_site]{ rvp.begin(0) }
    cell.axon[cell.extstim_site]{ rvp.end(0.5) }
    sn = new String()
    sprint(sn.s,".*axon[%d]",cell.extstim_site) // cell.junction_site)
    sec_list = new SectionList()
    rvp.list(sec_list)
    avg_diam = 0
    totalL = 0
    forsec sec_list{
        if(issection(sn.s)){ // in extstim_site, only go as far as stim
            for(i=1;i<n3d() && arc3d(i)<0.5*L;i=i+1){
                arcd = (diam3d(i-1)+diam3d(i))/2
                arcl = arc3d(i)-arc3d(i-1)
                avg_diam = avg_diam + arcd*arcl
            }
            //arcl = arc3d(i)-arc3d(i-1)
            //l1 = (0.5*L-arc3d(i-1))/2
            //l2 = arcl-l1
            //avg_diam = avg_diam + ((diam3d(i-1)*l1+diam3d(i)*l2)/arcl)*(l1*2)
            //totalL = totalL + l1*2
            totalL = totalL + L/2
        }else{ // between extstim and junction_site, use full length
            for(i=1;i<n3d();i=i+1){
                arcd = (diam3d(i-1)+diam3d(i))/2
                arcl = arc3d(i)-arc3d(i-1)
                avg_diam = avg_diam + arcd*arcl
            }
            totalL = totalL + L
        }
    }
    avg_diam = avg_diam/totalL
    print "avg_diam=", avg_diam
    if(kinetics==1){
        cell.stim_length = elength(avg_diam)/8
    }else{
        cell.stim_length = elength(avg_diam)/3
    }
    print "elength=", elength(avg_diam)
    if(cell.stim_length>totalL-cell.stim_length){
        print "totalL = ", totalL, ", stim_length = ", cell.stim_length, ", extstim_length = ", totalL-cell.stim_length
        return
    }
    curL = 0
    forsec sec_list{
        if(cell.stim_loc == 0 && curL+L>cell.stim_length){
            cell.stim_loc = (cell.stim_length-curL)/L
            cell.stim_sr = new SectionRef()
        }
        if(cell.extstimrec_loc == 0 && curL+L>totalL-cell.stim_length){
            cell.extstimrec_loc = (totalL-cell.stim_length-curL)/L
            cell.extstim_sr = new SectionRef()
        }
        curL = curL + L
    }
}

begintemplate SimpleAxon
public axon,extstim_site,stim_sr,stim_loc,junction_site,prop_site,IS,se,nsection,stim_length,extstimrec_loc,extstim_sr
external insert_Traub,insert_Jonas,insert_Yu,insert_conglomerate,insert_modTraub,tstop,choose_stim_site

create axon[31]
objref this,se,stim_sr,extstim_sr

proc init(){
    xopen("C270999B-P2axgeom.hoc")
    nsection = 31
    extstim_site = 8
    junction_site = 5
    prop_site = 10
    IS = 4
    
    forall {
        nseg = int(L/$2) + 1
    }
    
    kinetics = $1
    if(kinetics==0){
        insert_Traub()
    }else if(kinetics==1){
        insert_Jonas()
    }else if(kinetics==2){
        insert_Yu()
    }else if(kinetics==3){
        insert_conglomerate()
    }else{
        insert_modTraub(this,IS)
    }
    forall{
        insert extr
    }
    stim_length = 0
    stim_loc = 0
    extstimrec_loc = 0
  choose_stim_site(this,kinetics)
  print "in SimpleAxon init: axon[junction_site].diam=", axon[junction_site].diam

  axon[IS]{
    se = new SEClamp(0)
    se.dur1 = tstop
    se.amp1 = axon[extstim_site].e_pas
    somaL = 0.01
    somar = 16.136
    rho_soma = Ra*somaL/(PI*somar^2)*10^(-2)
    rho_IS = Ra*L/(PI*(diam/2)^2)*10^(-2)
    se.rs = (rho_soma+rho_IS)/2
  }
}

endtemplate SimpleAxon


begintemplate TypicalAxon
public axon,extstim_site,stim_sr,stim_loc,junction_site,prop_site,IS,se,nsection,stim_length,extstim_sr,extstimrec_loc
external insert_Traub,insert_Jonas,insert_Yu,insert_conglomerate,insert_modTraub,tstop,choose_stim_site

create axon[65]
objref this,se,stim_sr,extstim_sr

proc init(){
    xopen("C040896A-P3axgeom.hoc")
    nsection = 65
    extstim_site = 43
    //stim_site = 42
    junction_site = 40 
    //axon[junction_site]{ distance() }
    //axon[stim_site]{
    //stim_loc = (35-distance(0))/L
    //}
    prop_site = 2
    IS = 0
    
    forall {
        nseg = int(L/$2) + 1
    }
    
    kinetics = $1
    if(kinetics==0){
        insert_Traub()
    }else if(kinetics==1){
        insert_Jonas()
    }else if(kinetics==2){
        insert_Yu()
    }else if(kinetics==3){
        insert_conglomerate()
    }else{
        insert_modTraub(this,IS)
    }
    forall{
        insert extr
    }
    stim_length = 0
    stim_loc = 0
    extstimrec_loc = 0
    choose_stim_site(this,kinetics)
    print "axon[junction_site]=", axon[junction_site].diam
    
    axon[IS]{
        se = new SEClamp(0)
        se.dur1 = tstop
        se.amp1 = axon[extstim_site].e_pas
        somaL = 0.01
        somar = 23.13
        rho_soma = Ra*somaL/(PI*somar^2)*10^(-2)
        rho_IS = Ra*L/(PI*(diam/2)^2)*10^(-2)
        se.rs = (rho_soma+rho_IS)/2
    }
}

endtemplate TypicalAxon


begintemplate ComplexAxon
public axon,extstim_site,stim_loc,prop_site,junction_site,IS,se,nsection,stim_sr,stim_length,extstim_sr,extstimrec_loc
external insert_Traub,insert_Jonas,insert_Yu,insert_conglomerate,insert_modTraub,tstop,choose_stim_site

create axon[75]
objref this,se,stim_sr,extstim_sr

proc init(){local stim_dist,avg_diam,totalL localobj seclist
    xopen("C290500C-P1axgeom.hoc")
    nsection = 75
    extstim_site = 12
    prop_site = 39
    IS = 0
    
    forall {
        nseg = int(L/$2) + 1
    }
    
    kinetics = $1
    if(kinetics==0){
        insert_Traub()
    }else if(kinetics==1){
        insert_Jonas()
    }else if(kinetics==2){
        insert_Yu()
    }else if(kinetics==3){
        insert_conglomerate()
    }else{
        insert_modTraub(this,IS)
    }
    forall{
        insert extr
    }
    junction_site = 10
    stim_length = 0
    stim_loc = 0
    extstimrec_loc = 0
    choose_stim_site(this,kinetics)
  //stim_site = 10
  //stim_loc = elength(axon[stim_site].diam)/(4*axon[stim_site].L)

  axon[IS]{
    se = new SEClamp(0)
    se.dur1 = tstop
    se.amp1 = axon[extstim_site].e_pas
    somaL = 5.3779
    somar = 2.49
    rho_soma = Ra*somaL/(PI*somar^2)*10^(-2)
    rho_IS = Ra*L/(PI*(diam/2)^2)*10^(-2)
    se.rs = (rho_soma+rho_IS)/2
}

}

endtemplate ComplexAxon