Fit to complex-valued functions

Function DoComplexFitMC(ndat,a,b,sigma)          // driver to demonstrate fit to complex func
    variable ndat,a,b,sigma                 // uses Monte Carlo (MC) to create fake data
 
    make /c/d/free/n=(ndat) yw
    make /d/free/n=(2*ndat) sigw
    make /d/o/n=(2*ndat) ww2=0,yw2          // real waves for real and imag parts of yw, ww
    variable wmax=30
    ww2[0,ndat-1] = p/(ndat-1)*wmax         // linear span of frequencies  (often would be log...)
    ww2[ndat, ] = ww2[p-ndat]               // these values are not used (dummy placeholders)
    
    sigw = sigma                                    // wave of standard deviations
    make /d/free pars={a,b}                     // pack true parameters into wave, use for fit, too
    yw = ComplexFunc(pars,ww2)                  // noiseless "core"  for fit using true a,b
    yw2[0,ndat-1] = real(yw)
    yw2[ndat, ] = imag(yw[p-ndat])              // yw20 concatenates real & imaginary responses
    yw2 += gnoise(sigma)                            // add noise for fake data (to be fit)
    
    FuncFit/NTHR=0/W=2 ComplexFitAllAtOnce pars  yw2 /X=ww2 /W=sigw /I=1 /D       // does the fit
    
    //_______________   graphics stuff__________________________
    make /d/o/n=(ndat) yr,yi,w0                 //  waves for nice plotting
    w0 = ww2[p]                                 // get first ndat points
    yr = yw2[p];        yi = yw2[p+ndat]            // get real and imag parts
    
    make /c/d/free/n=1000 fit_yc
    make /d/o/n=1000 fit_yr, fit_yi
    SetScale/I x 0,wmax,"", fit_yi,fit_yr,fit_yc
    fit_yc = ComplexFunc(pars, x)
    fit_yr=real(fit_yc);        fit_yi = imag(fit_yc)
End

Function ComplexFitAllAtOnce(pars,y2,w2) : FitFunc  // kluge for complex fits using FuncFit
    wave pars,y2,w2
    variable ndat = numpnts(w2)/2                       // break between real and imag data
    make /c/d/free /n=(ndat) temp = ComplexFunc(pars,w2)                // call to complex fit func
    y2[0,ndat-1] = real(temp);  y2[ndat, ] = imag(temp[p-ndat])     // concatenates real & imag
End
 
Function /c ComplexFunc(pars,w)
    wave pars;  variable w
    variable a=pars[0], b=pars[1]                       // unpack parameter wave
    variable /c i=cmplx(0,1)                            // i = sqrt(-1)
    return exp(i*w*a) / (1-i*w*b)                   // this is the function we fit!
end