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

Last change on this file since 42 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files

File size: 7.4 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
73        //[1] Disc Radius (A)
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
100        Variable answer
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
153      answer = (vb-va)/2.0*summ
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
163        answer /= vcyl
164
165//Convert to [cm-1]
166        answer *= 1.0e8
167
168//Scale
169        answer *= scale
170       
171// add in the background
172        answer += bkg
173
174        Return (answer)
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 is all there is to the smeared calculation!
234Function SmearedStackedDiscs(w,x) :FitFunc
235        Wave w
236        Variable x
237       
238        Variable ans
239        SVAR sq = gSig_Q
240        SVAR qb = gQ_bar
241        SVAR sh = gShadow
242        SVAR gQ = gQVals
243       
244        //the name of your unsmeared model is the first argument
245        ans = Smear_Model_20(StackedDiscs,$sq,$qb,$sh,$gQ,w,x)
246
247        return(ans)
248End
249
Note: See TracBrowser for help on using the repository browser.