// Author: Etay Hay, 2011
//    Models of Neocortical Layer 5b Pyramidal Cells Capturing a Wide Range of
//    Dendritic and Perisomatic Active Properties
//    (Hay et al., PLoS Computational Biology, 2011) 
//
// Model of L5 Pyramidal Cell, constrained both for BAC firing and Current Step Firing


// Copied and adapted some needed functions to work with LFPy from file
// models/L5PCtemplate.hoc
//
// 30.01.2012 ehagen@umb.no

// Fixed problem of looping over the section lists
//
// 08.06.2012 ehagen@umb.no



objref this


proc geom_nseg() {local nSec, L1, L2, D1, D2, nSeg1, nSeg2
  soma area(.5) // make sure diam reflects 3d points
  nSec = 0
  forsec all {
    nseg = 1 + 2*int(L/40)
    nSec = nSec + 1
  }

  nSecAll = nSec
  nSec = 0
  soma { nSec = nSec + 1}
  nSecSoma	= 	nSec
  nSec = 0
  apic { nSec = nSec + 1}
  nSecApical= 	nSec
  nSec = 0
  dend { nSec = nSec + 1}
  nSecBasal	= 	nSec
  nSec = 0
  axon { nSec = nSec + 1}
  nSecAxonalOrig = nSecAxonal	= 	nSec
}


proc distribute_channels()	{local dist,val,base,maxLength
	base = $8
	soma distance()
	maxLength = getLongestBranch($s1)

	forsec $s1		{
		if(0==strcmp($s2,"Ra")){
			Ra = $8
		} else {
			for(x) {
				if ($3==3) {
					dist = distance(x)
				} else {
					dist = distance(x)/maxLength
				}
				val = calculate_distribution($3,dist,$4,$5,$6,$7,$8)
				sprint(tstr,"%s(%-5.10f) = %-5.10f",$s2,x,val)
				execute(tstr)
			}
		}
	}
}


// $s1 section
func getLongestBranch(){local maxL,d localobj distallist,sref
    sprint(tstr,"%s distance()",$s1)
    execute(tstr,this)    
    
  	if(0==strcmp($s1,"axon")){
            sprint(tstr,"%s[0] distance(1)",$s1)
            execute(tstr,this)    
  	}
		maxL = 0
		d = 0
		distallist = new SectionList()
		forsec $s1 {
			sref = new SectionRef()
			if (sref.nchild==0) distallist.append()
		}
		forsec distallist{
			d = distance(1)
			if(maxL<d) maxL = d
		}
		// for the soma case
		if (maxL == 0) {
      $s1 {
        maxL = L
      }
    }
		return maxL
}

// $1 is the distribution type:
//     0 linear, 1 sigmoid, 2 exponential
//     3 step for absolute distance (in microns)
func calculate_distribution()	{local value
	if ($1==0)	{value = $3 + $2*$4}
	if ($1==1) {value = $3 + ($4/(1+exp(($2-$5)/$6)))}
  	if ($1==2) {value = $3 + $6*exp($4*($2-$5))}
	if ($1==3) {
		if (($2 > $5) && ($2 < $6)) {
			value = $3
		} else {
			value = $4
		}
	}
	value = value*$7
	return value
}

// deleting axon, keeping only first 60 micrometers
proc delete_axon(){
    axon {delete_section()}
    create axon[2]
    access axon[0]{
      L= 30
      diam = 1
      nseg = 1+2*int(L/40)
    }
    access axon[1]{
      L= 30
      diam = 1
      nseg = 1+2*int(L/40)
    }

  nSecAxonal = 2
  connect axon(0), soma(0.5)
  connect axon[1](0), axon[0](1) 
  access soma
}