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

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

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

24 JUL 2007 SRK

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