source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Cylinder_PolyLength_v40.ipf @ 510

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 8.0 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","parameters_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","smear_parameters_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.