#-*- awk -*-


# Change N where N is the desired number of runs
BEGIN {
  N=399;
  # how many times do we step through the model
  TIMESTEPS=200;
  # how fine should the numeric appr be
  UNITSTEP=5;
  ENOUGH=10;
  found=0; 
}
m4_include(includes/base_awk.m4);

# Model:
#    dv/dt = 0.04v^2 + 5v + 140 - u + I
#    du/dt = a(bv - u)
#    if v>=30 then v=c and u=u+d
#    v has scale of mV and t has ms scale

# model parameters
Choice(getaaaVal,=,n){
  return linear(getaaaVal_min,getaaaVal_max,n); }
Choice(getbbbVal,=,n){
  return linear(getbbbVal_min,getbbbVal_max,n); }
Choice(getcccVal,=,n){
  return linear(getcccVal_min,getcccVal_max,n); }
Choice(getdddVal,=,n){
  return linear(getdddVal_min,getdddVal_max,n); }

# return in the range 0 to .1, typical cited value ".02"
function getaVal() {
  #return .02;
  return getaaaVal()/10.0;
}
# return in range 0 to 1, typical cited value ".2"
function getbVal(){
  #return .25;
  return getbbbVal();
}
# range -80 to -40
function getcVal(){
  return -(getcccVal()*40 + 40);
}
# range 0 to 10
function getdVal(){
  return getdddVal()*10;
}

function dv(v,u,current) {
  return .04*v*v + 5*v + 140 - u + current;
}
function du(v,u) {
  return getaVal()*(getbVal()*v - u);
}

# do time stepping in main function
function main(warmup,   n,delv,delu,vt,ut,current) {
  getaaaVal(n);
  getbbbVal(n);
  getcccVal(n);getdddVal(n);
  reset=0;
  # output model parameters 
  if (!warmup)
    print getaVal() " " getbVal() " " getcVal() " " \
      getdVal() >> "parameters.dat";

  # don't randomize the initial conditions
  #vt=getv0Val();ut=getu0Val();
  vt= -64;
  ut= getbVal()*vt;
  sumworth=0.0;
  # this comment is wrong, even though it came from izhikevich!
  # pulse in current is needed to stimulate "neuron"
  current = 0.0;

  if (!warmup)
    print vt " " ut " " current >> "points.dat"

  # variables for detecting phasic spiking

  for (time=0; time<TIMESTEPS; time++) {
    for (deltat=(time==0?1:0); deltat<UNITSTEP; deltat++) {
      if (time>1)
	current=.5;

      delv = dv(vt,ut,current)/UNITSTEP;
      delu = du(vt,ut)/UNITSTEP;
      vt += delv;
      ut += delu;

      # 30 mV is hard max
      if (vt > 30) {vt = 30;reset++;}

      if (!warmup)
	print vt " " ut " " current >> "points.dat"

      # reset equations, contains "spied" variables
      if (vt >= 30) {
	vt = getcVal();
	ut = ut + getdVal();
      }
    }
  }
  if (reset==1) {
    found++;
    return 1;
  } else {
    return 0;
  }
}

function gotEnough() {
  if(found>ENOUGH) return 1;
  else return 0;
}
#############################################