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

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

Typo in dialog in every model...

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 PlotCoreShellCylinderForm(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 PlotSmearedCSCylinderForm(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 := fSmearedCoreShellCylinderForm(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 SmearedCoreShellCylinderForm(s) :FitFunc
226        Struct ResSmearAAOStruct &s
227
228////the name of your unsmeared model is the first argument
229        s.yW = Smear_Model_20(CoreShellCylinder,s.coefW,s.xW,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 fSmearedCoreShellCylinderForm(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 = SmearedCoreShellCylinderForm(fs)
256       
257        return (0)
258End
Note: See TracBrowser for help on using the repository browser.