source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/VSANS/V_GaussSpheres_WB.ipf @ 1095

#pragma TextEncoding = "MacRoman"
#pragma rtGlobals=3             // Use modern global access method and strict wave access.


// turn the White Beam resolution smearing into a fitting function
// so that the wavelength smeared function can then be smeared by
// a gaussian resolution function that has geometry only.
//
// The geometry only resolution is generated by passing dl/l=0 to the resolution
// calculation.
//
// This representation uses the "middle" of the distribution

//
//

//#include "sphere_v40"
// plots the form factor of  spheres with a Gaussian radius distribution
//
// also can plot the distribution itself, based on the current model parameters
//
// integration is currently done using 20-pt quadrature, but may benefit from 
//switching to an adaptive integration.
//

Proc PlotGaussSpheresWB(num,qmin,qmax)
        Variable num=128,qmin=0.001,qmax=0.7
        Prompt num "Enter number of data points for model: "
        Prompt qmin "Enter minimum q-value (A^-1) for model: "
        Prompt qmax "Enter maximum q-value (A^-1) for model: "
        
        Make/O/D/N=(num) xwave_pgsWB,ywave_pgsWB
        xwave_pgsWB = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
        Make/O/D coef_pgsWB = {0.01,60,0.2,1e-6,3e-6,0.001}
        make/O/T parameters_pgsWB = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
        Edit parameters_pgsWB,coef_pgsWB
        
        Variable/G root:g_pgsWB
        g_pgsWB := GaussSpheresWB(coef_pgsWB,ywave_pgsWB,xwave_pgsWB)
        Display ywave_pgsWB vs xwave_pgsWB
        ModifyGraph log=1,marker=29,msize=2,mode=4
        Label bottom "q (A\\S-1\\M)"
        Label left "Intensity (cm\\S-1\\M)"
        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
        
        AddModelToStrings("GaussSpheresWB","coef_pgsWB","parameters_pgsWB","pgsWB")
End

// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
Proc PlotSmearedGaussSpheresWB(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_pgsWB = {0.01,60,0.2,1e-6,3e-6,0.001}                                       
        make/o/t smear_parameters_pgsWB = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}   
        Edit smear_parameters_pgsWB,smear_coef_pgsWB                                    
        
        // 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_pgsWB,smeared_qvals                             
        SetScale d,0,0,"1/cm",smeared_pgsWB                                                     
                                        
        Variable/G gs_pgsWB=0
        gs_pgsWB := fSmearedGaussSpheresWB(smear_coef_pgsWB,smeared_pgsWB,smeared_qvals)        //this wrapper fills the STRUCT
        
        Display smeared_pgsWB vs smeared_qvals                                                                  
        ModifyGraph log=1,marker=29,msize=2,mode=4
        Label bottom "q (A\\S-1\\M)"
        Label left "Intensity (cm\\S-1\\M)"
        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
        
        SetDataFolder root:
        AddModelToStrings("SmearedGaussSpheresWB","smear_coef_pgsWB","smear_parameters_pgsWB","pgsWB")
End
        



//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 GaussSpheresWB(cw,yw,xw) : FitFunc
        Wave cw,yw,xw
        
#if exists("GaussSpheresX")
//      MultiThread yw = GaussSpheresX(cw,xw)
        yw = V_fGaussSpheresWB(cw,xw)

#else
//      yw = fGaussSpheresWB(cw,xw)
        yw = 1
#endif
        return(0)
End

Function V_fGaussSpheresWB(w,xx) : FitFunc
        wave w
        variable xx
        
        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho,inten,loLim,upLim
        
        //the coefficient values
//      scale=w[0]
//      rad=w[1]
//      pd=w[2]
//      sig=pd*rad
//      rho=w[3]
//      rhos=w[4]
//      delrho=rho-rhos
//      bkg=w[5]
        
        
        // define limits based on lo/mean, hi/mean of the wavelength distribution
        // using the empirical definition, "middle" of the peaks
        loLim = 3.37/5.3
        upLim = 8.37/5.3
        
        inten = V_IntegrGaussSphereWB_mid(w,loLim,upLim,xx)

// why do I need this? Is this because this is defined as the mean of the distribution
//  and is needed to normalize the integral? verify this on paper.      
        inten *= 5.3

// normalize the integral       
        inten /= 19933          // "middle"  of peaks

// additional normalization???
        inten /= 1.05           // 
        
        Return(inten)
End


// the trick here is that declaring the last qVal wave as a variable
// since this is implicitly called N times in the wave assignment of the answer wave
Function V_IntegrGaussSphereWB_mid(cw,loLim,upLim,qVal)
        Wave cw
        Variable loLim,upLim
        Variable qVal
        
        Variable/G root:qq = qval
        Variable ans
        
//      ans = Integrate1D(V_intgrnd_top,lolim,uplim,2,0,cw)             //adaptive quadrature
        ans = Integrate1D(V_integrand_pgsWB,lolim,uplim,1,0,cw)         // Romberg integration
        
        return ans
end

Function V_integrand_pgsWB(cw,dum)
        Wave cw
        Variable dum            // the dummy of the integration

        Variable val
        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
//      SVAR funcStr = root:gFunctionString
//      FUNCREF SANSModel_proto func = $funcStr

        val = V_WhiteBeamDist_mid(dum*5.3)*GaussSpheresX(cw,qq/dum)
        
        return (val)
End

//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 fSmearedGaussSpheresWB(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 = SmearedGaussSpheresWB(fs)
        
        return (0)
End

// this is all there is to the smeared calculation!
Function SmearedGaussSpheresWB(s) :FitFunc
        Struct ResSmearAAOStruct &s

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

        return(0)
End