// *** some simple tools 

// globals

objref rect, recv, recapt, recapv, recsp1t, recsp1v, recsp2t, recsp2v, recrbt, recrbv
objref recsrt, recsrv, recmrt, recmrv, recfrt, recfrv, cdf
strdef fnm, cmd

rect    = new Vector(32766)
recv    = new Vector(32766)
recsp1t = new Vector(32766)
recsp1v = new Vector(32766)
recsp2t = new Vector(32766)
recsp2v = new Vector(32766)
recrbt  = new Vector(32766)
recrbv  = new Vector(32766)
recsrt  = new Vector(32766)
recsrv  = new Vector(32766)
recmrt  = new Vector(32766)
recmrv  = new Vector(32766)
recfrt  = new Vector(32766)
recfrv  = new Vector(32766)
recapt = new Vector(1024)
recapv = new Vector(1024)   

// apply TTX: at present we only target Na and NaL 
//            also, only const in soma/dend...

proc applyTTX() {local Namodif,NaLmodif
  // $1 is a percentage of how much TTX influence to apply.
  
  Namodif  = 1.0     // TTX influences this channel 100%
  NaLmodif = 1.0     // TTX also influences NaL
  
  forall(gna_Na=(100.0-($1*Namodif))*default_gNa_dend/100.0)
  SThcells[0].soma.gna_Na = (100.0-($1*Namodif))*default_gNa_soma/100.0
  
  forall(gna_NaL=(100.0-($1*NaLmodif))*default_gNaL_dend/100.0)
  SThcells[0].soma.gna_NaL = (100.0-($1*NaLmodif))*default_gNaL_soma/100.0
  
  printf(" [applying TTX]\n")
}

// wash TTX: return to original values (or near)

proc washTTX() {local washmodif
  // $1 is a percentage of how much TTX gets washed out
  washmodif = ($1/100.0)
  
  forall(gna_Na=washmodif*default_gNa_dend)
  SThcells[0].soma.gna_Na = washmodif*default_gNa_soma
  
  forall(gna_NaL=washmodif*default_gNaL_dend)
  SThcells[0].soma.gna_NaL = washmodif*default_gNaL_soma
  
  printf(" [washing out TTX]\n")
}

// apply Apamin: sKCa antagonist

proc applyApamin() {local sKCamodif,new_gsKCa_soma,new_gsKCa_dl,new_gsKCa_db
  printf(" [applying Apamin]\n")
  
  cset(0,"gk_sKCa","-apamin0.9")
}

// wash Apamin: return to original values

proc washApamin() {local washmodif,new_gsKCa_soma,new_gsKCa_dl,new_gsKCa_db
  printf(" [washing out Apamin]\n")
  
  cset(0,"gk_sKCa","")
}


// isolate a single Action potential (this can then be passed to an AP shape analyser)

func findAP() {local i, c, inAP, nAP, starti, endi, vmin, APOK
 // $1 is time to find first AP after...
 // $2 is rest voltage
 // $o3 is voltage vector
 // $o4 is time vector
 // $o5 is voltage vector to be created
 // $o6 is time vector to be created
 
 len = $o3.size()
 
 // find three APs after $1, and focus on the middle one
 
 inAP=0
 nAP=0
 starti= -1
 APOK=0
 vmin=50
 
 // seach for the starting index of time $1
 for i=0,len-1 {
   if ($o4.x[i]>=$1) {
     // in the "start looking" zone
     if ($o3.x[i]>$2) {
       if (!inAP) {
	 //a new AP
	 inAP=1
	 nAP=nAP+1
	 vmin=50
	 // stop if we are at rise of third AP
	 if (nAP==3) {
	   endi=i
	   APOK=1
	 }	 
       } else {
	 //still above threshold in the AP
       } 
     } else {
       if (inAP) {
	 //we were in an AP, but not now, only reset inAP wen we hit rock bottom
	 if ($o3.x[i]<vmin) {
	   vmin=$o3.x[i]
	 } else {
	   inAP=0
	   if (nAP==1) starti=i
	 }
       }       
     }
   }
   if (APOK) break
 }
 
 // create new vectors just to hold this AP 
 $o5.resize(endi-starti+1)
 $o6.resize(endi-starti+1) 
 
 c=0
 for i=starti,endi { 
   $o5.x[c]=$o3.x[i]
   $o6.x[c]=$o4.x[i]
   c=c+1
 }
 
 return APOK
}


