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

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

(1) - cursors can now be used to select a subrange of USANS data to fit. This is done by th fit wrapper, assigning a subrange of resW to the struct

(2) all of the smeared model functions are now in the latest form of Smear_Model_N() that is NOT a pointwise calculation anymore, since the USANS matrix smearing in inherently not so.

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