source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/RectPolySpheres.ipf @ 56

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

initial checkin of Analysis v.3.00 files

File size: 5.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////
4// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
5// Adopting these into the experiment will insure that they are always present
6////////////////////////////////////////////////
7// this example is for the form factor of polydisperse spheres
8// no interparticle interactions are included
9// the polydispersity in radius is a rectangular distribution
10//
11// 13 JAN 99 SRK
12////////////////////////////////////////////////
13
14Proc PlotPolyRectSpheres(num,qmin,qmax)
15        Variable num=128,qmin=0.001,qmax=0.7
16        Prompt num "Enter number of data points for model: "
17        Prompt qmin "Enter minimum q-value (^-1) for model: "
18        Prompt qmax "Enter maximum q-value (^-1) for model: "
19       
20        Make/O/D/n=(num) xwave_rect,ywave_rect
21        xwave_rect =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
22        Make/O/D coef_rect = {1,60,0.12,3.0e-6,0.}
23        make/o/t parameters_rect = {"scale","Radius (A)","polydispersity","contrast (A^-2)","background (cm^-1)"}
24        Edit parameters_rect,coef_rect
25        ywave_rect := PolyRectSpheres(coef_rect,xwave_rect)
26        Display ywave_rect vs xwave_rect
27        ModifyGraph log=1,marker=29,msize=2,mode=4
28        Label bottom "q (\\S-1\\M)"
29        Label left "Intensity (cm\\S-1\\M)"
30//      DoAlert 0,"The form facor is not properly normalized with the polydisperse volume"
31        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
32End
33
34
35///////////////////////////////////////////////////////////
36
37Proc PlotSmearedPolyRectSpheres()               
38        //no input parameters necessary, it MUST use the experimental q-values
39        // from the experimental data read in from an AVE/QSIG data file
40       
41        // if no gQvals wave, data must not have been loaded => abort
42        if(ResolutionWavesMissing())
43                Abort
44        endif
45       
46        // Setup parameter table for model function
47        Make/O/D smear_coef_rect = {1,60,0.12,3.0e-6,0.}
48        make/o/t smear_parameters_rect = {"scale","Radius (A)","polydispersity","contrast (A^-2)","background (cm^-1)"}
49        Edit smear_parameters_rect,smear_coef_rect
50       
51        // output smeared intensity wave, dimensions are identical to experimental QSIG values
52        // make extra copy of experimental q-values for easy plotting
53        Duplicate/O $gQvals smeared_rect,smeared_qvals
54        SetScale d,0,0,"1/cm",smeared_rect
55
56        smeared_rect := SmearedPolyRectSpheres(smear_coef_rect,$gQvals)
57        Display smeared_rect vs smeared_qvals
58        ModifyGraph log=1,marker=29,msize=2,mode=4
59        Label bottom "q (\\S-1\\M)"
60        Label left "Intensity (cm\\S-1\\M)"
61        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
62End
63
64
65///////////////////////////////////////////////////////////////
66// unsmeared model calculation
67///////////////////////////
68Function PolyRectSpheres(w,x) : FitFunc
69        Wave w                  // the coefficient wave
70        Variable x              // the x values, as a variable
71         
72        //* reassign names to the variable set */
73         Variable scale,rad,pd,cont,bkg
74         
75        scale = w[0]
76        rad = w[1]              // radius (A)
77        pd = w[2]               //polydispersity of rectangular distribution
78        cont = w[3]             // contrast (A^-2)
79        bkg = w[4]              // background (1/cm)
80
81// local variables
82        Variable inten,h1,qw,qr,width,sig,averad3,Vavg,Rg2
83
84        // as usual, poly = sig/ravg
85        // for the rectangular distribution, sig = width/sqrt(3)
86        // width is the HALF- WIDTH of the rectangular distribution
87       
88        sig = pd*rad
89        width = sqrt(3)*sig
90       
91        //x is the q-value
92        qw = x*width
93        qr = x*rad
94       
95        // as for the low QR "crud", the function is calculating the sines and cosines just fine
96        // - the problem seems to be that the
97        // leading terms nearly cancel with the last term (the -6*qr... term), to within machine
98        // precision - the difference is on the order of 10^-20
99        // so just use the limiting Guiner value
100        if(qr<0.1)
101                h1 = scale*cont*cont*1e8*4*pi/3*rad^3
102                h1 *=   (1 + 15*pd^2 + 27*pd^4 +27/7*pd^6)                              //6th moment
103                h1 /= (1+3*pd^2)                //3rd moment
104                Rg2 = 3/5*rad*rad*(1+28*pd^2+126*pd^4+108*pd^6+27*pd^8)
105                Rg2 /= (1+15*pd^2+27*pd^4+27/7*pd^6)
106               
107                h1 *= exp(-1/3*Rg2*x*x)
108                h1 += bkg
109                return(h1)
110        endif
111       
112       
113        //normal calculation
114        h1 = -0.5*qw + qr*qr*qw + (qw^3)/3
115        h1 -= 5/2*cos(2*qr)*sin(qw)*cos(qw)
116        h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw)
117        h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw)
118        h1 += qw*qr*sin(2*qr)*cos(2*qw)
119        h1 += 3*qw*(cos(qr)*cos(qw))^2
120        h1 += 3*qw*(sin(qr)*sin(qw))^2
121       
122        h1 -= 6*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw)
123       
124        // calculate P(q) = <f^2>
125        inten = 8*Pi*Pi*cont*cont/width/x^7*h1
126       
127// beta(q) would be calculated as 2/width/x/h1*h2*h2
128// with
129// h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
130
131        // normalize to the average volume
132        // <R^3> = ravg^3*(1+3*pd^2)
133        // or... "zf"  = (1 + 3*p^2), which will be greater than one
134       
135        averad3 =  rad^3*(1+3*pd^2)
136        inten /= 4*pi/3*averad3
137        //resacle to 1/cm
138        inten *= 1.0e8
139        //scale the result
140        inten *= scale
141        // then add in the background
142        inten += bkg
143       
144        return (inten)
145
146End   // end of PolyRectSpheres()
147
148// this is all there is to the smeared calculation!
149Function SmearedPolyRectSpheres(w,x) :FitFunc
150        Wave w
151        Variable x
152       
153        Variable ans
154        SVAR sq = gSig_Q
155        SVAR qb = gQ_bar
156        SVAR sh = gShadow
157        SVAR gQ = gQVals
158       
159        //the name of your unsmeared model is the first argument
160        ans = Smear_Model_20(PolyRectSpheres,$sq,$qb,$sh,$gQ,w,x)
161
162        return(ans)
163End
Note: See TracBrowser for help on using the repository browser.