// The axon is the venerable Mainen 1996 axon.  Not perfect, yet very stable.  

create a_soma
n_axon_seg = 2
create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]

proc create_axon() {
  n_axon_seg = 2
  create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]

  a_soma {
    equiv_diam = sqrt(area(.5)/(4*PI))
  }

  iseg {               // initial segment between hillock + myelin
     L = 15
     nseg = 5
     diam = equiv_diam/10        // see Sloper and Powell 1982, Fig.71
  }

  hill {                
    L = 20
    nseg = 5
    diam(0:1) = 2*iseg.diam:iseg.diam
  }

  // construct myelinated axon with nodes of ranvier

  for i=0,n_axon_seg-1 {
    myelin[i] {         // myelin element
      nseg = 2
      L = 100
      diam = iseg.diam
    }
    node[i] {           // nodes of Ranvier
      nseg = 1
      L = 1.0
      diam = iseg.diam*.75       // nodes are thinner than axon
    }
  }

  a_soma connect hill(0), 0.5
  hill connect iseg(0), 1
  iseg connect node[0](0), 1
  node[0] connect myelin[0](0), 1

  for i=0,n_axon_seg-2  { 
      myelin[i] connect node[i+1](0), 1
      node[i+1] connect myelin[i+1](0), 1 
  }
}

// --------
// Spines
// --------
proc add_spines() { local a
  forsec "dend" {
    cm*=$1
    g_pas*=$1
  }
  forsec "apic" {
    cm*=$1
    g_pas*=$1
  }
	
}


proc load_3dcell() {

// $s1 filename

  forall delete_section()
  xopen($s1)
  access a_soma
  
  // make sure no compartments exceed 50 uM length
  forall {
    diam_save = diam
    n = int(L/20)    
    if (n==0) n=1
    nseg = n 
    if (n3d() == 0) diam = diam_save
  }    
  create_axon()
  init_cell()
}


proc init_cell() {

  // passive
  forall {
    insert pas
    Ra = ra 
    cm = c_m 
    g_pas = 1/rm
    e_pas = epas_sim
  }
    // exceptions along the axon
     forsec "myelin" cm = cm_myelin
     forsec "node" g_pas = g_pas_node
  forall {
  	insert iA 
	insert kslow 
	insert na
	insert iH
	insert cah
	insert car
	insert cad
	insert bk
	insert sk
  }

  forall if(ismembrane("k_ion")) ek = Ek
  forall if(ismembrane("na_ion")) ena = Ena    
  forall eh=-33
  forall if(ismembrane("ca_ion")) {
 	vshift_ca = 0
  }
}

