save_can_from_fig9.hoc

back=0

//////// load needed templates////////////
if(!back){ load_file("nrngui.hoc")               }
if(!back){ load_file("template/ObliquePath.hoc") }
if(!back){ load_file("template/BasalPath.hoc"  ) }

objref econ,f1,f2,f3,ss,cvode
f1= new File()
f2= new File()
f3= new File()
ss=new SaveState()                                                                  
cvode = new CVode(1)
x=cvode.active(1)

restart=0
v_init=-70

strdef morphology_location, morpho_path, ObliqueTrunkSection, BasalTrunkSection
objref vRP, vAPEX

proc xopen_morphology(){
	sprint(morpho_path,"%s/%s",morphology_location,$s1)
	xopen(morpho_path)
}

// Carmen
morphology_location = "pc2b"
ObliqueTrunkSection = "trunk[17]"
BasalTrunkSection   = "trunk[7]"

xopen_morphology("cell.hoc")				//reads morpholofy from the file
xopen_morphology("cell-analysis-simple.hoc") //reads simplified version of CA 


			// --------------- Creating lists-----------------

xopen("lib/TP-lib.hoc")
Tip_sections(apical_non_trunk_list,apical_trunk_list,"Apical")		
objref apical_tip_list
apical_tip_list=TP_list							// Apical Tip list

print "apical_tip_list"
apical_tip_list.printnames()
print  "END apical_tip_list"

objref tmp_pl[num_tips],pl[num_tips],opl[num_tips],degree_apical_tip,peri_trunk_list
objref bl[num_tips],obl[num_tips],degree_basal_tip
xopen("lib/Oblique-lib.hoc")
oblique_sections(apical_tip_list,apical_trunk_list,num_tips)		// apical dendrite path lists and degree of tips

xopen("lib/vector-distance.hoc")

//-----------------------------------------------------------------------------------------------------
printf("Setting up cell\n")                                         // load cell-setup to
xopen("cell_setup_pc2b_CCh.hoc") 
//load_file("atolscalekd.ses")                                            // specify all mechanisms, membrane properties etc


///////////////////Load File//////////////////////

if(!back) { load_file("w_soma_cap.ses") } 


/*cvode_active(1)*/

//////   Spike counter//////////

objref apc, v1

proc insert_APC() {
	apc = new APCount(0.5)
	apc.thresh = $1
	v1 = new Vector()
	apc.record(v1)
}

proc init() {
	if(restart){
		f1.ropen("statekdsoma")
		ss.fread(f1)
		f1.close
		finitialize(v_init)
		ss.restore()
		t=0
		fcurrent()
		cvode.re_init()
	} else {
		finitialize(v_init)
		fcurrent()
	}
}

            
////////////////////Main///////////////////

nsyn=10
objref  s[nsyn], rsyn[nsyn], nc[nsyn]
objref rsynmda[nsyn], ncnmda[nsyn]

objref rect,recv,reci,savv,savi,savt
objref tvec,ampvec,iclamp
objref savtotcurr, savtotna16a, savtottrpm4, totcurr, totna16a, tottrpm4
strdef cmd, filename

