#pragma rtGlobals=1		// Use modern global access method.
print;timestamp()

Macro timestamp()
	print Secs2Date(DateTime,0),Secs2Time(DateTime,1)
End

// ==================== Fitting functions ======================================

function cf_ln(coefs,t)	//	elimination, one-compartment model, logarithms 19/7/2021
  wave coefs
  variable t
  return ln(coefs[0]*exp(-coefs[1]*t))
end function

function cf_bo(coefs,t)	//	elimination, one-compartment model, 17/7/2021
  wave coefs
  variable t
  return coefs[0]*exp(-coefs[1]*t)
end function

// models naming scheme: 1st letter is b, z, or f for bolus, zero-, or first-order absoprtion
// 2nd letter is o or w for One or tWo compartment models
// t or t1, t2, t3 for one or multiple termination times
// _k are extensions indicating modification of two-compartment models
// parameters are given as k10 and k12 instead of a and b (or l1 and l2), although the latter are the actual parameters used in the calculations

function cf_fo0(coefs,t)	//	classical absorption-elimination, one-compartment model, ka = ke 20/6/2021
  wave coefs		// coefs[0] = F D / Vd, coefs[1]=k
  variable t
  return coefs[0]*coefs[1]*t*exp(-coefs[1]*t)
end function

function cf_fo0t(coefs,t)	//	truncated classical absorption-elimination, one-compartment model, ka = ke 20/6/2021
  wave coefs		// coefs[0] = F D / Vd, coefs[1]=k, coefs[2]=tau
  variable t
  variable y
  if (t<=coefs[2])
    y = coefs[0]*coefs[1]*t*exp(-coefs[1]*t)
  else
    y = cf_fo0t(coefs,coefs[2])*exp(-coefs[1]*(t-coefs[2]))
  endif
  return y
end function

function cf_fo(coefs,t)	//	classical absorption-elimination, one-compartment model 7/4/2021
  wave coefs		// coefs[0] = F D / Vd, coefs[1]=ka, coefs[2]=kel
  variable t
  return coefs[0]*coefs[1]/(coefs[1]-coefs[2])*(exp(-coefs[2]*t)-exp(-coefs[1]*t))
end function

function cf_fot(coefs,t)	//	truncated classical absorption-elimination, one-compartment model 25/4/2021
  wave coefs		// coefs[0] = F D / Vd, coefs[1]=ka, coefs[2]=kel, coefs[3]=tau
  variable t
  variable y
  if (t<=coefs[3])
    y=coefs[0]*coefs[1]/(coefs[1]-coefs[2])*(exp(-coefs[2]*t)-exp(-coefs[1]*t))
  else
    y=cf_fot(coefs,coefs[3])*exp(-coefs[2]*(t-coefs[3]))
  endif
  return y
end function

function cf_zot1(coefs,t)		// finite, zero-order absorption, one-compartment elimination 16/3/2021
// changed to compact format on 7/4/2021
  wave coefs		// coefs[0] = F D / Vd, coefs[1]=kel, coefs[2]=tau
  variable t
  variable y
  if (t<=coefs[2])
    y= coefs[0]/coefs[1]/coefs[2]*(1-exp(-coefs[1]*t))	// /coefs[2] added 10/4/2021
  else
    y= cf_zot1(coefs,coefs[2])*exp(-coefs[1]*(t-coefs[2]))
  endif
  return y
end function

function cd_zot1(coefs,t)		// finite, zero-order absorption, one-compartment elimination, time derivative of cf_zot1 2022-01-16
  wave coefs		// coefs[0] = F D / Vd, coefs[1]=kel, coefs[2]=tau
  variable t
  variable y
  if (t<=coefs[2])
    y= coefs[0]/coefs[2]*exp(-coefs[1]*t)
  else
    y= -cf_zot1(coefs,coefs[2])*coefs[1]*exp(-coefs[1]*(t-coefs[2]))
  endif
  return y
end function

function cf_zot2(coefs,t)		// finite, 2-stage, zero-order absorption, one-compartment elimination 5/4/2021
// 2021-04-07 compact syntax for previous segments
// corrected expression for elimination of first term 7/4/2021
  wave coefs
  variable t
  variable y
  if (t<=coefs[1])
    y= coefs[0]/coefs[4]/coefs[1]*(1-exp(-coefs[4]*t))	// /coefs[1] added 10/4/2021
  else
    if (t<=coefs[1]+coefs[3])
      y= cf_zot2(coefs,coefs[1])*exp(-coefs[4]*(t-coefs[1])) + coefs[2]/coefs[4]/coefs[3]*(1-exp(-coefs[4]*(t-coefs[1])))	// /coefs[3] added 10/4/2021
    else
      y= cf_zot2(coefs,coefs[1]+coefs[3])*exp(-coefs[4]*(t-coefs[1]-coefs[3]))
    endif
  endif
  return y
end function

function cf_zot3(coefs,t)		// finite, 3-stage, zero-order absorption, one-compartment elimination 7/4/2021
  wave coefs
  variable t
  variable y
  if (t<=coefs[1])
    y= coefs[0]/coefs[6]/coefs[1]*(1-exp(-coefs[6]*t))	// /coefs[1] added 10/4/2021
  else
    if (t<=coefs[1]+coefs[3])
      y= cf_zot3(coefs,coefs[1])*exp(-coefs[6]*(t-coefs[1])) + coefs[2]/coefs[6]/coefs[3]*(1-exp(-coefs[6]*(t-coefs[1])))	// /coefs[3] added 10/4/2021
    else
      if (t<=coefs[1]+coefs[3]+coefs[5])
        y= cf_zot3(coefs,coefs[1]+coefs[3])*exp(-coefs[6]*(t-coefs[1]-coefs[3])) + coefs[4]/coefs[6]/coefs[5]*(1-exp(-coefs[6]*(t-coefs[1]-coefs[3])))	// /coefs[5] added 10/4/2021
      else
        y= cf_zot3(coefs,coefs[1]+coefs[3]+coefs[5])*exp(-coefs[6]*(t-coefs[1]-coefs[3]-coefs[5]))
      endif
    endif
  endif
  return y
end function

function cd_zot3(coefs,t)		// finite, 3-stage, zero-order absorption, one-compartment elimination time derivative 2022-01-16
  wave coefs
  variable t
  variable y
  if (t<=coefs[1])
    y= coefs[0]/coefs[1]*exp(-coefs[6]*t)
  else
    if (t<=coefs[1]+coefs[3])
      y= -cf_zot3(coefs,coefs[1])*coefs[6]*exp(-coefs[6]*(t-coefs[1])) + coefs[2]/coefs[3]*exp(-coefs[6]*(t-coefs[1]))
    else
      if (t<=coefs[1]+coefs[3]+coefs[5])
        y= -cf_zot3(coefs,coefs[1]+coefs[3])*coefs[6]*exp(-coefs[6]*(t-coefs[1]-coefs[3])) + coefs[4]/coefs[5]*exp(-coefs[6]*(t-coefs[1]-coefs[3]))
      else
        y= -cf_zot3(coefs,coefs[1]+coefs[3]+coefs[5])*coefs[6]*exp(-coefs[6]*(t-coefs[1]-coefs[3]-coefs[5]))
      endif
    endif
  endif
  return y
end function

function cf_zot4(coefs,t)		// finite, 4-stage, zero-order absorption, one-compartment elimination 6/9/2022
  wave coefs
  variable t
  variable y
  if (t<=coefs[1])
    y= coefs[0]/coefs[8]/coefs[1]*(1-exp(-coefs[8]*t))
  else
    if (t<=coefs[1]+coefs[3])
      y= cf_zot4(coefs,coefs[1])*exp(-coefs[8]*(t-coefs[1])) + coefs[2]/coefs[8]/coefs[3]*(1-exp(-coefs[8]*(t-coefs[1])))
    else
      if (t<=coefs[1]+coefs[3]+coefs[5])
        y= cf_zot4(coefs,coefs[1]+coefs[3])*exp(-coefs[8]*(t-coefs[1]-coefs[3])) + coefs[4]/coefs[8]/coefs[5]*(1-exp(-coefs[8]*(t-coefs[1]-coefs[3])))
      else
        if (t<=coefs[1]+coefs[3]+coefs[5]+coefs[7])
          y= cf_zot4(coefs,coefs[1]+coefs[3]+coefs[5])*exp(-coefs[8]*(t-coefs[1]-coefs[3]-coefs[5])) + coefs[6]/coefs[8]/coefs[7]*(1-exp(-coefs[8]*(t-coefs[1]-coefs[3]-coefs[5])))
        else
          y= cf_zot4(coefs,coefs[1]+coefs[3]+coefs[5]+coefs[7])*exp(-coefs[8]*(t-coefs[1]-coefs[3]-coefs[5]-coefs[7]))
        endif
      endif
    endif
  endif
  return y
end function

function cf_bwe(coefs,t)	// double exponential decay 2022-01-27
  wave coefs
  variable t
  variable y
  y= coefs[0]*exp(-coefs[2]*t)+coefs[1]*exp(-coefs[3]*t)
  return y
end function

function cf_bwet(coefs,t)	// double exponential decay 2022-01-28
  wave coefs
  variable t
  variable y
  y= coefs[0]*exp(-coefs[2]*(t-coefs[4]))+coefs[1]*exp(-coefs[3]*(t-coefs[4]))
  return y
end function

function cf_bw(coefs,t)	// two-compartment elimination with peripheral calculation 2022-05-31 [improvement of 2021-08-18 version]
  wave coefs
  variable t
  return real(cf_bwcp(coefs,t))
end

function/c cf_bwcp(coefs,t)	// two-compartment elimination with peripheral calculation 2022-05-31 [improvement of 2021-08-18 version]
  wave coefs
  variable t
  variable k12=coefs[2]+coefs[3]-coefs[1]-coefs[2]*coefs[3]/coefs[1]
  variable y,z
  y= coefs[0]/(coefs[2]-coefs[3])*((coefs[1]-coefs[3])*exp(-coefs[3]*t) - (coefs[1]-coefs[2])*exp(-coefs[2]*t))
  z= coefs[0]/(coefs[2]-coefs[3])*k12*(exp(-coefs[3]*t) - exp(-coefs[2]*t))
  return cmplx(y,z)
end function

function cf_bw_k(coefs,t)	// two-compartment elimination with peripheral calculation with k12 and k10 (not a and b) 2022-05-31 [improvement of 2021-08-18 version]
  wave coefs
  variable t
  return real(cf_bw_kcp(coefs,t))
end

function/c cf_bw_kcp(coefs,t)	// two-compartment elimination with peripheral calculation with k12 and k10 (not a and b) 2022-05-31 [improvement of 2021-08-18 version]
  wave coefs
  variable t
  variable k12=coefs[2],k10=coefs[3],k21=coefs[1]
  variable b=k21+k12+k10,c=k21*k10,d=sqrt(b*b-4*c),l1=0.5*(b-d),l2=0.5*(b+d)		// b,c,d nomenclature from trinomial
  variable y,z
  y= coefs[0]/(l1-l2)*((k21-l2)*exp(-l2*t) - (k21-l1)*exp(-l1*t))
  z= coefs[0]/(l1-l2)*k12*(exp(-l2*t) - exp(-l1*t))
  return cmplx(y,z)
end function

function cf_zwt1(coefs,t)	// finite, 1-stage, zero-order absorption, two-compartment elimination 6/12/2021
// old format without peripheral compartment calculation 9/4/2021, modified 20/7/2021
  wave coefs
  variable t
  return real(cf_zwt1cp(coefs,t))
end

function/c cf_zwt1cp(coefs,t)		// finite, 1-stage, zero-order absorption, two-compartment elimination 17/11/2021
// xxxcp functions are used for two-compartment models.
// they return both central (c) and peripheral (p) compartment concentrations in complex number format
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable l1=coefs[3],l2=coefs[4],k21=coefs[2],k10=l1*l2/k21,k12=l1+l2-k21-k10
  if (t<=coefs[1])
    y= coefs[0]/coefs[1]/(coefs[4]-coefs[3])*((coefs[2]-coefs[3])/coefs[3]*(1-exp(-coefs[3]*t)) - (coefs[2]-coefs[4])/coefs[4]*(1-exp(-coefs[4]*t)))
    z= coefs[0]*k12/coefs[1]/(coefs[4]-coefs[3])*(1/coefs[3]*(1-exp(-coefs[3]*t)) - 1/coefs[4]*(1-exp(-coefs[4]*t)))
  else
    yt=real(cf_zwt1cp(coefs,coefs[1]))
    zt=imag(cf_zwt1cp(coefs,coefs[1]))
    tt=t-coefs[1]
    y=1/(coefs[4]-coefs[3])*(coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))-yt*(coefs[3]*exp(-coefs[3]*tt)-coefs[4]*exp(-coefs[4]*tt)))
    z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
    z*=k12/(coefs[4]-coefs[3])
  endif
  return cmplx(y,z)
end function

function cf_zwt1_k(coefs,t)	//	finite, 1-stage, zero-order absorption, two-compartment elimination with ks 2021-12-22
// old format without peripheral compartment calcuation: finite, 1-stage, zero-order absorption, two-compartment elimination (k12,k10, not l1,l2) 16/8/2021
  wave coefs
  variable t
  return real(cf_zwt1_kcp(coefs,t))
end

function/c cf_zwt1_kcp(coefs,t)	//	finite, 1-stage, zero-order absorption, two-compartment elimination with k12, k10 2021-12-22
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable k21=coefs[2],k12=coefs[3],k10=coefs[4],t1=coefs[1]
  variable b=coefs[2]+coefs[3]+coefs[4],c=coefs[2]*coefs[4],d=sqrt(b*b-4*c),l1=0.5*(b-d),l2=0.5*(b+d)		// b,c,d nomenclature from trinomial
  if (t<=t1)
    y= coefs[0]/t1/(l2-l1)*((k21-l1)/l1*(1-exp(-l1*t)) - (k21-l2)/l2*(1-exp(-l2*t)))
    z= coefs[0]*k12/t1/(l2-l1)*(1/l1*(1-exp(-l1*t)) - 1/l2*(1-exp(-l2*t)))
  else
    yt=real(cf_zwt1_kcp(coefs,t1))
    zt=imag(cf_zwt1_kcp(coefs,t1))
    tt=t-t1
    y=1/(l2-l1)*(k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))-yt*(l1*exp(-l1*tt)-l2*exp(-l2*tt)))
    z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(coefs[2]-l2)*zt)*exp(-l2*tt)
    z*=k12/(l2-l1)
  endif
  return cmplx(y,z)
