#pragma rtGlobals=3		// Use modern global access method and strict wave access.

Function/S DoubleSFG_PeakFuncInfo(InfoDesired)
	Variable InfoDesired
	
	String info=""
	
	Switch (InfoDEsired)
		case PeakFuncInfo_ParamNames:
			info = "Center1;Amplitude1;Width1;Center2;Amplitude2;Width2" // Center1 = peak center, Amp1 = peak's amplitude,  
			break;							//Width1 = the peak's FWHM which is also refered to as "Gama". 2's refer to the second peak. Name the parameters as you like
		case PeakFuncInfo_PeakFName:        // BUT you must maintain their order of apperance in the parameter w[] wave!
			info = "DoubleSFG_Peak"
			break;
		case PeakFuncInfo_GaussConvFName:
			info = "SimpleLorentzianGuess"
			break;
		case PeakFuncInfo_ParameterFunc:
			info = "DoubleSFG_PeakParams"
			break;
		case PeakFuncInfo_DerivedParamNames:
			info = "PeakCenter1;Amplitude1;Area1;PeakWidth1;PeakCenter2;Amplitude2;Area2;PeakWidth2" // PeakWidth = FWHM = gama...
			break;
		default:
			break;
	endSwitch

	return info
end

Function SimpleLorentzianGuess(w) // Let's assume that a simple Lorentzian is a good guess to start with...
	Wave w
	//Note that the wave is already pre-filled with parameters from the automatic or "by-hand" guess from the side of MPF
	//assuming a possibly asymmetric Gaussian.
	//Parameters w[0],w[1],w[2], are corresponding to position, width, and height respectively,
	
	Redimension/N=6 w
	// w[0] = Position1==>Center1				w[3] = Center2
	// w[1] = Width1==>Width1				w[4] = Amplitude2
	// w[2] = Height1==>Amplitude1			w[5] = Width2
	
	// 1st peak: width of Lorenzian needs to be modified from the estimated Gaussian width
	w[2] = w[1]*w[2]*sqrt(pi)
	
	// 2nd peak: width of Lorenzian needs to be modified from the estimated Gaussian width
	w[3] = w[4]*w[5]*sqrt(pi)
	return 0
	
	
	
end



Function DoubleSFG_Peak(w, yw, xw)  // The w wave is a 6 object wave that has the fit parameters. Use w[0] to w[6], i.e.;
									 // w[0] = Center1
									 // w[1] = Amp1
								       // w[2] = Width1
									 // w[3] = Center2
									 // w[4] = Amp2
									 // w[5] = Width2
	
	Wave w, yw, xw
	
	Variable/C sqrtJ = cmplx(0,1)
	yw = ((abs(w[1]/((xw-w[0])+w[2]*sqrtJ) + w[4]/(xw-w[3])+w[5]*sqrtJ)))^2
end

Function DoubleSFG_PeakParams(cw, sw, outWave) //cw is the wave parameter from the fit, 
	Wave cw, sw, outWave						 //sw is the standard deviations for the parameters - both are 8 wave parameters as defined in "PeakFuncInfo_DerivedParamNames"
	
	
	//Location (central position) 1st peak - treated as a simple Lorenztian
	outWave[0][0] = cw[0]			// This is actually the peak center of a Lorentian.
	outWave[0][1] = sqrt(sw[0][0])  // The error is treated accordingly.
	
	//Amplitude 1st peak - treated as a simple Lorenztian
	outWave[1][0] = cw[1]
	outWave[2][1] =  NaN // sqrt( (outWave[1][0]^2)*((sw[2][2]/cw[2]^2) + (sw[1][1]/cw[1]^2) - 2*sw[1][2]/(cw[1]*cw[2])) )
	
	//Area 1st peak = Intensity or oscillator strength
	outWave[2][0] =  cw[1]/sqrt(cw[2])	// Amplitude/sqrt(Gama)
	outWave[2][1] = NaN
	
	//FWHM 1st peak - this is "Gama", i.e., the width of the peak...
	outWave[3][0] = cw[2]
	outWave[3][1] = NaN
	
	//***********************************Second Peak*************************************//
	
	//Location 2nd peak
	outWave[4][0] = cw[4]
	outWave[4][1] = sqrt(sw[4][4])
	
	//Amplitude 2nd peak
	outWave[5][0] = cw[5] // 2*cw[5]/(pi*cw[5])
	outWave[5][1] = NaN //sqrt( (outWave[5][0]^2)*((sw[5][5]/cw[5]^2) + (sw[4][4]/cw[4]^2) - 2*sw[4][5]/(cw[4]*cw[5])) )
	
	//Area 2nd peak
	outWave[6][0] = cw[5]/sqrt(cw[6])
	outWave[6][1] = NaN
	
	//FWHM 2nd peak 
	outWave[3][0] = cw[6]
	outWave[3][1] = NaN // sqrt(sw[4][4])
	
end