source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/StackedDiscs_v40.ipf @ 381

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

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

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       
41        AddModelToStrings("StackedDiscs","coef_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","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       
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.