// setup_wrote_vecs_standard.hoc
// basic version of declarations and functions to record v and cai from
// the model and write these to folders passed to the setup_vecs and
// write_vecs procs.
/*

Contains the create and write vector functions and associated
declarations modified from param3.hoc (which calculated simulations
for changes in the 3 parameters mshift (in ca current), ECl, gbar_ca)

(Modified from param3 documentation:)

For the following send traces

I) control Ca traces (protocol bAP)
II) Ca with inhib (bAPdendinhib)
III) calculate dCa for control and inhibited and their ratio for each
compartment:
dCa(inhibited)/dCa(control) = ([Ca2+]peak,inh - B)/([Ca2+]peak,ctrl - B)
calculated from I), II):
dCa_ratio.txt
IV) record also the dCa(inhibited), dCa(control) values separately
dCa_inhib.txt
dCa_ctrl.txt

The calcium and voltage trace files are stored in "proto_dir/file" 
where protocols directory proto_dir can be upto

inhibspine, 1xinhibdend, 10xinhibdend, bAP, bAPinhibspine, bAP10xinhibdend,
bAP1xinhibdend

and files can be as many as spine0_ca, spine1_ca, spine2_ca, dend_ca, spine0_v,
spine1_v, spine2, dend_v

however here unlike in FiveConditions.hoc, the protocols and files are
limited to X/bAP X/bAPspineinhib and spine1_ca, spine0_ca where spine1
is the middle spine head, and spine0 is the proximal spine head
(receives a smaller bAP typically because of the slight increase in
bAP amplitude as a function of distance along the dendrite.
*/
// "denominator" divides Na K max conductances so is bAP height controller 
denominator = 1 // will be denominators of 1/1, 1/2, 1/4, 1/8, 1/16
strdef proto_dir, file_suffix
file_suffix = ""
strdef filename
objref f // temporary file object
f = new File()

objref dend_v_vec, spine0_v_vec, spine1_v_vec, spine2_v_vec
objref dend_cai_vec, spine0_cai_vec, spine1_cai_vec, spine2_cai_vec
// tstop= 1 // test condition
// below are graphs (the 5 conditions turned into 6 with the addition
// of protocol bAP + 1x inhibition in dendrite
num_of_conditions = 10 // ctrl and inhib for 1,1/2,1/4,1/8,1/16 // 18
objref  graphs_v[num_of_conditions], graphs_ca[num_of_conditions], graphs_gates[num_of_conditions]
objref graph_bAP_envelopes
// add the activation and inactivation curves for the ca chan in spines and dend
objref m_ca_dend_vec, h_ca_dend_vec, m_ca_spine1_vec, h_ca_spine1_vec
objref  m_ca_spine0_vec, h_ca_spine0_vec
objref bAP_peak_v_vec, bAP_peak_v_x_vec, bAP_profile_vec

for i=0, num_of_conditions-1 {
  graphs_v[i] = new Graph(0)
  graphs_ca[i] = new Graph(0)
  graphs_gates[i] = new Graph(0)
}
graph_bAP_envelopes = new Graph()

proto_num=0  // index into above protocol graphs, iterates over num_conditions

