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

	Christina Weaver
	August 2011

    auxiliary procedures for loading data, and model fitting.

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




print "top of aux_procs.hoc"




/*
 * 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 "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(100, 0.1)  // nseg according to frequency
  forall {
    nseg *= 5
  }

  surface_area_one_spine = $2
  spine_dens = $3

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

  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.  
    //
    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)
}