source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models_v3.00/NewModels_2006/Cylinder_PolyLength.ipf @ 338

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

merging structural changes from Dev/trunk into Release/trunk

File size: 6.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3#include "CylinderForm"
4
5// calculates the form factor of a cylinder with polydispersity of length
6// the length distribution is a Schulz distribution, and any normalized distribution
7// could be used, as the average is performed numerically
8//
9// since the cylinder form factor is already a numerical integration, the size average is a
10// second integral, and significantly slows the calculation, and smearing adds a third integration.
11//
12//CORRECTED 12/5/2000 - Invariant is now correct vs. monodisperse cylinders
13// + upper limit of integration has been changed to account for skew of
14//Schulz distribution at high (>0.5) polydispersity
15//Requires 20 gauss points for integration of the radius (5 is not enough)
16//Requires either CylinderFit XOP (MacOSX only) or the normal CylinderForm Function
17//
18Proc PlotCyl_PolyLength(num,qmin,qmax)
19        Variable num=100,qmin=0.001,qmax=0.7
20        Prompt num "Enter number of data points for model: "
21        Prompt qmin "Enter minimum q-value (^-1) for model: "
22        Prompt qmax "Enter maximum q-value (^-1) for model: "
23       
24        make/o/d/n=(num) xwave_cypl,ywave_cypl
25        xwave_cypl = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
26        make/o/d coef_cypl = {1.,20.,1000,0.2,3.0e-6,0.01}
27        make/o/t parameters_cypl = {"scale","radius (A)","length (A)","polydispersity of Length","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
28        Edit parameters_cypl,coef_cypl
29        ywave_cypl := Cyl_PolyLength(coef_cypl,xwave_cypl)
30        Display ywave_cypl vs xwave_cypl
31        ModifyGraph log=1,marker=29,msize=2,mode=4
32        Label bottom "q (\\S-1\\M)"
33        Label left "Intensity (cm\\S-1\\M)"
34        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
35End
36
37// plot the smeared version - quite slow - use only for final fit
38Proc PlotSmearedCyl_PolyLength()       
39        //no input parameters necessary, it MUST use the experimental q-values
40        // from the experimental data read in from an AVE/QSIG data file
41       
42        // if no gQvals wave, data must not have been loaded => abort
43        if(ResolutionWavesMissing())
44                Abort
45        endif
46       
47        // Setup parameter table for model function
48        make/o/D smear_coef_cypl = {1.,20.,1000,0.2,3.0e-6,0.01}
49        make/o/t smear_parameters_cypl = {"scale","radius (A)","length (A)","polydispersity of Length","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
50        Edit smear_parameters_cypl,smear_coef_cypl
51       
52        // output smeared intensity wave, dimensions are identical to experimental QSIG values
53        // make extra copy of experimental q-values for easy plotting
54        Duplicate/O $gQvals smeared_cypl,smeared_qvals
55        SetScale d,0,0,"1/cm",smeared_cypl     
56
57        smeared_cypl := SmearedCyl_PolyLength(smear_coef_cypl,$gQvals)
58        Display smeared_cypl vs smeared_qvals
59        ModifyGraph log=1,marker=29,msize=2,mode=4
60        Label bottom "q (\\S-1\\M)"
61        Label left "Intensity (cm\\S-1\\M)"
62        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
63End
64
65//calculate the form factor averaged over the size distribution
66// both integrals are done using quadrature, although both may benefit from an
67// adaptive integration
68Function Cyl_PolyLength(w,x)  : FitFunc
69        Wave w
70        Variable x
71
72        //The input variables are (and output)
73        //[0] scale
74        //[1] avg RADIUS (A)
75        //[2] Length (A)
76        //[3] polydispersity (0<p<1)
77        //[4] contrast (A^-2)
78        //[5] background (cm^-1)
79        Variable scale,radius,pd,delrho,bkg,zz,length
80        scale = w[0]
81        radius = w[1]
82        length = w[2]
83        pd = w[3]
84        delrho = w[4]
85        bkg = w[5]
86       
87        zz = (1/pd)^2-1
88//
89// the OUTPUT form factor is <f^2>/Vavg [cm-1]
90//
91// local variables
92        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
93        Variable answer,zp1,zp2,zp3,vpoly
94        String weightStr,zStr
95       
96//      nord = 5       
97//      weightStr = "gauss5wt"
98//      zStr = "gauss5z"
99        nord = 20
100        weightStr = "gauss20wt"
101        zStr = "gauss20z"
102
103// 5 Gauss points (not enough for cylinder radius = high q oscillations)
104// use 20 Gauss points
105        if (WaveExists($weightStr) == 0) // wave reference is not valid,
106                Make/D/N=(nord) $weightStr,$zStr
107                Wave wtGau = $weightStr
108                Wave zGau = $zStr               
109                Make20GaussPoints(wtGau,zGau)   
110                //Make5GaussPoints(wtGau,zGau) 
111        else
112                if(exists(weightStr) > 1)
113                         Abort "wave name is already in use"   
114                endif
115                Wave wtGau = $weightStr
116                Wave zGau = $zStr
117        endif
118
119// set up the integration
120// end points and weights
121// limits are technically 0-inf, but wisely choose non-zero region of distribution
122        Variable range=3.4              //multiples of the std. dev. fom the mean
123        a = length*(1-range*pd)
124        if (a<0)
125                a=0             //otherwise numerical error when pd >= 0.3, making a<0
126        endif
127        If(pd>0.3)
128                range = 3.4 + (pd-0.3)*18
129        Endif
130        b = length*(1+range*pd) // is this far enough past avg length?
131//      printf "a,b,len_avg = %g %g %g\r", a,b,length
132        va =a
133        vb =b
134
135        qq = x          //current x point is the q-value for evaluation
136        summ = 0.0              // initialize integral
137   ii=0
138   do
139   //printf "top of nord loop, i = %g\r",i
140        // Using 5 Gauss points         
141                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
142                yyy = wtGau[ii] * len_kernel(qq,radius,length,zz,delrho,zi)
143                summ = yyy + summ
144                ii+=1
145        while (ii<nord)                         // end of loop over quadrature points
146//   
147// calculate value of integral to return
148   answer = (vb-va)/2.0*summ
149     
150//  contrast^2 is included in integration rad_kernel
151//      answer *= delrho*delrho
152//normalize by polydisperse volume
153// now volume depends on polydisperse Length - so normalize by the FIRST moment
154// 1st moment = volume!
155        vpoly = Pi*(radius)^2*length
156//Divide by vol, since volume has been "un-normalized" out
157        answer /= vpoly
158//convert to [cm-1]
159        answer *= 1.0e8
160//scale
161        answer *= scale
162// add in the background
163        answer += bkg
164
165        Return (answer)
166End             //End of function PolyRadCylForm()
167
168Function len_kernel(qw,rad,len_avg,zz,delrho,len)
169        Variable qw,rad,len_avg,zz,delrho,len
170       
171        Variable Pq,vcyl,dl
172       
173        //calculate the orientationally averaged P(q) for the input rad
174        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170.
175        Make/O/n=5 kernpar
176        Wave kp = kernpar
177        kp[0] = 1               //scale fixed at 1
178        kp[1] = rad
179        kp[2] = len
180        kp[3] = delrho
181        kp[4] = 0               //bkg fixed at 0
182       
183        Pq = CylinderForm(kp,qw)
184        //Pq = CylinderFit(kp,qw)       //from the XOP
185       
186        // undo the normalization that CylinderForm does
187        //CylinderForm returns P(q)/V, we want P(q)
188        vcyl=Pi*rad*rad*len
189        Pq *= vcyl
190        //un-convert from [cm-1]
191        Pq /= 1.0e8
192       
193        // calculate normalized distribution at len value
194        dl = Schulz_Point_pollen(len,len_avg,zz)
195       
196        return (Pq*dl) 
197End
198
199Function Schulz_Point_pollen(x,avg,zz)
200        Variable x,avg,zz
201       
202        Variable dr
203       
204        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
205       
206        return (exp(dr))
207End
208
209// this is all there is to the smeared calculation!
210Function SmearedCyl_PolyLength(w,x) :FitFunc
211        Wave w
212        Variable x
213       
214        Variable ans
215        SVAR sq = gSig_Q
216        SVAR qb = gQ_bar
217        SVAR sh = gShadow
218        SVAR gQ = gQVals
219       
220        //the name of your unsmeared model is the first argument
221        ans = Smear_Model_20(Cyl_PolyLength,$sq,$qb,$sh,$gQ,w,x)
222
223        return(ans)
224End
Note: See TracBrowser for help on using the repository browser.