// $Id: control_spiketuft.hoc,v 1.10 2007/03/17 01:39:06 ted Exp $

// based on control_probetuft.hoc

// doit() runs a simulation, then plots max v in tuft vs. distance from soma


objref paramfil, areafil, resultfil, vpeakfil
strdef paramfilname, areafilname, resultfilname, vpeakfilname
paramfilname="parameters.dat"
areafilname="areas.dat"
resultfilname="results.dat"
vpeakfilname="vpeak.dat"

proc write_paramfil() {  local tally
  // open, write, and close paramfil
  paramfil = new File()
  paramfil.wopen(paramfilname)

  paramfil.printf("Parameters\n")
  paramfil.printf("Model cell\n")
  paramfil.printf("Ra\t cm\t e_pas\t g_pas\n")
  paramfil.printf("%5.1f\t %5.3g\t %5.1g\t %g\n", Ra, cm, e_pas, g_pas)
  tally = 0
  forsec tuft tally+=nseg
  paramfil.printf("%d segments in tuft\n", tally)

  paramfil.printf("Tuft origin clamped to follow a spike waveform\n")
  paramfil = new File()
}

proc write_areafil() { local ii
  // open, write, and close areafil
  areafil = new File()
  areafil.wopen(areafilname)
  areafil.printf("Segment areas are tab separated in same order as segments in tuft SectionList\n")
  areavec.printf(areafil, "\t %9.4f")
  areafil = new File()
}


proc output_to_files() { local ii, jj
print "in output_to_files"

  write_paramfil()
  write_areafil()
  resultfil = new File()
  resultfil.wopen(resultfilname)
  resultfil.printf("Simulation results\n")
  resultfil.printf("\t  mean\t  min\t  max\t  var\t  stdev\n")
  vpeakfil = new File()
  vpeakfil.wopen(vpeakfilname)
  vpeakfil.printf("Peak depolarization\n")
  vpeakfil.printf("Values are tab separated in same sequence as segments in tuft SectionList\n")

// we have already analyzed the results!  no need to do it again!
//  analyze()

  resultvec.printf(resultfil, "\t %7.4f")
  // no point printing run number as first item on the line
  vpeakvec.printf(vpeakfil, "\t %7.4f")

  resultfil = new File()
  vpeakfil = new File()
}

///////////////////////////

// nice controls

BOTHCONTROL = 0
CMx2 = 1
RAx2 = 2
BOTHx2 = 3

status = BOTHCONTROL

proc restoreRa() {
  forall Ra/=2
}

proc restorecm() {
  forall for (x,0) cm(x)/=2
}

proc doubleRa() {
  forall Ra*=2
}

proc doublecm() {
  forall for (x,0) cm(x)*=2
}


proc setparams() {
  if ($1==BOTHCONTROL) {
    // if (status==BOTHCONTROL) 
    if (status==CMx2) restorecm()
    if (status==RAx2) restoreRa()
    if (status==BOTHx2) { restorecm()  restoreRa() }
  }
  if ($1==CMx2) {
    if (status==BOTHCONTROL) doublecm()
    // if (status==CMx2)
    if (status==RAx2) { doublecm()  restoreRa() }
    if (status==BOTHx2) restoreRa()
  }
  if ($1==RAx2) {
    if (status==BOTHCONTROL) doubleRa()
    if (status==CMx2) { doubleRa()  restorecm() }
    // if (status==RAx2)
    if (status==BOTHx2) restorecm()
  }
  if ($1==BOTHx2) {
    if (status==BOTHCONTROL) { doublecm()  doubleRa() }
    if (status==CMx2) doubleRa()
    if (status==RAx2) doublecm()
    // if (status==BOTHx2)
  }
  status = $1
}


proc doit() {
  run()
  plotvmax()
  plotvmaxnorm(1)

  analyze()
  plotpercentiles()

  output_to_files()
}


{
xpanel("Simulation control", 0)
xlabel("Parameters (default is control)")
xradiobutton("control", "setparams(BOTHCONTROL)")
xradiobutton("cm = 2*control", "setparams(CMx2)")
xradiobutton("Ra = 2*control", "setparams(RAx2)")
xradiobutton("both cm and Ra = 2*control", "setparams(BOTHx2)")
xbutton("Execute a run","doit()")
xpanel(2,507)
}