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

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

Second round of changes to ipf files to (hopefully) complete the changes to bring all of the model functions to be:

  • this set of changes is for the "2006" set of models
  • AAO
  • new dependencies
  • smearing using structures
  • #if conditional compile for XOPS

-this needs to be FULLY tested for correctness and bugs
-packages are certainly still broken (data folders)

SRK 25 JUL 07

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 conatining 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        s.yW = Smear_Model_20(Cyl_PolyLength,s.coefW,s.xW,s.resW)
268
269        return(0)
270End
271       
Note: See TracBrowser for help on using the repository browser.