end function

function cf_zwt2(coefs,t)	//	finite, 2-stage, zero-order absorption, two-compartment elimination 5/12/2021
  wave coefs
  variable t
  return real(cf_zwt2cp(coefs,t))
end

function/c cf_zwt2cp(coefs,t)	//	finite, 2-stage, zero-order absorption, two-compartment elimination 5/12/2021
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable l1=coefs[3],l2=coefs[4],k21=coefs[2],k10=l1*l2/k21,k12=l1+l2-k21-k10
  if (t<=coefs[1])
    y= coefs[0]/coefs[1]/(coefs[4]-coefs[3])*((coefs[2]-coefs[3])/coefs[3]*(1-exp(-coefs[3]*t)) - (coefs[2]-coefs[4])/coefs[4]*(1-exp(-coefs[4]*t)))
    z= coefs[0]*k12/coefs[1]/(coefs[4]-coefs[3])*(1/coefs[3]*(1-exp(-coefs[3]*t)) - 1/coefs[4]*(1-exp(-coefs[4]*t)))
  else
    if (t<=coefs[1]+coefs[6])
      yt=real(cf_zwt2cp(coefs,coefs[1]))
      zt=imag(cf_zwt2cp(coefs,coefs[1]))
      tt=t-coefs[1]
      y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
      y+=coefs[5]/coefs[6]*((coefs[2]/coefs[3]-1)*(1-exp(-coefs[3]*tt))-(coefs[2]/coefs[4]-1)*(1-exp(-coefs[4]*tt)))
      y/=coefs[4]-coefs[3]
      z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
      z+=coefs[5]/coefs[6]*((1-exp(-coefs[3]*tt))/coefs[3]-(1-exp(-coefs[4]*tt))/coefs[4])
      z*=k12/(coefs[4]-coefs[3])
    else
      yt=real(cf_zwt2cp(coefs,coefs[1]+coefs[6]))
      zt=imag(cf_zwt2cp(coefs,coefs[1]+coefs[6]))
      tt=t-coefs[1]-coefs[6]
      y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
      y/=coefs[4]-coefs[3]
      z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
      z*=k12/(coefs[4]-coefs[3])
    endif
  endif
  return cmplx(y,z)
end function

function cf_zwt2_k(coefs,t)	//	finite, 2-stage, zero-order absorption, two-compartment elimination 5/12/2021
  wave coefs
  variable t
  return real(cf_zwt2_kcp(coefs,t))
end

function/c cf_zwt2_kcp(coefs,t)	//	finite, 2-stage, zero-order absorption, two-compartment elimination 5/12/2021, with ks 2021-12-22
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable k21=coefs[2],k12=coefs[3],k10=coefs[4],t1=coefs[1],t2=coefs[6]
  variable b=coefs[2]+coefs[3]+coefs[4],c=coefs[2]*coefs[4],d=sqrt(b*b-4*c),l1=0.5*(b-d),l2=0.5*(b+d)		// b,c,d nomenclature from trinomial
  if (t<=t1)
    y= coefs[0]/t1/(l2-l1)*((k21-l1)/l1*(1-exp(-l1*t)) - (k21-l2)/l2*(1-exp(-l2*t)))
    z= coefs[0]*k12/t1/(l2-l1)*(1/l1*(1-exp(-l1*t)) - 1/l2*(1-exp(-l2*t)))
  else
    if (t<=t1+t2)
      yt=real(cf_zwt2_kcp(coefs,t1))
      zt=imag(cf_zwt2_kcp(coefs,t1))
      tt=t-t1
      y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
      y+=coefs[5]/t2*((k21/l1-1)*(1-exp(-l1*tt))-(k21/l2-1)*(1-exp(-l2*tt)))
      y/=l2-l1
      z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
      z+=coefs[5]/t2*((1-exp(-l1*tt))/l1-(1-exp(-l2*tt))/l2)
      z*=k12/(l2-l1)
    else
      yt=real(cf_zwt2_kcp(coefs,t1+t2))
      zt=imag(cf_zwt2_kcp(coefs,t1+t2))
      tt=t-t1-t2
      y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
      y/=l2-l1
      z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
      z*=k12/(l2-l1)
    endif
  endif
  return cmplx(y,z)
end function

function cf_zwt3(coefs,t)	//	finite, 3-stage, zero-order absorption, two-compartment elimination 5/12/2021
  wave coefs
  variable t
  return real(cf_zwt3cp(coefs,t))
end

function/c cf_zwt3cp(coefs,t)		//	finite, 3-stage, zero-order absorption, two-compartment elimination 5/12/2021
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable l1=coefs[3],l2=coefs[4],k21=coefs[2],k10=l1*l2/k21,k12=l1+l2-k21-k10
  if (t<=coefs[1])
    y= coefs[0]/coefs[1]/(coefs[4]-coefs[3])*((coefs[2]-coefs[3])/coefs[3]*(1-exp(-coefs[3]*t)) - (coefs[2]-coefs[4])/coefs[4]*(1-exp(-coefs[4]*t)))
    z= coefs[0]*k12/coefs[1]/(coefs[4]-coefs[3])*(1/coefs[3]*(1-exp(-coefs[3]*t)) - 1/coefs[4]*(1-exp(-coefs[4]*t)))
  else
    if (t<=coefs[1]+coefs[6])
      yt=real(cf_zwt3cp(coefs,coefs[1]))
      zt=imag(cf_zwt3cp(coefs,coefs[1]))
      tt=t-coefs[1]
      y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
      y+=coefs[5]/coefs[6]*((coefs[2]/coefs[3]-1)*(1-exp(-coefs[3]*tt))-(coefs[2]/coefs[4]-1)*(1-exp(-coefs[4]*tt)))
      y/=coefs[4]-coefs[3]
      z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
      z+=coefs[5]/coefs[6]*((1-exp(-coefs[3]*tt))/coefs[3]-(1-exp(-coefs[4]*tt))/coefs[4])
      z*=k12/(coefs[4]-coefs[3])
    else
      if (t<=coefs[1]+coefs[6]+coefs[8])
        yt=real(cf_zwt3cp(coefs,coefs[1]+coefs[6]))
        zt=imag(cf_zwt3cp(coefs,coefs[1]+coefs[6]))
        tt=t-coefs[1]-coefs[6]
        y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
        y+=coefs[7]/coefs[8]*((coefs[2]/coefs[3]-1)*(1-exp(-coefs[3]*tt))-(coefs[2]/coefs[4]-1)*(1-exp(-coefs[4]*tt)))
        y/=coefs[4]-coefs[3]
        z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
        z+=coefs[7]/coefs[8]*((1-exp(-coefs[3]*tt))/coefs[3]-(1-exp(-coefs[4]*tt))/coefs[4])
        z*=k12/(coefs[4]-coefs[3])
      else
        yt=real(cf_zwt3cp(coefs,coefs[1]+coefs[6]+coefs[8]))
        zt=imag(cf_zwt3cp(coefs,coefs[1]+coefs[6]+coefs[8]))
        tt=t-coefs[1]-coefs[6]-coefs[8]
        y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
        y/=coefs[4]-coefs[3]
        z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
        z*=k12/(coefs[4]-coefs[3])
      endif
    endif
  endif
  return cmplx(y,z)
end function

function cf_zwt3_k(coefs,t)	//	finite, 3-stage, zero-order absorption, two-compartment elimination 5/12/2021 with ks 2021-12-22
  wave coefs
  variable t
  return real(cf_zwt3_kcp(coefs,t))
end

function/c cf_zwt3_kcp(coefs,t)	//	finite, 3-stage, zero-order absorption, two-compartment elimination 5/12/2021 with ks 2021-12-22
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable k21=coefs[2],k12=coefs[3],k10=coefs[4],t1=coefs[1],t2=coefs[6],t3=coefs[8]
  variable b=coefs[2]+coefs[3]+coefs[4],c=coefs[2]*coefs[4],d=sqrt(b*b-4*c),l1=0.5*(b-d),l2=0.5*(b+d)		// b,c,d nomenclature from trinomial
  if (t<=t1)
    y= coefs[0]/t1/(l2-l1)*((k21-l1)/l1*(1-exp(-l1*t)) - (k21-l2)/l2*(1-exp(-l2*t)))
    z= coefs[0]*k12/t1/(l2-l1)*(1/l1*(1-exp(-l1*t)) - 1/l2*(1-exp(-l2*t)))
  else
    if (t<=t1+t2)
      yt=real(cf_zwt3_kcp(coefs,t1))
      zt=imag(cf_zwt3_kcp(coefs,t1))
      tt=t-t1
      y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
      y+=coefs[5]/t2*((k21/l1-1)*(1-exp(-l1*tt))-(k21/l2-1)*(1-exp(-l2*tt)))
      y/=l2-l1
      z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
      z+=coefs[5]/t2*((1-exp(-l1*tt))/l1-(1-exp(-l2*tt))/l2)
      z*=k12/(l2-l1)
    else
      if (t<=t1+t2+t3)
        yt=real(cf_zwt3_kcp(coefs,t1+t2))
        zt=imag(cf_zwt3_kcp(coefs,t1+t2))
        tt=t-t1-t2
        y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
        y+=coefs[7]/t3*((k21/l1-1)*(1-exp(-l1*tt))-(k21/l2-1)*(1-exp(-l2*tt)))
        y/=l2-l1
        z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
        z+=coefs[7]/t3*((1-exp(-l1*tt))/l1-(1-exp(-l2*tt))/l2)
        z*=k12/(l2-l1)
      else
        yt=real(cf_zwt3_kcp(coefs,t1+t2+t3))
        zt=imag(cf_zwt3_kcp(coefs,t1+t2+t3))
        tt=t-t1-t2-t3
        y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
        y/=l2-l1
        z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
        z*=k12/(l2-l1)
      endif
    endif
  endif
  return cmplx(y,z)
end function

function cf_zwt4(coefs,t)	//	finite, 4-stage, zero-order absorption, two-compartment elimination 6/9/2022
  wave coefs
  variable t
  return real(cf_zwt4cp(coefs,t))
end

function/c cf_zwt4cp(coefs,t)		//	finite, 4-stage, zero-order absorption, two-compartment elimination 6/9/2022
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable l1=coefs[3],l2=coefs[4],k21=coefs[2],k10=l1*l2/k21,k12=l1+l2-k21-k10
  if (t<=coefs[1])
    y= coefs[0]/coefs[1]/(coefs[4]-coefs[3])*((coefs[2]-coefs[3])/coefs[3]*(1-exp(-coefs[3]*t)) - (coefs[2]-coefs[4])/coefs[4]*(1-exp(-coefs[4]*t)))
    z= coefs[0]*k12/coefs[1]/(coefs[4]-coefs[3])*(1/coefs[3]*(1-exp(-coefs[3]*t)) - 1/coefs[4]*(1-exp(-coefs[4]*t)))
  else
    if (t<=coefs[1]+coefs[6])
      yt=real(cf_zwt4cp(coefs,coefs[1]))
      zt=imag(cf_zwt4cp(coefs,coefs[1]))
      tt=t-coefs[1]
      y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
      y+=coefs[5]/coefs[6]*((coefs[2]/coefs[3]-1)*(1-exp(-coefs[3]*tt))-(coefs[2]/coefs[4]-1)*(1-exp(-coefs[4]*tt)))
      y/=coefs[4]-coefs[3]
      z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
      z+=coefs[5]/coefs[6]*((1-exp(-coefs[3]*tt))/coefs[3]-(1-exp(-coefs[4]*tt))/coefs[4])
      z*=k12/(coefs[4]-coefs[3])
    else
      if (t<=coefs[1]+coefs[6]+coefs[8])
        yt=real(cf_zwt4cp(coefs,coefs[1]+coefs[6]))
        zt=imag(cf_zwt4cp(coefs,coefs[1]+coefs[6]))
        tt=t-coefs[1]-coefs[6]
        y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
        y+=coefs[7]/coefs[8]*((coefs[2]/coefs[3]-1)*(1-exp(-coefs[3]*tt))-(coefs[2]/coefs[4]-1)*(1-exp(-coefs[4]*tt)))
        y/=coefs[4]-coefs[3]
        z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
        z+=coefs[7]/coefs[8]*((1-exp(-coefs[3]*tt))/coefs[3]-(1-exp(-coefs[4]*tt))/coefs[4])
        z*=k12/(coefs[4]-coefs[3])
      else
        if (t<=coefs[1]+coefs[6]+coefs[8]+coefs[10])
          yt=real(cf_zwt4cp(coefs,coefs[1]+coefs[6]+coefs[8]))
          zt=imag(cf_zwt4cp(coefs,coefs[1]+coefs[6]+coefs[8]))
          tt=t-coefs[1]-coefs[6]-coefs[8]
          y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
          y+=coefs[9]/coefs[10]*((coefs[2]/coefs[3]-1)*(1-exp(-coefs[3]*tt))-(coefs[2]/coefs[4]-1)*(1-exp(-coefs[4]*tt)))
          y/=coefs[4]-coefs[3]
          z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
          z+=coefs[9]/coefs[10]*((1-exp(-coefs[3]*tt))/coefs[3]-(1-exp(-coefs[4]*tt))/coefs[4])
          z*=k12/(coefs[4]-coefs[3])
        else
          yt=real(cf_zwt4cp(coefs,coefs[1]+coefs[6]+coefs[8]+coefs[10]))
          zt=imag(cf_zwt4cp(coefs,coefs[1]+coefs[6]+coefs[8]+coefs[10]))
          tt=t-coefs[1]-coefs[6]-coefs[8]-coefs[10]
          y=yt*(coefs[4]*exp(-coefs[4]*tt)-coefs[3]*exp(-coefs[3]*tt))+coefs[2]*(yt+zt)*(exp(-coefs[3]*tt)-exp(-coefs[4]*tt))
          y/=coefs[4]-coefs[3]
          z=(yt+coefs[2]/(coefs[2]-coefs[3])*zt)*exp(-coefs[3]*tt)-(yt+coefs[2]/(coefs[2]-coefs[4])*zt)*exp(-coefs[4]*tt)
          z*=k12/(coefs[4]-coefs[3])
        endif
      endif
    endif
  endif
  return cmplx(y,z)
