Making a fitfunction multithreaded

Function myfitfunc_wrapper(w, yy, xx) : fitfunc
Wave w, yy, xx  
    Multithread yy = myfitfunc_kernel(w, xx)
End
 
Threadsafe function myfitfunc_kernel(w, xx)
Wave w
variable xx
    //your code goes here
    return w[0] + sqrt(xx) * w[1] * sin(xx) * cos(x) * ln(x)
End

Function myfitfunc_AAO_wrapper(w, yy, xx) : fitfunc
Wave w, yy, xx
 
    variable nt = ThreadProcessorCount
    variable tgID = ThreadGroupCreate(nt)
    variable ii, startP, perT, remainingP
 
    perT = floor(dimsize(yy, 0) / nt)
    remainingP = dimsize(yy, 0)
 
    make/n=(nt)/free/wave theDataY, theDataX
    for(ii = 0 ; ii < nt ; ii+=1)
        if(ii < nt -1)
            theDataY[ii] = newfreewave(0x04, perT)
            theDataX[ii] = newfreewave(0x04, perT)
        else
            theDataY[ii] = newfreewave(0x04, remainingP)
            theDataX[ii] = newfreewave(0x04, remainingP)        
        endif
        Wave yyy = theDataY[ii]
        Wave xxx = theDataX[ii]
 
        yyy = yy[p+startP]
                xxx = xx[p+startP]
 
        startP += perT
        remainingP -= perT
        ThreadStart tgID, ii, myfitfunc_AAO_kernel(w, theDataY[ii], theDataX[ii])
        WaveClear yyy, xxx
    endfor
    variable finished = ThreadGroupWait(tgID, inf)
    finished = threadgrouprelease(tgID)
    
    //now copy back into structure wave 
    startP = 0
    perT = floor(dimsize(yy, 0) / nt)
    remainingP = dimsize(yy, 0)
    for(ii = 0 ; ii < nt ; ii+=1)
        Wave yyy = theDataY[ii]
        yy[startP, startP+numpnts(yyy)] = yyy[p-startP]
        Waveclear yyy
 
        startP += perT
        remainingP -= perT
    endfor
End
 
Threadsafe function myfitfunc_AAO_kernel(w, yy, xx)
Wave w, yy, xx
    //your code goes here
    yy = w[0] + w[1] * sqrt(xx) * sin(xx) * cos(x) * ln(x)
End

Structure fitfuncStruct   
    Wave w
    wave y
    wave x[50]
 
    int16 numVarMD
    wave ffsWaves[50]
    wave ffsTextWaves[10]
    variable ffsvar[5]
    string ffsstr[5]
    nvar ffsnvars[5]
    svar ffssvars[5]
    funcref allatoncefitfunction ffsfuncrefs[10]
    uint32 ffsversion    // Structure version. 
    EndStructure 
 
Threadsafe Function allatoncefitfunction(w,y,x)
    Wave w,y,x
    //you need this to act as a signature for the fitfuncStruct
    //it needs to be threadsafe
End
 
Threadsafe function myfitfunc_AAO_kernel(w, yy, xx)
Wave w, yy, xx
    //your code goes here
    yy = w[0] + w[1] * sqrt(xx) * sin(xx) * cos(x) * ln(x)
End
 
Function myfitfunc_AAO_wrapper(s) : fitfunc
    Struct fitfuncStruct &s
 
    //work out how many points per thread
    variable nt =  ThreadProcessorCount
    variable tgID = ThreadGroupCreate(nt)
    variable ii, startP, perT, remainingP, finished
    perT = floor(dimsize(s.y, 0) / nt)
    remainingP = dimsize(s.y, 0)
        
    //create waveref waves to hold references to the y and x waves to be distributed amongst each thread.
    //you could create more than y and x if your function took more arguments
    make/n=(nt)/free/wave theDataY, theDataX
    
    //create the coefficients wave to be sent to all the threads, as you can't pass in s.w to a thread
    Wave theCoefs = s.w
 
    //this is the function we are going to pass the coefficients, and y and x waves onto. 
     //-It doesn't have to be an allatoncefitfunction it can be any function type you want (so long as it doesn't take structures).
    //- It needs to be threadsafe!
    //- the function can take as many parameters as you want, but you'll have to create a free wave for each of them (like theCoefs) to 
    //    be passed into the thread.
    Funcref allatoncefitfunction theFitFunction = s.ffsfuncrefs[0]
    
    for(ii = 0 ; ii < nt ; ii+=1)
        if(ii < nt -1)
            theDataY[ii] = newfreewave(0x04, perT)
            theDataX[ii] = newfreewave(0x04, perT)
        else
            theDataY[ii] = newfreewave(0x04, remainingP)
            theDataX[ii] = newfreewave(0x04, remainingP)        
        endif
        Wave yyy = theDataY[ii]
        Wave xxx = theDataX[ii]
        
        yyy = s.y[p+startP]
        xxx = s.x[0][p+startP]
        
        startP += perT
        remainingP -= perT
 
        //NOTE that you can't pass in a structure parameter in a threadstart operation. See:
        //DisplayHelptopic "Threadstart"
        // for more information
        ThreadStart tgID, ii, myfitfunc_AAO_kernel(theCoefs, yyy, xxx)
        WaveClear yyy, xxx
    endfor
     finished = ThreadGroupWait(tgID, inf)
    finished = threadgrouprelease(tgID)
 
    //now copy back into structure wave 
    startP = 0
    perT = floor(dimsize(s.y, 0) / nt)
    remainingP = dimsize(s.y, 0)
    for(ii = 0 ; ii < nt ; ii+=1)
        Wave yyy = theDataY[ii]
        s.y[startP, startP+numpnts(yyy)] = yyy[p-startP]
        Waveclear yyy
 
        startP += perT
        remainingP -= perT
    endfor
End