# source:sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/Cylinder_PolyLength_v40.ipf@379

Last change on this file since 379 was 379, checked in by srkline, 14 years ago

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

File size: 7.9 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "Cylinder_v40"
5
6// calculates the form factor of a cylinder with polydispersity of length
7// the length distribution is a Schulz distribution, and any normalized distribution
8// could be used, as the average is performed numerically
9//
10// since the cylinder form factor is already a numerical integration, the size average is a
11// second integral, and significantly slows the calculation, and smearing adds a third integration.
12//
13//CORRECTED 12/5/2000 - Invariant is now correct vs. monodisperse cylinders
14// + upper limit of integration has been changed to account for skew of
15//Schulz distribution at high (>0.5) polydispersity
16//Requires 20 gauss points for integration of the radius (5 is not enough)
17//Requires either CylinderFit XOP (MacOSX only) or the normal CylinderForm Function
18//
19Proc PlotCyl_PolyLength(num,qmin,qmax)
20        Variable num=100,qmin=0.001,qmax=0.7
21        Prompt num "Enter number of data points for model: "
22        Prompt qmin "Enter minimum q-value (A^-1) for model: "
23        Prompt qmax "Enter maximum q-value (A^-1) for model: "
24
25        make/o/d/n=(num) xwave_cypl,ywave_cypl
26        xwave_cypl = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
27        make/o/d coef_cypl = {1.,20.,1000,0.2,1e-6,6.3e-6,0.01}
28        make/o/t parameters_cypl = {"scale","radius (A)","length (A)","polydispersity of Length","SLD cylinder (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
29        Edit parameters_cypl,coef_cypl
30
31        Variable/G root:g_cypl
32        g_cypl := Cyl_PolyLength(coef_cypl,ywave_cypl,xwave_cypl)
33        Display ywave_cypl vs xwave_cypl
34        ModifyGraph log=1,marker=29,msize=2,mode=4
35        Label bottom "q (A\\S-1\\M)"
36        Label left "Intensity (cm\\S-1\\M)"
37        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
38
39        AddModelToStrings("Cyl_PolyLength","coef_cypl","cypl")
40End
41
42// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
43Proc PlotSmearedCyl_PolyLength(str)
44        String str
45        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
46
47        // if any of the resolution waves are missing => abort
48        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
49                Abort
50        endif
51
52        SetDataFolder \$("root:"+str)
53
54        // Setup parameter table for model function
55        make/o/D smear_coef_cypl = {1.,20.,1000,0.2,1e-6,6.3e-6,0.01}
56        make/o/t smear_parameters_cypl = {"scale","radius (A)","length (A)","polydispersity of Length","SLD cylinder (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
57        Edit smear_parameters_cypl,smear_coef_cypl
58
59        // output smeared intensity wave, dimensions are identical to experimental QSIG values
60        // make extra copy of experimental q-values for easy plotting
61        Duplicate/O \$(str+"_q") smeared_cypl,smeared_qvals
62        SetScale d,0,0,"1/cm",smeared_cypl
63
64        Variable/G gs_cypl=0
65        gs_cypl := fSmearedCyl_PolyLength(smear_coef_cypl,smeared_cypl,smeared_qvals)   //this wrapper fills the STRUCT
66
67        Display smeared_cypl vs smeared_qvals
68        ModifyGraph log=1,marker=29,msize=2,mode=4
69        Label bottom "q (A\\S-1\\M)"
70        Label left "Intensity (cm\\S-1\\M)"
71        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
72
73        SetDataFolder root:
74        AddModelToStrings("SmearedCyl_PolyLength","smear_coef_cypl","cypl")
75End
76
77
78
79//AAO version, uses XOP if available
80// simply calls the original single point calculation with
81// a wave assignment (this will behave nicely if given point ranges)
82Function Cyl_PolyLength(cw,yw,xw) : FitFunc
83        Wave cw,yw,xw
84
85#if exists("Cyl_PolyLengthX")
86        yw = Cyl_PolyLengthX(cw,xw)
87#else
88        yw = fCyl_PolyLength(cw,xw)
89#endif
90        return(0)
91End
92
93//calculate the form factor averaged over the size distribution
94// both integrals are done using quadrature, although both may benefit from an
95// adaptive integration
96Function fCyl_PolyLength(w,x)  : FitFunc
97        Wave w
98        Variable x
99
100        //The input variables are (and output)
101        //[0] scale
102        //[1] avg RADIUS (A)
103        //[2] Length (A)
104        //[3] polydispersity (0<p<1)
105        //[4] contrast (A^-2)
106        //[5] background (cm^-1)
107        Variable scale,radius,pd,delrho,bkg,zz,length,sldc,slds
108        scale = w[0]
109        radius = w[1]
110        length = w[2]
111        pd = w[3]
112        sldc = w[4]
113        slds = w[5]
114        delrho = sldc - slds
115        bkg = w[6]
116
117        zz = (1/pd)^2-1
118//
119// the OUTPUT form factor is <f^2>/Vavg [cm-1]
120//
121// local variables
122        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
123        Variable answer,zp1,zp2,zp3,vpoly
124        String weightStr,zStr
125
126//      nord = 5
127//      weightStr = "gauss5wt"
128//      zStr = "gauss5z"
129        nord = 20
130        weightStr = "gauss20wt"
131        zStr = "gauss20z"
132
133// 5 Gauss points (not enough for cylinder radius = high q oscillations)
134// use 20 Gauss points
135        if (WaveExists(\$weightStr) == 0) // wave reference is not valid,
136                Make/D/N=(nord) \$weightStr,\$zStr
137                Wave wtGau = \$weightStr
138                Wave zGau = \$zStr
139                Make20GaussPoints(wtGau,zGau)
140                //Make5GaussPoints(wtGau,zGau)
141        else
142                if(exists(weightStr) > 1)
143                         Abort "wave name is already in use"
144                endif
145                Wave wtGau = \$weightStr
146                Wave zGau = \$zStr
147        endif
148
149// set up the integration
150// end points and weights
151// limits are technically 0-inf, but wisely choose non-zero region of distribution
152        Variable range=3.4              //multiples of the std. dev. fom the mean
153        a = length*(1-range*pd)
154        if (a<0)
155                a=0             //otherwise numerical error when pd >= 0.3, making a<0
156        endif
157        If(pd>0.3)
158                range = 3.4 + (pd-0.3)*18
159        Endif
160        b = length*(1+range*pd) // is this far enough past avg length?
161//      printf "a,b,len_avg = %g %g %g\r", a,b,length
162        va =a
163        vb =b
164
165        qq = x          //current x point is the q-value for evaluation
166        summ = 0.0              // initialize integral
167   ii=0
168   do
169   //printf "top of nord loop, i = %g\r",i
170        // Using 5 Gauss points
171                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0
172                yyy = wtGau[ii] * len_kernel(qq,radius,length,zz,sldc,slds,zi)
173                summ = yyy + summ
174                ii+=1
175        while (ii<nord)                         // end of loop over quadrature points
176//
177// calculate value of integral to return
178   answer = (vb-va)/2.0*summ
179
180//  contrast^2 is included in integration rad_kernel
181//      answer *= delrho*delrho
182//normalize by polydisperse volume
183// now volume depends on polydisperse Length - so normalize by the FIRST moment
184// 1st moment = volume!
185        vpoly = Pi*(radius)^2*length
186//Divide by vol, since volume has been "un-normalized" out
187        answer /= vpoly
188//convert to [cm-1]
189        answer *= 1.0e8
190//scale
191        answer *= scale
192// add in the background
193        answer += bkg
194
195        Return (answer)
196End             //End of function PolyRadCylForm()
197
198Function len_kernel(qw,rad,len_avg,zz,sldc,slds,len)
199        Variable qw,rad,len_avg,zz,sldc,slds,len
200
201        Variable Pq,vcyl,dl
202
203        //calculate the orientationally averaged P(q) for the input rad
204        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170.
205        Make/O/n=6 kernpar
206        Wave kp = kernpar
207        kp[0] = 1               //scale fixed at 1
208        kp[1] = rad
209        kp[2] = len
210        kp[3] = sldc
211        kp[4] = slds
212        kp[5] = 0               //bkg fixed at 0
213
214#if exists("CylinderFormX")
215        Pq = CylinderFormX(kp,qw)
216#else
217        Pq = fCylinderForm(kp,qw)
218#endif
219
220        // undo the normalization that CylinderForm does
221        //CylinderForm returns P(q)/V, we want P(q)
222        vcyl=Pi*rad*rad*len
223        Pq *= vcyl
224        //un-convert from [cm-1]
225        Pq /= 1.0e8
226
227        // calculate normalized distribution at len value
228        dl = Schulz_Point_pollen(len,len_avg,zz)
229
230        return (Pq*dl)
231End
232
233Function Schulz_Point_pollen(x,avg,zz)
234        Variable x,avg,zz
235
236        Variable dr
237
238        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
239
240        return (exp(dr))
241End
242
243//wrapper to calculate the smeared model as an AAO-Struct
244// fills the struct and calls the ususal function with the STRUCT parameter
245//
246// used only for the dependency, not for fitting
247//
248Function fSmearedCyl_PolyLength(coefW,yW,xW)
249        Wave coefW,yW,xW
250
251        String str = getWavesDataFolder(yW,0)
252        String DF="root:"+str+":"
253
254        WAVE resW = \$(DF+str+"_res")
255
256        STRUCT ResSmearAAOStruct fs
257        WAVE fs.coefW = coefW
258        WAVE fs.yW = yW
259        WAVE fs.xW = xW
260        WAVE fs.resW = resW
261
262        Variable err
263        err = SmearedCyl_PolyLength(fs)
264
265        return (0)
266End
267
268// this is all there is to the smeared calculation!
269Function SmearedCyl_PolyLength(s) :FitFunc
270        Struct ResSmearAAOStruct &s
271
272//      the name of your unsmeared model (AAO) is the first argument
273        Smear_Model_20(Cyl_PolyLength,s.coefW,s.xW,s.yW,s.resW)
274
275        return(0)
276End
277
Note: See TracBrowser for help on using the repository browser.