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

Last change on this file since 127 was 127, checked in by srkline, 16 years ago

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

24 JUL 2007 SRK

File size: 6.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        ywave_cypl := Cyl_PolyLength(coef_cypl,xwave_cypl)
31        Display ywave_cypl vs xwave_cypl
32        ModifyGraph log=1,marker=29,msize=2,mode=4
33        Label bottom "q (\\S-1\\M)"
34        Label left "Intensity (cm\\S-1\\M)"
35        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
36End
37
38// plot the smeared version - quite slow - use only for final fit
39Proc PlotSmearedCyl_PolyLength()       
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_cypl = {1.,20.,1000,0.2,3.0e-6,0.01}
50        make/o/t smear_parameters_cypl = {"scale","radius (A)","length (A)","polydispersity of Length","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
51        Edit smear_parameters_cypl,smear_coef_cypl
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_cypl,smeared_qvals
56        SetScale d,0,0,"1/cm",smeared_cypl     
57
58        smeared_cypl := SmearedCyl_PolyLength(smear_coef_cypl,$gQvals)
59        Display smeared_cypl vs smeared_qvals
60        ModifyGraph log=1,marker=29,msize=2,mode=4
61        Label bottom "q (\\S-1\\M)"
62        Label left "Intensity (cm\\S-1\\M)"
63        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
64End
65
66//calculate the form factor averaged over the size distribution
67// both integrals are done using quadrature, although both may benefit from an
68// adaptive integration
69Function Cyl_PolyLength(w,x)  : FitFunc
70        Wave w
71        Variable x
72
73        //The input variables are (and output)
74        //[0] scale
75        //[1] avg RADIUS (A)
76        //[2] Length (A)
77        //[3] polydispersity (0<p<1)
78        //[4] contrast (A^-2)
79        //[5] background (cm^-1)
80        Variable scale,radius,pd,delrho,bkg,zz,length
81        scale = w[0]
82        radius = w[1]
83        length = w[2]
84        pd = w[3]
85        delrho = w[4]
86        bkg = w[5]
87       
88        zz = (1/pd)^2-1
89//
90// the OUTPUT form factor is <f^2>/Vavg [cm-1]
91//
92// local variables
93        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
94        Variable answer,zp1,zp2,zp3,vpoly
95        String weightStr,zStr
96       
97//      nord = 5       
98//      weightStr = "gauss5wt"
99//      zStr = "gauss5z"
100        nord = 20
101        weightStr = "gauss20wt"
102        zStr = "gauss20z"
103
104// 5 Gauss points (not enough for cylinder radius = high q oscillations)
105// use 20 Gauss points
106        if (WaveExists($weightStr) == 0) // wave reference is not valid,
107                Make/D/N=(nord) $weightStr,$zStr
108                Wave wtGau = $weightStr
109                Wave zGau = $zStr               
110                Make20GaussPoints(wtGau,zGau)   
111                //Make5GaussPoints(wtGau,zGau) 
112        else
113                if(exists(weightStr) > 1)
114                         Abort "wave name is already in use"   
115                endif
116                Wave wtGau = $weightStr
117                Wave zGau = $zStr
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 = length*(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 = length*(1+range*pd) // is this far enough past avg length?
132//      printf "a,b,len_avg = %g %g %g\r", a,b,length
133        va =a
134        vb =b
135
136        qq = x          //current x point is the q-value for evaluation
137        summ = 0.0              // initialize integral
138   ii=0
139   do
140   //printf "top of nord loop, i = %g\r",i
141        // Using 5 Gauss points         
142                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
143                yyy = wtGau[ii] * len_kernel(qq,radius,length,zz,delrho,zi)
144                summ = yyy + summ
145                ii+=1
146        while (ii<nord)                         // end of loop over quadrature points
147//   
148// calculate value of integral to return
149   answer = (vb-va)/2.0*summ
150     
151//  contrast^2 is included in integration rad_kernel
152//      answer *= delrho*delrho
153//normalize by polydisperse volume
154// now volume depends on polydisperse Length - so normalize by the FIRST moment
155// 1st moment = volume!
156        vpoly = Pi*(radius)^2*length
157//Divide by vol, since volume has been "un-normalized" out
158        answer /= vpoly
159//convert to [cm-1]
160        answer *= 1.0e8
161//scale
162        answer *= scale
163// add in the background
164        answer += bkg
165
166        Return (answer)
167End             //End of function PolyRadCylForm()
168
169Function len_kernel(qw,rad,len_avg,zz,delrho,len)
170        Variable qw,rad,len_avg,zz,delrho,len
171       
172        Variable Pq,vcyl,dl
173       
174        //calculate the orientationally averaged P(q) for the input rad
175        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170.
176        Make/O/n=5 kernpar
177        Wave kp = kernpar
178        kp[0] = 1               //scale fixed at 1
179        kp[1] = rad
180        kp[2] = len
181        kp[3] = delrho
182        kp[4] = 0               //bkg fixed at 0
183       
184        Pq = CylinderForm(kp,qw)
185        //Pq = CylinderFit(kp,qw)       //from the XOP
186       
187        // undo the normalization that CylinderForm does
188        //CylinderForm returns P(q)/V, we want P(q)
189        vcyl=Pi*rad*rad*len
190        Pq *= vcyl
191        //un-convert from [cm-1]
192        Pq /= 1.0e8
193       
194        // calculate normalized distribution at len value
195        dl = Schulz_Point_pollen(len,len_avg,zz)
196       
197        return (Pq*dl) 
198End
199
200Function Schulz_Point_pollen(x,avg,zz)
201        Variable x,avg,zz
202       
203        Variable dr
204       
205        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
206       
207        return (exp(dr))
208End
209
210// this is all there is to the smeared calculation!
211Function SmearedCyl_PolyLength(w,x) :FitFunc
212        Wave w
213        Variable x
214       
215        Variable ans
216        SVAR sq = gSig_Q
217        SVAR qb = gQ_bar
218        SVAR sh = gShadow
219        SVAR gQ = gQVals
220       
221        //the name of your unsmeared model is the first argument
222        ans = Smear_Model_20(Cyl_PolyLength,$sq,$qb,$sh,$gQ,w,x)
223
224        return(ans)
225End
Note: See TracBrowser for help on using the repository browser.