/* Parses a cell specification file implemented as a hoc template.
   It is assumed that dendritic sections are named "dend[n]".  If the global
   flag_spines is set to 1, then any explicit spines (whose names match the
   pattern "spine[n]") are counted on each section, and then used to adjust the
   section's length and diameter according to the normalization procedure
   discussed in Guy Major's PhD thesis.
     Regardless of the setting of flag_spines, all explicit spines are then
   deleted, and tree properties set according to the passive properties
   defined at the top of this file.

   Arguments:
     $1: the index of the cell to read.
 */

/* Passive properties.  These are mapped to NEURON's definitions of the
   same properties below.  */
/*  The values we used in the study */
CM=0.966 // average of Hanbing's WTD1 result.
RM=1/G_PAS
RA=100	// used by default for Wolf, Evans MSN studies

/* Traub values 
 * CM = 0.9
 * RM = 50000
 * RA = 250
 */

G_PAS=1/RM

ApicalHeadDiam = .47
ApicalHeadLen = .71
ApicalNeckDiam = .19
ApicalNeckLen = .44
BasalHeadDiam = .56
BasalHeadLen = .82
BasalNeckDiam = .16
BasalNeckLen = .54

SurfaceAreaOneApicalSpine = (ApicalNeckDiam * PI * ApicalNeckLen + \
                             ApicalHeadDiam * PI * ApicalHeadLen)
SurfaceAreaOneBasalSpine = (BasalNeckDiam * PI * BasalNeckLen + \
                            BasalHeadDiam * PI * BasalHeadLen)
VolumeOneApicalSpine = \
    PI * (ApicalNeckDiam/2.0) * (ApicalNeckDiam/2.0) * ApicalNeckLen + \
    PI * (ApicalHeadDiam/2.0) * (ApicalHeadDiam/2.0) * ApicalHeadLen
VolumeOneBasalSpine = \
    PI * (BasalNeckDiam/2.0) * (BasalNeckDiam/2.0) * BasalNeckLen + \
    PI * (BasalHeadDiam/2.0) * (BasalHeadDiam/2.0) * BasalHeadLen

/*
 * Applies spines to a cell on all dendrites matching the provided forsec
 * pattern.  The
 * global variable flag_spines is ignored, since this method only makes sense
 * to call when spine processing is desired.
 *
 * Arguments:
 * $1: "forsec" pattern describing the dendrites in the tree.
 */
proc applySubtreeSpecificSpines() { local total_surface_area, dend_surface_area, \
    surface_area_one_spine, surface_area_all_spines \
    localobj dendrite_pattern, each_section
  dendrite_pattern = new String($s1)
  surface_area_one_spine = $2
  dendrite_count = 0
  total_surface_area = 0
  forsec dendrite_pattern.s {
    each_section = new SectionRef()

    dendrite_count = dendrite_count + 1
    temp = area(0.5)
    num_spines = 0
    for i = 0, (each_section.nchild - 1) each_section.child[i] {
      if (issection("spine.*")) {
        num_spines = num_spines + 1
      }
    }

    dend_surface_area = 0
    for (x) {
      dend_surface_area = dend_surface_area + area(x)
    }
    total_surface_area = total_surface_area + dend_surface_area

    if (dend_surface_area > 0 && num_spines > 0) {
      surface_area_all_spines = (surface_area_one_spine * num_spines)
      factor = (dend_surface_area + surface_area_all_spines) / dend_surface_area
      L = L * (factor^(2/3))
      diam = diam * (factor^(1/3))
    }
  }
  printf("Dendrite_count: %d\n", dendrite_count)
  printf("total surface area before spine correction: %f\n", total_surface_area)
}


/*
 * Loads a cell, replaces the spines with the Major spine correction factor.
 *
 * Arguments:
 * $1: Path to neuron
 * $2: 1 if apical spine parameters should be assumed, 2 if basal spines,
 *     3 if both should be present, depending on whether the sections are
 *     named dend_apical* or dend_basal*.
 */
