source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models_v3.00/NewModels_2006/Cylinder_PolyRadius.ipf @ 381

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

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

File size: 6.7 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 radius
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//CORRECT! 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_PolyRadius(num,qmin,qmax)
19        Variable num=128,qmin=0.001,qmax=0.7
20        Prompt num "Enter number of data points for model: "
21        Prompt qmin "Enter minimum q-value (A^-1) for model: "
22        Prompt qmax "Enter maximum q-value (A^-1) for model: "
23       
24        make/o/d/n=(num) xwave_cypr,ywave_cypr
25        xwave_cypr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
26        make/o/d coef_cypr = {1.,20.,400,0.2,3.0e-6,0.01}
27        make/o/t parameters_cypr = {"scale","radius (A)","length (A)","polydispersity of Radius","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
28        Edit parameters_cypr,coef_cypr
29        ywave_cypr := Cyl_PolyRadius(coef_cypr,xwave_cypr)
30        Display ywave_cypr vs xwave_cypr
31        ModifyGraph log=1,marker=29,msize=2,mode=4
32        Label bottom "q (A\\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///////////////////////////////////////////////////////////
38
39Proc PlotSmearedCyl_PolyRadius()       
40        //no input parameters necessary, it MUST use the experimental q-values
41        // from the experimental data read in from an AVE/QSIG data file
42       
43        // if no gQvals wave, data must not have been loaded => abort
44        if(ResolutionWavesMissing())
45                Abort
46        endif
47       
48        // Setup parameter table for model function
49        make/o/D smear_coef_cypr = {1.,20.,400,0.2,3.0e-6,0.01}
50        make/o/t smear_parameters_cypr = {"scale","radius (A)","length (A)","polydispersity of Radius","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
51        Edit smear_parameters_cypr,smear_coef_cypr
52       
53        // output smeared intensity wave, dimensions are identical to experimental QSIG values
54        // make extra copy of experimental q-values for easy plotting
55        Duplicate/O $gQvals smeared_cypr,smeared_qvals
56        SetScale d,0,0,"1/cm",smeared_cypr     
57
58        smeared_cypr := SmearedCyl_PolyRadius(smear_coef_cypr,$gQvals)
59        Display smeared_cypr vs smeared_qvals
60        ModifyGraph log=1,marker=29,msize=2,mode=4
61        Label bottom "q (A\\S-1\\M)"
62        Label left "Intensity (cm\\S-1\\M)"
63        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
64End
65
66Function Cyl_PolyRadius(w,x) :FitFunc
67        Wave w
68        Variable x
69
70        //The input variables are (and output)
71        //[0] scale
72        //[1] avg RADIUS (A)
73        //[2] Length (A)
74        //[3] polydispersity (0<p<1)
75        //[4] contrast (A^-2)
76        //[5] background (cm^-1)
77        Variable scale,radius,pd,delrho,bkg,zz,length
78        scale = w[0]
79        radius = w[1]
80        length = w[2]
81        pd = w[3]
82        delrho = w[4]
83        bkg = w[5]
84       
85        zz = (1/pd)^2-1
86//
87// the OUTPUT form factor is <f^2>/Vavg [cm-1]
88//
89// local variables
90        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
91        Variable answer,zp1,zp2,zp3,vpoly
92        String weightStr,zStr
93       
94//      nord = 5       
95//      weightStr = "gauss5wt"
96//      zStr = "gauss5z"
97        nord = 20
98        weightStr = "gauss20wt"
99        zStr = "gauss20z"
100
101//      if wt,z waves don't exist, create them
102// 5 Gauss points (not enough for cylinder radius = high q oscillations)
103// use 20 Gauss points
104        if (WaveExists($weightStr) == 0) // wave reference is not valid,
105                Make/D/N=(nord) $weightStr,$zStr
106                Wave wtGau = $weightStr
107                Wave zGau = $zStr               // wave references to pass
108                Make20GaussPoints(wtGau,zGau)   
109                //Make5GaussPoints(wtGau,zGau) 
110//      //                  printf "w[0],z[0] = %g %g\r", wtGau[0],zGau[0]
111        else
112                if(exists(weightStr) > 1)
113                         Abort "wave name is already in use"    // execute if condition is false
114                endif
115                Wave wtGau = $weightStr
116                Wave zGau = $zStr
117//      //          printf "w[0],z[0] = %g %g\r", wtGau[0],zGau[0]     
118        endif
119
120// set up the integration
121// end points and weights
122// limits are technically 0-inf, but wisely choose non-zero region of distribution
123        Variable range=3.4              //multiples of the std. dev. fom the mean
124        a = radius*(1-range*pd)
125        if (a<0)
126                a=0             //otherwise numerical error when pd >= 0.3, making a<0
127        endif
128        If(pd>0.3)
129                range = 3.4 + (pd-0.3)*18
130        Endif
131        b = radius*(1+range*pd) // is this far enough past avg radius?
132//      printf "a,b,ravg = %g %g %g\r", a,b,radius
133        va =a
134        vb =b
135
136// evaluate at Gauss points
137        // remember to index from 0,size-1     
138        qq = x          //current x point is the q-value for evaluation
139        summ = 0.0              // initialize integral
140   ii=0
141   do
142   //printf "top of nord loop, i = %g\r",i
143        // Using 5 Gauss points         
144                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
145                yyy = wtGau[ii] * rad_kernel(qq,radius,length,zz,delrho,zi)
146                summ = yyy + summ
147                ii+=1
148        while (ii<nord)                         // end of loop over quadrature points
149//   
150// calculate value of integral to return
151   answer = (vb-va)/2.0*summ
152     
153//  contrast^2 is included in integration rad_kernel
154//      answer *= delrho*delrho
155//normalize by polydisperse volume
156// now volume depends on polydisperse RADIUS - so normalize by the second moment
157// 2nd moment = (zz+2)/(zz+1)
158        vpoly = Pi*(radius)^2*length*(zz+2)/(zz+1)
159//Divide by vol, since volume has been "un-normalized" out
160        answer /= vpoly
161//convert to [cm-1]
162        answer *= 1.0e8
163//scale
164        answer *= scale
165// add in the background
166        answer += bkg
167
168        Return (answer)
169End             //End of function PolyRadCylForm()
170
171Function rad_kernel(qw,ravg,len,zz,delrho,rad)
172        Variable qw,ravg,len,zz,delrho,rad
173       
174        Variable Pq,vcyl,dr
175       
176        //calculate the orientationally averaged P(q) for the input rad
177        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170.
178        Make/O/D/n=5 kernpar
179        Wave kp = kernpar
180        kp[0] = 1               //scale fixed at 1
181        kp[1] = rad
182        kp[2] = len
183        kp[3] = delrho
184        kp[4] = 0               //bkg fixed at 0
185       
186        Pq = CylinderForm(kp,qw)
187//      Pq = CylinderFormX(kp,qw)       //from the XOP
188       
189        // undo the normalization that CylinderForm does
190        vcyl=Pi*rad*rad*len
191        Pq *= vcyl
192        //un-convert from [cm-1]
193        Pq /= 1.0e8
194       
195        // calculate normalized distribution at len value
196        dr = Schulz_Point_pr(rad,ravg,zz)
197       
198        return (Pq*dr) 
199End
200
201Function Schulz_Point_pr(x,avg,zz)
202        Variable x,avg,zz
203       
204        Variable dr
205       
206        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
207       
208        return (exp(dr))
209End
210
211// this is all there is to the smeared calculation!
212Function SmearedCyl_PolyRadius(w,x) :FitFunc
213        Wave w
214        Variable x
215       
216        Variable ans
217        SVAR sq = gSig_Q
218        SVAR qb = gQ_bar
219        SVAR sh = gShadow
220        SVAR gQ = gQVals
221       
222        //the name of your unsmeared model is the first argument
223        ans = Smear_Model_20(Cyl_PolyRadius,$sq,$qb,$sh,$gQ,w,x)
224
225        return(ans)
226End
Note: See TracBrowser for help on using the repository browser.