source: sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/HollowCylinders_v40.ipf @ 379

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

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

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