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