end function

function cf_zwt4_k(coefs,t)	//	finite, 4-stage, zero-order absorption, two-compartment elimination with ks 2022-09-06
  wave coefs
  variable t
  return real(cf_zwt4_kcp(coefs,t))
end

function/c cf_zwt4_kcp(coefs,t)	//	finite, 4-stage, zero-order absorption, two-compartment elimination with ks 2022-09-06
  wave coefs
  variable t
  variable y,z,yt,zt,tt
  variable k21=coefs[2],k12=coefs[3],k10=coefs[4],t1=coefs[1],t2=coefs[6],t3=coefs[8],t4=coefs[10]
  variable b=coefs[2]+coefs[3]+coefs[4],c=coefs[2]*coefs[4],d=sqrt(b*b-4*c),l1=0.5*(b-d),l2=0.5*(b+d)		// b,c,d nomenclature from trinomial
  if (t<=t1)
    y= coefs[0]/t1/(l2-l1)*((k21-l1)/l1*(1-exp(-l1*t)) - (k21-l2)/l2*(1-exp(-l2*t)))
    z= coefs[0]*k12/t1/(l2-l1)*(1/l1*(1-exp(-l1*t)) - 1/l2*(1-exp(-l2*t)))
  else
    if (t<=t1+t2)
      yt=real(cf_zwt4_kcp(coefs,t1))
      zt=imag(cf_zwt4_kcp(coefs,t1))
      tt=t-t1
      y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
      y+=coefs[5]/t2*((k21/l1-1)*(1-exp(-l1*tt))-(k21/l2-1)*(1-exp(-l2*tt)))
      y/=l2-l1
      z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
      z+=coefs[5]/t2*((1-exp(-l1*tt))/l1-(1-exp(-l2*tt))/l2)
      z*=k12/(l2-l1)
    else
      if (t<=t1+t2+t3)
        yt=real(cf_zwt4_kcp(coefs,t1+t2))
        zt=imag(cf_zwt4_kcp(coefs,t1+t2))
        tt=t-t1-t2
        y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
        y+=coefs[7]/t3*((k21/l1-1)*(1-exp(-l1*tt))-(k21/l2-1)*(1-exp(-l2*tt)))
        y/=l2-l1
        z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
        z+=coefs[7]/t3*((1-exp(-l1*tt))/l1-(1-exp(-l2*tt))/l2)
        z*=k12/(l2-l1)
      else
        if (t<=t1+t2+t3+t4)
          yt=real(cf_zwt4_kcp(coefs,t1+t2+t3))
          zt=imag(cf_zwt4_kcp(coefs,t1+t2+t3))
          tt=t-t1-t2-t3
          y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
          y+=coefs[9]/t4*((k21/l1-1)*(1-exp(-l1*tt))-(k21/l2-1)*(1-exp(-l2*tt)))
          y/=l2-l1
          z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
          z+=coefs[9]/t4*((1-exp(-l1*tt))/l1-(1-exp(-l2*tt))/l2)
          z*=k12/(l2-l1)
        else
          yt=real(cf_zwt4_kcp(coefs,t1+t2+t3+t4))
          zt=imag(cf_zwt4_kcp(coefs,t1+t2+t3+t4))
          tt=t-t1-t2-t3-t4
          y=yt*(l2*exp(-l2*tt)-l1*exp(-l1*tt))+k21*(yt+zt)*(exp(-l1*tt)-exp(-l2*tt))
          y/=l2-l1
          z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
          z*=k12/(l2-l1)
        endif
      endif
    endif
  endif
  return cmplx(y,z)
end function

function cf_fw(coefs,t)	//	1-stage, first-order absorption, two-compartment elimination 8/6/2021
  wave coefs
  variable t
  return real(cf_fwcp(coefs,t))
end function

function cf_fw_k(coefs,t)	//	1-stage, first-order absorption, two-compartment elimination with k12 and k10 22/8/2021
  wave coefs
  variable t
  return real(cf_fw_kcp(coefs,t))
end function

function/c cf_fwcp(coefs,t)	//	1-stage, first-order absorption, two-compartment elimination 10/6/2021 with calculation for peripheral compartment
  wave coefs
  variable t
  variable y,z,B1,B2,B3,l1=coefs[3],l2=coefs[4],ka=coefs[1],k21=coefs[2],k10=l1*l2/k21,k12=l1+l2-k21-k10
  B1=coefs[0]*ka*(k21-l1)/(l2-l1)/(ka-l1)
  B2=coefs[0]*ka*(k21-l2)/(ka-l2)/(l1-l2)
  B3=coefs[0]*ka*(k21-ka)/(l1-ka)/(l2-ka)
  y=B1*exp(-l1*t)+B2*exp(-l2*t)+B3*exp(-ka*t)
  z=(-(l1*B1*exp(-l1*t)+l2*B2*exp(-l2*t)+ka*B3*exp(-ka*t))-ka*coefs[0]*exp(-ka*t)+(k12+k10)*y)/k21
  return cmplx(y,z)
end function

function/c cf_fw_kcp(coefs,t)	//	1-stage, first-order absorption, two-compartment elimination with k12 and k10 22/8/2021
  wave coefs
  variable t
  variable y,z,B1,B2,B3,ka=coefs[1],k21=coefs[2],k12=coefs[3],k10=coefs[4]
  variable b=k12+k21+k10,c=k21*k10,d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
  B1=coefs[0]*ka*(k21-l1)/(l2-l1)/(ka-l1)
  B2=coefs[0]*ka*(k21-l2)/(ka-l2)/(l1-l2)
  B3=coefs[0]*ka*(k21-ka)/(l1-ka)/(l2-ka)
  y=B1*exp(-l1*t)+B2*exp(-l2*t)+B3*exp(-ka*t)
  z=(-(l1*B1*exp(-l1*t)+l2*B2*exp(-l2*t)+ka*B3*exp(-ka*t))-ka*coefs[0]*exp(-ka*t)+(k12+k10)*y)/k21
  return cmplx(y,z)
end function

function cf_fwt(coefs,t)	//	finite, 1-stage, first-order absorption, two-compartment elimination 21/7/2021
  wave coefs
  variable t
  return real(cf_fwtcp(coefs,t))
end function

function/c cf_fwtcp(coefs,t)	//	finite, 1-stage, first-order absorption, two-compartment elimination 21/7/2021 with calculation for peripheral compartment
  wave coefs
  variable t
  variable y,yt,zt,z,tt,B1,B2,B3,l1=coefs[3],l2=coefs[4],ka=coefs[1],k21=coefs[2],k10=l1*l2/k21,k12=l1+l2-k21-k10,A1,A2
  if (t<=coefs[5])
    B1=coefs[0]*ka*(k21-l1)/(l2-l1)/(ka-l1)
    B2=coefs[0]*ka*(k21-l2)/(ka-l2)/(l1-l2)
    B3=coefs[0]*ka*(k21-ka)/(l1-ka)/(l2-ka)
    y=B1*exp(-l1*t)+B2*exp(-l2*t)+B3*exp(-ka*t)
    z=(-(l1*B1*exp(-l1*t)+l2*B2*exp(-l2*t)+ka*B3*exp(-ka*t))-ka*coefs[0]*exp(-ka*t)+(k12+k10)*y)/k21
  else
    yt=real(cf_fwtcp(coefs,coefs[5]))
    zt=imag(cf_fwtcp(coefs,coefs[5]))
    tt=t-coefs[5]
    A1=((k21-l1)*yt+k21*zt)/(l2-l1)
    A2=((k21-l2)*yt+k21*zt)/(l1-l2)
    y=A1*exp(-l1*tt)+A2*exp(-l2*tt)
    z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
    z*=k12/(l2-l1)
  endif
  return cmplx(y,z)
end function

function cf_fwt_k(coefs,t)	//	finite, 1-stage, first-order absorption, two-compartment elimination with k12 and k10 (not l1,l2) 16/8/2021
  wave coefs
  variable t
  return real(cf_fwt_kcp(coefs,t))
end function

function/c cf_fwt_kcp(coefs,t)	//	finite, 1-stage, first-order absorption, two-compartment elimination with k12 and k10 (not l1, l2) 16/8/2021 with calculation for peripheral compartment
  wave coefs
  variable t
  variable y,z,yt,zt,tt,B1,B2,B3,A1,A2,ka=coefs[1],k21=coefs[2],k12=coefs[3],k10=coefs[4]
  variable b=k12+k21+k10,c=k21*k10,d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
  if (t<=coefs[5])
    B1=coefs[0]*ka*(k21-l1)/(l2-l1)/(ka-l1)
    B2=coefs[0]*ka*(k21-l2)/(ka-l2)/(l1-l2)
    B3=coefs[0]*ka*(k21-ka)/(l1-ka)/(l2-ka)
    y=B1*exp(-l1*t)+B2*exp(-l2*t)+B3*exp(-ka*t)
    z=(-(l1*B1*exp(-l1*t)+l2*B2*exp(-l2*t)+ka*B3*exp(-ka*t))-ka*coefs[0]*exp(-ka*t)+(k12+k10)*y)/k21
  else
    yt=real(cf_fwt_kcp(coefs,coefs[5]))
    zt=imag(cf_fwt_kcp(coefs,coefs[5]))
    tt=t-coefs[5]
    A1=((k21-l1)*yt+k21*zt)/(l2-l1)
    A2=((k21-l2)*yt+k21*zt)/(l1-l2)
    y=A1*exp(-l1*tt)+A2*exp(-l2*tt)
    z=(yt+k21/(k21-l1)*zt)*exp(-l1*tt)-(yt+k21/(k21-l2)*zt)*exp(-l2*tt)
    z*=k12/(l2-l1)
  endif
  return cmplx(y,z)
end function

function cf_zzt2(coefs,t)	//	zero order absorption and elimination for duration tau1 and tau2
  wave coefs
  variable t
  variable y
  if (t<=coefs[1])
    y=coefs[0]*t/coefs[1]
  else
    if (t<=coefs[1]+coefs[2])
      y=coefs[0]*(1-(t-coefs[1])/coefs[2])
    else
      y=0
    endif
  endif
  return y
end function

// +++++++++++++++ fitting tools +++++++++++++++++++++++

Function chisq(x,y,i1,i2)		//  2021-12-29
  wave x,y
  variable i1,i2
  return real(chisq_corrcoef(x,y,i1,i2))
end

Function corrcoef(x,y,i1,i2)		// correlation coefficient for fit (improvised) 2021-04-09
  wave x,y
  variable i1,i2
  return imag(chisq_corrcoef(x,y,i1,i2))
end

Function/c chisq_corrcoef(x,y,i1,i2)		// correlation coefficient for fit (improvised) 2021-04-09, turned complex 2021-12-29 to include ss
  wave x,y
  variable i1,i2
  variable i=i1,ss=0,sxy=0,sx2=0,sy2=0,xm=mean(x,i1,i2),ym=mean(y,i1,i2)
  do
    ss+=(x[i]-y[i])^2
    sxy+=(x[i]-xm)*(y[i]-ym)
    sx2+=(x[i]-xm)^2
    sy2+=(y[i]-ym)^2
    i+=1
  while (i<i2)
  return cmplx(ss,sxy/sqrt(sx2*sy2))
End

// ================== Fitting macros ============================

Proc fitp_ln()			//2021-06-12
  string lbl=data_label,fitlabel		//	2021-08-17
  DoWindow/F Gcf_ln
  TextBox/C/N=text0 lbl
  duplicate/o c_data lnc_data
  lnc_data=ln(c_data)
  string t_unit=StringFromList(1,note(c_data))
  label bottom "\f02t\f00 ("+t_unit+")"
  variable i=0,c_max,t_max,i_max
  do
    if (c_data[i]>c_max)
      c_max=c_data[i]
      t_max=t_data[i]
      i_max=i
    endif
    i+=1
  while (i<numpnts(c_data))
  CurveFit/Q/M=2/NTHR=0 line  lnc_data[i_max,] /X=t_data /D /I=1	// fit lnc = f(t) without weights
  SetScale/P x 0,t_data[numpnts(t_data)-1],"" lnc_line
  lnc_line[0]=poly(W_coef,0)
  lnc_line[1]=poly(W_coef,t_data[numpnts(t_data)-1])
  duplicate/o c_data cs_lndata
  cs_lndata=poly(W_coef,t_data)
  cc=corrcoef(cs_lndata,lnc_data,i_max,numpnts(c_data)-1)
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_lnc_data fitlabel
  printf  "%s ln fit:	tmax = %1.2f "+t_unit+"	kel = %.3f (±%.3f) "+t_unit+"^-1	R^2 = %1.5f\r",data_label,t_max,-W_coef[1],W_sigma[1],cc
//  print lbl,":", i_max,", C_max =",c_max," @ t_max =",t_max,"h"
End

