# source:sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/ProlateCoreShell_v40.ipf@633

Last change on this file since 633 was 633, checked in by srkline, 13 years ago

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.

File size: 7.9 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4////////////////////////////////////////////////
5// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
6// Adopting these into the experiment will insure that they are always present
7////////////////////////////////////////////////
8// this function is for the form factor of a prolate ellipsoid
9//
10// 06 NOV 98 SRK
11////////////////////////////////////////////////
12
13Proc PlotProlateForm(num,qmin,qmax)
14        Variable num=128,qmin=0.001,qmax=0.7
15        Prompt num "Enter number of data points for model: "
16        Prompt qmin "Enter minimum q-value (A^-1) for model: "
17        Prompt qmax "Enter maximum q-value (A^-1) for model: "
18
19        Make/O/D/n=(num) xwave_pef,ywave_pef
20        xwave_pef =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
21        Make/O/D coef_pef = {1.,100,50,110,60,1e-6,2e-6,6.3e-6,0.001}
22        make/o/t parameters_pef = {"scale","major core radius (A)","minor core radius (A)","major shell radius (A)","minor shell radius (A)","SLD core (A-2)","SLD shell (A-2)","SLD solvent (A-2)","bkg (cm-1)"}
23        Edit parameters_pef,coef_pef
24        Variable/G root:g_pef
25        g_pef := ProlateForm(coef_pef,ywave_pef,xwave_pef)
26//      ywave_pef := ProlateForm(coef_pef,xwave_pef)
27        Display ywave_pef vs xwave_pef
28        ModifyGraph log=1,marker=29,msize=2,mode=4
29        Label bottom "q (A\\S-1\\M)"
30        Label left "Intensity (cm\\S-1\\M)"
31        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
32
33        AddModelToStrings("ProlateForm","coef_pef","parameters_pef","pef")
34End
35
36///////////////////////////////////////////////////////////
37// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
38Proc PlotSmearedProlateForm(str)
39        String str
40        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
41
42        // if any of the resolution waves are missing => abort
43        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
44                Abort
45        endif
46
47        SetDataFolder \$("root:"+str)
48
49        // Setup parameter table for model function
50        Make/O/D smear_coef_pef = {1.,100,50,110,60,1e-6,2e-6,6.3e-6,0.001}
51        make/o/t smear_parameters_pef = {"scale","major core radius (A)","minor core radius (A)","major shell radius (A)","minor shell radius (A)","SLD core (A-2)","SLD shell (A-2)","SLD solvent (A-2)","bkg (cm-1)"}
52        Edit smear_parameters_pef,smear_coef_pef
53
54        // output smeared intensity wave, dimensions are identical to experimental QSIG values
55        // make extra copy of experimental q-values for easy plotting
56        Duplicate/O \$(str+"_q") smeared_pef,smeared_qvals                               //**** mod
57        SetScale d,0,0,"1/cm",smeared_pef
58
59        Variable/G gs_pef=0
60        gs_pef := fSmearedProlateForm(smear_coef_pef,smeared_pef,smeared_qvals) //this wrapper fills the STRUCT
61
62        Display smeared_pef vs smeared_qvals
63        ModifyGraph log=1,marker=29,msize=2,mode=4
64        Label bottom "q (A\\S-1\\M)"
65        Label left "Intensity (cm\\S-1\\M)"
66        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
67
68        SetDataFolder root:
69        AddModelToStrings("SmearedProlateForm","smear_coef_pef","smear_parameters_pef","pef")
70End
71
72//AAO version
73Function ProlateForm(cw,yw,xw) : FitFunc
74        Wave cw,yw,xw
75
76#if exists("ProlateFormX")
77        MultiThread yw = ProlateFormX(cw,xw)
78#else
79        yw = fProlateForm(cw,xw)
80#endif
81        return(0)
82End
83
84///////////////////////////////////////////////////////////////
85// unsmeared model calculation
86///////////////////////////
87Function fProlateForm(w,x) : FitFunc
88        Wave w
89        Variable x
90
91//The input variables are (and output)
92        //[0] scale
93        //[1] crmaj, major radius of core       [A]
94        //[2] crmin, minor radius of core
95        //[3] trmaj, overall major radius
96        //[4] trmin, overall minor radius
97        //[5] sld core, [A^-2]
98        //[6] sld shell,
99        //[7] sld solvent
100        //[8] bkg [cm-1]
101        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg,sldc,slds,sld
102        scale = w[0]
103        crmaj = w[1]
104        crmin = w[2]
105        trmaj = w[3]
106        trmin = w[4]
107        sldc = w[5]
108        slds = w[6]
109        sld = w[7]
110        bkg = w[8]
111
112        delpc = sldc - slds                     //core - shell
113        delps = slds - sld                      //shell - solvent
114// local variables
115        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,prolatevol
116        String weightStr,zStr
117
118        weightStr = "gauss76wt"
119        zStr = "gauss76z"
120
121//      if wt,z waves don't exist, create them
122
123        if (WaveExists(\$weightStr) == 0) // wave reference is not valid,
124                Make/D/N=76 \$weightStr,\$zStr
125                Wave w76 = \$weightStr
126                Wave z76 = \$zStr                // wave references to pass
127                Make76GaussPoints(w76,z76)
128        else
129                if(exists(weightStr) > 1)
130                         Abort "wave name is already in use"    // execute if condition is false
131                endif
132                Wave w76 = \$weightStr
133                Wave z76 = \$zStr
134        endif
135
136// set up the integration
137        // end points and weights
138        nord = 76
139        nfn = 2         //only <f^2> is calculated
140        npro = 1        // PROLATE ELLIPSOIDS
141        va =0
142        vb =1
143//move this zi(i) evaluation inside other nord loop, since I don't have an array
144//      i=0
145//      do
146//       zi[i] = ( z76[i]*(vb-va) + vb + va )/2.0
147 //       i +=1
148 //     while (i<nord)
149//
150// evaluate at Gauss points
151        // remember to index from 0,size-1
152
153        qq = x          //current x point is the q-value for evaluation
154        summ = 0.0
155        ii=0
156        do
157                //printf "top of nord loop, i = %g\r",i
158                if(nfn ==1) //then              // "f1" required for beta factor
159                        if(npro ==1) //then     // prolate
160                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
161//           yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
162                        Endif
163//
164                        if(npro ==0) //then     // oblate
165                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
166//            yyy = w76[i]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
167                        Endif
168                Endif           //nfn = 1
169          //
170                if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
171                        if(npro ==1) //then     //prolate
172                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
173                                yyy = w76[ii]*gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
174                                //printf "yyy = %g\r",yyy
175                        Endif
176//
177                        if(npro ==0) //then     //oblate
178                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
179//              yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
180                        Endif
181                Endif           //nfn <>1
182
183                summ = yyy + summ               // get running total of integral
184                ii+=1
185        while (ii<nord)                         // end of loop over quadrature points
186//
187// calculate value of integral to return
188
189        answer = (vb-va)/2.0*summ
190
191        //normailze by particle volume
192        prolatevol = 4*Pi/3*trmaj*trmin*trmin
193        answer /= prolatevol
194
195        // rescale from 1/A to 1/cm
196        answer *= 1.0e8
197        //scale (arb)
198        answer *= scale
199        ////then add in background
200        answer += bkg
201
202        Return (answer)
203End     //prolate form factor
204
205//
206//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
207//                       BY (53) AND (56,57) IN CHEN AND
208//                       KOTLARCHYK REFERENCE
209//
210//     <PROLATE ELLIPSOIDS>
211//
212Function gfn2(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
213        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
214        // local variables
215        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn
216        Variable siq,sit
217
218        PI43=4.0/3.0*PI
219        aa = crmaj
220        bb = crmin
221        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx))
222        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx))
223        uq = sqrt(u2)*qq
224        ut= sqrt(ut2)*qq
225        vc = PI43*aa*bb*bb
226        vt = PI43*trmaj*trmin*trmin
227
228        if(uq == 0.0)
229                siq = 1.0/3.0
230   else
231                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq
232        endif
233        if(ut == 0.0)
234                sit = 1.0/3.0
235        else
236                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut
237        endif
238
239        gfnc = 3.0*siq*vc*delpc
240        gfnt = 3.0*sit*vt*delps
241        gfn = gfnc+gfnt
242        gfn2 = gfn*gfn
243
244        return gfn2
245End             //function gfn2 for prolate ellipsoids
246
247
248// this is all there is to the smeared calculation!
249Function SmearedProlateForm(s) :FitFunc
250        Struct ResSmearAAOStruct &s
251
252////the name of your unsmeared model is the first argument
253        Smear_Model_20(ProlateForm,s.coefW,s.xW,s.yW,s.resW)
254
255        return(0)
256End
257
258//wrapper to calculate the smeared model as an AAO-Struct
259// fills the struct and calls the ususal function with the STRUCT parameter
260//
261// used only for the dependency, not for fitting
262//
263Function fSmearedProlateForm(coefW,yW,xW)
264        Wave coefW,yW,xW
265
266        String str = getWavesDataFolder(yW,0)
267        String DF="root:"+str+":"
268
269        WAVE resW = \$(DF+str+"_res")
270
271        STRUCT ResSmearAAOStruct fs
272        WAVE fs.coefW = coefW
273        WAVE fs.yW = yW
274        WAVE fs.xW = xW
275        WAVE fs.resW = resW
276
277        Variable err
278        err = SmearedProlateForm(fs)
279
280        return (0)
281End
Note: See TracBrowser for help on using the repository browser.