load_file("gNa_gK_gCa.hoc")
stim.amp = 0.2 // to ensure threshold is reached

objectvar APfile, measuredAP, measuredAPfitrange
objectvar modelAP, modelAPfitrange
objectvar outfile[3]
double parameters[2]

APfile = new File()
APfile.ropen("tAP1_2us.txt") /* peak position: 200  (counting from 0) */
measuredAP = new Vector()
measuredAP.scanf(APfile)
APfile.close()
measuredAPfitrange = measuredAP.c(0, 550) /* from 0.4 ms before the peak to 0.7 ms after */

forall insert pk
mfb[1] distance(0, 0.5)
mfb[3] dist13 = distance(0.5) /* microns */


func distancefunction() {
gNaMFB = abs($&2[0])   /* mS/cm^2 */
gNaaxon = abs($&2[0])  /* mS/cm^2 */
gKMFB = abs($&2[1])    /* mS/cm^2 */
gKaxon = abs($&2[1])   /* mS/cm^2 */
print "threshold = ", th, " gNaMFB = ", gNaMFB, " gNaaxon = ", gNaaxon, " gKMFB = ", gKMFB, " gKaxon = ", gKaxon

/* gNa */
forall vthreshold_GNa = th
forsec "mfb" gnabar_GNa = gNaMFB/1000    /* S/cm^2 */
forsec "filext" gnabar_GNa = gNaMFB/1000 /* S/cm^2 */
forsec "axon" gnabar_GNa = gNaaxon/1000  /* S/cm^2 */

/* gK */
forall vthreshold_GK = th
forsec "mfb" gkbar_GK = gKMFB/1000    /* S/cm^2*/
forsec "filext" gkbar_GK = gKMFB/1000 /* S/cm^2*/
forsec "axon" gkbar_GK = gKaxon/1000  /* S/cm^2*/

/* gCa */
forsec "mfb" vthreshold_GCa = th
forsec "filext" vthreshold_GCa = th

modelAP = new Vector()
modelAP.record(&mfb[2].v(0.5))

run()

mfb[1] t1 = tpeak_pk(0.5)
mfb[2] t2 = tpeak_pk(0.5)
mfb[3] t3 = tpeak_pk(0.5)

modelAPpeak = int(t2/dt + 1e-6)
if (modelAPpeak > 1150) modelAPpeak = 600 // do not stop the optimization if the AP occurs too late, but generate large APmeansquareerror if that happens
modelAPfitrange = modelAP.c(modelAPpeak - 200, modelAPpeak + 350)

APmeansquareerror = measuredAPfitrange.meansqerr(modelAPfitrange)
print "Mean square error of AP waveform = ", APmeansquareerror

t13 = 1000*(t3 - t1) /* microseconds */
if (t13 > 0) {
	velocity = dist13/t13 /* m/s */
} else {
	velocity = 0 // failed AP
}
print "velocity = ", velocity

return APmeansquareerror
}

proc guess() {
	parameters[0] = 50
	parameters[1] = 7
}

proc fit() {
	attr_praxis(0.00001, 5, 3)
	mindistance = fit_praxis(2, "distancefunction", &parameters[0])
	mindistance = distancefunction(2, &parameters[0])
}

for i=0,2 outfile[i] = new File()

proc output() {
	outfile[0].wopen("parameters1.txt")
	outfile[1].wopen("measuredAP1.txt")
	outfile[2].wopen("modelAP1.txt")

	for (th = -62; th <= -50; th += 1) {
		guess()
		fit()

		outfile[0].printf("th = %g\ngNaMFB = %g\ngNaaxon = %g\ngKMFB = %g\ngKaxon = %g\n", th, abs(parameters[0]), abs(parameters[0]), abs(parameters[1]), abs(parameters[1]))
		outfile[0].printf("mean square error of AP waveform = %g\n", APmeansquareerror)
		outfile[0].printf("velocity = %g\n", velocity)
		outfile[0].printf("\n")

		measuredAPfitrange.printf(outfile[1], "%g\n")
		modelAPfitrange.printf(outfile[2], "%g\n")
	}
	for i=0,2 outfile[i].close()
}

output()
quit()