load_file("nrngui.hoc")
load_file("vplaym.hoc")
load_file("about.hoc")


//soma
create soma
access soma

soma{
 nseg=1
 diam=10  //um
 L=10/PI  //um
 cm=16  // 1uF/cm2
 //Ra, the axial resistance used to connect sections together [34.5 ohm-cm]
 //100 um2 area (i_vc graphs on same scale as ina)
 
 insert CPR
 insert Kv
 insert h
}



//Place and Hidden the PrintWindowManager
{pwman_place(0,0,0)}

{
xpanel("soma(0 - 1) (Parameters)", 0)
xlabel("soma(0 - 1) (Parameters)")
xlabel("nseg = 1")
soma.L = 3.1831
xvalue("L","soma.L", 1,"", 0, 0 )
soma.diam = 10
xvalue("diam","soma.diam", 1,"", 0, 0 )
soma.cm = 16
xvalue("cm","soma.cm", 1,"", 0, 0 )
soma.gCabar_CPR = 4.92
xvalue("gCabar_CPR","soma.gCabar_CPR", 1,"", 0, 0 )
soma.eCa_CPR = 40
xvalue("eCa_CPR","soma.eCa_CPR", 1,"", 0, 0 )
soma.aoCa_CPR = 0.0031
xvalue("aoCa_CPR","soma.aoCa_CPR", 1,"", 0, 0 )
soma.VhalfCa_CPR = -16.6
xvalue("VhalfCa_CPR","soma.VhalfCa_CPR", 1,"", 0, 0 )
soma.SCa_CPR = 5.7
xvalue("SCa_CPR","soma.SCa_CPR", 1,"", 0, 0 )
//soma.ghbar_CPR = 4.9
//xvalue("ghbar_CPR","soma.ghbar_CPR", 1,"", 0, 0 )
//soma.eh_CPR = -32
//xvalue("eh_CPR","soma.eh_CPR", 1,"", 0, 0 )
//soma.aoh_CPR = 0.00164
//xvalue("aoh_CPR","soma.aoh_CPR", 1,"", 0, 0 )
//soma.Vhalfh_CPR = -56
//xvalue("Vhalfh_CPR","soma.Vhalfh_CPR", 1,"", 0, 0 )
//soma.Sh_CPR = -5.33
//xvalue("Sh_CPR","soma.Sh_CPR", 1,"", 0, 0 )
soma.eCl_CPR = -45
xvalue("eCl_CPR","soma.eCl_CPR", 1,"", 0, 0 )
soma.gClbar_CPR = 6.5
xvalue("gClbar_CPR","soma.gClbar_CPR", 1,"", 0, 0 )
soma.SCl_CPR = 0.09
xvalue("SCl_CPR","soma.SCl_CPR", 1,"", 0, 0 )
soma.FactorCaI_CPR = 0.45
xvalue("FactorCaI_CPR","soma.FactorCaI_CPR", 1,"", 0, 0 )
soma.eKca_CPR = -80
xvalue("eKca_CPR","soma.eKca_CPR", 1,"", 0, 0 )
soma.gKcabar_CPR = 0.5
xvalue("gKcabar_CPR","soma.gKcabar_CPR", 1,"", 0, 0 )
soma.gl_CPR = 0.01
xvalue("gl_CPR","soma.gl_CPR", 1,"", 0, 0 )
soma.el_CPR = 0
xvalue("el_CPR","soma.el_CPR", 1,"", 0, 0 )
soma.gCGMP_CPR = 0
xvalue("gCGMP_CPR","soma.gCGMP_CPR", 1,"", 0, 0 )
soma.eCGMP_CPR = 0.8
xvalue("eCGMP_CPR","soma.eCGMP_CPR", 1,"", 0, 0 )
soma.gKvbar_Kv = 2
xvalue("gKvbar_Kv","soma.gKvbar_Kv", 1,"", 0, 0 )
soma.eKv_Kv = -80
xvalue("eKv_Kv","soma.eKv_Kv", 1,"", 0, 0 )
xpanel(107,58)
}

//Control Panel
{
xpanel("RunControl", 0)
v_init = -39.6 //mV
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run"," myrun()")
t = 0
xvalue("t","t", 2 )
tstop = 1024*5-1
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 10
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 0.1
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
realtime = 0
xvalue("Real Time","realtime", 0,"", 0, 1 )
I=20
xvalue("I (pA)","I",1,"", 0, 1)
//xvalue("scale (pA)","scale",1,"", 0, 1)
xpanel(440,120)
}

//Current Injection
objectvar stim
soma stim = new IClamp(0.5)
stim.del=0
stim.dur=tstop
stim.amp=10  // in nA units


//plot the curves

objectvar save_window_, gObj
//the Display Window
{
 gObj= new VBox()
 gObj.intercept(1)
 save_window_ = new Graph(0)
 save_window_.size(0,tstop,-70,-30)
 save_window_.view(0,-70,tstop,40, 820,0,265, 250)
 graphList[0].append(save_window_)
 save_window_.addexpr("v(.5)", 2, 1)

 
 save_window_ = new Graph(0)
 save_window_.size(0,tstop,-70,-30)
 save_window_.view(0,0,tstop,40, 820,800, 265, 250)
 graphList[1].append(save_window_)
 save_window_.addexpr("soma.nCa_CPR( 0.5 )", 2,1)
 gObj.intercept(0)
 gObj.map("Graphs",820,0,0,0)
}


