# source:sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/StackedDiscs_v40.ipf@379

Last change on this file since 379 was 379, checked in by srkline, 14 years ago

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

File size: 8.7 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//
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
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:
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
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
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
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
194
195//Convert to [cm-1]
197
198//Scale
200
203
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
224        dr1 = rhoc-rhosolv
225        dr2 = rhol-rhosolv
226        area = Pi*rcore*rcore
227        totald=2*(thick+length)
228
229        besarg1 = qq*rcore*sin(dum)
230        besarg2 = qq*rcore*sin(dum)
231
232        sinarg1 = qq*length*cos(dum)
233        sinarg2 = qq*(length+thick)*cos(dum)
234
235        t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(bessJ(1,besarg1)/besarg1)
236        t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(bessJ(1,besarg2)/besarg2)
237
238        retval =((t1+t2)^2)*sin(dum)
239
240        // loop for the structure facture S(q)
241
242             kk=1
243             sqq=0.0
244      do
245                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0
246                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt)
247
248                kk+=1
249        while (kk<N)
250
251        // end of loop for S(q)
252
253        sqq=1.0+2.0*sqq/N
254
255        retval *= sqq
256
257    return retval
258
259End             //Function Stackdisc()
260
261///////////////////////////////////////////////////////////////
262
263// this model needs 76 Gauss points for a proper smearing calculation
264// since there can be sharp interference fringes that develop from the stacking
265Function SmearedStackedDiscs(s) :FitFunc
266        Struct ResSmearAAOStruct &s
267
268////the name of your unsmeared model is the first argument
269        Smear_Model_76(StackedDiscs,s.coefW,s.xW,s.yW,s.resW)
270
271        return(0)
272End
273
274
275//wrapper to calculate the smeared model as an AAO-Struct
276// fills the struct and calls the ususal function with the STRUCT parameter
277//
278// used only for the dependency, not for fitting
279//
280Function fSmearedStackedDiscs(coefW,yW,xW)
281        Wave coefW,yW,xW
282
283        String str = getWavesDataFolder(yW,0)
284        String DF="root:"+str+":"
285
286        WAVE resW = \$(DF+str+"_res")
287
288        STRUCT ResSmearAAOStruct fs
289        WAVE fs.coefW = coefW
290        WAVE fs.yW = yW
291        WAVE fs.xW = xW
292        WAVE fs.resW = resW
293
294        Variable err
295        err = SmearedStackedDiscs(fs)
296
297        return (0)
298End
Note: See TracBrowser for help on using the repository browser.