// This function loads all morphologyrelated files and templates that will
// most probably be needed in the experiments to follow.
// based in Terrence Brannon and modified Yiota Poirazi version, July 2001, poirazi@LNC.usc.edu
// Writen by Jose Gomez, October 2008, jfcgomez@ull.es
// I calculated the center of the soma, and unitary vector, NOW it is not
//used the adjustment (=0)
objref vRP, vAPEX
proc cell_analysis() {
xopen_filehoc("../../template","ObliquePath")
xopen_filehoc("../../template","BasalPath")
forall insert d3 // mod file to enable 3D mapping of each point along the cell
$o1.defvar("Distance Calculation", "adjustment", "0", "This adjustment factor is supplied to the vector distance function so that distance calculations are measured at the cell body.")
$o1.xopen_library("Terrence","vectordistance")
$o1.xopen_geometry_dependent("somalist") // It's the same for every pruned cell
$o1.xopen_geometry_dependent("axonseclist") // It's the same for every pruned cell
$o1.xopen_geometry_dependent("basaltreelist") // It's the same name but it's the file written by Jose Gomez
$o1.xopen_geometry_dependent("apicaltrunklist") // It's the same name but it's the file written by Jose Gomez
$o1.xopen_geometry_dependent("apicalnontrunklist") // It's the same name but it's the file written by Jose Gomez
vRP=new Vector(3)
vAPEX=new Vector(3)
CenterOfMass()
myunitvector(apical_trunk_list,apical_non_trunk_list)
}
// Calculate the Center of Mass of the soma
proc CenterOfMass(){
xcg=0
ycg=0
zcg=0
ncp=0
Sum_diam=0
strdef temp
forsec "soma" {
print secname()," Area= ",area(0.5)
n=n3d()
nseg=n
ncp+=n
for i=0,n1 {
print "# of segment",i,"\tDiam",diam3d(i),"\tx=",x3d(i),"\ty=",y3d(i),"\tz=",z3d(i)
xcg+=x3d(i)*diam3d(i)
ycg+=y3d(i)*diam3d(i)
zcg+=z3d(i)*diam3d(i)
Sum_diam+=diam3d(i)
}
}
xcg/=Sum_diam
ycg/=Sum_diam
zcg/=Sum_diam
print "Center of soma: ","Xcg= ", xcg ,"\tYcg= ", ycg,"\tZcg= ", zcg
vRP.x[0]=xcg
vRP.x[1]=ycg
vRP.x[2]=zcg
}
// End Center of Mass
// myunitvector() 
//Inputs: $o1 is the SectionList called apical_trunk_list
// $o2 is the SectionLIst called apical_non_trunk_list
proc myunitvector(){local loop, leght
vector_ux=0
vector_uy=0
vector_uz=0
forsec $o1{
sr=new SectionRef()
lenght=sr.sec.L //I use the lenght as a weight to give more
//importance to longer dendrites
vector_ux+=lenght*(xcg+x3d(0))/sqrt((xcgx3d(0))^2+(ycgy3d(0))^2+(zcgz3d(0))^2)
vector_uy+=lenght*(ycg+y3d(0))/sqrt((xcgx3d(0))^2+(ycgy3d(0))^2+(zcgz3d(0))^2)
vector_uz+=lenght*(zcg+z3d(0))/sqrt((xcgx3d(0))^2+(ycgy3d(0))^2+(zcgz3d(0))^2)
i+=1
}
forsec $o2{
sr=new SectionRef()
lenght=sr.sec.L
vector_ux+=lenght*(xcg+x3d(0))/sqrt((xcgx3d(0))^2+(ycgy3d(0))^2+(zcgz3d(0))^2)
vector_uy+=lenght*(ycg+y3d(0))/sqrt((xcgx3d(0))^2+(ycgy3d(0))^2+(zcgz3d(0))^2)
vector_uz+=lenght*(zcg+z3d(0))/sqrt((xcgx3d(0))^2+(ycgy3d(0))^2+(zcgz3d(0))^2)
i+=1
}
modulevector=sqrt((vector_ux)^2+(vector_uy)^2+(vector_uz)^2)
vector_ux=vector_ux/modulevector
vector_uy=vector_uy/modulevector
vector_uz=vector_uz/modulevector
vAPEX.x[0]=vector_ux
vAPEX.x[1]=vector_uy
vAPEX.x[2]=vector_uz
print "FINAL vector normal = ",(vector_ux)^2+(vector_uy)^2+(vector_uz)^2
print "FINAL vector coord = ",vector_ux,vector_uy,vector_uz
}
//END myunitvector()
