/* access dend1[7]  -- somatic recording site at 0.5    */
/* access dend1[134] -- dendritic recording site at 0.01 */

strdef somaoutname, dendoutname
somaoutname = "soma.out"
dendoutname = "dend.out"

objectvar elecstim[2]
access dend1[7]
elecstim[0] = new IClamp(0.5)
elecstim[1] = new IClamp(0.5)

func somaticv() {
	return dend1[7].v(0.5)
}

func dendriticv() {
	return dend1[134].v(0.01)
}

proc read_expt_data() {
	if (long>0) {
	  if (qblock<1) {
		somareadfile="ri22scon.asc"
		dendreadfile="ri22dcon.asc"
	  }
	  if (qblock>=1) {
		somareadfile="ri22scs.asc"
		dendreadfile="ri22dcs.asc"
	  }
	}
	if (short>0) {
		somareadfile="ri22scs2.asc"
		dendreadfile="ri22dcs2.asc"
	}
	vecdt = 0.126
	vect = 400
	veclength = int(vect/vecdt)
	print "Reading real somatic data for fit..."
	somarawfile.ropen(somareadfile)
	somarawvec.scanf(somarawfile,veclength,2,2)
	somarawfile.close()
	print "Reading real dendritic data for fit..."
	dendrawfile.ropen(dendreadfile)
	dendrawvec.scanf(dendrawfile,veclength,2,2)
	dendrawfile.close()
	print "Done\n"
}


proc create_axes() {
  xmin=$1
  xmax=$2
  xticks=$3
  ymin=$4
  ymax=$5
  yticks=$6
  plt(-4)
  plt(-3)
  setcolor(1)
  axis(xmin, xmax, xticks, ymin, ymax, yticks)
  axis()
}

proc somastep_cc() {
  dend1[7] elecstim.loc(0.5)
  elecstim.del=$2
  elecstim.dur=$3
  elecstim.amp=$1
  elecstim.i=0
  tmax=$4
  dt=$5
  steps_per_ms=1/dt
  maxvsoma=0
  maxvdend=0
  ssvsoma=0
  ssvdend=0
  peakvdendt=0
  peakvdend=0
  peakvsoma=0
  run()
  rn=abs(ssvsoma/elecstim.amp)
}

proc somastep_cc2() {
  dend1[7] elecstim.loc(0.5)
  elecstim.del=$2
  elecstim.dur=$3
  elecstim.amp=$1
  elecstim.i=0
  tmax=$4
  dt=$5
  steps_per_ms=1/dt
  maxvsoma=0
  maxvdend=0
  ssvsoma=0
  ssvdend=0
  peakvdendt=0
  peakvdend=0
  peakvsoma=0
 /* run() */
 /* rn=abs(ssvsoma/elecstim.amp) */
}

proc doublestep_cc() {
  dend1[7] elecstim[0].loc(0.5)
  elecstim[0].del=$2
  elecstim[0].dur=$3
  elecstim[0].amp=$1
  elecstim[0].i=0
  dend1[7] elecstim[1].loc(0.5)
  elecstim[1].del=$5
  elecstim[1].dur=$6
  elecstim[1].amp=$4
  elecstim[1].i=0
  tmax=$7
  dt=$8
  steps_per_ms=1/dt
  maxvsoma=0
  maxvdend=0
  ssvsoma=0
  ssvdend=0
}

proc stabilize_cc() {
  fstim(0)
  elecstim.amp=0
  otstop=tstop
  tstop=$1
  dt=$2
  steps_per_ms=1/dt
  run()
  tstop=otstop
}

proc findsf() {
  target=$1
  adjust=(rn-target)/target
  oldareascale=areascale
  areascale=areascale+(adjust*areascale)
  printf("rn = %f, old areascale = %f, new areascale = %f\n",rn, oldareascale, areascale)
}

proc findmax() {
    if (abs(dend1[7].v(0.5)-Vrestsoma)>maxvsoma) {
		maxvsoma=abs(dend1[7].v(0.5)-Vrestsoma)
    }
    if (abs(dendriticv()-Vrestdend)>maxvdend) {
		maxvdend=abs(dendriticv()-Vrestdend)
		if (short>0) {
			peakvdendt=t-delay
			peakvsoma=abs(dend1[7].v(0.5)-Vrestsoma)
			peakvdend=abs(dendriticv()-Vrestdend)
		}
    }
    for i=0,599 {
    	if (dend1[i].v(0.5)-vrest[i]>peakv[i]) {
			peakv[i]=dend1[i].v(0.5)-vrest[i]
    	}
    	if (dend1[i].v(0.5)-vrest[i]>0) {
			synint[i]+=(dend1[i].v(0.5)-vrest[i])*dt
    	}
	}
}

