Custom peak function in MPF2 - what happened to the fit wave?

Function/S DS_PeakFuncInfo(InfoDesired)
    Variable InfoDesired

    String info=""
       
    Switch (InfoDesired)
        case PeakFuncInfo_ParamNames:
            info = "Location;width;asymmetry;Height;"
            break;
        case PeakFuncInfo_PeakFName:
            info = "MPFDS"
            break;
        case PeakFuncInfo_GaussConvFName:
            info = "GausstoDSGuess"
            break;
        case PeakFuncInfo_ParameterFunc:
            info = "DSParams"
            break;
        case PeakFuncInfo_DerivedParamNames:
            info = "Location;Height;FWHM;Asymmetry;"
            break;
        default:
            break;
    endSwitch
   
    return info
end

Function MPFDS(w,yw,xw)
  wave w
  wave yw, xw

  yw = w[3]*cos(Pi*w[2]/2+(1-w[2])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[2])/2)
end

Function GaussToDSGuess(w)
  wave w
  variable height=w[2]
  Redimension/N=4 w
  w[2]=0.05 //just put in a small asymmetry to start the fit with.
  w[3]=height
  return 0
end

Function DSParams(cw, sw, outWave)
    Wave cw, sw, outWave
    Variable width = cw[1]
    Variable widthSigma = sqrt(sw[1][1])
    Variable asym=cw[2]
    Variable asymsigma=sqrt(sw[2][2])
    Variable amp = cw[3]
    Variable ampSigma = sqrt(sw[3][3])
   
    // location
    outWave[0][0] = cw[0]
    outWave[0][1] = sqrt(sw[0][0])
       
    // amplitude
    outWave[1][0] = amp
    outWave[1][1] = ampSigma
   
    // area    ..... haven't what to return here, the integral of DS diverges at large values of E, so the area is a bit meaningless.
   
    outWave[2][0] = nan
    outWave[2][1] = nan
   
    // FWHM
    outWave[3][0] = width*2*sqrt(ln(2))
    outWave[3][1] = widthSigma*2*sqrt(ln(2))
end

Function MPFDS(w,yw,xw)

yw = w[3]*cos(Pi*w[2]/2+(1-w[2])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[2])/2) //this is wrong
yw = w[3]*cos(Pi*w[2]/2+(1-w[2])*atan((xw-w[0])/w[1]))/((xw-w[0])^2+w[1]^2)^((1-w[2])/2) //this is right!!

...
Variable Resolution = 10                 // the sample resolution (higher = slower)
Variable Xdelta     = abs(xw[0] - xw[NumPnts(xw)-1]) / (Resolution * NumPnts(xw))    // calculate an high-res xdelta for scaling
Variable yPoints    = Resolution * NumPnts(yw)   // make the wave X times bigger than the original y wave
Make/FREE/D/N=(2*yPoints+1) yWave            // create a wave two times as big to have space before and after the peak
SetScale/P x xw[0]-yPoints/2*Xdelta, Xdelta, yWave   // reserve 1/2 size of space before the x wave starts (1/2 after)

yWave = w[2]*cos(Pi*w[3]/2+(1-w[3])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[3])/2)

if (w[4] != 0)                       // convolve with gaussian if width is set
    Variable width  = w[4]/Xdelta            // width normalized by x scaling
    yPoints         = 10*width       // create enough points to fit in the full gaussian
    Make/FREE/D/N=(yPoints + 1) GaussWave        // create a gaussian for the convolution
    GaussWave = Gauss(x,yPoints/2,width)
    Convolve/A GaussWave yWave           // convolve with the prepared gaussian
endif
yw = yWave(xw[p])                    // write the calculated values into the resultwave
...