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