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

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

File size: 7.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// this function is for the form factor of a right circular cylinder with core/shell scattering length density profile
9//
10// the core dimensions are given and a constant shell thickness is added to the radius and to dach and of the length
11// this way, the scattering amplitude is simply the difference between two cylinders of different dimensions
12//
13// 06 NOV 98 SRK
14////////////////////////////////////////////////
15
16Proc PlotCoreShellCylinder(num,qmin,qmax)
17        Variable num=128,qmin=0.001,qmax=0.7
18        Prompt num "Enter number of data points for model: "
19        Prompt qmin "Enter minimum q-value (^-1) for model: "
20        Prompt qmax "Enter maximum q-value (^-1) for model: "
21       
22        make/o/d/n=(num) xwave_cscyl,ywave_cscyl
23        xwave_cscyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
24        make/o/d coef_cscyl = {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01}
25        make/o/t parameters_cscyl = {"scale","core radius (A)","shell THICKNESS (A)","CORE length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
26        Edit parameters_cscyl,coef_cscyl
27        Variable/G root:g_cscyl
28        g_cscyl := CoreShellCylinder(coef_cscyl,ywave_cscyl,xwave_cscyl)
29//      ywave_cscyl := CoreShellCylinder(coef_cscyl,xwave_cscyl)
30        Display ywave_cscyl vs xwave_cscyl
31        ModifyGraph log=1,marker=29,msize=2,mode=4
32        Label bottom "q (\\S-1\\M)"
33        Label left "Intensity (cm\\S-1\\M)"
34        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
35End
36
37///////////////////////////////////////////////////////////
38// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
39Proc PlotSmearedCoreShellCylinder(str)                                                         
40        String str
41        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
42       
43        // if any of the resolution waves are missing => abort
44        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
45                Abort
46        endif
47       
48        SetDataFolder $("root:"+str)
49       
50        // Setup parameter table for model function
51        make/o/d smear_coef_cscyl =  {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01}
52        make/o/t smear_parameters_cscyl = {"scale","core radius (A)","shell radius (A)","length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
53        Edit smear_parameters_cscyl,smear_coef_cscyl
54       
55        // output smeared intensity wave, dimensions are identical to experimental QSIG values
56        // make extra copy of experimental q-values for easy plotting
57        Duplicate/O $(str+"_q") smeared_cscyl,smeared_qvals
58        SetScale d,0,0,"1/cm",smeared_cscyl     
59                                       
60        Variable/G gs_cscyl=0
61        gs_cscyl := fSmearedCoreShellCylinder(smear_coef_cscyl,smeared_cscyl,smeared_qvals)     //this wrapper fills the STRUCT
62       
63        Display smeared_cscyl vs smeared_qvals
64        ModifyGraph log=1,marker=29,msize=2,mode=4
65        Label bottom "q (\\S-1\\M)"
66        Label left "Intensity (cm\\S-1\\M)"
67        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
68       
69        SetDataFolder root:
70End
71       
72
73//AAO version
74Function CoreShellCylinder(cw,yw,xw) : FitFunc
75        Wave cw,yw,xw
76
77#if exists("CoreShellCylinderX")
78        yw = CoreShellCylinderX(cw,xw)
79#else
80        yw = fCoreShellCylinder(cw,xw)
81#endif
82        return(0)
83End
84
85///////////////////////////////////////////////////////////////
86// unsmeared model calculation
87///////////////////////////
88Function fCoreShellCylinder(w,x) : FitFunc
89        Wave w
90        Variable x
91
92//The input variables are (and output)
93        //[0] scale
94        //[1] cylinder CORE RADIUS (A)
95        //[2] shell Thickness (A)
96        //[3]  cylinder CORE LENGTH (A)
97        //[4] core SLD (A^-2)
98        //[5] shell SLD (A^-2)
99        //[6] solvent SLD (A^-2)
100        //[7] background (cm^-1)
101        Variable scale,length,delrho,bkg,rcore,thick,rhoc,rhos,rhosolv
102        scale = w[0]
103        rcore = w[1]
104        thick = w[2]
105        length = w[3]
106        rhoc = w[4]
107        rhos = w[5]
108        rhosolv = w[6]
109        bkg = w[7]
110//
111// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
112//
113
114// local variables
115        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
116        Variable answer
117        String weightStr,zStr
118       
119        weightStr = "gauss76wt"
120        zStr = "gauss76z"
121
122       
123//      if wt,z waves don't exist, create them
124// 20 Gauss points is not enough for cylinder calculation
125       
126        if (WaveExists($weightStr) == 0) // wave reference is not valid,
127                Make/D/N=76 $weightStr,$zStr
128                Wave w76 = $weightStr
129                Wave z76 = $zStr                // wave references to pass
130                Make76GaussPoints(w76,z76)     
131        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
132        else
133                if(exists(weightStr) > 1)
134                         Abort "wave name is already in use"    // execute if condition is false
135                endif
136                Wave w76 = $weightStr
137                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
138        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
139        endif
140
141
142// set up the integration
143        // end points and weights
144        nord = 76
145        va = 0
146        vb = Pi/2
147      halfheight = length/2.0
148
149// evaluate at Gauss points
150        // remember to index from 0,size-1
151
152        qq = x          //current x point is the q-value for evaluation
153      summ = 0.0                // initialize integral
154      ii=0
155      do
156                // Using 76 Gauss points
157                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
158                yyy = w76[ii] * CoreShellcyl(qq, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi)
159                summ += yyy
160
161                ii+=1
162        while (ii<nord)                         // end of loop over quadrature points
163//   
164// calculate value of integral to return
165
166      answer = (vb-va)/2.0*summ
167     
168// contrast is now explicitly included in the core-shell calculation
169
170//normalize by cylinder volume
171//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
172//calculate TOTAL volume
173// length is the total core length
174        vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2*thick)
175        answer /= vcyl
176//convert to [cm-1]
177        answer *= 1.0e8
178//Scale
179        answer *= scale
180// add in the background
181        answer += bkg
182
183        Return (answer)
184       
185End             //End of function CoreShellCylinderForm()
186
187///////////////////////////////////////////////////////////////
188// F(qq, rcore, thick, rhoc,rhos,rhosolv, length, zi)
189//
190Function CoreShellcyl(qq, rcore, thick, rhoc,rhos,rhosolv, length, dum)
191        Variable qq, rcore, thick, rhoc,rhos,rhosolv, length, dum
192       
193// qq is the q-value for the calculation (1/A)
194// rcore is the core radius of the cylinder (A)
195//thick is the uniform thickness
196// rho(n) are the respective SLD's
197
198// length is the *Half* CORE-LENGTH of the cylinder = L (A)
199
200// dum is the dummy variable for the integration (x in Feigin's notation)
201
202   //Local variables
203        Variable dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval
204       
205        dr1 = rhoc-rhos
206        dr2 = rhos-rhosolv
207        vol1 = Pi*rcore*rcore*(2*length)
208        vol2 = Pi*(rcore+thick)*(rcore+thick)*(2*length+2*thick)
209       
210        besarg1 = qq*rcore*sin(dum)
211        besarg2 = qq*(rcore+thick)*sin(dum)
212        sinarg1 = qq*length*cos(dum)
213        sinarg2 = qq*(length+thick)*cos(dum)
214       
215        t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1
216        t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2
217       
218        retval = ((t1+t2)^2)*sin(dum)
219       
220    return retval
221   
222End     //Function CoreShellcyl()
223
224// this is all there is to the smeared calculation!
225Function SmearedCoreShellCylinder(s) :FitFunc
226        Struct ResSmearAAOStruct &s
227
228////the name of your unsmeared model is the first argument
229        Smear_Model_20(CoreShellCylinder,s.coefW,s.xW,s.yW,s.resW)
230
231        return(0)
232End
233
234
235//wrapper to calculate the smeared model as an AAO-Struct
236// fills the struct and calls the ususal function with the STRUCT parameter
237//
238// used only for the dependency, not for fitting
239//
240Function fSmearedCoreShellCylinder(coefW,yW,xW)
241        Wave coefW,yW,xW
242       
243        String str = getWavesDataFolder(yW,0)
244        String DF="root:"+str+":"
245       
246        WAVE resW = $(DF+str+"_res")
247       
248        STRUCT ResSmearAAOStruct fs
249        WAVE fs.coefW = coefW   
250        WAVE fs.yW = yW
251        WAVE fs.xW = xW
252        WAVE fs.resW = resW
253       
254        Variable err
255        err = SmearedCoreShellCylinder(fs)
256       
257        return (0)
258End
Note: See TracBrowser for help on using the repository browser.