// Channel analysis.hoc executes analysis of a model current generated by channel generator.hoc.
// I-V and normalized g-V plot measured at steady point and maximum value (G/Gmax)(Figure 5B) are formed.
// Risetime to half-maximum activation (Figures 5C,D)and time constant of activation phase 
// are ploted for K+ channel. Representative trances of K+ channels are shown in Figure 5A.
// These analysis help fitting procedure against experimental data.
//
// Takaki Watanabe
// wtakaki@m.u-tokyo.ac.jp

objref vdat4, ikdat, inadat, ik1dat,tdat2,tdat3, iktest, tdat4,kconduc,pkconduc
objref TEVC1, TEVC2, TEVC3, TEVC4, TEVC5, TEVC6, TEVC7, TEVC8, TEVC9, TEVC10, TEVC11, TEVC12, TEVC13, TEVC14, ikdat1

psize=10000

vdat4 = new Vector(psize,0)
ikdat = new Vector(psize,0)
inadat = new Vector(psize,0)
ik1dat = new Vector(psize,0)
tdat2 = new Vector(psize,0)
tdat3 = new Vector(5,0)
iktest = new Vector(psize,0)
tdat4 = new Vector(5,0)
kconduc= new Vector(100,0)
pkconduc= new Vector(100,0)
ikdat1 = new Vector(psize,0)
TEVC1 = new File()
TEVC2 = new File()
TEVC3 = new File()
TEVC4 = new File()
TEVC5 = new File()
TEVC6 = new File()
TEVC7 = new File()
TEVC8 = new File()
TEVC9 = new File()
TEVC10 = new File()
TEVC11 = new File()
TEVC12 = new File()
TEVC13 = new File()
TEVC14 = new File()

vdat4.record(&a.v(0.5))
ikdat.record(&a.ik(0.5))
inadat.record(&a.ina(0.5))
tdat2.record(&t)
ikdat1.record(&a.ik(0.5))

tpoint1=xvdur
tpoint0=(xvpredur+tpoint1)/dt
color6=1

w[11]= new VBox()
w[11].intercept(1)

w[12]= new HBox()
w[12].intercept(1)
yt.x[11]=10
yb.x[11]=-1

ana[12]= new Graph()
ana[12].size(vlow,vhigh+10,yt.x[11],yb.x[11])
ana[12].label(0.05,0.85,"IV@steady")

yt.x[13]=10
yb.x[13]=-1
ana[14]= new Graph()
ana[14].size(vlow,vhigh+10,yt.x[13],yb.x[13])
ana[14].label(0.05,0.85,"IV@Imax")
	
xpanel("")
xbutton("Clear","Clearrecord2()")
xvalue("color","color6",1,"",0,0)
xbutton("New traces","num2=num2+1 num4=num4+1 newtraces2() ")
xpanel()
w[12].intercept(0)
w[12].map()

w[14]= new HBox()
w[14].intercept(1)
yt.x[15]=1
yb.x[15]=0
ana[16]= new Graph()
ana[16].size(vlow,vhigh+10,yt.x[15],yb.x[15])
ana[16].label(0.05,0.85,"Ngk@steady")

yt.x[16]=1
yb.x[16]=0
ana[17]= new Graph()
ana[17].size(vlow,vhigh+10,yt.x[16],yb.x[16])
ana[17].label(0.05,0.85,"Ngk@Gmax(Fig.5B)")
	
w[14].intercept(0)
w[14].map()

w[13]= new HBox()
w[13].intercept(1)
yt.x[12]=10
yb.x[12]=-1
ana[13]= new Graph()
ana[13].size(vlow,vhigh+10,yt.x[12],yb.x[12])
ana[13].label(0.05,0.85,"half rise(Fig.5C,D)")

yt.x[14]=10
yb.x[14]=-1
ana[15]= new Graph()
ana[15].size(vlow,vhigh+10,yt.x[14],yb.x[14])
ana[15].label(0.05,0.85,"tau rise")
	
w[13].intercept(0)
w[13].map()

w[11].intercept(0)
w[11].map("Channel analysis",1000,400,400,400)

num3=0