proc newgraph() {local i, ivi
 // $o1 is destination time vector
 // $o2 is destination voltage vector
 // $o3 is the graph object
 // $4  is the start point
 // $5  screen x pos 
 // $6  screen y pos
 // $7  cut and create or just plot?
 // $s8 string label for graph
 
  if ($7>0) {
    ivi= -1
    for i=0,rect.size()-1 {
      if (rect.x[i]>=$4) {
	if (ivi<0) ivi=i
	break
      }
    }     
    
    $o1.copy(rect,ivi,rect.size()-1)
    $o2.copy(recv,ivi,recv.size()-1)
    $o1.resize(rect.size()-ivi)
    $o2.resize(recv.size()-ivi)
  }
  $o3 = new Graph(0)
  $o3.view($4,-75.0-5,20+10,90+10, $5, $6, 325, 215)
  $o3.erase_all()
  $o3.size($4,$o1.x[$o1.size()-1]+5,$o2.min()-5,$o2.max()+5)
  $o2.plot($o3,$o1) 
  $o3.label(0.04, 0.93, $s8)
}


proc set_aCSF() {local req
  // $1 is the requested environment:
  //   =0 NEURON defaults
  //   =3 Beurrier et al (1999)
  //   =4 Bevan & Wilson (1999)
  
  req=$1
  
  // WARNING: if changing cai0_ca_ion, must hand edit Cacum.mod!
  
  // Beurrier et al (1999)
  //  Calculated aCSF
  if (req==3) {
    printf("Setting in vitro parameters based on Beurrier et al (1999)\n")
    
    nai0_na_ion = 15
    nao0_na_ion = 150
    
    ki0_k_ion = 140
    ko0_k_ion = 3.6
    
    cai0_ca_ion = 1e-04
    cao0_ca_ion = 2.4
    
    cli0_cl_ion = 4
    clo0_cl_ion = 135
  }    
  
  // Bevan & Wilson (1999)
  //  Calculated aCSF (same as Bevan et al. 2002, Hallworth et al. 2003)
  if (req==4) {
    printf("Setting in vitro parameters based on Bevan & Wilson (1999)\n")
    
    nai0_na_ion = 15
    nao0_na_ion = 128.5
    
    ki0_k_ion = 140
    ko0_k_ion = 2.5
    
    cai0_ca_ion = 1e-04
    cao0_ca_ion = 2.0
    
    cli0_cl_ion = 4
    clo0_cl_ion = 132.5
  }

  // NEURON's defaults
  if (req==0) {
    printf("WARNING: Using NEURON defaults for in vitro parameters\n")
   
    nai0_na_ion = 10
    nao0_na_ion = 140
    
    ki0_k_ion = 54
    ko0_k_ion = 2.5
    
    cai0_ca_ion = 5e-05
    cao0_ca_ion = 2
    
    cli0_cl_ion = 0
    clo0_cl_ion = 0
  }  
}


proc cset() { local tree, sec, ws, val, tst
  // $1 is the cell
  // $s2 is the string name of the conductance to set
  // $s3 is the string file modifier
  
  if (name_declared($s2)>0) {
    
    sprint(fnm,"sth-data/cell-%s%s",$s2,$s3)
    cdf = new File()
    cdf.ropen(fnm)
 
    while (!cdf.eof()) {
      tree = cdf.scanvar()
      sec = cdf.scanvar()
      ws = cdf.scanvar()
      val = cdf.scanvar()
      
      if (tree==-1) {
	sprint(cmd,"SThcells[%d].soma.%s = %.9f",$1,$s2,val)
      }
      if (tree==0) {
	sprint(cmd,"SThcells[%d].dend0[%d].%s(%.2f) = %.9f",$1,sec,$s2,ws,val)
      }       
      if (tree==1) {
	sprint(cmd,"SThcells[%d].dend1[%d].%s(%.2f) = %.9f",$1,sec,$s2,ws,val)
      }     
      execute(cmd)
    }
    cdf.close()
  }   
}