source: sans/Dev/trunk/NCNR_SANS_Utilities/Single_Model_Template/My_Model_v40.ipf @ 740

#pragma rtGlobals=1             // Use modern global access method.
#pragma IgorVersion=6.0

////////////////////////////////////////////////////
//
// A template for writing your own scattering functions
//
// FIVE things to do:
// 1) identify the equation and the free variables before coding
//
// --then you may start coding ---
//
// 2) change coef_MyModel to have the coefficients (defined in w[])
//    in the exact order you have specified
//    - this change is one line, search for "CH#1"
// 3) change parameters_MyModel to has parameter names, that match
//    the coefficients - this line follows the coefficients of step 1
//          - this is "CH#2"
// 4) change the function MyModel() to calculate what you want
//    w[N] must contain all the fitting parameters - search for "CH#3"
// 5) repeat steps 2 and 3 (copy/paste) changing the RHS of the same
//    lines in the macro for the smeared model calculation
//    - search for "CH#4", and the following line
//
// - then, all you need to do is (re)plot the macro from the Macros menu
//
// As-is, the model calculates a form factor for a uniform sphere 
//
// v 2.00 -- 01 MAY 2006 SRK
// v 2.1 -- Jul 2007 SRK
//                              converted to use structures, AAO, like the rest of the package
// v 2.2 -- Nov 2009 SRK
//                              updated AddModelToStrings() for the new release
//
////////////////////////////////////////////////////

//this macro sets up all the necessary parameters and waves that are
//needed to calculate the model function.
//
Macro PlotMyModel(num,qmin,qmax)
        Variable num=200, qmin=0.001, qmax=0.7
        Prompt num "Enter number of data points for model: "
        Prompt qmin "Enter minimum q-value (^-1) for model: " 
        Prompt qmax "Enter maximum q-value (^-1) for model: "
//
        Make/O/D/n=(num) xwave_MyModel, ywave_MyModel
        xwave_MyModel =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
        Make/O/D coef_MyModel = {1.,60,1e-6,6.3e-6,0.01}                        //CH#1
        make/o/t parameters_MyModel = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"}       //CH#2
        Edit parameters_MyModel, coef_MyModel
        
        Variable/G root:g_MyModel=0
        root:g_MyModel := MyModel(coef_MyModel, ywave_MyModel, xwave_MyModel)
        Display ywave_MyModel vs xwave_MyModel
        ModifyGraph marker=29, msize=2, mode=4
        ModifyGraph log=1,grid=1,mirror=2
        Label bottom "q (\\S-1\\M) "
        Label left "I(q) (cm\\S-1\\M)"
        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
        
        AddModelToStrings("MyModel","coef_MyModel","parameters_MyModel","MyModel")
//
End

//
//this macro sets up all the necessary parameters and waves that are
//needed to calculate the  smeared model function.
//
//no input parameters are necessary, it MUST use the experimental q-values
// from the experimental data read in from an AVE/QSIG data file
////////////////////////////////////////////////////
// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
Macro PlotSmearedMyModel(str)                                                           
        String str
        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
        
        // if any of the resolution waves are missing => abort
        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
                Abort
        endif
        
        SetDataFolder $("root:"+str)
        
        // Setup parameter table for model function
        Make/O/D smear_coef_MyModel = {1.,60,1e-6,6.3e-6,0.01}          //CH#4 change this line and the following line
        make/o/t smear_parameters_MyModel = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"}
        Edit smear_parameters_MyModel,smear_coef_MyModel                                        //display parameters in a table
        
        // output smeared intensity wave, dimensions are identical to experimental QSIG values
        // make extra copy of experimental q-values for easy plotting
        Duplicate/O $(str+"_q") smeared_MyModel,smeared_qvals
        SetScale d,0,0,"1/cm",smeared_MyModel
                                        
        Variable/G gs_MyModel=0
        gs_MyModel := fSmearedMyModel(smear_coef_MyModel,smeared_MyModel,smeared_qvals) //this wrapper fills the STRUCT
        
        Display smeared_MyModel vs smeared_qvals
        ModifyGraph log=1,marker=29,msize=2,mode=4
        Label bottom "q (\\S-1\\M)"
        Label left "I(q) (cm\\S-1\\M)"
        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
        
        SetDataFolder root:
        AddModelToStrings("SmearedMyModel","smear_coef_MyModel","smear_parameters_MyModel","MyModel")
End


// nothing to change here
//
//AAO version, uses XOP if available
// simply calls the original single point calculation with
// a wave assignment (this will behave nicely if given point ranges)
Function MyModel(cw,yw,xw) : FitFunc
        Wave cw,yw,xw
        
#if exists("MyModelX")
        yw = MyModelX(cw,xw)
#else
        yw = fMyModel(cw,xw)
#endif
        return(0)
End


//CH#3
// you should write your function to calculate the intensity
// for a single q-value (that's the input parameter x)
// based on the wave (array) of parameters that you send it (w)
//
// unsmeared model calculation
//
Function fMyModel(w,x) : FitFunc
        Wave w
        Variable x
        
        // variables are:                                                       
        //[0] scale
        //[1] radius ()
        //[2] sld sphere (-2)
        //[3] sld solv
        //[4] background (cm-1)
        
        Variable scale,radius,delrho,bkg,sldSph,sldSolv,qval            
        scale = w[0]
        radius = w[1]
        sldSph = w[2]
        sldSolv = w[3]
        bkg = w[4]
        
        delrho = sldSph - sldSolv       
        
        qval = x                //rename the input q-value, purely for readability
        
        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
        // and is rescaled to give [=] cm^-1
        
        Variable bes,f,vol,f2
        //
        //handle q==0 separately to avoid divide-by-zero in the Bessel func
        If(qval==0)
                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
                return(f)
        Endif
        
        bes = 3*(sin(qval*radius)-qval*radius*cos(qval*radius))/qval^3/radius^3
        vol = 4*pi/3*radius^3
        f = vol*bes*delrho              // [=] 
        // normalize to single particle volume and convert to 1/cm
        f2 = f * f / vol * 1e8          // [=] 1/cm
        
        return (scale*f2+bkg)   // Scale, then add in the background
        
End

//CH#4  
///////////////////////////////////////////////////////////////
// smeared model calculation
//
// you don't need to do anything with this function, as long as
// your MyModel works correctly, you get the resolution-smeared
// version for free.
//
// this is all there is to the smeared model calculation!
Function SmearedMyModel(s) : FitFunc
        Struct ResSmearAAOStruct &s

//      the name of your unsmeared model (AAO) is the first argument
        Smear_Model_20(MyModel,s.coefW,s.xW,s.yW,s.resW)

        return(0)
End


///////////////////////////////////////////////////////////////


// nothing to change here
//
//wrapper to calculate the smeared model as an AAO-Struct
// fills the struct and calls the ususal function with the STRUCT parameter
//
// used only for the dependency, not for fitting
//
Function fSmearedMyModel(coefW,yW,xW)
        Wave coefW,yW,xW
        
        String str = getWavesDataFolder(yW,0)
        String DF="root:"+str+":"
        
        WAVE resW = $(DF+str+"_res")
        
        STRUCT ResSmearAAOStruct fs
        WAVE fs.coefW = coefW   
        WAVE fs.yW = yW
        WAVE fs.xW = xW
        WAVE fs.resW = resW
        
        Variable err
        err = SmearedMyModel(fs)
        
        return (0)
End