i=0
forall{
	for(x){
	
	sprint(cmd,"%s%d","objref totcurrentpart",i)
	execute(cmd)
	
	
	if(ismembrane("ican")){
			sprint(cmd,"%s%d","objref tottrpm4apart",i)
			execute(cmd)

			sprint(cmd,"%s%d","objref savtrpm4apart",i)
			execute(cmd)

			}
	
	if(ismembrane("na16a")){
			sprint(cmd,"%s%d","objref totna16apart",i)
			execute(cmd)
			}		
	
	i=i+1
	}
}
proc main(){

	recv =new Vector()
	rect =new Vector()
	reci =new Vector()
	
	i=0
	forall{
		for(x){
		
		sprint(cmd,"%s%d%s","totcurrentpart",i,"=new Vector()")
		execute(cmd)
		
		
		if(ismembrane("ican")){
				sprint(cmd,"%s%d%s","tottrpm4apart",i,"=new Vector()")
				execute(cmd)
				}
		
		if(ismembrane("na16a")){
				sprint(cmd,"%s%d%s","totna16apart",i,"=new Vector()")
				execute(cmd)
				}		

		
		i=i+1
		}
	}

	//tol=cvode.atolscale(&soma.v(0.5),1e-5)
	period=95.75
	phase0=67.8234
	dur=3
	amp2=0.0
	tstop= 16000

	in=0

	th=-14
	access soma
	insert_APC(th)

	f3.wopen("spike.dat")
	current=0.0 // 0.18
	del=50
	npulse = 10
	pulsdur = 1000
	pulsamp =  0.18//
		starttime = 13000
        curbase = 0.09

	soma {
		iclamp = new IClamp( 0.5 )
		tvec	= new Vector(5 )
		ampvec = new Vector(5 )
	}

	iclamp.del = 0
	iclamp.dur = 1e9
		tvec.x[0 ] = 0
		tvec.x[1 ] = starttime
		tvec.x[2 ] = starttime + pulsdur
		tvec.x[3 ] = starttime + 2*pulsdur
		tvec.x[4 ] = tstop
		ampvec.x[0  ] = curbase
		ampvec.x[1  ] = curbase
		ampvec.x[2  ] = curbase + pulsamp
		ampvec.x[3  ] = curbase
		ampvec.x[4  ] = curbase
	ampvec.play(&iclamp.amp,tvec,1)
	
	access soma[0]
	recv.record(&v(0.5))
	rect.record(&t)
	reci.record(&iclamp.i)
	
	
	i=0
	forall{
		for(x){
		
		sprint(cmd,"%s%d%s%s%s%f%s","totcurrentpart",i,".record(&",secname(),".i_cap(",x,"))")
		execute(cmd)
		
		
		if(ismembrane("ican")){
				sprint(cmd,"%s%d%s%s%s%f%s","tottrpm4apart",i,".record(&",secname(),".can_ican(",x,"))")
				execute(cmd)

				}
		
		if(ismembrane("na16a")){
				sprint(cmd,"%s%d%s%s%s%f%s","totna16apart",i,".record(&",secname(),".ina_na16a(",x,"))")
				execute(cmd)
				}		

		
		i=i+1
		}
	}
	

	//////////////run//////////////////////// 
	run()
	/////////////////////////////////////

	savv= new File()
	savt= new File()
	savi = new File()

	savv.wopen("v.txt")
	savt.wopen("time.txt")
	savi.wopen("i.txt")
	
	recv.printf(savv)
	rect.printf(savt)
	reci.printf(savi)
	
	savv.close
	savt.close
	savi.close
	
	totcurr=new Vector(rect.size())
	totna16a=new Vector(rect.size())
	//tottrpm4=new Vector(rect.size())

	
	i=0
	forall{
		for(x){
		
		sprint(cmd,"%s%d%s%f%s","totcurrentpart",i,".mul(area(",x,")*1e-8)")
		execute(cmd)
		sprint(cmd,"%s%d%s","totcurr.add(totcurrentpart",i,")")
		execute(cmd)
		
		if(ismembrane("ican")){
		
				sprint(cmd,"%s%d%s","savtrpm4apart",i,"= new File()")
				execute(cmd)
				
				sprint(filename,"%s%d%s","nanodomain_data/can",i,".txt")
				sprint(cmd,"%s%d%s","savtrpm4apart",i,".wopen(filename)")
				execute(cmd)
				
				sprint(cmd,"%s%d%s%d%s","tottrpm4apart",i,".printf(savtrpm4apart",i,")")
				print cmd
				execute(cmd)

				sprint(cmd,"%s%d%s","savtrpm4apart",i,".close")
				execute(cmd)
				
				
				}
		
		if(ismembrane("na16a")){
				sprint(cmd,"%s%d%s%f%s","totna16apart",i,".mul(area(",x,")*1e-8)")
				execute(cmd)
				sprint(cmd,"%s%d%s","totna16a.add(totna16apart",i,")")
				execute(cmd)
				}		

		
		i=i+1
		

		}
	}
	
	savtotcurr= new File()
	savtotna16a= new File()
	//savtottrpm4 = new File()

	savtotcurr.wopen("totcurr.txt")
	savtotna16a.wopen("totna16a.txt")
	//savtottrpm4.wopen("tottrpm4.txt")
	
	totcurr.printf(savtotcurr)
	totna16a.printf(savtotna16a)
	//tottrpm4.printf(savtottrpm4)
	
	savtotcurr.close
	savtotna16a.close
	//savtottrpm4.close
	
	
	// print the spike number 
	if(!back) printf("\n Current: %.4f nA | # Spikes: %d\n", current,  apc.n)
	if(!back) printf("___________________\n")
	v1.printf(f3)

	f3.close
	f2.wopen("state.new")
	ss.save
	ss.fwrite(f2)
	f2.close
}


main()
/********    end file    ******/