objectvar save_window_1, gObj1
//the Display Window
{

 gObj1= new VBox()
 
 gObj1.intercept(1)
 save_window_ = new Graph(0)
 save_window_1 = new Graph(0)
 save_window_1.size(0,tstop,-70,-30)
 save_window_1.view(0,0,tstop,40, 820,800, 265, 250)
 graphList[2].append(save_window_1)
 save_window_1.addexpr("soma.iCl( 0.5 )", 2,1)
 save_window_1.addexpr("soma.ih_h( 0.5 )", 1,1)
save_window_1.addexpr("soma.iKv( 0.5 )", 3,1)
save_window_1.addexpr("soma.iCa( 0.5 )", 4,1)
save_window_1.addexpr("soma.iKca( 0.5 )", 5,1)
save_window_1.addexpr("soma.il_CPR( 0.5 )", 6,1)
save_window_1.addexpr("soma.iCGMP_CPR( 0.5 )", 7,1)

 gObj1.intercept(0)
 
 gObj1.map("Graphs",720,0,0,0)
}

objectvar save_window_2, gObj2
{
 gObj2= new VBox()
 
 gObj2.intercept(1)
 save_window_ = new Graph(0)
 save_window_2 = new Graph(0)
 save_window_2.size(0,tstop,0.68,1.66)
 save_window_2.view(0, 0.68, tstop, 0.98, 303, 308, 299.2, 199.2)
 graphList[2].append(save_window_2)
 save_window_2.addexpr("Cas_CPR", 1,1,0.8, 0.9, 2)
 
 gObj2.intercept(0)
 
 gObj2.map("Graphs",720,0,0,0)

}






totalsteps=int(tstop/dt)+1

//Begin VectorPlay[1]
objectvar ocbox_
ocbox_=new VectorPlay(1)
ocbox_.vy = new Vector(totalsteps)
ocbox_.vx = new Vector(totalsteps)


//the photocurrent for the light response of cones

func pcurrent_cone() { local photocurrent
              tc=$1
scale=-I
a1	=0.0100
tau1	=60.0000
a2	=0.0600
s2	=120.0000
t2	=555.0000
a3	=0.7000
s3	=200.0000
t3	=1500.0000
s4	=80.0000
t4	=1200.0000
if (tc<=2000) photocurrent=I
if (tc>2000) {
	tc=tc-2000
	iexp1=(1+a1)*(1-exp(-tc/tau1))^3
	isigma2=-(a1+a2)/(1+exp(-(tc-t2)/s2))
	isigma3=-a3/((1+exp(-(tc-t3)/s3)))^1.3
	isigma4=-(1-a2-a3)/((1+exp(-(tc-t4)/s4)))^2
	inormal=iexp1+isigma2+isigma3+isigma4
	photocurrent=I+scale*inormal
	}
  return photocurrent*0.001  //the result is in nA
}  
  











//Apply Current
proc PlayInto() { local totalt,stept
totalt=$1
stept=$2

 
stim.dur=totalt
// resize the vx and vy vector according to the total length of the waveform
ocbox_.vy.resize(totalsteps)
ocbox_.vx.resize(totalsteps)

for (i=0;i<totalsteps;i=i+1) { ocbox_.vx.x[i]=i*stept  
ocbox_.vy.x[i]=pcurrent_cone(i*stept) } // time in ms
ocbox_.g.erase()
ocbox_.g.unmap()
//input into
ocbox_.vy.plot(ocbox_.g, ocbox_.vx)
ocbox_.sname = "IClamp[0].amp"
ocbox_.con1(1)
ocbox_.g.size(0,tstop,0,40)
}
//End VectorPlay[1]

object_push(ocbox_)
have_name = 1   //this is generated by save session, but how to realize it by
//handmade codes is still unsolved
//which should be solved, push/pop is to access into the private variables of
//the class
object_pop()


objref Vm,gVm, tempV
tempV= new Vector()
Vm = new Vector() // vector to record the voltage response
gVm= new Graph() // the graph to plot Vm
gVm.unmap()


 print tstop
// run the first time to construct the window for the figures
PlayInto(tstop,dt)
ocbox_.b.map("VectorPlay[1]", 424, 519, 306, 283)
run()

    
proc myrun() { local ii 
 totalsteps=int(tstop/dt)+1  //0:dt:tstop
 Vm.record(&v(0.5))
  
 PlayInto(tstop,dt)
 run()

 Vm.resize(totalsteps-2)

 gVm.erase()
 gVm.unmap()
 Vm.plot(gVm)
 gVm.size(0,Vm.size(),-39.9,-35.2)
 gVm.view(0, -39.9, Vm.size(), 4.7, 35, 138, 264.75, 250)

 doNotify()
}