Proc fitp_dc()		// 2021-06-21
  string lbl=data_label,fitlabel,printlbl=lbl,templbl
  variable i=0,i_max,t_max,c_max,i_start_dc
  do
    if (c_data[i]>c_max)
      c_max=c_data[i]
      t_max=t_data[i]
      i_max=i
    endif
    i+=1
  while (i<numpnts(c_data))
  string c_unit=StringFromList(0,note(c_data))
  string t_unit=StringFromList(1,note(c_data))
  DoWindow/F Gcf_dC
  TextBox/C/N=text0 lbl
  label bottom "\f02C\f00 (\E"+c_unit+")"
  label left "d\f02C\f00/d\f02t\f00 (\E"+c_unit+" "+t_unit+")"
  duplicate/o c_data, dc_data
  Differentiate dc_data/X=t_data
  if (strlen(CsrInfo(B)) > 0)
    i_start_dc=pcsr(B)
  else
    i_start_dc=i_max+2
  endif
  CurveFit/Q/M=2/W=0 line, dc_data[i_start_dc,]/X=c_data/D	// fit dc/dt = f(c)
  duplicate/o dc_data dcs_data
  wavestats/q dc_data
  setaxis left V_min*1.1,0
  dcs_data=poly(W_coef,c_data)
  cc=corrcoef(dcs_data,dc_data,i_start_dc,numpnts(c_data)-1)
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
  fitlabel+="\rd\\f02C\\f00/d\f02t\\f00 = "+num2str(W_coef[0])+" ± "+num2str(W_sigma[0])+" "+c_unit+" "+t_unit
  fitlabel+="\r\\f02?\\f00\\S2\\M = "+num2str(V_chisq)+",  \\f02R\\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_dc_data fitlabel
//  print W_coef[0]," ±", W_sigma[0],W_coef[1]," ±", W_sigma[1], lbl
//  sprintf templbl " kel = %.3f (±%.3f) "+t_unit+"^-1 %s",-W_coef[1],W_sigma[1],lbl
End

Proc fitp(mdl,flag)		// 2021-08-17 main fitting procedure
	string mdl			// fo,fot,fo0,fo0t,zot1,zot2,zot3,zwt1_k,bo ,bw, bw_k, etc.
	variable flag			// 0, 1, 2
	string graph="Gcf_"+mdl,table="Tcf_"+mdl,fit="fit_"+mdl,res="res_"+mdl,H="H_"+mdl,cf="cf_"+mdl
	string par="par_"+mdl,yt="yt_"+mdl,parl="parl_"+mdl,parh="parh_"+mdl
	variable tau
	string cflag			// 2022-01-06 for showing in output the use of fitting constraints
	string c_unit=StringFromList(0,note(c_data))
	string t_unit=StringFromList(1,note(c_data))
	DoWindow/F $graph
	label left "\f02C\f00 (\E"+c_unit+")"
	label bottom "\f02t\f00 ("+t_unit+")"
	TextBox/C/N=text0 data_label
	variable i=0
	holdflags=""
	do
		holdflags+=num2str($H[i])
		i+=1
	while (i<numpnts($H)-2)
	if (flag==1)			// just simulation, no optimization
		holdflags=ReplaceString("0",holdflags,"1")
		cflag="s"
	endif
	variable j=numpnts($H)-1, startp=$H[j-1], endp=($H[j]==0)*(numpnts(c_data)-1)+($H[j]>0)*$H[j]
	if ((flag==2) && (exists(parh)==1))			// fitting with constraints, if high parameter limits are defined
		make/o/t/n=(2*(numpnts($par)-sum($H))) T_par_constr	// 2021-12-30
		i=0
		j=0
		do
			if ($H[j]==0)						// 2022-01-05 to constrain parameters than are not held constant
				T_par_constr[i]="K"+num2str(j)+" > 0"
				T_par_constr[i+numpnts($par)-sum($H)]="K"+num2str(j)+" < "+num2str($parh[j])
				i+=1
			endif
			j+=1
		while (j<numpnts($par))
		cflag="c"
	else
		flag=0
	endif
	RemoveFromGraph/Z $res
	killwaves/Z $res
	Make/o/n=500 $fit
	SetScale/I x t_data[0],t_data[numpnts(t_data)-1],"", $fit
	if (flag==2)			// fitting with constraints
		FuncFit/Q/H=holdflags/NTHR=0/M=2/TBOX=768 $cf $par c_data[startp,endp] /X=t_data /I=1 /R /C=T_par_constr
	else
		FuncFit/Q/H=holdflags/NTHR=0/M=2/TBOX=768 $cf $par c_data[startp,endp] /X=t_data /I=1 /R
	endif
	AddCorrelation()
	$fit=$cf($par,x)
	duplicate/o $fit auc_fit
	auc_fit=area($fit,0,x)
	if (ItemsInList(wavelist(fit,";","WIN:"))==0)
		AppendToGraph $fit
	endif
	rename res_c_data,$res
	modifyGraph mode(c_data)=3,mode($res)=4,msize=0
	duplicate/o c_data,sim_data
	sim_data=$cf($par,t_data)
	cc=corrcoef(sim_data,c_data,startp,endp)
	string fitlabel			// 2022-06-13 added units in tau, k and FD/Vd or FD/Vc; labels still need trimming for publication purposes
	if (cmpstr(mdl,"fo0")==0)
		fitlabel ="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00 = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" h\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"fo0t")==0)
		fitlabel ="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00 = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02o\f00 = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[2]
	endif
	if (cmpstr(mdl,"fo")==0)
		fitlabel="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\\Ba\\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"fot")==0)
		fitlabel="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Ba\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02o\f00 = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[3]
	endif
	if (cmpstr(mdl,"zot1")==0)
		fitlabel="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Bel\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02o\f00 = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[2]
	endif
	if (cmpstr(mdl,"zot2")==0)
		fitlabel="\f02F\f00\B1\MD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\MD/V\f00\Bd\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+c_unit
		fitlabel+="\r\f02o\f00\B2\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit
		fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[3]
	endif
	if (cmpstr(mdl,"zot3")==0)
		fitlabel="\f02F\f00\B1\MD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+=", \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\MD/V\f00\Bd\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+c_unit
		fitlabel+=", \f02o\f00\B2\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit
		fitlabel+="\r\f02F\f00\B3\MD/V\f00\Bd\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+c_unit
		fitlabel+=", \f02o\f00\B3\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+t_unit
		fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[3]+$par[5]
	endif
	if (cmpstr(mdl,"zot4")==0)
		fitlabel="\f02F\f00\B1\MD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+=", \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\MD/V\f00\Bd\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+c_unit
		fitlabel+=", \f02o\f00\B2\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit
		fitlabel+="\r\f02F\f00\B3\MD/V\f00\Bd\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+c_unit
		fitlabel+=", \f02o\f00\B3\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+t_unit
		fitlabel+="\r\f02F\f00\B4\MD/V\f00\Bd\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+c_unit
		fitlabel+=", \f02o\f00\B4\M = "+num2str($par[7])+" ± "+num2str(W_sigma[7])+" "+t_unit
		fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str($par[8])+" ± "+num2str(W_sigma[8])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[3]+$par[5]+$par[7]
	endif
	if (cmpstr(mdl,"zwt1")==0)
		fitlabel="\f02FD/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02o\f00 = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]
	endif
	if (cmpstr(mdl,"zwt1_k")==0)
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02FD/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02o\f00 = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]
	endif
	if (cmpstr(mdl,"zwt2")==0)
		fitlabel="\f02F\f00\B1\M\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="  \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\M\f02D/V\f00\Bc\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+c_unit
		fitlabel+="  \f02o\f00\B2\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[6]
	endif
	if (cmpstr(mdl,"zwt2_k")==0)
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02F\f00\B1\M\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+=",  \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\M\f02D/V\f00\Bc\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+c_unit
		fitlabel+=",  \f02o\f00\B2\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[6]
	endif
	if (cmpstr(mdl,"zwt3")==0)
		fitlabel="\f02F\f00\B1\M\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="  \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\M\f02D/V\f00\Bc\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+c_unit
		fitlabel+="  \f02o\f00\B2\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit
		fitlabel+="\r\f02F\f00\B3\M\f02D/V\f00\Bc\M = "+num2str($par[7])+" ± "+num2str(W_sigma[7])+" "+c_unit
		fitlabel+="  \f02o\f00\B3\M = "+num2str($par[8])+" ± "+num2str(W_sigma[8])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[6]+$par[8]
	endif
	if (cmpstr(mdl,"zwt3_k")==0)
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02F\f00\B1\M\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+=",  \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\M\f02D/V\f00\Bc\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+c_unit
		fitlabel+=",  \f02o\f00\B2\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit
		fitlabel+="\r\f02F\f00\B3\M\f02D/V\f00\Bc\M = "+num2str($par[7])+" ± "+num2str(W_sigma[7])+" "+c_unit
		fitlabel+="  \f02o\f00\B3\M = "+num2str($par[8])+" ± "+num2str(W_sigma[8])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[6]+$par[8]
	endif
	if (cmpstr(mdl,"zwt4")==0)
		fitlabel="\f02F\f00\B1\M\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="  \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\M\f02D/V\f00\Bc\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+c_unit
		fitlabel+="  \f02o\f00\B2\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit
		fitlabel+="\r\f02F\f00\B3\M\f02D/V\f00\Bc\M = "+num2str($par[7])+" ± "+num2str(W_sigma[7])+" "+c_unit
		fitlabel+="  \f02o\f00\B3\M = "+num2str($par[8])+" ± "+num2str(W_sigma[8])+" "+t_unit
		fitlabel+="\r\f02F\f00\B4\M\f02D/V\f00\Bc\M = "+num2str($par[9])+" ± "+num2str(W_sigma[9])+" "+c_unit
		fitlabel+="  \f02o\f00\B4\M = "+num2str($par[10])+" ± "+num2str(W_sigma[10])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[6]+$par[8]+$par[10]
	endif
	if (cmpstr(mdl,"zwt4_k")==0)
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02F\f00\B1\M\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+=",  \f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02F\f00\B2\M\f02D/V\f00\Bc\M = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+c_unit
		fitlabel+=",  \f02o\f00\B2\M = "+num2str($par[6])+" ± "+num2str(W_sigma[6])+" "+t_unit
		fitlabel+="\r\f02F\f00\B3\M\f02D/V\f00\Bc\M = "+num2str($par[7])+" ± "+num2str(W_sigma[7])+" "+c_unit
		fitlabel+="  \f02o\f00\B3\M = "+num2str($par[8])+" ± "+num2str(W_sigma[8])+" "+t_unit
		fitlabel+="\r\f02F\f00\B4\M\f02D/V\f00\Bc\M = "+num2str($par[9])+" ± "+num2str(W_sigma[9])+" "+c_unit
		fitlabel+="  \f02o\f00\B4\M = "+num2str($par[10])+" ± "+num2str(W_sigma[10])+" "+t_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[1]+$par[6]+$par[8]+$par[10]
	endif
	if (cmpstr(mdl,"bo")==0)
		fitlabel="\f02D/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Bel\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"bwe")==0)
		fitlabel  =  "\f02A\f00 = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02B\f00 = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+c_unit
		fitlabel+="\r\f02a\f00 = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00 = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		print "AUC =",$par[0]/$par[2]+$par[1]/$par[3]
	endif
	if (cmpstr(mdl,"bwet")==0)
		fitlabel="\f02?\f00 = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02I\f00 = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+c_unit
		fitlabel+="\r\f02a\f00 = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00 = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02o\f00 = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"bw")==0)
		fitlabel="\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit+", \f02A\f00 = "+num2str($par[0]*($par[1]-$par[2])/($par[3]-$par[2]))+" "+c_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"+", \f02B\f00 = "+num2str($par[0]*($par[1]-$par[3])/($par[2]-$par[3]))+" "+c_unit
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[2]+$par[3]-$par[1]-$par[2]*$par[3]/$par[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[2]*$par[3]/$par[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"bw_k")==0)
		variable b=$par[2]+$par[1]+$par[3],c=$par[1]*$par[3],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02D/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2]) +" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3]) +" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"fw")==0)
		fitlabel="\f02FD/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Ba\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"fw_k")==0)
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02FD/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Ba\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M" + ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M" + ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
	endif
	if (cmpstr(mdl,"fwt")==0)
		fitlabel="\f02FD/V\f00\Bc\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Ba\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3])+" "+t_unit+"\S-1\M"+", \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02a\f00\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M"+", \f02k\\f00\\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02o\f00 = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[5]
	endif
	if (cmpstr(mdl,"fwt_k")==0)
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		fitlabel="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02k\f00\Ba\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B21\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\f00\B12\\M = "+num2str($par[3])+" ± "+num2str(W_sigma[3]) +" "+t_unit+"\S-1\M"+ ", \\f02a\\f00 = "+num2str(l1)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02k\\f00\\B10\\M = "+num2str($par[4])+" ± "+num2str(W_sigma[4])+" "+t_unit+"\S-1\M" + ", \\f02a\\f00 = "+num2str(l2)+" "+t_unit+"\S-1\M"
		fitlabel+="\r\f02o\f00 = "+num2str($par[5])+" ± "+num2str(W_sigma[5])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[5]
	endif
	if (cmpstr(mdl,"zzt2")==0)
		fitlabel="\f02FD/V\f00\Bd\M = "+num2str($par[0])+" ± "+num2str(W_sigma[0])+" "+c_unit
		fitlabel+="\r\f02o\f00\B1\M = "+num2str($par[1])+" ± "+num2str(W_sigma[1])+" "+t_unit
		fitlabel+="\r\f02o\f00\B2\M = "+num2str($par[2])+" ± "+num2str(W_sigma[2])+" "+t_unit
		fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
		tau=$par[5]
	endif
	if (flag==2)
		fitlabel=replacestring("\f02?\f00",fitlabel,cflag+"\f02?\f00")
	endif
	if (strlen(fitlabel)>0)
		TextBox/C/N=CF_c_data fitlabel
	endif
	if(tau>0)
		make/o/n=1 $yt=$cf($par,tau)
		SetScale/P x tau,1,"", $yt
		if (ItemsInList(wavelist(yt,";","WIN:"))==0)
			AppendToGraph $yt
			ModifyGraph mode($yt)=3,marker($yt)=17,rgb($yt)=(0,0,0),msize=4
		endif
	endif
	print $par
	string templbl, printlabel=data_label+"\t"+mdl+"\t"+cflag
	sprintf templbl "\t%.6g\t%.5g",V_chisq,cc
	printlabel+=templbl
	i=0
	do
		sprintf templbl "\t%.5g \t± %.3g",$par[i],W_sigma[i]
		printlabel+=templbl
		i+=1
	while (i<numpnts($par))
	sprintf templbl "\tAUC = %.5g",auc_fit[numpnts(auc_fit)-1]
	printlabel+=templbl
	print printlabel
	DoWindow/F $table
