source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/HollowCylinders.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: 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 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 (^-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 PlotSmearedHollowCylinder(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.