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