Normalized Derivative Peak Functions

#pragma rtGlobals=1     // Use modern global access method.
 
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Normalized peak functions (Gauss, Lorentz, pseudo-Voigt) in absorption and dispersion mode.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  by W. Harneit (Berlin), November 2009
//      August 2010: Updated comments, included note on multi-peak fitting
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Each function takes three parameters (in wave w): intensity, position, and width.
//
//  For the ...Amp functions, the intensity parameter means amplitude (from zero). 
//      The area is not normalized in this case.
//  For the ...Area functions, intensity means integrated area (or double integral for Diff...).
//      The intensity parameter does not correspond to the peak amplitude in this case.
//
//  For the absorption functions (without Diff...), width means full width at half maximum (FWHM).
//  The Diff... functions are the derivatives (dispersion mode). Width is peak-to-peak width (p.p. width).
//      p.p. width (dispersion) and FWHM (absorption) are linearly related, depending on the peak form.
//
//  The Pseudo-Voigt functions are a simple linear interpolation between Gauss and Lorentz form, with 
//  _equal_ width. They take a fourth parameter 0 <= w[3] <= 1, which is the weight of the Lorentzian 
//  (the weight of the Gaussian is 1 - w[3]). Strange things will happen if w[3] > 1.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Examples: 
//
//  GaussAmp is the absorptive Gaussian, with amplitude w[0], position w[1], fwhm w[2].
//
//  DiffLorArea is the dispersive Lorentzian, with double integral w[0], position w[1], p.p. width w[2].
// 
//  To do multi-peak fitting, follow the steps below. 
//  (see also DisplayHelpTopic "Fitting Sums of Fit Functions")
//
//      We assume that a two-line absorption spectrum is present in spec_INT. The position parameters
//      are the most sensitive parameters, followed by the width. Intensities (here: areas) may be way off.
//      Tip: You can use "Decommentize" (Edit menu) on the following lines.
//
//      // Step 1: provide guesses for the two peaks
//      Make/D/O low={4e-10,0.33707,2e-5}, high={4e-10,0.33818,2e-5, 0.77}
//      // Step 2: provide hold strings (if any)
//      String hstr = "0001"    // to hold the shape parameter p=0.77 for the "high" peak
//      // Step 3: do the fit
//      FuncFit/M=0/L=4096/W=0/X {{LorArea,low}, {PseudoVoigtArea,high,hold=hstr}}, spec_INT /D
//      // /M=0 => no cov. matrix; /L= numpnts; /W=0 => don't wait; /X => full X range; /D => auto dest.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  code note:  Numerical constants are given in the comments with 5 digits precision. 
//              If speed is an issue, one could use global constants "constant kln2 = 0.69315" etc.
//
function GaussAmp(w, x)
wave w
variable x
    variable xi = 2 * (x - w[1]) / w[2], ln2 = ln(2)    // ln2 = 0.69315
    return w[0] * exp(-ln2*xi*xi)
end
 
function GaussArea(w, x)
wave w
variable x
    return GaussAmp(w, x) / ( sqrt( pi / ln(2) ) / 2 * w[2] )   // divisor = 1.0645 * w[2]
end
 
function LorAmp(w, x)
wave w
variable x
    variable xi = 2 * (x - w[1]) / w[2]
    return w[0] / (1 + xi*xi)
end
 
function LorArea(w, x)
wave w
variable x
    return LorAmp(w, x) / ( pi / 2 * w[2] ) // divisor = 1.5708 * w[2]
end
 
function DiffGaussAmp(w, x)
wave w
variable x
    variable xi = 2 * (x - w[1]) / w[2], ln2 = ln(2)    // ln2 = 0.69315
    return -xi * w[0] * exp(-0.5*(xi*xi-1))
end
 
function DiffGaussArea(w, x)
wave w
variable x
    return DiffGaussAmp(w, x) / ( sqrt( pi * e / 8 ) * w[2] * w[2] )    // divisor = 1.0332 * w[2]^2
end
 
function DiffLorAmp(w, x)
wave w
variable x
    variable xi = 2 * (x - w[1]) / w[2]
    return -xi * w[0] * ( 4 / (3 + xi*xi) )^2
end
 
function DiffLorArea(w, x)
wave w
variable x
    return DiffLorAmp(w, x) / ( sqrt( 4 / 3 ) * w[2] * w[2] )   // divisor = 1.1547 * w[2]^2
end
 
function PseudoVoigtAmp(w, x)
wave w
variable x
    return w[3] * LorAmp(w, x) + (1 - w[3]) * GaussAmp(w, x)
end
 
function PseudoVoigtArea(w, x)
wave w
variable x
    return w[3] * LorArea(w, x) + (1 - w[3]) * GaussArea(w, x)
end
 
function DiffPseudoVoigtAmp(w, x)
wave w
variable x
    return w[3] * DiffLorAmp(w, x) + (1 - w[3]) * DiffGaussAmp(w, x)
end
 
function DiffPseudoVoigtArea(w, x)
wave w
variable x
    return w[3] * DiffLorArea(w, x) + (1 - w[3]) * DiffGaussArea(w, x)
end