//IVplot, G/Gmax, activation curve, inactivation curve,
proc run4() {

if(k==-80){
TEVC1.wopen("TEVC1.dat")
ikdat1.printf(TEVC1,"%f\n")
}else if(k==-70){
TEVC2.wopen("TEVC2.dat")
ikdat1.printf(TEVC2,"%f\n")
}else if(k==-60){
TEVC3.wopen("TEVC3.dat")
ikdat1.printf(TEVC3,"%f\n")
}else if(k==-50){
TEVC4.wopen("TEVC4.dat")
ikdat1.printf(TEVC4,"%f\n")
}else if(k==-40){
TEVC5.wopen("TEVC5.dat")
ikdat1.printf(TEVC5,"%f\n")
}else if(k==-30){
TEVC6.wopen("TEVC6.dat")
ikdat1.printf(TEVC6,"%f\n")
}else if(k==-20){
TEVC7.wopen("TEVC7.dat")
ikdat1.printf(TEVC7,"%f\n")
}else if(k==-10){
TEVC8.wopen("TEVC8.dat")
ikdat1.printf(TEVC8,"%f\n")
}else if(k==0){
TEVC9.wopen("TEVC9.dat")
ikdat1.printf(TEVC9,"%f\n")
}else if(k==10){
TEVC10.wopen("TEVC10.dat")
ikdat1.printf(TEVC10,"%f\n")
}else if(k==20){
TEVC11.wopen("TEVC11.dat")
ikdat1.printf(TEVC11,"%f\n")
}else if(k==30){
TEVC12.wopen("TEVC12.dat")
ikdat1.printf(TEVC12,"%f\n")
}else if(k==40){
TEVC13.wopen("TEVC13.dat")
ikdat1.printf(TEVC13,"%f\n")
}else if(k==50){
TEVC14.wopen("TEVC14.dat")
ikdat1.printf(TEVC14,"%f\n")
}

    ana[num2].beginline(color,1)
    ikdat.c.line(ana[num2],tdat2.c,color6,1)
    ana[num2].size(0,tstop,yt.x[11],yb.x[11])
	ana[num2].label(0.9,num3,channelname,2,1,1,1,color6)
    ana[num2].flush
    ana[num2].exec_menu("View = plot")

iktest = new Vector(psize,0)
    tpoint1=xvdur
	tpoint0=((xvpredur+tpoint1)/dt)

	ana[12].mark(k,ikdat.c.x[tpoint0]*area(0.5)*1e-2,"O",3,color6,1)
    ana[12].line(k,ikdat.c.x[tpoint0]*area(0.5)*1e-2)
	ana[12].size(vlow-10,vhigh+10,yt.x[11],yb.x[11])

	kconduc.x[m]= ikdat.c.x[tpoint0]*area(0.5)*1e-2/(k-ek)
	iktest.copy(ikdat,xvpredur/dt,(xvpredur+xvdur)/dt)
	
	pkconduc.x[m]= iktest.c.max()*area(0.5)*1e-2/(k-ek)
	m=m+1
	tdat3 =iktest.c.indvwhere(iktest,">=",iktest.max()*0.5)
	tdat4 =iktest.c.indvwhere(iktest,">=",iktest.max()*0.63)
	ana[14].mark(k,iktest.max()*area(0.5)*1e-2,"O",3,color6,1)
    ana[14].line(k,iktest.max()*area(0.5)*1e-2)
	ana[14].size(vlow,vhigh+10,yt.x[13],yb.x[13])
if(tdat4.x[0]*dt>=xvdur){
	ana[13].mark(k,0,"O",3,color6,1)
    ana[13].line(k,0)
	ana[13].size(vlow,vhigh+10,yt.x[12],yb.x[12])
	ana[15].mark(k,0,"O",3,color6,1)
    ana[15].line(k,0)
	ana[15].size(vlow,vhigh+10,yt.x[14],yb.x[14])
	}else if(tdat4.x[0]*dt<xvdur){
		ana[13].mark(k,(tdat3.x[0]*dt),"O",3,color6,1)
    ana[13].line(k,(tdat3.x[0]*dt))
	ana[13].size(vlow,vhigh+10,yt.x[12],yb.x[12])
	//print k, tdat3.x[0]*dt
	ana[15].mark(k,(tdat4.x[0]*dt),"O",3,color6,1)
    ana[15].line(k,(tdat4.x[0]*dt))
	ana[15].size(vlow,vhigh+10,yt.x[14],yb.x[14])
	}
	if (k == vhigh) { 
	for n=0,(vhigh-vlow)/10{
	ana[16].mark(vlow+10*n,kconduc.c.x[n]/kconduc.c.max(),"O",3,color6,1)
	ana[16].line(vlow+10*n,kconduc.c.x[n]/kconduc.c.max(),"O",3,color6,1)

	ana[17].mark(vlow+10*n,pkconduc.c.x[n]/pkconduc.c.max(),"O",3,color6,1)
	ana[17].line(vlow+10*n,pkconduc.c.x[n]/pkconduc.c.max(),"O",3,color6,1)
	//	print pkconduc.c.x[n]/pkconduc.c.max()
	}
    }
	for p = 12,17 {
    ana[p].exec_menu("View = plot")
	ana[p].exec_menu("10% Zoom out")
	}
}

proc Clearrecord2(){ 
ana[12].erase_all()
ana[12].yaxis(0)
ana[12].size(vlow-10,vhigh+10,yt.x[11],yb.x[11])
ana[13].erase_all()
ana[13].yaxis(0)
ana[13].size(vlow-10,vhigh+10,yt.x[12],yb.x[12])
ana[14].erase_all()
ana[14].yaxis(0)
ana[14].size(vlow-10,vhigh+10,yt.x[13],yb.x[13])
ana[15].erase_all()
ana[15].yaxis(0)
ana[15].size(vlow-10,vhigh+10,yt.x[14],yb.x[14])
ana[16].erase_all()
ana[16].yaxis(0)
ana[16].size(vlow-10,vhigh+10,yt.x[15],yb.x[15])
ana[17].erase_all()
ana[17].yaxis(0)
ana[17].size(vlow-10,vhigh+10,yt.x[16],yb.x[16])
	ana[12].label(0.1,0.9,"IV@steady")
	ana[13].label(0.1,0.9,"half rise ")
	ana[14].label(0.1,0.9,"IV@peak")
	ana[15].label(0.1,0.9,"tau rise")
	ana[16].label(0.1,0.9,"Ngk@steady")
	ana[17].label(0.1,0.9,"Ngk@peak")
}

     newtraces2()