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

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

Updated StackedDiscs? model to Make double precision waves for XOP

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