#pragma rtGlobals=1		// Use modern global access method.
Menu "macros"
	"TripleEfit"
	"RatioCa"
End

// INSTRUCTIONS   ( last updated by Haijiang Cai on 12/13/06 )
//
// These functions perform triple-component exponential fit to the Capacitance wave(Cm), on which you put Cursor A.
// Some of the calculations used here are adapted from MPI TripleExpFit (http://www.mpibpc.gwdg.de/abteilungen/140/software/index.html)
//
// TripleEfit ( )
//  Operation:
//	Find the point that capacitance starts to increase and put Cursor A on this point. 
// 	This point is not the actual start time of Cm increase. It is the initial guess of the start time of Cm increase.
//	If the fit does not meet the criterions, the function will use the point before or after this point as a new initial guess, and perform another fit.
//	This function will search this initial guess until the criterions are satisfied or the guess point is before the start of flash.
//  Criterions for fit: 
//	tau_fast  < 100 msec, tau_slow < 0.5 sec
//	Fast and slow amplitudes are larger than 0 
//   Notes:
// 	(1)	t0, t_start and fit_t0 are the initial guess of the time that Cm starts to increase. 
//		If you put Cursor A at the returned value of fit_t0, the macro should achieve a good fit with one attempt.
//		But if the returned fit_t0 is 0.499 (that is, one point before Flash_start), this means the attempted fit failed to meet the criterions.
//		In this case, you have to check the fitted curve (black), if it looks ok, we call it "trace fit". 
//		Only the value of Burstsize and Sustainedsize are trustable, the kinetic values are not reliable.
//		Or you can put your Cursor A on a larger x value (the bigger t0 above the command line) and try again.
//	(2)	Sometimes, you can find two fit_t0 points that can give satisfied fitting, this should not cause much problem. 
//		The results I checked so far are not significantly different.
//
// RatioCa ( )
//	This function convert the 350/380 ratio into the Ca concentration according to the standard Ratio-Calcium calibration
//	Before using this function, you have to make a ratio wave and a time wave first.

Function TripleEfit( )
	string CmName = "cm_1", fitCmName = "fit_cm_1"
	variable Length = 5, Flash_start = 0.5                                  // Usually we fit the Cm segment 5sec after the flash, which is triggered at 0.5sec
	prompt CmName, "Name the wave after baseline subtraction"
	prompt Length, "Set the time length of the fit (sec)"
	prompt Flash_start, "At what time(sec) the flash is triggered?"
	doprompt "3-Component Exponential Fit of the Cm Wave after Cursor(A)", CmName, Length, Flash_start
	
	string  cmWave	
	cmWave = CsrWave (A)							// Name of  the Cm wave that the cursor is on
	Duplicate/O $cmWave, wavetouse					// Never modify the original wave
	Variable Cm_base, Cm_start, Cm_end, dCm_preflash, t_end
	t_end = Flash_start + Length
	Cm_base = mean(wavetouse, 0.05, 0.1) 			// Cm value at the beginning, do not use the first point which is always strange
	Cm_start = mean(wavetouse, Flash_start - 0.1, Flash_start - 0.05) 		// Cm value preflash
	Cm_end = mean(wavetouse, t_end-0.05, t_end)            //  Cm value at the end of the fitting
	dCm_preflash = Cm_start - Cm_base				// This value is used to test how much secretion occurs before flash
	Wavetouse -= Cm_start							// Subtract the baseline
	
	Variable t_start, t_test1, t_test2, swtch, cnt = 0
	t_start = xCsr(A)									//  Important, initial guess of the start time of release
	t_test1 = t_start
	t_test2 = t_start
	swtch = 0