proc setup_vectors_to_record() {
// record voltages in the dendrite and the Spine[X].head X=0,1,2 (corresponds to paper spines 1 2 3)
  dend_v_vec = new Vector()
  dend_cai_vec= new Vector()

  spine0_v_vec= new Vector()
  spine1_v_vec= new Vector()
  spine2_v_vec= new Vector()
  spine0_cai_vec = new Vector()
  spine1_cai_vec = new Vector()
  spine2_cai_vec = new Vector()

  m_ca_dend_vec = new Vector()
  h_ca_dend_vec = new Vector()
  m_ca_spine1_vec = new Vector()
  h_ca_spine1_vec = new Vector()
  m_ca_spine0_vec = new Vector()
  h_ca_spine0_vec = new Vector()
  bAP_peak_v_vec = new Vector()
  bAP_peak_v_x_vec = new Vector()
  t_vec = new Vector()
  t_vec.indgen(0, tstop, dt)

  spinelocation=adjacent_shaft_x_loc // record the dendrite where Spine[1] is attached
  sprint(tmpstr, "dendrite.v(%f)", adjacent_shaft_x_loc) // used to be 180/dendrite.L
  dend_v_vec.label(tmpstr)
  spine0_v_vec.label("Spine[0].head.v(0.5)")
  spine1_v_vec.label("Spine[1].head.v(0.5)")
  spine2_v_vec.label("Spine[2].head.v(0.5)")

  sprint(tmpstr,"dendrite.cai(%g)", spinelocation) // dend. comp. adj. to spine
  dend_cai_vec.label(tmpstr)
  spine0_cai_vec.label("Spine[0].head.cai(0.5)")
  spine1_cai_vec.label("Spine[1].head.cai(0.5)")
  spine2_cai_vec.label("Spine[2].head.cai(0.5)")

  sprint(tmpstr,"dendrite.m_ca(%g)", spinelocation) // dend. comp. adj. to spine
  m_ca_dend_vec.label(tmpstr)
  sprint(tmpstr,"dendrite.h_ca(%g)", spinelocation) // dend. comp. adj. to spine
  h_ca_dend_vec.label(tmpstr)
  m_ca_spine1_vec.label("Spine[1].head.m_ca(0.5)")
  h_ca_spine1_vec.label("Spine[1].head.h_ca(0.5)")
  m_ca_spine0_vec.label("Spine[0].head.m_ca(0.5)")
  h_ca_spine0_vec.label("Spine[0].head.h_ca(0.5)")

  dend_v_vec.record(&dendrite.v(spinelocation), t_vec)
  dend_cai_vec.record(&dendrite.cai(spinelocation), t_vec)

  spine0_v_vec.record(&Spine[0].head.v(0.5), t_vec)
  spine1_v_vec.record(&Spine[1].head.v(0.5), t_vec)
  spine2_v_vec.record(&Spine[2].head.v(0.5), t_vec)

  spine0_cai_vec.record(&Spine[0].head.cai(0.5), t_vec)
  spine1_cai_vec.record(&Spine[1].head.cai(0.5), t_vec)
  spine2_cai_vec.record(&Spine[2].head.cai(0.5), t_vec)

  m_ca_dend_vec.record(&dendrite.m_ca(spinelocation), t_vec)
  h_ca_dend_vec.record(&dendrite.h_ca(spinelocation), t_vec)
  m_ca_spine1_vec.record(&Spine[1].head.m_ca(0.5), t_vec)
  h_ca_spine1_vec.record(&Spine[1].head.h_ca(0.5), t_vec)
  m_ca_spine0_vec.record(&Spine[0].head.m_ca(0.5), t_vec)
  h_ca_spine0_vec.record(&Spine[0].head.h_ca(0.5), t_vec)
}