//	DoWindow/H/F
End

// ======================== Graph Button functions =======================

Function ButtonRunFit(ctrlName) : ButtonControl			// for "Fit" and "Simu" buttons and (2022-01-04) "cFit" button in Gcf_xxx graphs
  String ctrlName
  string mdl=ctrlName[4,strlen(ctrlName)-1],comm="fitp(\""+mdl+"\""
  if ((cmpstr(ctrlName,"fit_ln")==0) || (cmpstr(ctrlName,"fit_dc")==0))
    comm=ctrlName+"()"
    comm=ReplaceString("t_",comm,"tp_")
  else
    if (cmpstr(ctrlName[0,2],"Fit")==0)
      comm+=",0)"
    else
      if (cmpstr(ctrlName[0,2],"Cit")==0)
        comm+=",2)"
      else
        comm+=",1)"
      endif
    endif
  endif
  execute comm
End

Function ButtonShowParams(ctrlName) : ButtonControl		// for "Params" buttons
  String ctrlName
  string comm=ReplaceString("Show_coefs_",ctrlName,"Tcf_")
  DoWindow/F $comm
End

Function Update_data_ButtonProc(ctrlName) : ButtonControl		// for the "Run" button
  String ctrlName
  string comm="Update_data()"
  execute comm
End

Macro Update_data(lbl,c_unit,t_unit)			//2021-06-12 for the "Run" button in G_data, 2022-08-08 added t_units
  string lbl,c_unit,t_unit
  prompt lbl "Data label", popup labellist
  prompt c_unit "Concentration units", popup unitslist
  prompt t_unit "Time units", popup tunitslist
  
  data_label=lbl
  labellist=RemoveListItem(0,labellist)		// 2022-01-04
  labellist=AddListItem(lbl,labellist)
  unitslist=RemoveListItem(0,unitslist)
  unitslist=AddListItem(c_unit,unitslist)
  tunitslist=RemoveListItem(0,tunitslist)
  tunitslist=AddListItem(t_unit,tunitslist)
  DoWindow/F G_data
  TextBox/C/N=text0 lbl
  label left "\f02C\f00 (\E"+c_unit+")"
  label bottom "\f02t\f00 (\E"+t_unit+")"
  note/K c_data c_unit+";"+ t_unit				// 2021-08-17
  duplicate/o c_data, dc_data
  Differentiate dc_data/X=t_data
  duplicate/o c_data lnc_data
  lnc_data=ln(c_data)
  variable i=0,c_max,t_max,i_max,AUC_tot,AUC_tmax
  do
    if (c_data[i]>c_max)
      c_max=c_data[i]
      t_max=t_data[i]
      i_max=i
    endif
    i+=1
  while (i<numpnts(c_data))
  AUC_tot=areaXY(t_data,c_data)			// 2021-12-30
  AUC_tmax=areaXY(t_data,c_data,0,t_data[i_max])		// corrected 2022-04-07
  duplicate/o t_data auc_data				// 2022-04-02
  auc_data=areaXY(t_data,c_data,0,t_data[x])			// corrected 2022-04-07
  lbl+=": i_max = "+num2str(i_max)+" C_max = "+num2str(c_max)+" "+c_unit+" @ t_max = "+num2str(t_max)+" "+t_unit
  lbl+=",  AUC[0-t_max] = "+num2str(AUC_tmax)+", AUC = "+num2str(AUC_tot)+", AUC ratio = "+num2str(AUC_tmax/AUC_tot)
  print lbl
End

//  ------------------------------------------------=================== Simulation Panel ================== -----------------------------

Window P_sim() : Panel		// 2021-07-17,18,21,22; help added 2022-05-24
	DoWindow/F T_models
	PauseUpdate; Silent 1		// building window...
	NewPanel /W=(1550,76,1820,312)
	ModifyPanel cbRGB=(49152,65280,32768)
	SetVariable duration,pos={1,2},size={130,24},title="Duration"
	SetVariable duration,labelBack=(0,65280,65280),fSize=16,format="%g h"
	SetVariable duration,limits={1,inf,0},value= duration
	SetVariable duration help={"Define simulation duration"}
	PopupMenu model,pos={136,2},size={130,24},bodyWidth=80,proc=model_select,title="Model"
	string pval="\""
	variable i= 0
	do
		pval+=models[i]+";"
		i+=1
	while (i<numpnts(models))
	pval+="\""
	PopupMenu model,fSize=16
	PopupMenu model,mode=1,popvalue=model_name,value= #pval
	PopupMenu model help={"Select model for simulations and fits. For model descriptions see 'Show models'"}
	Button Params,pos={1,32},size={130,24},proc=ShowWindowProc,title="Show parameters"
	Button Params,fSize=16,fColor=(52224,0,20736)
	Button Params help={"Bring to front parameter table for selected model"}
	Button Show_models,pos={136,32},size={130,24},proc=ShowModelsProc,title="Show models"
	Button Show_models,fSize=16
	Button Show_models help={"Bring to front table with models and descriptions"}
	Button Show_sim,pos={1,62},size={130,24},proc=ShowWindowProc,title="Show simulation"
	Button Show_sim,fSize=16,fColor=(65280,65280,0)
	Button Show_sim help={"Bring to front simulation window for selected model"}
	Button Simulate,pos={136,62},size={70,24},proc=Simulate_proc,title="Simulate"
	Button Simulate,labelBack=(65535,65535,65535),fSize=16,fColor=(65280,43520,0)
	Button Simulate help={"Run simulation for selected model"}
	PopupMenu annot_simul,pos={206,62},size={60,24},bodyWidth=60,proc=annotate_simulation_select,title=""
	PopupMenu annot_simul,fSize=16
	PopupMenu annot_simul,mode=1,popvalue=annotate_simulation,value= #"\"yes;fold;no\""
	PopupMenu annot_simul help={"Select position of simulation parameters in simulation plot"}
	Button ShowDataTable,pos={1,92},size={130,24},proc=ShowDataTableProc,title="Data Table"
	Button ShowDataTable,fSize=16,fColor=(0,12800,52224)
	Button ShowDataTable help={"Bring to front table to paste or edit data (and error bars)"}
	Button AddNoise,pos={136,92},size={80,24},proc=AddNoiseProc,title="Add noise"
	Button AddNoise,fSize=16,fColor=(65280,43520,0)
	Button AddNoise help={"Add proportional noise to simulated data and discretize them"}
	SetVariable nslevel,pos={226,92},size={40,24}
	SetVariable nslevel,labelBack=(0,65280,65280),fSize=16,title=" "
	SetVariable nslevel,limits={0,inf,0},value= nslevel
	SetVariable nslevel help={"Proportional noise level for simulation"}
	Button ShowData,pos={1,122},size={130,24},proc=ShowDataProc,title="Show data graph"
	Button ShowData,fSize=16,fColor=(0,39168,0)
	Button ShowData help={"Bring to front window with data and simple operations buttons"}
	Button Show_fit,pos={136,122},size={98,24},proc=ShowWindowProc,title="Show fitting"
	Button Show_fit,fSize=16
	Button Show_fit help={"Bring to front window for fitting selected model"}
	Button Map,pos={226,122},size={40,24},proc=map_chisq,title="Map"
	Button Map,fSize=16,fColor=(33056,0,0)
	Button Map help={"Do multiple calculations of chi squared for ranges of parameters defined in parameter table"}
	Button CopySimulated,pos={1,152},size={130,24},proc=CopySimProc,title="Copy simulated"
	Button CopySimulated,fSize=16
	Button CopySimulated help={"Transfer simulated data to Data table"}
	Button Many_fits,pos={136,152},size={80,24},proc=NSimFit_proc,title="Many fits"
	Button Many_fits,fSize=16,fColor=(33056,33056,0)
	Button Many_fits help={"Perform multiple fits on simulated data"}
	SetVariable nfits,pos={226,152},size={40,24}
	SetVariable nfits,fsize=16,title=" "
	SetVariable nfits,limits={1,inf,0},value= n_fits
	Button Many_fits help={"Select number for multiple fits"}
	Button ReDrawPanel,pos={1,182},size={130,24},proc=ReDrawPanel,title="Refrech Panel"
	Button ReDrawPanel,fSize=16,fColor=(13056,0,0)
	Button ReDrawPanel help={"Update panel for possible new models"}
	Button AddDataLabels,pos={136,182},size={130,24},proc=AddDataLabelsProc,title="Add data labels"
	Button AddDataLabels,fSize=16
	Button AddDataLabels help={"Expand listing of names of data sets"}
	Button HideWindows,pos={1,212},size={130,24},proc=HideWindowsProc,title="Hide windows"
	Button HideWindows,fSize=16,fColor=(13056,0,0)
	Button HideWindows help={"Keep only essential graphs and tables"}
	Button Show_correlations,pos={136,212},size={130,24},proc=ShowCovarProc,title="Show correlations"
	Button Show_correlations,fSize=16,fColor=(65056,0,65000)
	Button Show_correlations help={"Bring to front table with covariance and correlation matrix"}
EndMacro

