# source:sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/StackedDiscs.ipf@178

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

Changed StackedDiscs? and MultiShell? to use 76 point quadrature. I though this was done long ago, but the RT ticket was still open. Now it's done. The branch also has these changes.

File size: 7.5 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////
4// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile.
5// Adopting these into the experiment will insure that they are always present.
6////////////////////////////////////////////////
7//
8// This function calculates the total coherent scattered intensity from stacked discs (tactoids) with a core/layer
9// structure.  Assuming the next neighbor distance (d-spacing) in a stack of parallel discs obeys a Gaussian
10// distribution, a strcture factor S(q) proposed by Kratky and Porod in 1949 is used in this function.
11//
12// 04 JUL 01   DLH
13////////////////////////////////////////////////
14
15Proc PlotStackedDiscs(num,qmin,qmax)
16        Variable num=500,qmin=0.001,qmax=1.0
17        Prompt num "Enter number of data points for model: "
18        Prompt qmin "Enter minimum q-value (^-1) for model: "
19        Prompt qmax "Enter maximum q-value (^-1) for model: "
20
21        make/o/n=(num) xwave_scyl,ywave_scyl
22        xwave_scyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
23        make/o coef_scyl = {0.01,3000.,10.,15.,4.0e-6,-4.0e-7,5.0e-6,1,0,1.0e-3}
24        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)"}
25        Edit parameters_scyl,coef_scyl
26        ywave_scyl := StackedDiscs(coef_scyl,xwave_scyl)
27        Display ywave_scyl vs xwave_scyl
28        ModifyGraph log=1,marker=29,msize=2,mode=4
29        Label bottom "q (\\S-1\\M)"
30        Label left "Intensity (cm\\S-1\\M)"
31        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
32End
33///////////////////////////////////////////////////////////
34
35///////////////////////////////////////////////////////////
36Proc PlotSmearedStackedDiscs()
37        //no input parameters necessary, it MUST use the experimental q-values
38        // from the experimental data read in from an AVE/QSIG data file
39
40        // if no gQvals wave, data must not have been loaded => abort
41        if(ResolutionWavesMissing())
42                Abort
43        endif
44
45        // Setup parameter table for model function
46        make/o smear_coef_scyl = {0.01,3000.,10.,15.,4.0e-6,-4.0e-7,5.0e-6,1,0,1.0e-3}
47        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)"}
48        Edit smear_parameters_scyl,smear_coef_scyl
49
50        // output smeared intensity wave, dimensions are identical to experimental QSIG values
51        // make extra copy of experimental q-values for easy plotting
52        Duplicate/O \$gQvals smeared_scyl,smeared_qvals
53        SetScale d,0,0,"1/cm",smeared_scyl
54
55        smeared_scyl := SmearedStackedDiscs(smear_coef_scyl,\$gQvals)
56        Display smeared_scyl vs smeared_qvals
57        ModifyGraph log=1,marker=29,msize=2,mode=4
58        Label bottom "q (\\S-1\\M)"
59        Label left "Intensity (cm\\S-1\\M)"
60        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
61End
62///////////////////////////////////////////////////////////////
63
64///////////////////////////////////////////////////////////////
65// unsmeared model calculation
66///////////////////////////
67Function StackedDiscs(w,x) : FitFunc
68        Wave w
69        Variable x
70
71//The input variables are (and output)
72        //[0] Scale
74        //[2] Disc Core Thickness (A)
75        //[3] Disc Layer Thickness (A)
76        //[4] Core SLD (A^-2)
77        //[5] Layer SLD (A^-2)
78        //[6] Solvent SLD (A^-2)
79        //[7] Number of Discs Stacked
80        //[8] Gaussian Standrad Deviation of d-Spacing
81        //[9] background (cm^-1)
82
83        Variable scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd
84        scale = w[0]
85        rcore = w[1]
86        length = w[2]
87        thick = w[3]
88        rhoc = w[4]
89        rhol = w[5]
90        rhosolv = w[6]
91        N = w[7]
92        gsd = w[8]
93        bkg = w[9]
94//
95// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
96//
97
98// local variables
99        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight,kk,sqq,dexpt,d
101
102        d=2*thick+length
103
104        String weightStr,zStr
105
106        weightStr = "gauss76wt"
107        zStr = "gauss76z"
108
109
110//      if wt,z waves don't exist, create them
111// 20 Gauss points is not enough for cylinder calculation
112
113        if (WaveExists(\$weightStr) == 0) // wave reference is not valid,
114                Make/D/N=76 \$weightStr,\$zStr
115                Wave w76 = \$weightStr
116                Wave z76 = \$zStr                // wave references to pass
117                Make76GaussPoints(w76,z76)
118        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
119        else
120                if(exists(weightStr) > 1)
121                         Abort "wave name is already in use"    // execute if condition is false
122                endif
123                Wave w76 = \$weightStr
124                Wave z76 = \$zStr                // Not sure why this has to be "declared" twice
125        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
126        endif
127
128
129// set up the integration
130        // end points and weights
131        nord = 76
132        va = 0
133        vb = Pi/2
134      halfheight = length/2.0
135
136// evaluate at Gauss points
137        // remember to index from 0,size-1
138
139      qq = x            //current x point is the q-value for evaluation
140      summ = 0.0                // initialize integral
141      ii=0
142      do
143                // Using 76 Gauss points
144                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
145                yyy = w76[ii] * Stackdisc_kern(qq, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N)
146                summ += yyy
147
148                ii+=1
149        while (ii<nord)                         // end of loop over quadrature points
150//
151// calculate value of integral to return
152
154
155// contrast is now explicitly included in the core-shell calculation
156
157//Normalize by total disc volume
158//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
159//Calculate TOTAL volume
160// length is the total core thickness
161
162        vcyl=Pi*rcore*rcore*(2*thick+length)*N
164
165//Convert to [cm-1]
167
168//Scale
170
173
175
176End             //End of function StackDiscs()
177
178///////////////////////////////////////////////////////////////
179
180// F(qq, rcore, rhoc,rhosolv, length, zi)
181//
182Function Stackdisc_kern(qq, rcore, rhoc,rhol,rhosolv, length,thick,dum,gsd,d,N)
183        Variable qq, rcore, rhoc,rhol,rhosolv, length,thick,dum,gsd,d,N
184
185// qq is the q-value for the calculation (1/A)
186// rcore is the core radius of the cylinder (A)
187// rho(n) are the respective SLD's
188// length is the *Half* CORE-LENGTH of the cylinder = L (A)
189// dum is the dummy variable for the integration (x in Feigin's notation)
190
191   //Local variables
192        Variable totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,kk,sqq,dexpt
193
194        dr1 = rhoc-rhosolv
195        dr2 = rhol-rhosolv
196        area = Pi*rcore*rcore
197        totald=2*(thick+length)
198
199        besarg1 = qq*rcore*sin(dum)
200        besarg2 = qq*rcore*sin(dum)
201
202        sinarg1 = qq*length*cos(dum)
203        sinarg2 = qq*(length+thick)*cos(dum)
204
205        t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(bessJ(1,besarg1)/besarg1)
206        t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(bessJ(1,besarg2)/besarg2)
207
208        retval =((t1+t2)^2)*sin(dum)
209
210        // loop for the structure facture S(q)
211
212             kk=1
213             sqq=0.0
214      do
215                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0
216                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt)
217
218                kk+=1
219        while (kk<N)
220
221        // end of loop for S(q)
222
223        sqq=1.0+2.0*sqq/N
224
225        retval *= sqq
226
227    return retval
228
229End             //Function Stackdisc()
230
231///////////////////////////////////////////////////////////////
232
233// this model needs 76 Gauss points for a proper smearing calculation
234// since there can be sharp interference fringes that develop from the stacking
235Function SmearedStackedDiscs(w,x) :FitFunc
236        Wave w
237        Variable x
238
239        Variable ans
240        SVAR sq = gSig_Q
241        SVAR qb = gQ_bar