proc  write_vectors() {
//previous versions looked like:
//  sprint(tmpstr,"NaKgs/1/%d/%s/spine0_cai.dat", denominator, proto_dir)
// however this version merely relies on proto_dir, file_suffix
  sprint(tmpstr,"%s/spine0_cai%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec, spine0_cai_vec)
  sprint(tmpstr,"%s/spine1_cai%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec, spine1_cai_vec)
  sprint(tmpstr,"%s/spine2_cai%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec, spine2_cai_vec)
  sprint(tmpstr,"%s/spine0_v%s.dat", proto_dir, file_suffix)

  write_vec(tmpstr,t_vec, spine0_v_vec)
  sprint(tmpstr,"%s/spine1_v%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec, spine1_v_vec)
  sprint(tmpstr,"%s/spine2_v%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec, spine2_v_vec)
  
  sprint(tmpstr,"%s/dend_v%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec,dend_v_vec)
  sprint(tmpstr,"%s/dend_cai%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr,t_vec, dend_cai_vec)
/* leave out ca activation inactivation for now
  Note will need to fix below by adding denominator folder if decide to include
  sprint(tmpstr,"m_ca_dend_vec.dat" , proto_dir)
  write_vec(tmpstr,t_vec, m_ca_dend_vec)
  
  sprint(tmpstr,"h_ca_dend_vec.dat" , proto_dir)
  write_vec(tmpstr,t_vec, h_ca_dend_vec)

  sprint(tmpstr,"m_ca_spine1_vec.dat" , proto_dir)
  write_vec(tmpstr,t_vec, m_ca_spine1_vec)

  sprint(tmpstr,"h_ca_spine1_vec.dat" , proto_dir)
  write_vec(tmpstr,t_vec, h_ca_spine1_vec)

  sprint(tmpstr,"m_ca_spine0_vec.dat" , proto_dir)
  write_vec(tmpstr,t_vec, m_ca_spine0_vec)

  sprint(tmpstr,"h_ca_spine0_vec.dat" , proto_dir)
  write_vec(tmpstr,t_vec, h_ca_spine0_vec)
*/
  spine0_cai_vec.line(graphs_ca[proto_num],t_vec,3,1)
  spine1_cai_vec.line(graphs_ca[proto_num],t_vec,1,1)
  spine2_cai_vec.line(graphs_ca[proto_num],t_vec,2,1)
  dend_cai_vec.line(graphs_ca[proto_num],t_vec,5,1)
  // graphs
  graphs_ca[proto_num].exec_menu("View = plot")
  sprint(tmpstr,"cai under protocol %s %s" , proto_dir, file_suffix)
  graphs_ca[proto_num].label(.5,.95,tmpstr)
  graphs_ca[proto_num].view(0, 0, 500 , 0.0011, 180, 280, 300, 200)
  
  spine0_v_vec.line(graphs_v[proto_num],t_vec,3,1)
  spine1_v_vec.line(graphs_v[proto_num],t_vec,1,1)
  spine2_v_vec.line(graphs_v[proto_num],t_vec,2,1)
  dend_v_vec.line(graphs_v[proto_num],t_vec,5,1)
  graphs_v[proto_num].exec_menu("View = plot")

  sprint(tmpstr,"v under protocol %s %s" , proto_dir, file_suffix)
  graphs_v[proto_num].label(.5,.95,tmpstr)
  graphs_v[proto_num].view(114, -70, 6 , 78, 200, 300, 300, 200)

  m_ca_spine0_vec.line(graphs_gates[proto_num], t_vec,1, 1)
  h_ca_spine0_vec.line(graphs_gates[proto_num], t_vec,2, 1)
  m_ca_spine1_vec.line(graphs_gates[proto_num], t_vec,3, 1)
  h_ca_spine1_vec.line(graphs_gates[proto_num], t_vec,4, 1)
  m_ca_dend_vec.line(graphs_gates[proto_num], t_vec,5, 1)
  h_ca_dend_vec.line(graphs_gates[proto_num], t_vec,6, 1)
  graphs_gates[proto_num].exec_menu("View = plot")
  sprint(tmpstr,"m_ca, h_ca under protocol %s" ,  proto_dir)
  graphs_gates[proto_num].label(.5,.95,tmpstr)
  // save the bAP peak v envelope (column 2) as a function of distance (column 1)
  bAP_peak_v_vec = new Vector()
  bAP_peak_v_x_vec = new Vector()

  dendrite for (x,0) {
    bAP_peak_v_vec.append(vmax_ds(x))
    bAP_peak_v_x_vec.append(x*L)
  }
  sprint(tmpstr,"%s/bAP_envelop%s.dat", proto_dir, file_suffix)
  write_vec(tmpstr, bAP_peak_v_x_vec, bAP_peak_v_vec)
  bAP_peak_v_vec.line(graph_bAP_envelopes, bAP_peak_v_x_vec, proto_num%2+1, 1) // proto_num colors avoids white
  graph_bAP_envelopes.exec_menu("View = plot")
  // xopen("hoc/annotate_bAP_envelope_graph.ses") // hoc only opens this once
}

proc run_model() {
  print "model running with protocol ", proto_dir
  setup_vectors_to_record()
  init()
  run()
  write_vectors()
}

// Start of top-level script to run model with X Conditions

// first switch to semi-high (0.05 is high) neck resistance state:
forsec "neck" {diam=0.07}
forsec "neck" {Ra=200}
// function is being broken up into 1)without and 2)with
// functions so that graphs can be made with the vectors while they are
// available.  These new functions immediately follow below

dend_choice = 0 // the special subscript 0 selects the NetCon for the
// dendrite compartment near Spine[1].  Note spine_choice (defined to
// be 2) gives the NetCon for the inhibited middle spine (Spine[1] of
// [0], [1], [2] spines)