source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/StackedDiscs_v40.ipf @ 633

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

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.

File size: 8.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
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//
9// This function calculates the total coherent scattered intensity from stacked discs (tactoids) with a core/layer
10// structure.  Assuming the next neighbor distance (d-spacing) in a stack of parallel discs obeys a Gaussian
11// distribution, a strcture factor S(q) proposed by Kratky and Porod in 1949 is used in this function.
12//
13// 04 JUL 01   DLH
14//
15// SRK - 2007
16// this model needs 76 Gauss points for a proper smearing calculation
17// since there can be sharp interference fringes that develop from the stacking
18////////////////////////////////////////////////
19
20Proc PlotStackedDiscs(num,qmin,qmax)
21        Variable num=500,qmin=0.001,qmax=1.0
22        Prompt num "Enter number of data points for model: "
23        Prompt qmin "Enter minimum q-value (A^-1) for model: "
24        Prompt qmax "Enter maximum q-value (A^-1) for model: "
25       
26        make/o/D/n=(num) xwave_scyl,ywave_scyl
27        xwave_scyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
28        make/o/D coef_scyl = {0.01,3000.,10.,15.,4.0e-6,-4.0e-7,5.0e-6,1,0,1.0e-3}
29        make/o/t parameters_scyl = {"scale","Disc Radius (A)","Core Thickness (A)","Layer Thickness (A)","Core SLD (A^-2)","Layer SLD (A^-2)","Solvent  SLD(A^-2)","# of Stacking","GSD of d-Spacing","incoh. bkg (cm^-1)"}
30        Edit parameters_scyl,coef_scyl
31       
32        Variable/G root:g_scyl
33        g_scyl := StackedDiscs(coef_scyl,ywave_scyl,xwave_scyl)
34//      ywave_scyl := StackedDiscs(coef_scyl,xwave_scyl)
35        Display ywave_scyl vs xwave_scyl
36        ModifyGraph log=1,marker=29,msize=2,mode=4
37        Label bottom "q (A\\S-1\\M)"
38        Label left "Intensity (cm\\S-1\\M)"
39        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
40       
41        AddModelToStrings("StackedDiscs","coef_scyl","parameters_scyl","scyl")
42End
43///////////////////////////////////////////////////////////
44
45///////////////////////////////////////////////////////////
46// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
47Proc PlotSmearedStackedDiscs(str)                                                               
48        String str
49        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
50       
51        // if any of the resolution waves are missing => abort
52        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
53                Abort
54        endif
55       
56        SetDataFolder $("root:"+str)
57       
58        // Setup parameter table for model function
59        make/o/D smear_coef_scyl = {0.01,3000.,10.,15.,4.0e-6,-4.0e-7,5.0e-6,1,0,1.0e-3}
60        make/o/t smear_parameters_scyl = {"scale","Disc Radius (A)","Core Thickness (A)","Layer Thickness (A)","Core SLD (A^-2)","Layer SLD (A^-2)","Solvent SLD (A^-2)","# of Stacking","GSD of d-Spacing","incoh. bkg (cm^-1)"}
61        Edit smear_parameters_scyl,smear_coef_scyl
62       
63        // output smeared intensity wave, dimensions are identical to experimental QSIG values
64        // make extra copy of experimental q-values for easy plotting
65        Duplicate/O $(str+"_q") smeared_scyl,smeared_qvals
66        SetScale d,0,0,"1/cm",smeared_scyl                             
67               
68        Variable/G gs_scyl=0
69        gs_scyl := fSmearedStackedDiscs(smear_coef_scyl,smeared_scyl,smeared_qvals)     //this wrapper fills the STRUCT
70       
71        Display smeared_scyl vs smeared_qvals
72        ModifyGraph log=1,marker=29,msize=2,mode=4
73        Label bottom "q (A\\S-1\\M)"
74        Label left "Intensity (cm\\S-1\\M)"
75        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
76       
77        SetDataFolder root:
78        AddModelToStrings("SmearedStackedDiscs","smear_coef_scyl","smear_parameters_scyl","scyl")
79End
80
81///////////////////////////////////////////////////////////////
82
83//AAO version
84Function StackedDiscs(cw,yw,xw) : FitFunc
85        Wave cw,yw,xw
86       
87#if exists("StackedDiscsX")
88        yw = StackedDiscsX(cw,xw)
89#else
90        yw = fStackedDiscs(cw,xw)
91#endif
92        return(0)
93End
94///////////////////////////////////////////////////////////////
95// unsmeared model calculation
96///////////////////////////
97Function fStackedDiscs(w,x) : FitFunc
98        Wave w
99        Variable x
100
101//The input variables are (and output)
102        //[0] Scale
103        //[1] Disc Radius (A)
104        //[2] Disc Core Thickness (A)
105        //[3] Disc Layer Thickness (A)
106        //[4] Core SLD (A^-2)
107        //[5] Layer SLD (A^-2)
108        //[6] Solvent SLD (A^-2)
109        //[7] Number of Discs Stacked
110        //[8] Gaussian Standrad Deviation of d-Spacing
111        //[9] background (cm^-1)
112       
113        Variable scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd
114        scale = w[0]
115        rcore = w[1]
116        length = w[2]
117        thick = w[3]
118        rhoc = w[4]
119        rhol = w[5]
120        rhosolv = w[6]
121        N = w[7]
122        gsd = w[8]
123        bkg = w[9]
124//
125// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
126//
127
128// local variables
129        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight,kk,sqq,dexpt,d
130        Variable answer
131       
132        d=2*thick+length
133       
134        String weightStr,zStr
135       
136        weightStr = "gauss76wt"
137        zStr = "gauss76z"
138
139       
140//      if wt,z waves don't exist, create them
141// 20 Gauss points is not enough for cylinder calculation
142       
143        if (WaveExists($weightStr) == 0) // wave reference is not valid,
144                Make/D/N=76 $weightStr,$zStr
145                Wave w76 = $weightStr
146                Wave z76 = $zStr                // wave references to pass
147                Make76GaussPoints(w76,z76)     
148        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
149        else
150                if(exists(weightStr) > 1)
151                         Abort "wave name is already in use"    // execute if condition is false
152                endif
153                Wave w76 = $weightStr
154                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
155        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
156        endif
157
158
159// set up the integration
160        // end points and weights
161        nord = 76
162        va = 0
163        vb = Pi/2
164      halfheight = length/2.0
165
166// evaluate at Gauss points
167        // remember to index from 0,size-1
168
169      qq = x            //current x point is the q-value for evaluation
170      summ = 0.0                // initialize integral
171      ii=0
172      do
173                // Using 76 Gauss points
174                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
175                yyy = w76[ii] * Stackdisc_kern(qq, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N)
176                summ += yyy
177
178                ii+=1
179        while (ii<nord)                         // end of loop over quadrature points
180//   
181// calculate value of integral to return
182
183      answer = (vb-va)/2.0*summ
184     
185// contrast is now explicitly included in the core-shell calculation
186
187//Normalize by total disc volume
188//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
189//Calculate TOTAL volume
190// length is the total core thickness
191
192        vcyl=Pi*rcore*rcore*(2*thick+length)*N
193        answer /= vcyl
194
195//Convert to [cm-1]
196        answer *= 1.0e8
197
198//Scale
199        answer *= scale
200       
201// add in the background
202        answer += bkg
203
204        Return (answer)
205       
206End             //End of function StackDiscs()
207
208///////////////////////////////////////////////////////////////
209
210// F(qq, rcore, rhoc,rhosolv, length, zi)
211//
212Function Stackdisc_kern(qq, rcore, rhoc,rhol,rhosolv, length,thick,dum,gsd,d,N)
213        Variable qq, rcore, rhoc,rhol,rhosolv, length,thick,dum,gsd,d,N
214               
215// qq is the q-value for the calculation (1/A)
216// rcore is the core radius of the cylinder (A)
217// rho(n) are the respective SLD's
218// length is the *Half* CORE-LENGTH of the cylinder = L (A)
219// dum is the dummy variable for the integration (x in Feigin's notation)
220
221   //Local variables
222        Variable totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,kk,sqq,dexpt
223        Variable be1,be2,si1,si2
224       
225        dr1 = rhoc-rhosolv
226        dr2 = rhol-rhosolv
227        area = Pi*rcore*rcore
228        totald=2*(thick+length)
229       
230        besarg1 = qq*rcore*sin(dum)
231        besarg2 = qq*rcore*sin(dum)
232       
233        sinarg1 = qq*length*cos(dum)
234        sinarg2 = qq*(length+thick)*cos(dum)
235       
236        if(besarg1 == 0)
237                be1 = 0.5
238        else
239                be1 = bessJ(1,besarg1)/besarg1
240        endif
241        if(besarg2 == 0)
242                be2 = 0.5
243        else
244                be2 = bessJ(1,besarg2)/besarg2
245        endif
246       
247        if(sinarg1 == 0)
248                si1 = 1
249        else
250                si1 = sin(sinarg1)/sinarg1
251        endif
252        if(sinarg2 == 0)
253                si2 = 1
254        else
255                si2 = sin(sinarg2)/sinarg2
256        endif
257       
258        t1 = 2*area*(2*length)*dr1*si1*be1
259        t2 = 2*area*dr2*(totald*si2-2*length*si1)*be2
260       
261        retval =((t1+t2)^2)*sin(dum)
262       
263        // loop for the structure facture S(q)
264       
265        kk=1
266        sqq=0.0
267   do
268                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0
269                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt)
270
271                kk+=1
272        while (kk<N)                           
273       
274        // end of loop for S(q)
275
276        sqq=1.0+2.0*sqq/N
277       
278        retval *= sqq
279   
280    return retval
281   
282End             //Function Stackdisc()
283
284///////////////////////////////////////////////////////////////
285
286// this model needs 76 Gauss points for a proper smearing calculation
287// since there can be sharp interference fringes that develop from the stacking
288Function SmearedStackedDiscs(s) :FitFunc
289        Struct ResSmearAAOStruct &s
290
291////the name of your unsmeared model is the first argument
292        Smear_Model_76(StackedDiscs,s.coefW,s.xW,s.yW,s.resW)
293
294        return(0)
295End
296
297
298//wrapper to calculate the smeared model as an AAO-Struct
299// fills the struct and calls the ususal function with the STRUCT parameter
300//
301// used only for the dependency, not for fitting
302//
303Function fSmearedStackedDiscs(coefW,yW,xW)
304        Wave coefW,yW,xW
305       
306        String str = getWavesDataFolder(yW,0)
307        String DF="root:"+str+":"
308       
309        WAVE resW = $(DF+str+"_res")
310       
311        STRUCT ResSmearAAOStruct fs
312        WAVE fs.coefW = coefW   
313        WAVE fs.yW = yW
314        WAVE fs.xW = xW
315        WAVE fs.resW = resW
316       
317        Variable err
318        err = SmearedStackedDiscs(fs)
319       
320        return (0)
321End
Note: See TracBrowser for help on using the repository browser.