source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/RectPolySpheres.ipf @ 146

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

Typo in dialog in every model...

File size: 6.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
6// Adopting these into the experiment will insure that they are always present
7////////////////////////////////////////////////
8// this example is for the form factor of polydisperse spheres
9// no interparticle interactions are included
10// the polydispersity in radius is a rectangular distribution
11//
12// 13 JAN 99 SRK
13////////////////////////////////////////////////
14
15Proc PlotPolyRectSpheres(num,qmin,qmax)
16        Variable num=128,qmin=0.001,qmax=0.7
17        Prompt num "Enter number of data points for model: "
18        Prompt qmin "Enter minimum q-value (^-1) for model: "
19        Prompt qmax "Enter maximum q-value (^-1) for model: "
20       
21        Make/O/D/n=(num) xwave_rect,ywave_rect
22        xwave_rect =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
23        Make/O/D coef_rect = {1,60,0.12,3.0e-6,0.}
24        make/o/t parameters_rect = {"scale","Radius (A)","polydispersity","contrast (A^-2)","background (cm^-1)"}
25        Edit parameters_rect,coef_rect
26        Variable/G root:g_rect
27        g_rect := PolyRectSpheres(coef_rect,ywave_rect,xwave_rect)
28//      ywave_rect := PolyRectSpheres(coef_rect,xwave_rect)
29        Display ywave_rect vs xwave_rect
30        ModifyGraph log=1,marker=29,msize=2,mode=4
31        Label bottom "q (\\S-1\\M)"
32        Label left "Intensity (cm\\S-1\\M)"
33//      DoAlert 0,"The form facor is not properly normalized with the polydisperse volume"
34        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
35End
36
37
38///////////////////////////////////////////////////////////
39// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
40Proc PlotSmearedPolyRectSpheres(str)                                                           
41        String str
42        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
43       
44        // if any of the resolution waves are missing => abort
45        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
46                Abort
47        endif
48       
49        SetDataFolder $("root:"+str)
50       
51        // Setup parameter table for model function
52        Make/O/D smear_coef_rect = {1,60,0.12,3.0e-6,0.}
53        make/o/t smear_parameters_rect = {"scale","Radius (A)","polydispersity","contrast (A^-2)","background (cm^-1)"}
54        Edit smear_parameters_rect,smear_coef_rect
55       
56        // output smeared intensity wave, dimensions are identical to experimental QSIG values
57        // make extra copy of experimental q-values for easy plotting
58        Duplicate/O $(str+"_q") smeared_rect,smeared_qvals
59        SetScale d,0,0,"1/cm",smeared_rect                     
60               
61        Variable/G gs_rect=0
62        gs_rect := fSmearedPolyRectSpheres(smear_coef_rect,smeared_rect,smeared_qvals)  //this wrapper fills the STRUCT
63
64        Display smeared_rect vs smeared_qvals
65        ModifyGraph log=1,marker=29,msize=2,mode=4
66        Label bottom "q (\\S-1\\M)"
67        Label left "Intensity (cm\\S-1\\M)"
68        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) 
69       
70        SetDataFolder root:
71End
72
73
74//AAO version
75Function PolyRectSpheres(cw,yw,xw) : FitFunc
76                Wave cw,yw,xw
77
78#if exists("PolyRectSpheresX")
79        yw = PolyRectSpheresX(cw,xw)
80#else
81        yw = fPolyRectSpheres(cw,xw)
82#endif
83        return(0)
84End
85///////////////////////////////////////////////////////////////
86// unsmeared model calculation
87///////////////////////////
88Function fPolyRectSpheres(w,x) : FitFunc
89        Wave w                  // the coefficient wave
90        Variable x              // the x values, as a variable
91         
92        //* reassign names to the variable set */
93         Variable scale,rad,pd,cont,bkg
94         
95        scale = w[0]
96        rad = w[1]              // radius (A)
97        pd = w[2]               //polydispersity of rectangular distribution
98        cont = w[3]             // contrast (A^-2)
99        bkg = w[4]              // background (1/cm)
100
101// local variables
102        Variable inten,h1,qw,qr,width,sig,averad3,Vavg,Rg2
103
104        // as usual, poly = sig/ravg
105        // for the rectangular distribution, sig = width/sqrt(3)
106        // width is the HALF- WIDTH of the rectangular distribution
107       
108        sig = pd*rad
109        width = sqrt(3)*sig
110       
111        //x is the q-value
112        qw = x*width
113        qr = x*rad
114       
115        // as for the low QR "crud", the function is calculating the sines and cosines just fine
116        // - the problem seems to be that the
117        // leading terms nearly cancel with the last term (the -6*qr... term), to within machine
118        // precision - the difference is on the order of 10^-20
119        // so just use the limiting Guiner value
120        if(qr<0.1)
121                h1 = scale*cont*cont*1e8*4*pi/3*rad^3
122                h1 *=   (1 + 15*pd^2 + 27*pd^4 +27/7*pd^6)                              //6th moment
123                h1 /= (1+3*pd^2)                //3rd moment
124                Rg2 = 3/5*rad*rad*(1+28*pd^2+126*pd^4+108*pd^6+27*pd^8)
125                Rg2 /= (1+15*pd^2+27*pd^4+27/7*pd^6)
126               
127                h1 *= exp(-1/3*Rg2*x*x)
128                h1 += bkg
129                return(h1)
130        endif
131       
132       
133        //normal calculation
134        h1 = -0.5*qw + qr*qr*qw + (qw^3)/3
135        h1 -= 5/2*cos(2*qr)*sin(qw)*cos(qw)
136        h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw)
137        h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw)
138        h1 += qw*qr*sin(2*qr)*cos(2*qw)
139        h1 += 3*qw*(cos(qr)*cos(qw))^2
140        h1 += 3*qw*(sin(qr)*sin(qw))^2
141       
142        h1 -= 6*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw)
143       
144        // calculate P(q) = <f^2>
145        inten = 8*Pi*Pi*cont*cont/width/x^7*h1
146       
147// beta(q) would be calculated as 2/width/x/h1*h2*h2
148// with
149// h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
150
151        // normalize to the average volume
152        // <R^3> = ravg^3*(1+3*pd^2)
153        // or... "zf"  = (1 + 3*p^2), which will be greater than one
154       
155        averad3 =  rad^3*(1+3*pd^2)
156        inten /= 4*pi/3*averad3
157        //resacle to 1/cm
158        inten *= 1.0e8
159        //scale the result
160        inten *= scale
161        // then add in the background
162        inten += bkg
163       
164        return (inten)
165
166End   // end of PolyRectSpheres()
167
168// this is all there is to the smeared calculation!
169Function SmearedPolyRectSpheres(s) :FitFunc
170        Struct ResSmearAAOStruct &s
171
172////the name of your unsmeared model is the first argument
173        s.yW = Smear_Model_20(PolyRectSpheres,s.coefW,s.xW,s.resW)
174
175        return(0)
176End
177
178//wrapper to calculate the smeared model as an AAO-Struct
179// fills the struct and calls the ususal function with the STRUCT parameter
180//
181// used only for the dependency, not for fitting
182//
183Function fSmearedPolyRectSpheres(coefW,yW,xW)
184        Wave coefW,yW,xW
185       
186        String str = getWavesDataFolder(yW,0)
187        String DF="root:"+str+":"
188       
189        WAVE resW = $(DF+str+"_res")
190       
191        STRUCT ResSmearAAOStruct fs
192        WAVE fs.coefW = coefW   
193        WAVE fs.yW = yW
194        WAVE fs.xW = xW
195        WAVE fs.resW = resW
196       
197        Variable err
198        err = SmearedPolyRectSpheres(fs)
199       
200        return (0)
201End
Note: See TracBrowser for help on using the repository browser.