Proc Simulate_proc(ctrlName) : ButtonControl		// 2021-07-17
	String ctrlName
	string data = "sim_"+model_name,par="par_"+model_name,parlabels="plb_"+model_name,w_name="Gsim_"+model_name, model="cf_"+model_name
	string dataP = "simP_"+model_name,yt="yt_"+model_name, modelCP="cf_"+model_name+"cp"
	make/o/n=1000 $data
	SetScale x 0,duration,"", $data
	$data=$model($par,x)
	wavestats/q $data
	if (cmpstr(annotate_simulation,"yes")==0)
		Display/W=(10,10,450,160) $data
		modify axisEnab(bottom)={0,0.6}
	else
		Display/W=(10,10,240,160) $data
	endif
	modify fsize=12, lblMargin(left)=5, lsize=2, standoff=0
	label left "\f02C"	// \f00 (\Eng/mL)"
	label bottom "\f02t\f00 (hour)"
	if (exists(modelCP)==6)
		make/o/n=1000 $dataP
		SetScale x 0,duration,"", $dataP
		$dataP=imag($modelCP($par,x))
		append $dataP
		modify lsize=2,rgb($dataP)=(0,0,0),  lStyle($dataP)=1
	endif
	DoWindow/K $w_name
	DoWindow/C $w_name
	string tagtxt
	variable i = 0
	variable tau
	do
		tagtxt+="\f02"+$parlabels[i]+"\f00 = "+num2str($par[i])+"\r"	// generic building of tagtxt
		i+=1
	while (i<numpnts($par))
	tagtxt+="\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
	tagtxt=replacestring("tau",tagtxt,"o")
	if (cmpstr(model_name,"bo")==0)		// 2021-08-18
		tagtxt =  "\\f02FD/V\\f00\\Bd\\M = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00\\Bel\\M = "+num2str($par[1])+" h\\S-1\\M"
	endif
	if (cmpstr(model_name,"bw")==0)		// 2021-08-18
		tagtxt = "\\f02A\\f00 = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[2])+" h\\S-1\\M, \\f02k\\f00\\B12\\M = "+num2str($par[2]+$par[3]-$par[1]-$par[2]*$par[3]/$par[1])+" h\S-1\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M, \\f02k\\f00\\B10\\M = "+num2str($par[2]*$par[3]/$par[1])+" h\S-1\M"
	endif
	if (cmpstr(model_name,"bw_k")==0)		// 2021-08-18
		variable b=$par[2]+$par[1]+$par[3],c=$par[1]*$par[3],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		tagtxt = "\\f02A\\f00 = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00\\B12\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\S-1\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\S-1\M"
	endif
	if (cmpstr(model_name,"fo0")==0)		// 2021-08-18
		tagtxt =  "\\f02FD/V\\f00\\Bd\\M = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00 = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02t\\f00\\Bmax\\M = "+num2str(V_maxloc)+" h"
	endif
	if (cmpstr(model_name,"fo0t")==0)		// 2021-08-18
		tagtxt =  "\\f02FD/V\\f00\\Bd\\M = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00 = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02o\\f00 = "+num2str($par[2])+" h"
		tagtxt+="\r\\f02t\\f00\\Bmax\\M = "+num2str(V_maxloc)+" h"
		tau=$par[2]
	endif
	if (cmpstr(model_name,"fo")==0)		// 2021-08-18
		tagtxt =  "\\f02FD/V\\f00\\Bd\\M = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00\\Ba\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\Bel\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02t\\f00\\Bmax\\M = "+num2str(V_maxloc)+" h"
	endif
	if (cmpstr(model_name,"fot")==0)		// 2021-08-18
		tagtxt =  "\\f02FD/V\\f00\\Bd\\M = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00\\Ba\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\Bel\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02o\\f00 = "+num2str($par[3])+" h"
		tagtxt+="\r\\f02t\\f00\\Bmax\\M = "+num2str(V_maxloc)+" h"
		tau=$par[3]
	endif
	if (cmpstr(model_name,"zot1")==0)		// 2021-08-18
		tagtxt =  "\\f02FD/V\\f00\\Bd\\M = "+num2str($par[0])
		tagtxt+="\r\\f02o\\f00 = "+num2str($par[2])+" h"
		tagtxt+="\r\\f02k\\f00\\Bel\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[2]
	endif
	if (cmpstr(model_name,"zot2")==0)		// 2021-08-16
		tagtxt =  "\\f02F\\f00\\B1\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[0])+", \\f02o\\B1\\M\\f00 = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\f00\\B2\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[2])+", \\f02o\\B2\\M\\f00 = "+num2str($par[3])+" h"
		tagtxt+= "\r\\f02k\\Bel\\M\\f00 = "+num2str($par[4])+" h\\S-1\\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[3]
	endif
	if (cmpstr(model_name,"zot3")==0)		// 2021-08-18
		tagtxt =  "\\f02F\\f00\\B1\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[0])+", \\f02o\\B1\\M\\f00 = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\f00\\B2\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[2])+", \\f02o\\B2\\M\\f00 = "+num2str($par[3])+" h"
		tagtxt+="\r\\f02F\\f00\\B3\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[4])+", \\f02o\\B3\\M\\f00 = "+num2str($par[5])+" h"
		tagtxt+= "\r\\f02k\\Bel\\M\\f00 = "+num2str($par[6])+" h\\S-1\\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[3]+$par[5]
	endif
	if (cmpstr(model_name,"zot4")==0)		// 2022-09-06
		tagtxt =  "\\f02F\\f00\\B1\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[0])+", \\f02o\\B1\\M\\f00 = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\f00\\B2\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[2])+", \\f02o\\B2\\M\\f00 = "+num2str($par[3])+" h"
		tagtxt+="\r\\f02F\\f00\\B3\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[4])+", \\f02o\\B3\\M\\f00 = "+num2str($par[5])+" h"
		tagtxt+="\r\\f02F\\f00\\B4\\M\\f02D/V\\f00\\Bd\\M = "+num2str($par[6])+", \\f02o\\B4\\M\\f00 = "+num2str($par[7])+" h"
		tagtxt+= "\r\\f02k\\Bel\\M\\f00 = "+num2str($par[8])+" h\\S-1\\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[3]+$par[5]+$par[7]
	endif
	if (cmpstr(model_name,"zwt1")==0)		// 2021-08-18, 2021-11-17
		tagtxt = "\\f02FD/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00 = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M, \\f02k\\f00\\B12\\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[4])+" h\\S-1\\M, \\f02k\\f00\\B10\\M = "+num2str($par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]
	endif
	if (cmpstr(model_name,"zwt1_k")==0)		// 2021-08-16,18
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		variable fa,fb,f,auc=$par[0]/$par[4]
		fa=l1^2*($par[2]-l2)*(exp(l2*$par[1])-1)
		fb=l2^2*($par[2]-l1)*(exp(l1*$par[1])-1)
		f=l1*l2*(l1-l2)*$par[2]*$par[1]/(fa-fb)
		tagtxt = "\\f02FD/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00 = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M, AUC = "+num2str(auc)
		tagtxt+="\r\\f02k\\f00\\B12\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\S-1\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[4])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h, F = "+num2str(f)
		tau=$par[1]
	endif
	if (cmpstr(model_name,"zwt2")==0) 		// 2021-12-06
		tagtxt = "\\f02F\\B1\\MD/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00\\B1\\M = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\B2\\MD/V\\f00\\Bc\\M = "+num2str($par[5])+",  \\f02o\\f00\\B2\\M = "+num2str($par[6])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M\t\\f02k\\f00\\B12\\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[4])+" h\\S-1\\M\t\\f02k\\f00\\B10\\M = "+num2str($par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[6]
	endif
	if (cmpstr(model_name,"zwt2_k")==0)		// 2021-12-22
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		variable fa,fb,f,auc=($par[0]+$par[5])/$par[4]
		fa=$par[0]/$par[1]*(1-exp(-l1*$par[1]))/l1*exp(-l1*$par[6])+$par[5]/$par[6]*(1-exp(-l1*$par[6])/l1)
		fb=$par[0]/$par[1]*(1-exp(-l2*$par[1]))/l2*exp(-l2*$par[6])+$par[5]/$par[6]*(1-exp(-l2*$par[6])/l2)
		f=auc*(l1-l2)/(($par[2]/l2-1)*fb*exp(l2*($par[1]+$par[6]))-($par[2]/l1-1)*fa*exp(l1*($par[1]+$par[6])))
		tagtxt = "\\f02F\f00\B1\M\f02D/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00\B1\M = "+num2str($par[1])+" h"
		tagtxt+= "\r\\f02F\f00\B2\M\f02D/V\\f00\\Bc\\M = "+num2str($par[5])+",  \\f02o\\f00\B2\M = "+num2str($par[6])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M, AUC = "+num2str(auc)
		tagtxt+="\r\\f02k\\f00\\B12\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\S-1\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[4])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h, \f02F\f00 = "+num2str(f)
		tau=$par[1]+$par[6]
	endif
	if (cmpstr(model_name,"zwt3")==0) 		// 2021-12-06
		tagtxt = "\\f02F\\B1\\MD/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00\\B1\\M = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\B2\\MD/V\\f00\\Bc\\M = "+num2str($par[5])+",  \\f02o\\f00\\B2\\M = "+num2str($par[6])+" h"
		tagtxt+="\r\\f02F\\B3\\MD/V\\f00\\Bc\\M = "+num2str($par[7])+",  \\f02o\\f00\\B3\\M = "+num2str($par[8])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M,  \\f02k\\f00\\B12\\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[4])+" h\\S-1\\M, \\f02k\\f00\\B10\\M = "+num2str($par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[6]+$par[8]
	endif
	if (cmpstr(model_name,"zwt3_k")==0)		// 2021-12-22
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		tagtxt = "\\f02F\f00\B1\M\f02D/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00\B1\M = "+num2str($par[1])+" h"
		tagtxt+= "\r\\f02F\f00\B2\M\f02D/V\\f00\\Bc\\M = "+num2str($par[5])+",  \\f02o\\f00\B2\M = "+num2str($par[6])+" h"
		tagtxt+= "\r\\f02F\f00\B3\M\f02D/V\\f00\\Bc\\M = "+num2str($par[7])+",  \\f02o\\f00\B3\M = "+num2str($par[8])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M\r\\f02k\\f00\\B12\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\S-1\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[4])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[6]+$par[8]
	endif
	if (cmpstr(model_name,"zwt4")==0) 		// 2022-09-06
		tagtxt = "\\f02F\\f00\\B1\\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00\\B1\\M = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\f00\\B2\\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[5])+",  \\f02o\\f00\\B2\\M = "+num2str($par[6])+" h"
		tagtxt+="\r\\f02F\\f00\\B3\\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[7])+",  \\f02o\\f00\\B3\\M = "+num2str($par[8])+" h"
		tagtxt+= "\r\\f02F\\f00\\B4\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[9])+",  \\f02o\\f00\B4\\M = "+num2str($par[10])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M,  \\f02k\\f00\\B12\\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[4])+" h\\S-1\\M, \\f02k\\f00\\B10\\M = "+num2str($par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[6]+$par[8]
	endif
	if (cmpstr(model_name,"zwt4_k")==0)		// 2022-09-06
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		tagtxt = "\\f02F\\f00\\B1\\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[0])+",  \\f02o\\f00\\B1\\M = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02F\\f00\\B2\\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[5])+",  \\f02o\\f00\\B2\\M = "+num2str($par[6])+" h"
		tagtxt+="\r\\f02F\\f00\\B3\\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[7])+",  \\f02o\\f00\\B3\\M = "+num2str($par[8])+" h"
		tagtxt+= "\r\\f02F\\f00\\B4\M\\f02D/V\\f00\\Bc\\M = "+num2str($par[9])+",  \\f02o\\f00\B4\\M = "+num2str($par[10])+" h"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M\r\\f02k\\f00\\B12\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\S-1\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[4])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]+$par[6]+$par[8]
	endif
	if (cmpstr(model_name,"fw")==0)		// 2021-08-17
		tagtxt =  "\\f02FD/V\\f00\\Bc\\M = "+num2str($par[0])
		tagtxt+="\r\\f02k\\f00\\Ba\\M = "+num2str($par[1])+" h\\S-1\\M\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M, \f02k\f00\B12\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[4])+" h\\S-1\\M, \f02k\f00\B10\M = "+num2str($par[3]*$par[4]/$par[2])+" h\S-1\M"
		tagtxt+="\r\f02t\f00\Bmax\M = "+num2str(V_maxloc)+" h"
	endif
	if (cmpstr(model_name,"fw_k")==0)		// 2021-08-22
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		tagtxt = "\\f02FD/V\\f00\\Bc\\M = "+num2str($par[0])+"\r\\f02k\\f00\\Ba\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B12\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[4])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\\S-1\\M"
		tagtxt+="\r\\f02t\\f00\\Bmax\\M = "+num2str(V_maxloc)+" h"
	endif
	if (cmpstr(model_name,"fwt")==0)
		tagtxt= "\\f02FD/V\\f00\Bc\M = "+num2str($par[0])+", \\f02o\\f00 = "+num2str($par[5])+" h"
		tagtxt+="\r\\f02k\\f00\\Ba\\M = "+num2str($par[1])+" h\\S-1\\M\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[3])+" h\\S-1\\M, \\f02k\\f00\B12\\M = "+num2str($par[3]+$par[4]-$par[2]-$par[3]*$par[4]/$par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02a\\f00 = "+num2str($par[4])+" h\\S-1\\M, \\f02k\\f00\B10\\M = "+num2str($par[3]*$par[4]/$par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02t\\f00\\Bmax\\M = "+num2str(V_maxloc)+" h"
		tau=$par[5]
	endif
	if (cmpstr(model_name,"fwt_k")==0)		// 2021-08-16
		variable b=$par[3]+$par[2]+$par[4],c=$par[2]*$par[4],d=sqrt(b*b-4*c),l1=0.5*(b+d),l2=0.5*(b-d)		// b,c,d nomenclature from trinomial
		tagtxt = "\\f02FD/V\\f00\Bc\M = "+num2str($par[0])+", \\f02o\\f00 = "+num2str($par[5])+" h"
		tagtxt+="\r\\f02k\\f00\\Ba\\M = "+num2str($par[1])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B21\\M = "+num2str($par[2])+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B12\\M = "+num2str($par[3])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l1)+" h\\S-1\\M"
		tagtxt+="\r\\f02k\\f00\\B10\\M = "+num2str($par[4])+" h\\S-1\\M, \\f02a\\f00 = "+num2str(l2)+" h\\S-1\\M"
		tagtxt+="\r\\f02t\\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[5]
	endif
	if (cmpstr(model_name,"zzt2")==0)		// 2021-08-16
		tagtxt = "\\f02FD/V\\f00\Bd\M = "+num2str($par[0])
		tagtxt+="\r\\f02o\\f00\B1\M = "+num2str($par[1])+" h"
		tagtxt+="\r\\f02o\\f00\B2\M = "+num2str($par[2])+" h"
		tagtxt+="\r\\f02t\\f00\Bmax\M = "+num2str(V_maxloc)+" h"
		tau=$par[1]
	endif
	print $par
	if (cmpstr(annotate_simulation,"no")==0)
	else
		TextBox/C/N=text0/F=0/B=1/A=LC/X=65/Y=0 tagtxt
	endif
	if(tau>0)
		make/o/n=1 $yt=$model($par,tau)
		SetScale/P x tau,1,"", $yt
		if (ItemsInList(wavelist(yt,";","WIN:"))==0)
			AppendToGraph $yt
			ModifyGraph mode($yt)=3,marker($yt)=17,rgb($yt)=(0,0,0),msize=3
		endif
	endif
End

Proc map_chisq(ctrlName) : ButtonControl		// 2021-12-29
	String ctrlName
	variable steps=10		// number of values to be tried for logarithmic partitioning of each parameter range
	variable i=0
	string par="par_"+model_name,parl="parl_"+model_name,parh="parh_"+model_name,parlabels="plb_"+model_name
	string model="cf_"+model_name,H="H_"+model_name
	do
		if (!exists(parl) || !exists(parh))
			print "Creating ranges for model", model_name
			make/n=(numpnts($par)) $parl=0.01,$parh=1000
			i=0
			do
				if (char2num($parlabels[i][0])==107)		// parameter label starts with "k"
					$parh[i]=5
				endif
				i+=1
			while (i<numpnts($parh))
		endif
		i=0
		do
			if (char2num($parlabels[i][0])==116)			// parameter label starts with "t" (for "tau")
				$parh[i]=min($parh[i],t_data[numpnts(t_data)-1])
			endif
			i+=1
		while (i<numpnts($parh))
		variable npar=numpnts($par),step=0,tot_steps=steps^npar,ipar,rem,div,cc,ss,startp,endp,j,ss_min=1e100,step_min,cc_min
		variable/c ss_cc
		j=numpnts($H)-1
		startp=$H[j-1]
		endp=($H[j]==0)*(numpnts(c_data)-1)+($H[j]>0)*$H[j]
		duplicate/o $par pars
		duplicate/o c_data sim_data
		do
			ipar=0
			div=step
			do
				rem=mod(div,steps)		// 2022-01-08 corrected calculation
				div=(div-rem)/steps
				pars[ipar]=$parl[ipar]*($parh[ipar]/$parl[ipar])^(rem/(steps-1))
//				print step,ipar,rem,div
				ipar+=1
			while (ipar<npar)
			sim_data=$model(pars,t_data)
			ss_cc=chisq_corrcoef(sim_data,c_data,startp,endp)
			ss=real(ss_cc)
			cc=imag(ss_cc)
			if (ss<ss_min)
				ss_min=ss
				cc_min=cc
				step_min=step
				printf "%.0f: par%s= {", step,model_name
				ipar=0
				do
					printf "%f, ",pars[ipar]
					ipar+=1
				while (ipar<npar)
				printf "} ss = %f, cc = %f\r",ss,cc
			endif
			step+=1
		while (step<tot_steps)
		div=step_min
		ipar=0
		do
			rem=mod(div,steps)
			div=(div-rem)/steps
			pars[ipar]=$parl[ipar]*($parh[ipar]/$parl[ipar])^(rem/(steps-1))
//			print step,ipar,rem,div
			ipar+=1
		while (ipar<npar)
		printf "%s %s %.0f: ss = %g, cc = %g", model_name,data_label,step_min,ss_min,cc_min
		ipar=0
		do
			printf ", %s = %.5g ",$parlabels[ipar], pars[ipar]
			ipar+=1
		while (ipar<npar)
		printf "\r"
