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

Last change on this file since 510 was 510, checked in by srkline, 14 years ago

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

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","parameters_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","smear_parameters_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.