// Following do-while loop will search the start point from the xCsr(A) to find a fit that meet the criterions	
//*****************************************************************************************************************************
	do
		Print " "
		print "t0 =", t_start            						// This is the initial gues of start time for the fit, see Notes(1)
		cnt +=1										// Count how many times of fitting are used
		variable dT, dCm
		dT = Length	
		t_end = t_start+dT
		dCm = 1.5* (Cm_end- Cm_start)				// The value 1.5 is used by MPI TripleExpFit
		
		Duplicate/O/R=(t_start, t_end) wavetouse, wavetouse1
		Smooth/B 3, wavetouse1						// The after-flash wave is smoothed first before fitting 
		duplicate/o wavetouse1 WeightWave 			// Weight wave should have the same number of points as the Cm wave
		WeightWave = 1 								// Make y value of weight 1
		WeightWave[0, numpnts(WeightWave)/5] = 5 	// This weight value(5) and length(one fifth) are used by MPI TripleExpFit
	
		Make/N=7/D/O coWave 				// Coefficient wave with 7 coefficients	
		// Initial guesses:
		coWave[0] = 0.3*dCm					// Amplitude of the first (fast) component
		coWave[1] = 0.01*dT					// Tau of the first (fast) component, initial guess is 50msec
		coWave[2] = 0.3*dCm					// Amplitude of the second (slow) component
		coWave[3] = 0.1*dT					// Tau of the second (slow) component, initial guess is 0.5sec
		coWave[4] = 0.4*dCm 				// Amplitude of the third component
		coWave[5] = 2*dT					// Tau of the third component
		coWave[6] = t_start					// Start time of the capacitance increase, see Notes(1)
		
		FuncFit/Q/H="0000000" TripleE coWave wavetouse1/W=WeightWave /D/R 	// Do not hold the start time of release	
		
		print "A_fast =", coWave[0]*1e15, "fF;                   A_slow =", coWave[2]*1e15,"fF;"
		print "Tau_fast =", coWave[1], "sec;      Tau_slow =", coWave[3],"sec;"
		
		if (coWave[1]<0.1 && coWave[3]<1)	// See Criterions for fit	
			if (coWave[0]>0 && coWave[3]>0) 	// See Criterions for fit	
				break						// This will exit the do-while loop
			endif
		else
			if (swtch == 0)
				t_test1 += 0.001
				t_start = t_test1				// If the previous fit is not satisfied, then change the initial guess of start time
				swtch += 1
			elseif (swtch == 1)
				t_test2 -= 0.001
				t_start = t_test2				// If the previous fit is not satisfied, then change the initial guess of start time
				swtch -= 1
			endif	
		endif
	while (t_start >= Flash_start)
	print "times used for fitting", cnt	
	
	If (t_start == Flash_start-0.001)							// see Notes(1)
		doAlert 0, "Kinetic Fit failed. Click OK to check Trace Fit  or put your Cursor A on a larger x value and try again"
	endif
//*****************************************************************************************************************************
	variable dCm_startTime
	dCm_startTime = coWave[6]				// This is the "fitted" start time where Cm begins to increase after flash
	
	duplicate/o fit_wavetouse1, wavefit	
	variable burstCm, sustainCm, sustainRate, sustain_startTime, Burst, Sustain
	burstCm = coWave[0]+coWave[2]
	sustainCm = wavefit(dCm_startTime + Length) - coWave[0] - coWave[2] 
	Findlevel/Q wavefit, burstCm
	sustain_startTime = V_LevelX
	sustainRate = sustainCm/(dCm_startTime + Length - sustain_startTime)		// See Sorensen 2003, cell, method
	Burst = wavefit(dCm_startTime + 1)					// Release in the first 1sec after dCm start to increase
	Sustain = wavefit(dCm_startTime + Length) - Burst	// Release from 1sec after dCm starts to increase to the end of fit

	Make/o/n=1 fit_t0, Cm0 dCm0 Delay A_fast tau_fast A_slow tau_slow sustain_Rate BurstSize SustainedSize
		fit_t0[0] = t_start								// sec, see Notes(1)
		Cm0[0] = Cm_start *1e12						// pF, this is the basal Cm value, represents the size of the cell
		dCm0[0] = dCm_preflash*1e15					// fF, this reflects if any secretion occurs before flash
		Delay[0] = (dCm_startTime - Flash_start)*1e3 	// msec, Flash_start is set by your recording trigger
		A_fast[0] = coWave[0]*1e15 					// fF
		tau_fast[0] = coWave[1]*1e3 					// msec
		A_slow[0] = coWave[2]*1e15 					// fF
		tau_slow[0] = coWave[3] 						// sec
		sustain_Rate[0] = sustainRate*1e15			// fF/sec
		BurstSize[0] = Burst*1e15						// fF
		SustainedSize[0] = Sustain*1e15				// fF
				
	duplicate/o wavetouse $CmName
	display $CmName fit_wavetouse1
	ModifyGraph rgb(fit_wavetouse1)=(0,0,0)
	ModifyGraph lsize(fit_wavetouse1)=1.5
	
	Edit fit_t0 Cm0 dCm0 Delay A_fast A_slow tau_fast tau_slow sustain_Rate BurstSize SustainedSize
	
	Display Res_wavetouse1							// This is the residue wave, you can use this wave to estimate how well your fit is.
	Print "------------------------------------------------------------------------------------------------------------------------------------------------"
End

Function/D TripleE (w, x)
	wave/D w
	Variable/D x
	Variable result
	
	if (w[1]==0)
		w[1] = 1e-4
	endif
	if (w[3]==0)
		w[3] = 1e-4
	endif
	if(w[5]==0)
		w[5] = 1e-4
	endif
	
	result = w[0]*(1-exp(-(x-w[6])/w[1]))
	result += w[2]*(1-exp(-(x-w[6])/w[3]))
	result += w[4]*(1-exp(-(x-w[6])/w[5]))			// This equation is defined in Sorensen, 2003, cell
	return result
end


Function RatioCa ( )
	string Ratio_wave = "ratio1", Time_wave = "time1"
	variable Camin=1e-8, Camax=1e-2, CaFactor=1.01, Time_step = 0.5, Flash_startpnt = 6
	prompt Ratio_wave, "Source Ratio Wave", popup WaveList("*", ";", "")
	prompt Time_wave, "Source Time Wave", popup WaveList("*", ";", "")
	prompt Flash_startpnt, "How many points of the original ratio wave are before flash"
	prompt Time_step, "How often to average the Ca value (sec)"
	doprompt " Ratio = R0 + (R1*Ca)/(Ca + K1) + (R2*Ca)/(Ca + K2) ", Ratio_wave, Time_wave, Flash_startpnt, Time_step