proc readcell() { localobj sref
  strdef neuron_name
  neuron_name = $s1
  printf("Loading neuron: %s\n", neuron_name)
  spine_type = $2

  // Before this point, a "create soma" call will have been needed to declare the
  // soma.  However, we don't want that soma to prevent the new soma from loading.
  // So delete it here.  (Yes, it's awkward, but the earlier create seems to be
  // required by NEURON.)
  forall { delete_section() }

  xopen(neuron_name)

  // Ensure that NEURON evaluates the cell in 3D mode when calling diam(), by
  // using a side effect of the area() call.  It doesn't matter which section
  // is used for the call, and the return value of area() can be discarded.
  forall {
    area(0.5)
  }

  // This used to be at the end of the function.  I'm trying to move it to the
  // top, where it makes more sense, since the for(x) construct gets used to
  // do the spine adjustment.
  geom_nseg(100, 0.1)  // nseg according to frequency
  forall {
    nseg *= 9
  }

  surface_area_one_spine = -1
  volume_one_spine = -1
  if (spine_type == 3) {
    applySubtreeSpecificSpines("dend_apical", SurfaceAreaOneApicalSpine)
    applySubtreeSpecificSpines("dend_basal", SurfaceAreaOneBasalSpine)
  } else {
    if (spine_type == 1) {  // Apical
      printf("using apical spines\n")
      surface_area_one_spine = SurfaceAreaOneApicalSpine
      volume_one_spine = VolumeOneApicalSpine
    } else if (spine_type == 2) {
      printf("using basal spines\n")
      surface_area_one_spine = SurfaceAreaOneBasalSpine
      volume_one_spine = VolumeOneBasalSpine
    }

    totalSurfaceArea = 0
    spineSurfaceArea = 0
    spineVolume = 0

    if (flag_spines == 1) {  // Global
      forsec "dend" {
        temp = area(0.5)
        sref = new SectionRef()
        num_spines = 0
        for j = 0, sref.nchild-1 sref.child[j] {
          if (issection("spine.*")) {
            num_spines = num_spines + 1
          }
        }

        SurfaceAreaDend = 0
        volumeDend = 0
        for (x) {
          SurfaceAreaDend = SurfaceAreaDend + area(x)
        }
        totalSurfaceArea = totalSurfaceArea + SurfaceAreaDend

        if (SurfaceAreaDend > 0 && num_spines > 0) {
          SurfaceAreaAllSpines = (surface_area_one_spine * num_spines)
          factor = (SurfaceAreaDend + SurfaceAreaAllSpines) / SurfaceAreaDend
          L = L * (factor^(2/3))
          diam = diam * (factor^(1/3))
        }

        spineSurfaceArea = spineSurfaceArea + SurfaceAreaAllSpines
        spineVolume = spineVolume + (volume_one_spine * num_spines)
      }
    }
    // Print some summary info for tables 1 and 2
    printf("\n")
    // Technically, this is the surface area *before* spine correction.
    printf("surface area: %g\n", totalSurfaceArea)
  }  // spine_type != 3

  forsec "soma" {
    insert hh

    /* Units are S/cm^2, values are taken from Traub 2003 */
    gnabar_hh = 0.1875
    gkbar_hh = 0.125
    /* As far as I can find, Traub 2003 doesn't set a value for the leak
       conductance, and the Traub-approved model on ModelDB doesn't include
       one. */
    gl_hh = 0
    el_hh = -70
  }

  forsec "spine" {  delete_section() }
  forall {
    Ra = RA
    cm = CM
    insert pas  g_pas = G_PAS e_pas=E_PAS
    //insert max
    v = E_PAS
  }

  celsius = 21
}




/*
 * Adds spines to a cell on all dendrites that are part of the specified SectionList.
 * pattern.  The
 * global variable flag_spines is ignored, since this method only makes sense
 * to call when spine processing is desired.
 *
 * Arguments:
 * $o1: SectionList to loop over
 * $2:  Surface area of a single spine
 * $3:  spine density for branches in the SectionList
 *
 *  written by Christina Weaver, Jan 2012
 */
