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

Last change on this file since 252 was 252, checked in by ajj, 15 years ago

Renaming files with _v40 - see trac ticket #108

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