//*****************************************************************************************************************************
// following lines are used to generate a standard Ratio-Ca curve
// R0, R1, Rmax, K1 and K2 are from in vivo Calcium calibration

	//variable R0=0.043, R1=0.0258, R2, Rmax=0.075, K1=0.0000664, K2=0.000000346   // calibration with Fura-4F/Furaptra on 12/28/2005
	variable R0=0.116, R1=0.0205, R2, Rmax=0.215, K1=0.000000162, K2=0.000167  // calibration with Fura-2/Furaptra on 6/10/2006
	R2 = Rmax - R0 - R1
	variable tableptr=0, stdN, stdCaV
	stdN=0
	stdCaV = Camin
	do
		stdCaV *= CaFactor
		stdN+=1
	while((stdCaV <= Camax))         // decide how many points needed

	make/o/n=(stdN) stdRatio stdCa
	stdCaV = Camin
	do
		stdRatio[tableptr] = R0 + (R1*stdCaV)/(stdCaV + K1) + (R2*stdCaV)/(stdCaV + K2)
		stdCa[tableptr]=stdCaV
		stdCaV *= CaFactor
		tableptr+=1
	while (stdCaV <= Camax)
	print "Total number of points on standard Ca-Ratio curve:", numpnts(stdRatio)
//	display stdRatio vs stdCa
//	ModifyGraph log(bottom)=1
//	ModifyGraph mode=3
//*****************************************************************************************************************************
//*****************************************************************************************************************************
// following lines are used to convert all ratio value to Ca value according to the standard Ratio-Ca wave above

	duplicate/o $Ratio_wave Ratio_wave1
	duplicate/o Ratio_wave1 Ca_wave
	variable tableptr1 = 0, tableptr2 = 0, Ratio_value = 0, Ratio_delta = 0, stdRatio_value = 0, stdNpnts = 0, ratioNpnts
	stdNpnts = numpnts(stdRatio)
	ratioNpnts = numpnts(Ratio_wave1)
	do 
		Ratio_value = Ratio_wave1[tableptr1]
		tableptr2=0
		do
			stdRatio_value = stdRatio[tableptr2]
			//print stdRatio_value
			if (stdRatio_value == Ratio_value)
				Ca_wave[tableptr1] = stdCa[tableptr2]*1e6						// make the final unit uM
				break		
			elseif (stdRatio_value > Ratio_value)
				Ca_wave[tableptr1] = (stdCa[tableptr2-1]+stdCa[tableptr2])/2*1e6	// make the final unit uM
				break
			endif
			tableptr2+=1
		while(1)				// exit is via break statement		
		tableptr1+=1
	while(tableptr1 < ratioNpnts)
//	display Ca_wave
//	ModifyGraph mode=3,marker=18
//*****************************************************************************************************************************
// set the time of the first point after flash as 0, make a new time wave start from 0 but with the same scale 
// convert the x value of the Ca wave according to the time wave
 		
	duplicate/o $Time_wave Time_wave1
	variable time_start
	time_start = Time_wave1[Flash_startpnt]
	Time_wave1-= time_start
	
	variable point_number
	point_number = numpnts(Time_wave1)-1
	duplicate/o/R=[Flash_startpnt, point_number] Time_wave1 Time_wave2
	duplicate/o/R=[Flash_startpnt, point_number] Ca_wave Ca_wave2
	variable tscale_start, tscale_end
	tscale_start = Time_wave2[0]
	tscale_end = Time_wave2[point_number]
	setscale x, tscale_start, tscale_end, Ca_wave2
	Ca_wave2 =  interp(x, Time_wave2, Ca_wave2)
	Smooth/B 3, Ca_wave2
//*****************************************************************************************************************************
// make a new Ca wave, every 0.5sec pick a Ca value
// insert two points at the beginning of the new ca wave, first point is the basal Ca, second point is the Ca before flash
	
	variable time2=0, tableptr0=0, Ca_point
	Ca_point = (tscale_end-tscale_start)/time_step
	make/o/n=(Ca_point) Avg_Ca

	do
		Avg_Ca[tableptr0]=Ca_wave2(time2)
		time2+=time_step
		tableptr0+=1
	while(time2<=tscale_end)

	Variable Basal_Ca, Preflash_Ca
	Basal_Ca = Ca_wave[0]						// this is the basal Ca concentration before any illumination
	Preflash_Ca = Ca_wave[Flash_startpnt-1]		// this is the Ca conc. after 350/380 adjustment but before flash photolization
	
	insertPoints 0,2, Avg_Ca
	Avg_Ca[0] = Basal_Ca
	Avg_Ca[1] = Preflash_Ca
	display Avg_Ca
	ModifyGraph mode=3,marker=18
	Edit $Ratio_wave Ca_wave Avg_Ca
	Print "------------------------------------------------------------------------------------------------------------------------------------------------"
end