A custom fitting equation that also calls a fit

#pragma rtGlobals=1 // Use modern global access method.
 
//This function works as a basic fit function for Diffusion's dependence on oscillating amplitude 
Function DiffusionVsOscAmplitude(w,OscAmp) : FitFunc
    Wave w
    Variable OscAmp
 
    //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
    //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
    //CurveFitDialog/ Equation:
    //CurveFitDialog/ This runs a simulation to the differential equation determine the fit.
    //CurveFitDialog/ End of Equation
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ OscAmp
    //CurveFitDialog/ Coefficients 2
    //CurveFitDialog/ w[0] = Ds
    //CurveFitDialog/ w[1] = invtm
print "Starting fitfunc"
 
    //Set global variables
    variable /G position_cm = 200e-7
    variable /G SimLength_cm =400e-7
    Variable /G dx_cm = 2e-7
    variable /G dt_s = .0002
    variable /G SimTime_s = 3
    
    //Make the sensitive volume
    SenSliceFull(OscAmp)
    print OscAmp
    print w[0]
    print w[1]
print "Calling runSim"  
    return RunSimulation(w[0],w[1])
End
 
 
//To run the simulation 
Function RunSimulation(Ds, invtm)
    Variable Ds     //diffusion constant
    Variable invtm  //Relaxation rate for T_1 or \tau_m
print "In RunSim"
    variable /G position_cm
    variable /G SimLength_cm
    Variable /G dx_cm
    variable /G dt_s
    variable /G SimTime_s
    Wave ExpFits
    Wave SenSlice
    Wave Spins
    
    variable Iterations = round(SimTime_s/dt_s) //Number of spatial points to simulate
    variable Positions = round(SimLength_cm/dx_cm)  //Number of time steps to take
    variable frac = cosh(dt_s*invtm)*exp(-1*dt_s*invtm)  //the fraction of spins that flip in time dt.
    variable i  //Will hold the time step
    variable j //will hold the position step
    variable k //used to update plots every now and then
    variable sum  //will hold the integration of Spin*SenSlice
    variable result  //The resultant measured inverse corrolation time of the force
    
    Make /O /N=(Positions) TempSpins  //This will hold the spin state at the next time step temporarly
    Make /O /N=(Iterations) OrigSpinState  //This will hold how much of the original spin state is left convolved with the sen slice vs time (for each dt)
    SetScale/P x 0,(dt_s),"", OrigSpinState  //each point repersents a single time step
    Make /O /N=2 TempCoef  //Since we will have a fit inside a fit we need to do house cleaning so the two W-coef waves do not screw with each other
    
    OrigSpinState = nan  //just so I don't have lingering values from a previous run
    for(i=0;i<Iterations;i=i+1)
        sum = 0
        for(j=0;j<Positions; j=j+1)
            sum = sum + SenSlice[j]*Spins[j]  //Integrate to find "force"
            
            //Spins diffusing
            if(j==0)
                TempSpins[j] = Spins[j] + Ds*dt_s*(Spins[j+1]-Spins[j])/(dx_cm*dx_cm)   //spins are blocked from diffusing past on the left side
                //TempSpins[j] = 0                                                  //Spins diffuse past the left side (and disapear)   
            elseif(j==Positions-1)
                TempSpins[j] = Spins[j] + Ds*dt_s*(Spins[j-1]-Spins[j])/(dx_cm*dx_cm)  //spins are blocked from diffusing past on the right side
                //TempSpins[j] = 0                                                  //Spins diffuse past the right side (and disapear)
            else
                TempSpins[j] = Spins[j] + Ds*dt_s*(Spins[j-1]+Spins[j+1]-2*Spins[j])/(dx_cm*dx_cm)  
            endif
        endfor
        
        OrigSpinState[i] = sum  //The "force" at this time period
        
        for(j=0;j<Positions; j=j+1)
            Spins[j] = frac*TempSpins[j]  //Get ready for next iteration
        endfor
        if(k<1)
            doUpdate
            k=Iterations/100
        endif
        k=k-1
    endfor
    print "Calling fit"
    CurveFit/NTHR=0 exp kwCWave=ExpFits,  OrigSpinState /D 
    print ExpFits[2]
    print "Returning"
    return ExpFits[2]
end
 
Function SenSliceFull(OscAmp)
    variable OscAmp
    variable /G position_cm
    variable /G SimLength_cm
    Variable /G dx_cm
    
    variable i
    
    Make /O /N=(round(SimLength_cm/dx_cm)) SenSlice
    Make /O /N=(round(SimLength_cm/dx_cm)) Spins
    for(i=0;i<round(SimLength_cm/dx_cm);i=i+1)
        SenSlice[i]=0
        Spins[i]=0
    endfor
    for(i=ceil(-1*(OscAmp/dx_cm));i<=floor(OscAmp/dx_cm);i=i+1)
        SenSlice[round(position_cm/dx_cm)+i]=Cos((i*dx_cm/OscAmp)*(pi/2))
        Spins[round(position_cm/dx_cm)+i]=1
    endfor
    SetScale/P x 0,(dx_cm),"", SenSlice
    SetScale/P x 0,(dx_cm),"", Spins
end

threadsafe Function DoMyFit(cwave, datawave)
    Wave cwave, datawave
    
    CurveFit/NTHR=1 exp kwCWave=cwave,  datawave
end
 
Function runsimulation(Ds, invtm)
    Variable Ds     //diffusion constant
    Variable invtm  //Relaxation rate for T_1 or \tau_m
    
    //*** lots of code here****
    
    Variable threadid = ThreadGroupCreate(1)        // just one thread to encapsulate the inner fit
    ThreadStart threadid, 0, DoMyFit(ExpFits,  OrigSpinState)
    do
        variable tgs= ThreadGroupWait(threadid,5000)        // five-second wait is excessive, but we don't want to spin the loop.
                                                // We just want to wait until the inner fit is done, which could be fairly lengthy.
    while( tgs != 0 )
 
    return ExpFits[2]
end