load_file("nrngui.hoc")
xopen("$(NEURONHOME)/lib/hoc/noload.hoc") // standard run tools

xopen("tabchannels.hoc")
load_file("MC_def.hoc")         

tstop   = 6000
celsius = 35 

dt = 0.01
steps_per_ms = 5  // DT = 0.2 ms
secondorder  = 2  // 2!!!

gnic_tuft = 0.0e-3  // mS/cm2

objref mit
mit = new Mitral(gnic_tuft)

load_file("MC_save.hoc") 
nrncontrolmenu()

// Current injection
T1  = 1500     // 
Dur = 4000
Ic1 = 0.2   //

// Noisy current injection
F0  = 0.0   // 
F1  = 0.0   // 
   

 objref Stim1, Stim2
 mit.soma Stim1 = new IClamp(0.5)
 Stim1.del = T1
 Stim1.dur = Dur		
 Stim1.amp = Ic1  
 
 mit.tuft Stim2 = new OdorInput(0.5)
 Stim2.del  = 1000
 Stim2.dur  = Dur		
 Stim2.torn = 3000
 Stim2.r    = 100
 Stim2.f0   = F0
 Stim2.f1   = F1 
 
 
// Exponential synapse
objref syn1, syn2
mit.tuft syn1 = new ExpSyn(0.5)
syn1.tau =  10
syn1.e   =  0   // -70

// Alpha synapses   
mit.dend syn2 = new AlphaSynapse(0.5)     // Tuft  
//mit.tuft syn2 = new AlphaSynapse(0.0)   // Tuft
syn2.onset = 3362    // 3255 for 0.2 nA
syn2.tau   = 3
syn2.gmax  = 0.0    // Dend: 0.003/0.01; Tuft:0.005/0.02
syn2.e     = -80	 
	 
	 
// Artifical Spiking Input
objref nc
nc = new NetStim(.5)
nc.number   = 100
nc.start    = 1000
nc.interval = 100
nc.noise    = 0

objref ns
ns = new NetCon(nc, syn1) // mit.AMPA
ns.threshold = 0
ns.delay     = 0
ns.weight    = 0e-3       // 30 uS   
 

proc advance() {
fadvance()
 
}

proc run() {
  running_ = 1
  stdinit()
  continuerun(tstop)
}

 
objref mitrate, spk
spk = new Vector()
mitrate = new Vector()

T = (tstop-T1)/1000

proc firing_rate() { local i, n
    n = 0
	if (mit.spiketimes.size() > 0) {
      spk = mit.spiketimes
	  ind = spk.indwhere(">", T1)
	  if (ind != -1) {
	   //print ind
	   n = spk.size()-ind
	 }
    }
	mitrate.append(n/T)
	mitrate.printf()
} 
  
objref g1, g2, g3, g4, g5, g6, g7, g8, g9
proc fig()  {


g3 = new Graph(0)
addplot(g3, 0)
g3.size(0,tstop,-70,60)
g3.view(0,-80,tstop,140, 0,150,400,150)
g3.addvar("Tuft.V", "mit.tuft.v(0.5)", 1, 1, 0.9, 0.9, 2)  

g4 = new Graph(0)
addplot(g4, 0)
g4.size(0,tstop,-70,60)
g4.view(0,-80,tstop,140, 0,420,400,150)
g4.addvar("Prim.V", "mit.prim.v(0.5)", 5, 1, 0.9, 0.9, 2)  

g1 = new Graph(0)
addplot(g1, 0)
g1.size(0,tstop,-70,60)
g1.view(0,-80,tstop,140, 0,680,400,150)
g1.addvar("Soma.V", "mit.soma.v(0.5)", 3, 1, 0.9, 0.9, 2)  

g2 = new Graph(0)
addplot(g2, 0)
g2.size(0,tstop,-70,60)
g2.view(0,-80,tstop,140, 0,920,400,150)
g2.addvar("Dend.V", "mit.dend.v(1.0)", 2, 1, 0.9, 0.9, 2)  //1: black; 2: red; 3: blue

/*
g9 = new Graph(0)
addplot(g9, 0)
g9.size(0,tstop,0,1)
g9.view(0,0,tstop,1, 450,200,400,150)
g9.addvar("Conductance", "syn2.i", 2, 1, 0.9, 0.9, 2)  //1: black; 2: red; 3: blue
*/



g6 = new Graph(0)
addplot(g6, 0)
g6.size(0,tstop,0,1)
g6.view(0,0,tstop,1, 0,150,400,150)
g6.addvar("Ks Activation", "mit.soma.m_IKs", 2, 1, 0.9, 0.9, 2)  
g6.addvar("Ks Inactivation", "mit.soma.h_IKs", 1, 1, 0.9, 0.9, 2)

/*
g5 = new Graph(0)
addplot(g5, 0)
g5.size(0,tstop,-2,2)
g5.view(0,-2,tstop,4, 450,200,400,150)
g5.addvar("Noisy input", "Stim2.i", 2, 1, 0.9, 0.9, 2)  //1: black; 2: red; 3: blue


g7 = new Graph(0)
addplot(g7, 0)
g7.size(0,tstop,0,0.0001)
g7.view(0,0,tstop,0.0001, 0,420,400,150)
g7.addvar("Cai", "mit.soma.cai", 2, 1, 0.9, 0.9, 2)


g8 = new Graph(0)
addplot(g8, 0)
g8.size(0,tstop,-0.0005,0.0005)
g8.view(0,-0.005,tstop,0.01, 450,200,400,150)
g8.addvar("IKCa", "mit.soma.ik_IKCa", 1, 1, 0.9, 0.9, 2)  //1: black; 2: red; 3: blue
*/

}


fig()
run()
//firing_rate()
save_data()