source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/SC_ParaCrystal_v40.ipf @ 515

Last change on this file since 515 was 515, checked in by srkline, 13 years ago

A bunch of random changes:
-SASCALC - added the 2.5 cm aperture NG3/7g/polarizer back in, but made the 5 cm the default
-Some fiddling with 2D functionality to enable 2D resolution smearing. Still a work in progress
-Added two new model functions, and added them to the list.
-Changes to the wrapper for the new release to generate kw=val strings as needed for each function
and its coefficients and parameters. This behavior has been changed in the new release, so this
fix should keep old analysis experiments compatible with the new version.

File size: 8.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.0
3
4//
5//
6// Simple Cubic paracrystal, powder average
7//
8// VERY slow, since the function is so ill-behaved and needs LOTS of quadrature
9// points. Adaptive methods were even slower and troublesom to converge,
10// although in theory they should be a better choice than blindly increasing the number of points.
11//
12// 150 points seems to give reasonable reproduction of the peak heights in the paper.
13// peak locations are correct
14// 76 points of quadrature for the smearing is only a guess, it's not been tested yet.
15//
16// Original implementation - Danilo Pozzo
17//              modified and modernized for more efficient integration SRK Nov 2008
18//
19//REFERENCE
20//Hideki Matsuoka etal. Physical Review B, Vol 36 Num 3, p1754 1987   ORIGINAL PAPER
21//Hideki Matsuoka etal. Physical Review B, Vol 41 Num 6, p3854 1990   CORRECTIONS TO PAPER
22//
23////////////////////////////////////////////////////
24
25
26
27Proc PlotSC_ParaCrystal(num,qmin,qmax)
28        Variable num=100, qmin=0.001, qmax=0.7
29        Prompt num "Enter number of data points for model: "
30        Prompt qmin "Enter minimum q-value (^-1) for model: "
31        Prompt qmax "Enter maximum q-value (^-1) for model: "
32//
33        Make/O/D/n=(num) xwave_SC_ParaCrystal, ywave_SC_ParaCrystal
34        xwave_SC_ParaCrystal =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
35        Make/O/D coef_SC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
36        make/o/t parameters_SC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)", "Background (cm-1)"}   
37        Edit parameters_SC_ParaCrystal, coef_SC_ParaCrystal
38       
39        Variable/G root:gNordSC=150
40       
41        Variable/G root:g_SC_ParaCrystal
42        g_SC_ParaCrystal := SC_ParaCrystal(coef_SC_ParaCrystal, ywave_SC_ParaCrystal, xwave_SC_ParaCrystal)
43        Display ywave_SC_ParaCrystal vs xwave_SC_ParaCrystal
44        ModifyGraph marker=29, msize=2, mode=4
45        ModifyGraph grid=1,mirror=2
46        ModifyGraph log=0
47        Label bottom "q (\\S-1\\M) "
48        Label left "I(q) (cm\\S-1\\M)"
49        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
50       
51        AddModelToStrings("SC_ParaCrystal","coef_SC_ParaCrystal","parameters_SC_ParaCrystal","SC_ParaCrystal")
52//
53End
54
55//
56//this macro sets up all the necessary parameters and waves that are
57//needed to calculate the  smeared model function.
58//
59//no input parameters are necessary, it MUST use the experimental q-values
60// from the experimental data read in from an AVE/QSIG data file
61////////////////////////////////////////////////////
62// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
63Proc PlotSmearedSC_ParaCrystal(str)                                                             
64        String str
65        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
66       
67        // if any of the resolution waves are missing => abort
68        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
69                Abort
70        endif
71       
72        SetDataFolder $("root:"+str)
73       
74        // Setup parameter table for model function
75        Make/O/D smear_coef_SC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
76        make/o/t smear_parameters_SC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)", "Background (cm-1)"}
77        Edit smear_parameters_SC_ParaCrystal,smear_coef_SC_ParaCrystal                                  //display parameters in a table
78       
79        // output smeared intensity wave, dimensions are identical to experimental QSIG values
80        // make extra copy of experimental q-values for easy plotting
81        Duplicate/O $(str+"_q") smeared_SC_ParaCrystal,smeared_qvals
82        SetScale d,0,0,"1/cm",smeared_SC_ParaCrystal
83       
84        Variable/G gNordSC = 150       
85        Variable/G gs_SC_ParaCrystal=0
86        gs_SC_ParaCrystal := fSmearedSC_ParaCrystal(smear_coef_SC_ParaCrystal,smeared_SC_ParaCrystal,smeared_qvals)     //this wrapper fills the STRUCT
87       
88        Display smeared_SC_ParaCrystal vs smeared_qvals
89        ModifyGraph marker=29,msize=2,mode=4
90        ModifyGraph log=0
91        Label bottom "q (\\S-1\\M)"
92        Label left "I(q) (cm\\S-1\\M)"
93        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
94       
95        SetDataFolder root:
96        AddModelToStrings("SmearedSC_ParaCrystal","smear_coef_SC_ParaCrystal","smear_parameters_SC_ParaCrystal","SC_ParaCrystal")
97End
98
99
100// nothing to change here
101//
102//AAO version, uses XOP if available
103// simply calls the original single point calculation with
104// a wave assignment (this will behave nicely if given point ranges)
105Function SC_ParaCrystal(cw,yw,xw) : FitFunc
106        Wave cw,yw,xw
107       
108#if exists("SC_ParaCrystalX")
109        yw = SC_ParaCrystalX(cw,xw)
110#else
111        yw = fSC_ParaCrystal(cw,xw)
112#endif
113        return(0)
114End
115
116
117//
118// unsmeared model calculation
119//
120Function fSC_ParaCrystal(w,x) : FitFunc
121        Wave w
122        Variable x
123       
124//       Input (fitting) variables are not used
125//      you would give them nice names
126        Variable integral,loLim,upLim
127        loLim = 0
128        upLim = 2*Pi
129       
130        Variable/G root:gDumY=0         //root:gDumX=0
131       
132
133        Variable scale,Dnn,gg,Rad,contrast,background,yy,latticeScale
134        scale = w[0]
135        Dnn = w[1] //Nearest neighbor distance A
136        gg = w[2] //Paracrystal distortion factor
137        Rad = w[3] //Sphere radius
138        contrast = w[4] - w[5] //SLD contrast
139        background = w[6]
140               
141// always calculate for type 0, SC
142        latticeScale = (4/3)*pi*(Rad^3)/(Dnn^3) //Volume fraction calculated from lattice symmetry and sphere radius
143
144        NVAR/Z nord=root:gNordSC
145        if(NVAR_Exists(nord)!=1)
146                nord=20
147        endif
148       
149        integral = IntegrateFn_N(Integrand_SC_Outer,loLim,upLim,w,x,nord)
150       
151       
152        integral *= SphereForm_SC(Rad,contrast,x)*scale*latticeScale
153        //integral *= scale             //testing, returns Z(q) only
154
155        integral += background 
156       
157        Return (integral)
158       
159End
160
161// the outer integral is also an integral
162Function Integrand_SC_Outer(w,x,dum)
163        Wave w
164        Variable x,dum
165               
166        NVAR yy = root:gDumY           
167        yy = dum                                        // save the current dummy yy for use in the inner loop
168        Variable retVal,loLim,upLim
169        //
170        loLim = 0
171        upLim = Pi
172
173        NVAR/Z nord=root:gNordSC
174        if(NVAR_Exists(nord)!=1)
175                nord=20
176        endif
177       
178        retVal = IntegrateFn_N(Integrand_SC_Inner,loLim,upLim,w,x,nord)
179       
180        return(retVal)
181End
182
183//returns the value of the integrand of the inner integral
184Function Integrand_SC_Inner(w,qq,dum)
185        Wave w
186        Variable qq,dum
187       
188        NVAR yy = root:gDumY            //use the yy value from the outer loop
189        Variable xx,retVal
190        xx = dum
191
192        retVal = SC_Integrand(w,qq,xx,yy)
193       
194        return(retVal)
195End
196
197Function SC_Integrand(w,qq,xx,yy)
198        Wave w
199        Variable qq,xx,yy
200       
201        Variable retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg
202        Dnn = w[1] //Nearest neighbor distance A
203        gg = w[2] //Paracrystal distortion factor
204        aa = Dnn
205        Da = gg*aa
206       
207        temp1 = qq*qq*Da*Da
208        temp2 = (1-exp(-1*temp1))^3
209        temp3 = qq*aa
210        temp4 = 2*exp(-0.5*temp1)
211        temp5 = exp(-1*temp1)
212       
213       
214        retVal = temp2*SCeval(xx,yy,temp3,temp4,temp5)
215        retVal /=4*Pi
216       
217        return(retVal)
218end
219
220Function SCeval(Theta,Phi,temp3,temp4,temp5) //Function to calculate integrand values for simple cubic structure
221        Variable Theta,Phi,temp3,temp4,temp5 //Phi and theta independent parts of the equation. These are passed to the funtion in order to take them off the loop and increase speed
222        Variable temp6,temp7,temp8,temp9 //Theta and phi dependent parts of the equation
223        Variable result
224       
225        temp6 = sin(Theta)
226        temp7 = -1*temp3*sin(Theta)*cos(Phi)
227        temp8 = temp3*sin(Theta)*sin(Phi)
228        temp9 = temp3*cos(Theta)
229        result = temp6/((1-temp4*cos((temp7))+temp5)*(1-temp4*cos((temp8))+temp5)*(1-temp4*cos((temp9))+temp5))
230       
231        return (result)
232end
233
234Function SphereForm_SC(radius,delrho,x)                                 
235        Variable radius,delrho,x
236       
237        // variables are:                                                       
238        //[2] radius ()
239        //[3] delrho (-2)
240        //[4] background (cm-1)
241       
242        // calculates scale * f^2/Vol where f=Vol*3*delrho*(sin(qr)-qrcos(qr))/qr^3
243        // and is rescaled to give [=] cm^-1
244       
245        Variable bes,f,vol,f2
246        ////handle q==0 separately
247        If(x==0)
248                f = 4/3*pi*radius^3*delrho*delrho*1e8
249                return(f)
250        Endif
251       
252       
253        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
254        vol = 4*pi/3*radius^3
255        f = vol*bes*delrho              // [=] 
256        // normalize to single particle volume, convert to 1/cm
257        f2 = f * f / vol * 1.0e8                // [=] 1/cm
258       
259        return (f2)     
260       
261End
262
263
264
265
266///////////////////////////////////////////////////////////////
267// smeared model calculation
268//
269Function SmearedSC_ParaCrystal(s) : FitFunc
270        Struct ResSmearAAOStruct &s
271
272//      the name of your unsmeared model (AAO) is the first argument
273        Smear_Model_76(SC_ParaCrystal,s.coefW,s.xW,s.yW,s.resW)
274
275        return(0)
276End
277
278
279/////////////////////////////////////////////////////////////////
280//wrapper to calculate the smeared model as an AAO-Struct
281// fills the struct and calls the ususal function with the STRUCT parameter
282//
283// used only for the dependency, not for fitting
284//
285Function fSmearedSC_ParaCrystal(coefW,yW,xW)
286        Wave coefW,yW,xW
287       
288        String str = getWavesDataFolder(yW,0)
289        String DF="root:"+str+":"
290       
291        WAVE resW = $(DF+str+"_res")
292       
293        STRUCT ResSmearAAOStruct fs
294        WAVE fs.coefW = coefW   
295        WAVE fs.yW = yW
296        WAVE fs.xW = xW
297        WAVE fs.resW = resW
298       
299        Variable err
300        err = SmearedSC_ParaCrystal(fs)
301       
302        return (0)
303End
Note: See TracBrowser for help on using the repository browser.