proc findss() {
    if (t<3*dt) {
	ssvsoma=0
	ssvdend=0
	ssflag=0
	ndt=0.1*duration/dt
	nt=0
	duration=elecstim.dur
	delay=elecstim.del
    }
    if ((t>delay+duration-ndt*dt) && (t<delay+duration)) {
	ssvsoma=ssvsoma+abs(dend1[7].v(0.5)-Vrestsoma)
	ssvdend=ssvdend+abs(dend1[134].v(0.01)-Vrestdend)
	nt=nt+1
    }
    if ((t>delay+duration) && (ssflag<1)) {
	if (nt>0) {
	  ssvsoma=ssvsoma/nt
	  ssvdend=ssvdend/nt
	}
	ssflag=1
   }
}

proc syn_cc() {
  integral1=0
  integral2=0
  access dend1[134]
  create_syn2(1)
  loc_syn2(0,$1) /* $1 defines location between 0 and 1 synapse */
  {gmax_syn2=$2 tau0_syn2=$3 tau1_syn2=$4 e_syn2=$5 onset_syn2=$6 }
  delay=$6
  tmax=$7
  dt=$8
  t=0
  n=0
  maxvsoma=0
  maxvdend=0
  while (t<=tmax) {
    if (abs(soma.v-Vrestsoma)>maxvsoma) maxvsoma=abs(soma.v-Vrestsoma)
    if (abs(dend1[134].v(0.01)-Vrestdend)>maxvdend) maxvdend=abs(dend1[134].v(0.01)-Vrestdend)
    integral1=integral1+((soma.v-Vrestsoma)*dt)
    integral2=integral2+((dend1[134].v(0.01)-Vrestdend)*dt)
    setcolor(1)
    plot(0,t,soma.v)
    setcolor(2)
    plot(0,t,dend1[7].v(0.5))
/*    save[n]=soma.v */
    n=n+1
/*    write_psp_asc(soma.v)*/
    fadvance(1)
  }
/*  halfwidth() */
}

proc dendtau_cc() {
  integral1=0
  integral2=0
  access dend1[134]
/*  objectvar waveform */
/*  waveform = new tau2(0) */
/*  tau2.loc(0) */
  /* tau2.loc(0,$1) /* $1 defines location between 0 and 1 synapse */
/*  {tau2.imax=$2 tau2.tau0=$3 tau2.tau1=$4 tau2.onset=$5 } */
  delay=$5
  tmax=$6
  dt=$7
  t=0
  n=0
  maxvsoma=0
  maxvdend=0
  while (t<=tmax) {
    if (abs(soma.v-Vrestsoma)>maxvsoma) maxvsoma=abs(soma.v-Vrestsoma)
    if (abs(dend1[134].v(0.01)-Vrestdend)>maxvdend) maxvdend=abs(dend1[134].v(0.01)-Vrestdend)
    integral1=integral1+((soma.v-Vrestsoma)*dt)
    integral2=integral2+((dend1[134].v(0.01)-Vrestdend)*dt)
    setcolor(1)
    plot(0,t,soma.v)
    setcolor(2)
    plot(0,t,dend1[134].v(0.01))
    setcolor(3)
    plot(0,t,(i_tau2*10)-20)
/*    save[n]=soma.v */
    n=n+1
/*    write_psp_asc(soma.v) */
    fadvance(1)
  }
/*  halfwidth() */
}

proc sum_error() {
	if (short>0) {
		tstarts=7
		tstartd=5
	}
	if (long>0) {
		tstarts=50
		tstartd=50
	}
	if (both>0) {
		tstarts=50
		tstartd=50
	}
	if (t>tstarts) {
		mserrsoma=mserrsoma+($1-somarawv)^2
	}
	if (t>tstartd) {
		mserrdend=mserrdend+($2-dendrawv)^2
	}
}

proc avg_error() {
	mserrsoma2=mserrsoma/((tstop-tstarts)/dt)
	mserrdend2=mserrdend/((tstop-tstartd)/dt)
	mserrsoma=sqrt(mserrsoma2)
	mserrdend=sqrt(mserrdend2)
	mserr=(mserrsoma+mserrdend)/2
}

proc write_data() {
	printf("Writing output files...\n")
	outfile.wopen(somaoutname)
     /*	somasimvec.printf(outfile,"%6.3f\n") */
	somasimvec.vwrite(outfile)
	outfile.close()
	outfile.wopen(dendoutname)
     /*	dendsimvec.printf(outfile,"%6.3f\n") */
	dendsimvec.vwrite(outfile)
	outfile.close()
}