//		print ss_min,step_min
	while (0)
End

Proc AddNoiseProc(ctrlName) : ButtonControl		// 2021-10-02
	String ctrlName
	string data = "sim_"+model_name,dataP = "simP_"+model_name,par = "par_"+model_name
	string w_name = "Gsim_"+model_name, model = "cf_"+model_name,yt = "yt_"+model_name
	DoWindow/F $w_name
	make/o/n=(numpnts(sample_times)) $data=$model($par,sample_times)*(1+gnoise(nslevel))
	Append $data vs sample_times
	RemoveFromGraph/Z $data, $dataP
	modify mode=3,marker($data)=19,msize($data)=3
	if (ItemsInList(wavelist(yt,";","WIN:"))==1)
		ReorderTraces $yt {$data}
	endif
	string nstag="\f02s\f00 = "+num2str(nslevel)
	TextBox/C/N=text_ns/F=0/B=1/A=LC/X=91/Y=50 nstag
End

Proc CopySimProc(ctrlName) : ButtonControl		// 2021-10-04
	String ctrlName
	string data = "sim_"+model_name,par="par_"+model_name,w_name="Gsim_"+model_name, wfit_name="Gcf_"+model_name
	data_label=data
	DoWindow/F G_data
	Duplicate/o $data c_data
	if (deltax(c_data) != 1)
		Duplicate/o $data t_data
		t_data=x
	else
		Duplicate/o sample_times t_data
	endif
	TextBox/C/N=text0 data_label
	Update_data(data_label,"ng/mL")
	DoWindow/F $wfit_name
End

Proc NSimFit_proc(ctrlName) : ButtonControl		// 2021-10-04
	String ctrlName
	string par="par_"+model_name, fit_res="FitRes_"+model_name,Tr="Tres_"+model_name,cf="cf_"+model_name,H="H_"+model_name
	variable tau
	variable i=0
	holdflags=""
	do
		holdflags+=num2str($H[i])
		i+=1
	while (i<numpnts($H)-2)
	variable j=numpnts($H)-1, startp=$H[j-1], endp=($H[j]==0)*(numpnts(c_data)-1)+($H[j]>0)*$H[j]
	duplicate/o $par pars
	variable npar=numpnts(pars)
	make/o/n=(n_fits+1,2*npar+2) $fit_res=NaN
	DoWindow/K $Tr
	edit/N=$Tr/W=(2,200,600,500) $fit_res
	ModifyTable width(Point)=20,format($fit_res)=3,digits($fit_res)=4,width($fit_res)=40
	$fit_res[0][0,npar-1]=$par[y]
	t_data=sample_times
	duplicate/o t_data,sim_data
	i=0
	do
		make/o/n=(numpnts(sample_times)) c_data=$cf($par,sample_times)*(1+gnoise(nslevel))
		pars=$par
 		FuncFit/Q/H=holdflags/NTHR=0/M=2 $cf pars  c_data[startp,endp] /X=t_data
		sim_data=$cf(pars,t_data)
		$fit_res[i+1][0,npar-1]=pars[y]
		$fit_res[i+1][npar,2*npar-1]=W_sigma[y-npar]
		$fit_res[i+1][2*npar]=V_chisq
		$fit_res[i+1][2*npar+1]=corrcoef(sim_data,c_data,startp,endp)
		i+=1
	while (i<n_fits)	
End

Macro AddToLabelList(newlabel)		// 2022-01-28 improved version
  string newlabel
  prompt newlabel, "Enter new data label"
  if (strlen(newlabel)>0)
    labellist=RemoveListItem(0,labellist)
    labellist=SortList(AddListItem(newlabel,labellist))
    labellist=AddListItem(newlabel,labellist)
//    redimension/n=(ItemsInList(labellist)) data_labels
//    data_labels=StringFromList(x,labellist)
//    labellist=
  endif
End

Function/s RebuildDataLabels()
  wave/t data_labels=data_labels
  variable n = numpnts(data_labels),i=0
  string out=""
  do
    out+=data_labels[i]+";"
    i+=1
  while (i<n)
  return out
End

Macro List2Twave(list,twave)		//2022-08-27
// the converse of TextWave2Str
  string list,twave
  make/t/o/n=(ItemsInList(list)) $twave
  $twave=StringFromList(x,list)
End

Function/s TextWave2Str(twave)		//2022-08-27
// the converse of twave=StringFromList(out,x) or List2Twave
  wave/t twave
  variable n = numpnts(twave),i=0
  string out=""
  do
    out+=twave[i]+";"
    i+=1
  while (i<n)
  return out
End

Function ShowWindowProc(ctrlName) : ButtonControl		// 2021-07-17,18,21
	String ctrlName
	svar mn=model_name
	string win
	if (char2num(ctrlName[7])==116)
		win="Gcf_"+mn
	else
		if (char2num(ctrlName[7])==109)
			win="Gsim_"+mn
		else
			win="Tcf_"+mn
		endif
	endif
	if (ItemsInList(WinList(win,";","WIN:3"))==1)
		DoWindow/F $win
	else
		print "Window ",win," does not exist."
		if (char2num(ctrlName[7])==116)
			wave t_Data=t_data,c_data=c_data,s_data=s_data
			svar data_label=data_label
			Display /W=(331.2,186.8,678.6,394.4) c_data vs t_data
			DoWindow/C $win
			ModifyGraph mode=4, marker(c_data)=19, rgb(c_data)=(0,39168,0), fSize=12
			string c_unit=StringFromList(0,note(c_data))
			string t_unit=StringFromList(1,note(c_data))
			Label left "\f02C\f00 ("+c_unit+")"
			Label bottom "\\f02t\\f00 ("+t_unit+")"
			ErrorBars c_data Y,wave=(s_data,s_data)
			Cursor/P A c_data 10
			ShowInfo
			TextBox/C/N=text0/A=MC/X=3.76/Y=-35.86 data_label
			string Fit_button="Fit_"+mn,Show_coefs_button="Show_coefs_"+mn,Sim_button="Sim_"+mn,Cit_button="Cit_"+mn
			Button $Fit_button,pos={1,2},size={30,20},proc=ButtonRunFit,title="Fit"
			Button $Show_coefs_button,pos={61,2},size={50,20},proc=ButtonShowParams,title="Params"
			Button $Show_coefs_button,fColor=(52225,0,20736)
			Button $Sim_button,pos={31,2},size={33,20},proc=ButtonRunFit,title="Simu"
			Button $Sim_button,fColor=(65280,43520,0)
			Button $Cit_button,pos={1,22},size={30,20},proc=ButtonRunFit,title="cFit"
		endif
	endif
End

Function model_select(ctrlName,popNum,popStr) : PopupMenuControl		// 2021-07-17
	String ctrlName
	Variable popNum
	String popStr
	svar model_name=model_name
	model_name=popStr
End

Function annotate_simulation_select(ctrlName,popNum,popStr) : PopupMenuControl		// 2021-12-19/22
	String ctrlName
	Variable popNum
	String popStr
	svar annotate_simulation=annotate_simulation
	annotate_simulation=popStr
End

Function ShowModelsProc(ctrlName) : ButtonControl		// 2021-07-22
	String ctrlName
	string winname="T_models"
	if (ItemsInList(winlist(winname,",",""))==1)
		DoWindow/F $winname
	else
		execute winname+"()"
	endif
End

Function ShowDataTableProc(ctrlName) : ButtonControl		// 2021-07-22
	String ctrlName
	DoWindow/F T_curve
End

Function AddDataLabelsProc(ctrlName) : ButtonControl		// 2021-07-22
	String ctrlName
	execute "AddToLabelList()"
End

Function ShowDataProc(ctrlName) : ButtonControl		// 2021-07-22
	String ctrlName
	string wn="G_data"
	if (ItemsInList(WinList(wn,";","")))
		DoWindow/F $wn
	else
		string comm = wn+"()"
		execute comm
	endif
End

Function ShowCovarProc(ctrlName) : ButtonControl		// 2022-07-10
	String ctrlName
	DoWindow/F T_Covar
End

Function ReDrawPanel(ctrlName) : ButtonControl
	String ctrlName
	DoWindow/K P_sim
	execute "P_sim()"
End

Function HideWindowsProc(ctrlName) : ButtonControl	// 2022-07-10
	String ctrlName
	string wl=winlist("*",";","WIN:3"),wn
	wl=removefromlist("T_curve;G_data;T_Covar",wl)
	variable n=itemsinlist(wl),i=0
	do
		wn=stringfromlist(i,wl)
		DoWindow/HIDE=1 $wn
		i+=1
	while (i<n)
End

Function ButtonProc(ctrlName) : ButtonControl
	String ctrlName
	DoWindow/F Gcf_ln
End

Function ButtonProc_1(ctrlName) : ButtonControl
	String ctrlName
	DoWindow/F Gcf_dc
End

// ======================= Combined fitting macros =============

Macro fits1(lbl)		// 2021-04-30
  string lbl
  variable i=0,i_max,t_max,c_max,i_start_dc
  string fitlabel,printlbl=lbl,templbl
  do
    if (c_data[i]>c_max)
      c_max=c_data[i]
      t_max=t_data[i]
      i_max=i
    endif
    i+=1
  while (i<numpnts(c_data))
  DoWindow/F Gcf_dC		// 1st analysis
  TextBox/C/N=text0 lbl
  duplicate/o c_data, dc_data
  Differentiate dc_data/X=t_data
  if (strlen(CsrInfo(B)) > 0)
    i_start_dc=pcsr(B)
  else
    i_start_dc=i_max+2
  endif
  CurveFit/Q/M=2/W=0 line, dc_data[i_start_dc,]/X=c_data/D	// fit dc/dt = f(c)
  duplicate/o dc_data dcs_data
  wavestats/q dc_data
  setaxis left V_min*1.1,0
  dcs_data=poly(W_coef,c_data)
  cc=corrcoef(dcs_data,dc_data,i_max+1,numpnts(c_data)-1)
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\rd\f02c\f00/d\f02t\f00\\Boo\\M = "+num2str(W_coef[0])+" ± "+num2str(W_sigma[0])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_dc_data fitlabel
//  print W_coef[0]," ±", W_sigma[0],W_coef[1]," ±", W_sigma[1], lbl
  sprintf templbl " kel = %.3f (±%.3f) h^-1 %s",-W_coef[1],W_sigma[1],lbl
  printlbl+=templbl
  DoWindow/F Gcf_ln			// 2nd analysis
  TextBox/C/N=text0 lbl
  duplicate/o c_data lnc_data
  lnc_data=ln(c_data)
  CurveFit/Q/M=2/NTHR=0 line  lnc_data[i_max,] /X=t_data /D /I=1	// fit lnc = f(t) without weights
  duplicate/o c_data cs_lndata
  cs_lndata=poly(W_coef,t_data)
  cc=corrcoef(cs_lndata,lnc_data,i_max,numpnts(c_data)-1)
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_lnc_data fitlabel
  sprintf templbl "	ln fit:	tmax = %1.2f h	kel = %.3f (±%.3f) h^-1	R^2 = %1.5f",t_max,-W_coef[1],W_sigma[1],cc
  printlbl+=templbl
  duplicate/o c_data s_lndata		// 3rd analysis
  s_lndata=c_data^-2
  CurveFit/Q/W=0/M=2 line, lnc_data[i_max,]/X=t_data[i_max,] /W=s_lndata/I=1	// fit lnc = f(t) to line with weights for ln(c) data
  duplicate/o c_data cs_data
  cs_data=exp(poly(W_coef,t_data))
  cc=corrcoef(cs_data,c_data,i_max,numpnts(c_data)-1)
