A fit function for QENS data.

Function fitUsingNLor(cw, yw, xw1) : FitFunc
    Wave cw, yw, xw1
    //An all at once fit function to fit QENS data.
    //function fits a linear background
    // with a delta function
    // and N lorentzians.
    // delta and lorentzians are centred at the same position (cw[2])
    // xw1 contains the x values
    // res contains the instrumental resolution function at the same x values.
    // the returned y wave consists of the sum of the delta and N lorentzian functions
    // convolved with the resolution function
 
        //the resolution function MUST be contained in the following wave
    Wave res = root:res
 
    //!!!!!!!!!!!!!!
    //ASSUMES:
        //1) THAT THE POINT SPACING IN XW1 IS LINEAR!!!!!!!!!!!!!!!!!
        //2) THE RESOLUTION WAVE MUST HAVE THE SAME NUMBER OF POINTS AS XW1 AND BE THE WAVE ROOT:RES
        //3) DO NOT TRY TO INTERPRET THE FUNCTION AT ANY OTHER POINTS THAN THE INPUT XW1 WAVE (IT DOESN'T MAKE SENSE).
    //!!!!!!!!!!!!!!
    
    //explanatory
    //cw[0] = a
    //cw[1] = b
    //cw[2] = PEAK POS
    //cw[3] = area - delta
    //cw[4] = FWHM - lor 0
    //cw[5] = AREA - lor 0
    //cw[2*n + 4] = FWHM - lor n
    //cw[2*n + 5]  = AREA - lor n
    
    variable timer
    variable numconts, ii, offset, area_res, conv_offset, rot_remaining, shift, thearea
    numconts = (numpnts(cw) - 4)/2
    yw = 0
    make/free/d/n=3 lor_cw = 0
    make/free/d/n=(numpnts(yw)) tempyw
    MarkPerfTestTime 0
    for(ii = 0 ; ii < numconts ; ii+=1)
        offset = 2*ii + 4
        lor_cw[0] = cw[2]
        lor_cw[1] = cw[offset]
        lor_cw[2] = cw[offset + 1]
        //multithread tempyw = 2*lor_cw[2]/pi * lor_cw[1]/(4*(xw1-lor_cw[0])^2 + lor_cw[1]^2)
        MPFXLorenzianPeak(lor_cw, tempyw, xw1)
        multithread yw += tempyw
    endfor
    
    MarkPerfTestTime 1
    //convolve and account for the scaling
    convolve/a res, yw
    
    MarkPerfTestTime 2
    
        fastop yw = (xw1[1] - xw1[0]) *yw
    theArea = areaXY(xw1, res)
    fastop yw =  (1/theArea) * yw
    
    MarkPerfTestTime 3
    //have to rotate because of the convolution
    //the rotate operation gets you most of the way there
    conv_offset = binarysearchinterp(xw1, 0) - trunc(numpnts(xw1) / 2)
    rotate -trunc(conv_offset), yw
    duplicate/free yw, tempyw 
    rot_remaining = conv_offset - trunc(conv_offset)
    yw = tempyw[p + rot_remaining]
    
    MarkPerfTestTime 4
    //add in scaling of the delta function
    duplicate/free res, tempyw
    shift = cw[2]/(xw1[1] - xw1[0])
    area_res = areaxy(W_energy, res)
    yw += tempyw[p-shift] * cw[3] / area_res
    
    //add in the background
    multithread yw += cw[0] + cw[1] * xw1
    MarkPerfTestTime 5
End