proc applySubtreeConstantSpineDensity() { local total_surface_area, dend_surface_area, \
    surface_area_one_spine, surface_area_all_spines, spine_dens, mean_diam


print "readcell:  top of Subtree proc."

  // Ensure that NEURON evaluates the cell in 3D mode when calling diam(), by
  // using a side effect of the area() call.  It doesn't matter which section
  // is used for the call, and the return value of area() can be discarded.
  forall {
    area(0.5)
  }

  // This used to be at the end of the function.  I'm trying to move it to the
  // top, where it makes more sense, since the for(x) construct gets used to
  // do the spine adjustment.
//  geom_nseg(100, 0.1)  // nseg according to frequency
  geom_nseg(500, 0.1)  // nseg according to frequency
  forall {
    nseg *= 9
  }

  surface_area_one_spine = $2
  spine_dens = $3

printf("Adding spines to %s, density %.3f, spine SA %g\n",$o1,$3,$2)

dendrite_count = 0
  total_surface_area = 0
  forsec $o1 {

printf("\t%s\n",secname() )
    dendrite_count = dendrite_count + 1
    temp = area(0.5)
    num_spines = L * spine_dens

    dend_surface_area = 0
    mean_diam = 0
    for (x) {
      dend_surface_area = dend_surface_area + area(x)
      if( x > 0 && x < 1 )  mean_diam += diam(x)
    }
    mean_diam /= nseg
    total_surface_area = total_surface_area + dend_surface_area

    // adjusted by Christina Weaver, 5 Jan 12.  Still some error, but better than using 
    // Patrick's method which sets the diam throughout the section to whatever it is in the 
    // middle of the section.  See SpineDensFix.hoc for other attempts to fix this.
    //
    if (dend_surface_area > 0 && num_spines > 0) {
      surface_area_all_spines = (surface_area_one_spine * num_spines)
      factor = (dend_surface_area + surface_area_all_spines) / dend_surface_area
      L = L * (factor^(2/3))
      diam = mean_diam * (factor^(1/3))
    }
  }
  printf("Dendrite_count: %d\n", dendrite_count)
  printf("total surface area before spine correction: %f\n", total_surface_area)
}


/************************************************************************************

	find_totLen_byDist()

	DEPRECATED.  Use find_totLen() below instead.

	May 2012, cmw

	Does NOT calculate the total length accurately.  It misses half a section 
	length, equal to L/(2*nseg).  The find_totLen() function below calculates it 
	properly.  This code is kept here as legacy.

************************************************************************************/
func find_totLen_byDist() {  local a_prLen, a_totLen, b_prLen, b_totLen, prxCut, lastX

    prxCut = $1

    printf("\tProximal < %g um\n",prxCut)
    soma { distance() }

    a_prLen = 0
    a_totLen = 0
    b_prLen = 0
    b_totLen = 0

//printf("APICAL\n")
    forsec apical {
        a_totLen += L
	if( distance(1) < prxCut ) {
	    a_prLen += L
	    //printf("%s: **** 0 (%g) - 1 (%g):  prL %g\n",secname(),distance(0),distance(1),a_prLen)
        }
	if( distance(0) < prxCut  && distance(1) > prxCut ) {
            lastX = 0
            for(x) {
	        if( distance(x) < prxCut ) lastX = x
            }
            a_prLen += distance(lastX)-distance(0)
	    //printf("%s: ---- 0 (%g) - %g (%g) - 1 (%g):  prL %g\n",secname(),distance(0),lastX,distance(lastX),distance(1),a_prLen)
        }
        
    }

//printf("BASAL\n")
    forsec basal {
        b_totLen += L
	if( distance(1) < prxCut ) {
	    b_prLen += L
	    //printf("%s: **** 0 (%g) - 1 (%g):  prL %g\n",secname(),distance(0),distance(1),b_prLen)
        }
	if( distance(0) < prxCut  && distance(1) > prxCut ) {
            lastX = 0
            for(x) {
	        if( distance(x) < prxCut ) lastX = x
            }
            b_prLen += distance(lastX)-distance(0)
	    printf("%s: ---- 0 (%g) - lastX = %g (%g) - 1 (%g):  prL %g\n",secname(),distance(0),lastX,distance(lastX),distance(1),b_prLen)
        }
        
    }

    printf("Soma\t%g\n",soma.L)
    printf("Apical\n\tprox\t%g\n\tdist\t%g\n\ttotl\t%g\n",a_prLen,a_totLen-a_prLen,a_totLen)
    printf("Basal\n\tprox\t%g\n\tdist\t%g\n\ttotl\t%g\n",b_prLen,b_totLen-b_prLen,b_totLen)

    return a_prLen + b_prLen + soma.L
}




