load_proc("nrnmainmenu")             // load main NEURON library
//cvode.active(0)


//----------------------	xopen_filehoc()	------------------ 
//Outputs:	Load Files.hoc https://correoweb.ccti.ull.es/horde/imp/mailbox.php?mailbox=INBOX&actionID=149
//Inputs:	$s1 is a string root to "file.hoc" from the actual root
//		$s2 is a string with the file name (not put .hoc)

proc xopen_filehoc() {
	strdef tmp_str
	
	sprint(tmp_str,"xopen(\"%s/%s.hoc\")",$s1,$s2)	
	execute1(tmp_str)
}

//------------------------------------------------------------------------------------------------


xopen_filehoc("../../template","ExperimentControl")

strdef accstr                       // not confuse experimental variable bindings with neurophysiological variable bindings
objref econ                         // Create an experiment object
show_errs=1
debug_lev=1
econ=new ExperimentControl(show_errs,debug_lev)
econ.self_define(econ) 							// points the object at itself
econ.morphology_dir = "../../morphology/nXXX"                             // Setup morphology directory
//econ.morphology_dir = "../../morphology/n128"                             // Setup morphology directory
//econ.morphology_dir = "../../morphology/n130"                             // Setup morphology directory
//econ.morphology_dir = "../../morphology/n404"                             // Setup morphology directory
//econ.morphology_dir = "../../morphology/n129" 
//econ.morphology_dir = "../../morphology/n125" 


econ.data_dir       = "data"                       // Define directory to save produced data 
sprint(econ.syscmd, "mkdir -p %s", econ.data_dir)  // make the data directory
system(econ.syscmd)

//------------------------------------------- Setup cell-----------------------------------------------------------

econ.xopen_geometry_dependent("n123")       		//n123_trunk_axon , load the raw cell morphology
//econ.xopen_geometry_dependent("n128_trunk_axon")		//n128_trunk_axon , load the raw cell morphology
//econ.xopen_geometry_dependent("n130_trunk")
//econ.xopen_geometry_dependent("n404")
//econ.xopen_geometry_dependent("n129")
//econ.xopen_geometry_dependent("n125")


xopen_filehoc("../../morphology/nXXX","apical-trunk-list")
xopen_filehoc("../../morphology/nXXX","apical-non-trunk-list")
xopen_filehoc("../../morphology/nXXX","soma-list")



objref sr,distance_element
objref fileref						//to save in file



//*********************  Calculate the Center of Mass of the soma

xcg=0
ycg=0
zcg=0
ncp=0
Sum_diam=0
strdef temp
//fileref=new File()						//to save in file
//sprint(temp, "%s/trunk_distanceN125p.dat", econ.data_dir)	//to save in file
//fileref.wopen(temp)

forsec "soma" {
		print secname()," Area= ",area(0.5)
		
		n=n3d()
		nseg=n
		ncp+=n		
		for i=0,n-1 {
			print "# of segment",i,"\tDiam",diam3d(i),"\tx=",x3d(i),"\ty=",y3d(i),"\tz=",z3d(i)

//	fileref.printf( "%g  %g  %g\n",x3d(i),y3d(i),z3d(i))
			xcg+=x3d(i)*diam3d(i)
			ycg+=y3d(i)*diam3d(i)
			zcg+=z3d(i)*diam3d(i)
			Sum_diam+=diam3d(i)
		}
}

//fileref.close()	
xcg/=Sum_diam
ycg/=Sum_diam
zcg/=Sum_diam

print "Center of soma: ","Xcg= ", xcg ,"\tYcg= ", ycg,"\tZcg= ", zcg

//******************** End Center of Mass

access apical_dendrite[67] // for the n123
print "apical_dendrite[67]", x3d(0),y3d(0),z3d(0)

//access trunk[10]   //for the n125

//access apical_dendrite[202]  //for the n130

vector_ux=(-xcg+x3d(0))/sqrt((xcg-x3d(0))^2+(ycg-y3d(0))^2+(zcg-z3d(0))^2)
vector_uy=(-ycg+y3d(0))/sqrt((xcg-x3d(0))^2+(ycg-y3d(0))^2+(zcg-z3d(0))^2)
vector_uz=(-zcg+z3d(0))/sqrt((xcg-x3d(0))^2+(ycg-y3d(0))^2+(zcg-z3d(0))^2)

print "vector normal = ",(vector_ux)^2+(vector_uy)^2+(vector_uz)^2



//----------------------------------   oblique_sections()     ----------------------------
//Inputs:	$o1 is the SectionList called apical_non_trunk_list
//		$o2 is the SectionLIst called apical_trunk_list
//		$3 is the number of Tips

proc mydistance(){local loop, leght

	
	i=0
	forsec $o1 {i+=1}						//to save in file
	distance_element=new Vector(i)
	fileref=new File()						//to save in file
	sprint(temp, "%s/soma_distanceN123old.dat", econ.data_dir)	//to save in file
	fileref.wopen(temp)						//to save in file
	fileref.printf("Apical Dendrite Dist.soma   Lenght    Diam.\n")
	i=0
	forsec $o1{	
		sr=new SectionRef()
		lenght=sr.sec.L
		distance_element.x[i]=x3d(0)*vector_ux+y3d(0)*vector_uy+z3d(0)*vector_uz

	
		temp=secname()
		fileref.printf( "%s  %g  %g   %g\n", temp,distance_element.x[i],lenght,diam)
		printf( "%s   %g   %g      %g\n", temp,distance_element.x[i],lenght,diam)	
	i+=1
	}
	
	fileref.close()			//to save in file
}

//------------------------------------------

mydistance(apical_trunk_list)
mydistance(apical_non_trunk_list)
mydistance(soma_list)

//------------------------------------------ Display NEURON--------------
	objref display
	display=new Shape()
	display.show(0)		// mode=0 diameters, =1 centroid, =2	schematic



create center_soma
center_soma{pt3dclear()
	pt3dadd(xcg-1, ycg, zcg, 1)	//-1 and + 1 are added to can see center_soma at display (shape)
	pt3dadd(xcg, ycg, zcg, 1)
	pt3dadd(xcg+1, ycg, zcg, 1)
}
access center_soma

print "Center of soma: ","Xcg= ", x3d(1) ,"\tYcg= ", y3d(1),"\tZcg= ", z3d(1)

create vector_cell
vector_cell{pt3dclear()
	pt3dadd(xcg, ycg, zcg, 1)	//-1 and + 1 are added to can see center_soma at display (shape)
	pt3dadd(vector_ux*500, vector_uy*500, vector_uz*500, 1)
	


}


//center_list is to be able to display the center_soma
objref center_list
center_list=new SectionList()
center_soma center_list.append()


//vector_list is to be able to display the unit vector
objref vector_list
vector_list=new SectionList()
vector_cell vector_list.append()


	

	display.color_list(center_list,1)
        display.color_list(vector_list,6)
	display.color_list(apical_trunk_list,6)
//	display.color_list(axon_sec_list,3)
	display.color_list(soma_list,2)
//	display.color_list(dendrite_sec_list,4)
//	display.color_list(apical_dendrite_sec_list,5)