proc density() {

  forall {
    Ra = ra 
    cm = c_m 
    g_pas = 1/rm
    e_pas = epas_sim
    vshiftm_na= na_shift1
    vshifth_na= na_shift2
    taum_scale_na=na_taum_scale
    tauh_scale_na=na_tauh_scale
    q10_iH=ih_q10
    shift_cah=cah_shift
    shifth_cah=cah_shifth

    shift_car=car_shift
    shifth_car=car_shifth
    qm_car=car_qm
  }


  add_spines(2)
   // exceptions along the axon
   forsec "myelin" cm = cm_myelin
   forsec "node" g_pas = g_pas_node
   forall vshift_na = -5
   forsec "myelin" {
  	gbar_na    = gna_soma
	gbar_kslow =gkslow_beta 
	gbar_iA    = gka_beta
   }

  hill {
  	gbar_na = gna_node
	vshiftm_na=7
	vshifth_na=3
	gbar_iA = gka_node
	gbar_kslow = gkslow_node
  }
	
  iseg {
  	gbar_na = gna_node
	vshiftm_na=7
	vshifth_na=3
	gbar_iA = gka_node
	gbar_kslow = gkslow_node
  }
	
  forsec "node" {
  	gbar_na = gna_node
	vshiftm_na=7
	vshifth_na=3
	gbar_iA = gka_node
	gbar_kslow = gkslow_node
  }

   // we put in all of teh dendrites the same values as in the a_soma and then overwrite the appical values	
  forsec "dend" {
    gbar_na = gna_soma
    gbar_kslow = gkslow_start+gkslow_beta
    gbar_iA = gka_start+gka_beta
    gbar_iH=gih_start 
    pbar_car=pcar_soma
    pbar_cah=pcah_soma
    gbar_sk=gsk_dend
    gbar_bk=gbk_dend
  }
  
  a_soma {
    gbar_na = gna_soma
    gbar_kslow = gkslow_start+gkslow_beta
    gbar_iA = gka_start+gka_beta
    gbar_iH=gih_start
    pbar_car=pcar_soma
    pbar_cah=pcah_soma
    gbar_sk=gsk_soma
    gbar_bk=gbk_soma
  }

  
  a_soma distance(0,0.5)
 
  forsec  ApicalDendSectionName {
	dist1=distance(0)
 	dist2=distance(1)

	
	gbar_na = gna_api
	if (dist1<dist_na){
           gbar_na(0:1) = (gna_soma+dist1*(gna_api-gna_soma)/dist_na):(gna_soma+dist2*(gna_api-gna_soma)/dist_na)
        }


	pbar_cah = pcah_api
        if (dist1<dist_cah){
           pbar_cah(0:1) = (pcah_soma+dist1*(pcah_api-pcah_soma)/dist_cah):(pcah_soma+dist2*(pcah_api-pcah_soma)/dist_cah) 
        }

        pbar_car = pcar_api
        if (dist1<dist_car){
           pbar_car(0:1) = (pcar_soma+dist1*(pcar_api-pcar_soma)/dist_car):(pcar_soma+dist2*(pcar_api-pcar_soma)/dist_car) 
	   //(pcar_soma+dist1*(pcar_api-pcar_soma)/dist_car):(pcar_soma+dist2*(pcar_api-pcar_soma)/dist_car)
        } 


	gbar_iA(0:1) = (gka_start+gka_beta*exp(gka_alpha*dist1)):(gka_start+gka_beta*exp(gka_alpha*dist2))

	gbar_kslow(0:1) = (gkslow_start+gkslow_beta*exp(gkslow_alpha*dist1)):(gkslow_start+gkslow_beta*exp(gkslow_alpha*dist2))
	
	gbar_iH(0:1) = (gih_start+gih_end/(1+exp(gih_alpha*(dist1-gih_x2)))):(gih_start+gih_end/(1+exp(gih_alpha*(dist2-gih_x2))))
	 
	gbar_bk = gbk_dend
        if (dist1<dist_bk){
           gbar_bk(0:1) = (gbk_soma+dist1*(gbk_dend-gbk_soma)/dist_bk):(gbk_soma+dist2*(gbk_dend-gbk_soma)/dist_bk)
        }

	gbar_sk = gsk_dend
	if (dist1<dist_sk){
           gbar_sk(0:1) = (gsk_soma+dist1*(gsk_dend-gsk_soma)/dist_sk):(gsk_soma+dist2*(gsk_dend-gsk_soma)/dist_sk)
        }

	
    }

        forsec  ApicalDendSectionName {

                if (pcah_soma<pcah_api){
                        if (pbar_cah>pcah_api){
                                pbar_cah=pcah_api
                        }
                }
                if (pcah_soma>pcah_api){
                        if (pbar_cah<pcah_api){
                                pbar_cah=pcah_api
                        }
                }

                if (pcar_soma<pcar_api){
                        if (pbar_car>pcar_api){
                                pbar_car=pcar_api
                        }
                }
                if (pcar_soma>pcar_api){
                        if (pbar_car<pcar_api){
                                pbar_car=pcar_api
                        }
                }

                if (gbar_na<gna_api) gbar_na=gna_api
		if (gbar_bk<gbk_dend) gbar_bk=gbk_dend
		if (gbar_sk<gsk_dend) gbar_sk=gsk_dend	
        }


}