source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/Cylinder_PolyRadius.ipf @ 151

Last change on this file since 151 was 151, checked in by srkline, 15 years ago

(1) - cursors can now be used to select a subrange of USANS data to fit. This is done by th fit wrapper, assigning a subrange of resW to the struct

(2) all of the smeared model functions are now in the latest form of Smear_Model_N() that is NOT a pointwise calculation anymore, since the USANS matrix smearing in inherently not so.

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