func find_totLen() {  local prLen, totLen, prxCut, lastX, dstChg, LChg, prDistLen

    prxCut = $1

    //printf("\tProximal < %g um\n",prxCut)
    soma { distance() }

    prLen = 0
    totLen = 0


    forsec "soma" { 
	for(x) {
	    if( x==0 || x == 1 ) continue
	    Lcnt += L / nseg

	}
    }
    prLen = Lcnt
    totLen = Lcnt

    forsec dendritic {
	lastX = 0
	LChg = 0
	for(x) {
	    if( x==0 || x==1 ) continue
	    totLen += L / nseg
	    if( distance(x) > prxCut ) { 
	        continue 
	    }
	    //lastX = x
	    LChg += (L / nseg)
    	}
	prLen += LChg

        /****
	if( distance(0) < prxCut  && distance(1) > prxCut ) {
 	    dstChg = distance(lastX)-distance(0)
	    printf("%s: [L %g / nseg %d = %g; half = %g] ---- 0 (%g) - last X = %g (%g) - 1 (%g):  dstL %g\tLChg %g\tDIFF %g\n",secname(),L,nseg, L/nseg,L/(2*nseg),distance(0),lastX,distance(lastX),distance(1),dstChg,LChg,dstChg-LChg)
	}
	****/
    }


    printf("Total length, including soma (computed with L/nseg):\n\tprox\t%g\ndist\t%g\ntotal\t%g\n",prLen,totLen-prLen,totLen)

    return prLen
}




func find_totApicLen() {  local prLen, totLen, prxCut, lastX, dstChg, LChg, prDistLen

    prxCut = $1

    soma { distance() }

    prLen = 0
    totLen = 0

    forsec "soma" { 
	for(x) {
	    if( x==0 || x == 1 ) continue
	    Lcnt += L / nseg

	}
    }
    prLen = Lcnt
    totLen = Lcnt

    forsec apical {
	lastX = 0
	LChg = 0
	for(x) {
	    if( x==0 || x==1 ) continue
	    totLen += L / nseg
	    if( distance(x) > prxCut ) { 
	        continue 
	    }
	    //lastX = x
	    LChg += (L / nseg)
    	}
	prLen += LChg

    }


    printf("Total apical length, including soma (computed with L/nseg):\n\tprox\t%g\ndist\t%g\ntotal\t%g\n",prLen,totLen-prLen,totLen)

    return prLen
}


func find_totBasLen() {  local prLen, totLen, prxCut, lastX, dstChg, LChg, prDistLen

    prxCut = $1

    soma { distance() }

    prLen = 0
    totLen = 0

    forsec "soma" { 
	for(x) {
	    if( x==0 || x == 1 ) continue
	    Lcnt += L / nseg

	}
    }
    prLen = Lcnt
    totLen = Lcnt

    forsec basal {
	lastX = 0
	LChg = 0
	for(x) {
	    if( x==0 || x==1 ) continue
	    totLen += L / nseg
	    if( distance(x) > prxCut ) { 
	        continue 
	    }
	    //lastX = x
	    LChg += (L / nseg)
    	}
	prLen += LChg

    }


    printf("Total basal length, including soma (computed with L/nseg):\n\tprox\t%g\ndist\t%g\ntotal\t%g\n",prLen,totLen-prLen,totLen)

    return prLen
}