source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/StackedDiscs_v40.ipf @ 510

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

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","parameters_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","smear_parameters_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.