//  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])
//  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_c_data fitlabel
  sprintf templbl "	w-ln fit:	kel = %.3f (±%.3f) h^-1	R^2 = %1.5f",-W_coef[1],W_sigma[1],cc
  printlbl+=templbl
  DoWindow/F Gcf_zot1		// 4th analysis
  TextBox/C/N=text0 lbl
  par_zot1[2]=t_max
  par_zot1[1]=-w_coef[1]
  par_zot1[0]=exp(W_coef[0])
  RemoveFromGraph/Z fit_zot1,res_zot1
  killwaves/Z fit_zot1,res_zot1
  FuncFit/Q/NTHR=0/M=2 cf_zot1 par_zot1  c_data /X=t_data  /I=1 /D /R		// fit non-lin 2021-04-14
  AddCorrelation()
  rename fit_c_data fit_zot1
  rename res_c_data,res_zot1
  duplicate/o c_data,sim_data
  sim_data=cf_zot1(par_zot1,t_data)
  cc=corrcoef(sim_data,c_data,0,numpnts(c_data)-1)
  fitlabel="\f02FD/V\f00\Bd\M = "+num2str(par_zot1[0])+" ±"+num2str(W_sigma[0])+"\r\f02o\f00 = "+num2str(par_zot1[2])+" ± "+num2str(W_sigma[2])
  fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str(par_zot1[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_c_data fitlabel
  sprintf templbl "	o = %1.2f (±%.2f) h	kel = %.3f (±%.3f) h^-1	R^2 = %1.5f",par_zot1[2],W_sigma[2],par_zot1[1],W_sigma[1],cc
  printlbl+=templbl
  par_fo[0]=par_zot1[0]
  par_fo[2]=par_zot1[1]
  fit_fo(lbl)			// 5th analysis
  sprintf templbl "	Bateman:	ka = %.3f (±%.3f) h^-1	kel = %.3f (±%.3f) h^-1	R^2 = %1.5f",par_fo[1],W_sigma[1],par_fo[2],W_sigma[2],cc
  printlbl+=templbl
  par_fot[0]=par_fo[0]
  par_fot[1]=par_fo[1]
  par_fot[2]=par_fo[2]
  par_fot[3]=0.5*(t_max+par_zot1[2])
  fit_fot(lbl)			// 6h analysis
  sprintf templbl "	trunc-Bateman:	ka = %.3f (±%.3f) h^-1	kel = %.3f (±%.3f) h^-1	o = %1.2f (±%.2f) h	R^2 = %1.5f",par_fot[1],W_sigma[1],par_fot[2],W_sigma[2],par_fot[3],W_sigma[3],cc
  printlbl+=templbl
  print printlbl
  DoWindow/H/F
End

Macro AUCs1(lbl)		// 2021-04-12
  string lbl
  variable i=0,i_max,t_max,c_max,AUC,AUCiv,AUCr,s_AUCr,s_AUCr2,cc,kt,ekt
  string printlbl,templbl,fitlabel
  DoWindow/F G_AUC
  TextBox/C/N=text0 lbl
  do
    if (c_data[i]>c_max)
      c_max=c_data[i]
      t_max=t_data[i]
      i_max=i
    endif
    i+=1
  while (i<numpnts(c_data))
  duplicate/o c_data lnc_data
  lnc_data=ln(c_data)
  duplicate/o c_data s_lndata
  s_lndata=c_data^-2
  CurveFit/Q/W=0/M=2 line, lnc_data[i_max,]/X=t_data[i_max,] /W=s_lndata/D/I=1	// fit lnc = f(t) to line with weights for ln(c) data
  AUC=areaXY(t_data,c_data)-c_data[200]/W_coef[1]	// adding tail by extrapolation 2021-04-13
  print W_coef[0],W_coef[1]
  duplicate/o c_data cs_data
  cs_data=exp(poly(W_coef,t_data))
  cc=corrcoef(cs_data,c_data,i_max,numpnts(c_data)-1)
  AUCiv=-exp(W_coef[0])/W_coef[1]
  AUCr=AUC/AUCiv
  s_AUCr=AUC*exp(-W_coef[0])*sqrt(W_sigma[1]^2+(W_sigma[0]*W_coef[1])^2-2*V_Rab*W_sigma[0]*W_sigma[1]*W_coef[1])
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  fitlabel+="\rAUC ratio = "+num2str(AUCr)+" ± "+num2str(s_AUCr)
  TextBox/C/N=CF_c_data fitlabel
  sprintf printlbl "tmax = %1.2f h kel = %.3f (±%.3f) h^-1 F = %.3f (±%.3f)",t_max,-W_coef[1],W_sigma[1],sqrt(AUCr),s_AUCr
  print printlbl
  if (0)
  DoWindow/F Gcf_zot1		// 2nd analysis
  TextBox/C/N=text0 lbl
  par_zot1[2]=t_max
  par_zot1[1]=-w_coef[1]
  RemoveFromGraph/Z fit_zot1,res_zot1
  killwaves/Z fit_zot1,res_zot1
  FuncFit/Q/NTHR=0/M=2 cf_zot1 par_zot1  c_data /X=t_data  /I=1 /D /R		// fit non-lin 2021-04-14
  AddCorrelation()
  rename fit_c_data fit_zot1
  rename res_c_data,res_zot1
  duplicate/o c_data,sim_data
  sim_data=cf_zot1(par_zot1,t_data)
  kt=par_zot1[1]*par_zot1[2]
  ekt=exp(kt)-1
  cc=corrcoef(sim_data,c_data,0,numpnts(c_data)-1)
  s_AUCr2=(1-kt/ekt)/ekt*sqrt(M_covar[1][1]*par_zot1[2]^2+M_covar[2][2]*par_zot1[1]^2+2*M_covar[1][2]*kt)
  fitlabel="\f02FD/V\f00\Bd\M = "+num2str(par_zot1[0])+" ±"+num2str(W_sigma[0])+"\r\f02o\f00 = "+num2str(par_zot1[2])+" ± "+num2str(W_sigma[2])
  fitlabel+="\r\f02k\f00\\Bel\\M = "+num2str(par_zot1[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
//  fitlabel+="\rAUC = "+num2str(par_1[0]/par_1[1])
//  fitlabel+=" ± "+num2str(sqrt(M_covar[0][0]+M_covar[1][1]*(par_1[0]/par_1[1])^2-2*M_covar[0][1]*par_1[0]/par_1[1])/par_1[1])
//  fitlabel+="\rAUC ratio = "+num2str(kt/ekt)+" ± "+num2str(s_AUCr2)
  TextBox/C/N=CF_c_data fitlabel
  sprintf templbl " o = %1.2f (±%.2f) h  kel = %.3f (±%.3f) h^-1 F = %.3f (±%.3f)",par_zot1[2],W_sigma[2],par_zot1[1],W_sigma[1],kt/ekt,s_AUCr2
  printlbl+=templbl
  DoWindow/F Gcf_ln		// nth analysis
  TextBox/C/N=text0 lbl
  CurveFit/Q/M=2/NTHR=0 line  lnc_data[i_max,] /X=t_data /D /I=1	// fit lnc = f(t) without weights
  duplicate/o c_data cs_lndata
  cs_lndata=poly(W_coef,t_data)
  cc=corrcoef(cs_lndata,c_data,i_max,numpnts(c_data)-1)
  AUCiv=-exp(W_coef[0])/W_coef[1]
  AUCr=AUC/AUCiv
  s_AUCr=AUC*exp(-W_coef[0])*sqrt(W_sigma[1]^2+(W_sigma[0]*W_coef[1])^2-2*V_Rab*W_sigma[0]*W_sigma[1]*W_coef[1])
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
//  fitlabel+="\rAUC ratio = "+num2str(AUCr)+" ± "+num2str(s_AUCr)
  TextBox/C/N=CF_lnc_data fitlabel
  sprintf templbl " kel = %.3f (±%.3f) h^-1 F = %.3f (±%.3f)",par_zot1[1],W_sigma[1],AUCr,s_AUCr
//  printlbl+=templbl
  DoWindow/F Gcf_dC		// (n+1)th analysis
  TextBox/C/N=text0 lbl
  duplicate/o c_data, dc_data
  Differentiate dc_data/X=t_data
  CurveFit/Q/M=2/W=0/TBOX=(0x300) line, dc_data[i_max+2,]/X=c_data/D	// fit dc/dt = f(c)
  duplicate/o dc_data dcs_data
  wavestats/q dc_data
  setaxis left V_min*1.1,0
  dcs_data=poly(W_coef,c_data)
  cc=corrcoef(dcs_data,dc_data,i_max+1,numpnts(c_data)-1)
  fitlabel="\f02k\f00\\Bel\\M = "+num2str(-W_coef[1])+" ± "+num2str(W_sigma[1])
  fitlabel+="\rd\f02c\f00/d\f02t\f00\\Boo\\M = "+num2str(W_coef[0])+" ± "+num2str(W_sigma[0])
  fitlabel+="\r\f02?\f00\\S2\\M = "+num2str(V_chisq)+",  \f02R\f00\S2\M = "+num2str(cc)
  TextBox/C/N=CF_dc_data fitlabel
//  print W_coef[0]," ±", W_sigma[0],W_coef[1]," ±", W_sigma[1], lbl
  sprintf templbl " kel = %.3f (±%.3f) h^-1 %s",-W_coef[1],W_sigma[1],lbl
  printlbl+=templbl
//  print printlbl
  DoWindow/H/F
  endif
End

// =================== Window reconstruction macros ===================

Window T_models() : Table
	PauseUpdate; Silent 1		// building window...
	Edit/W=(708,224,1038,609.8) models,model_description
	ModifyTable format(Point)=1,width(Point)=26,width(models)=47,width(model_description)=236
EndMacro

Window G_data() : Graph		// macro made 2022-05-24
	PauseUpdate; Silent 1		// building window...
do
	if (itemsinlist(winlist("G_data",";","")))
		DoWindow/F G_data
		break
	endif
	Display /W=(407.4,45.8,756,253.4) c_data vs t_data
	AppendToGraph/R auc_data vs t_data
	ModifyGraph mode(c_data)=4
	ModifyGraph marker(c_data)=19
	ModifyGraph rgb(c_data)=(0,39168,0)
	ModifyGraph fSize(left)=12,fSize(bottom)=12
	ModifyGraph lblPos(left)=53
	Label left "\\f02C\\f00 (\\Eig/mL)"
	Label bottom "\\f02t\\f00 (hours)"
	ErrorBars c_data Y,wave=(s_data,s_data)
	Cursor/P A c_data 7
	ShowInfo
	TextBox/C/N=text0/A=MC/X=25.88/Y=20.08 data_label		//2022-07-10
	Button button0,pos={1,2},size={32,20},proc=Update_data_ButtonProc,title="Run"
	Button button0 help={"Select appropriate data label and concentration and time units, then calculate AUC curve, ln and derivative"}
	Button button_lnc,pos={1,22},size={32,20},proc=ButtonProc,title="ln c"
	Button button_lnc help={"Go to to ln c graph"}
	Button button_dc,pos={1,43},size={32,19},proc=ButtonProc_1,title="dc/dt"
	Button button_dc help={"Go to dc/dt graph"}
while (0)
EndMacro

Window Gcf_dc() : Graph		// macro made 2022-05-24
	PauseUpdate; Silent 1		// building window...
	Display /W=(376.8,341,727.8,548) dc_data vs c_data
	AppendToGraph fit_dc_data
	ModifyGraph mode(dc_data)=3
	ModifyGraph marker(dc_data)=18
	ModifyGraph rgb(dc_data)=(0,0,52224)
	ModifyGraph fSize=12
	ModifyGraph notation(left)=1
	Label left "d\\f02C\\f00/d\\f02t\\f00 (\\E)"
	Label bottom "\\f02C\\Bb\\M"
	SetAxis left -0.673827970027924,0
	Cursor/P B dc_data 3
	ShowInfo
	TextBox/C/N=text0/X=9.09/Y=66.13 ""
	TextBox/C/N=CF_dc_data ""
	Button Fit_dc,pos={1,2},size={50,20},proc=ButtonRunFit,title="Fit"
	Button Fit_dc,help={"Run fit from selected point to the end of data"}
EndMacro

Window Gcf_ln() : Graph		// macro made 2022-05-24
	PauseUpdate; Silent 1		// building window...
	Display /W=(394.8,418.4,743.4,625.4) lnc_data vs t_data
	AppendToGraph fit_lnc_data,lnc_line
	ModifyGraph mode(lnc_data)=4
	ModifyGraph marker(lnc_data)=17
	ModifyGraph lSize(fit_lnc_data)=3
	ModifyGraph rgb(lnc_data)=(26112,17408,0),rgb(lnc_line)=(0,0,52224)
	ModifyGraph msize(lnc_data)=3
	ModifyGraph fSize=12
	ModifyGraph standoff=0
	Label left "ln\\f02C\\f00"
	Label bottom "\\f02t\\f00 (hours)"
	Cursor/P A lnc_data 2;Cursor/P B lnc_data 15
	ShowInfo
	TextBox/C/N=text0/X=60.85/Y=78.23 ""
	TextBox/C/N=CF_lnc_data ""
	Button Fit_ln,pos={1,2},size={50,20},proc=ButtonRunFit,title="Fit",help={"Run fit from Cmax to end"}
EndMacro

Macro Gcf_setup(mdl,n)		// for building tables and graphs for a new function/model
	string mdl		// (new) model to set up fitting window and parameters 2022-05-28
	variable n		// number of parameters in model
	string plb="plb_"+mdl,par="par_"+mdl, cf="cf_"+mdl,H="H_"+mdl, parl="parl_"+mdl, parh="parh_"+mdl
	string Tcf="Tcf_"+mdl,Gcf="Gcf_"+mdl
	do
	if (exists(cf)!=6)
		print "There is no "+mdl+" procedure. You need to write that first."
		break
	endif
	make/t/o/n=(n) $plb
	make/o/n=(n) $par,$parl,$parh
	make/o/n=(n+2) $H
	if ($parh[0]==0)
		$parl=0.01
		$parh=1000
	endif
	if (itemsinlist(winlist(Tcf,";",""))==1)
		DoWindow/F $Tcf
	else
		edit $plb, $par, W_sigma, $H, $parl, $parh
		DoWindow/C $Tcf
	endif
	if (itemsinlist(winlist(Gcf,";",""))==1)
		DoWindow/F $Gcf
	else
		PauseUpdate; Silent 1		// building window...
		Display /W=(263.4,287.6,613.8,496.4) c_data vs t_data
		DoWindow/C $Gcf
		ModifyGraph mode=3, marker=19, rgb=(0,39168,0), fSize=12
		ModifyGraph lblPos(left)=53
		Label left "\\f02C\\f00 (\\Eig/mL)"
		Label bottom "\\f02t\\f00 (hours)"
		ErrorBars c_data Y,wave=(s_data,s_data)
		Cursor/P A c_data 0
		ShowInfo
		TextBox/C/N=CF_c_data/X=-1.15/Y=16.47 ""
		TextBox/C/N=text0/X=60.46/Y=59.44 ""
		string Fit_model="Fit_"+mdl,Show_coefs_model="Show_coefs_"+mdl,Sim_model="Sim_"+mdl,Cit_model="Cit_"+mdl
		Button $Fit_model,pos={1,2},size={30,20},proc=ButtonRunFit,title="Fit"
		Button $Fit_model,help={"H wave in parameter table indicates parameters held plus range of points to be used in fit; 0 means all at that end"}
		Button $Cit_model,pos={1,22},size={30,20},proc=ButtonRunFit,title="cFit"
		Button $Cit_model,help={"Do fit with constrains on parameters indicated in parl(ow) and parh(igh) wavesin parameter table"}
		Button $Sim_model,pos={31,2},size={33,20},proc=ButtonRunFit,title="Simu"
		Button $Sim_model,fColor=(65280,43520,0)
		Button $Sim_model,help={"Do not adjust parameters, only compare data with curve and calculate chi squared"}
		Button $Show_coefs_model,pos={61,2},size={50,20},proc=ButtonShowParams,title="Params"
		Button $Show_coefs_model,fColor=(52225,0,20736)
		Button $Show_coefs_model,help={"Bring to front parameter table"}
	endif
	while (0)
EndMacro