source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/HollowCylinders.ipf @ 42

Last change on this file since 42 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files
File size: 5.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////
4// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
5// Adopting these into the experiment will insure that they are always present
6////////////////////////////////////////////////
7// this function is for the form factor of a right circular cylinder with uniform scattering length density
8//
9// 06 NOV 98 SRK
10////////////////////////////////////////////////
11
12Proc PlotHollowCylinderForm(num,qmin,qmax)
13        Variable num=128,qmin=0.001,qmax=0.7
14        Prompt num "Enter number of data points for model: "
15        Prompt qmin "Enter minimum q-value (^-1) for model: "
16        Prompt qmax "Enter maximum q-value (^-1) for model: "
17       
18        Make/O/D/n=(num) xwave_Hcyl,ywave_Hcyl
19        xwave_Hcyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
20        Make/O/D coef_Hcyl = {1.,20.,30.,400,3.0e-6,0.01}
21        make/o/t parameters_Hcyl = {"scale","core radius (A)","shell radius (A)","length (A)","contrast (A^-2)","incoh. bkg (cm^-1)"}
22        Edit parameters_Hcyl,coef_Hcyl
23        ywave_Hcyl := HollowCylinderForm(coef_Hcyl,xwave_Hcyl)
24        Display ywave_Hcyl vs xwave_Hcyl
25        ModifyGraph log=1,marker=29,msize=2,mode=4
26        Label bottom "q (\\S-1\\M)"
27        Label left "Intensity (cm\\S-1\\M)"
28        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
29End
30///////////////////////////////////////////////////////////
31
32Proc PlotSmearedHollowCylinderForm()   
33        //no input parameters necessary, it MUST use the experimental q-values
34        // from the experimental data read in from an AVE/QSIG data file
35       
36        // if no gQvals wave, data must not have been loaded => abort
37        if(ResolutionWavesMissing())
38                Abort
39        endif
40       
41        // Setup parameter table for model function
42        Make/O/D smear_coef_Hcyl = {1.,20.,30.,400,3.0e-6,0.01}
43        make/o/t smear_parameters_Hcyl = {"scale","core radius (A)","shell radius (A)","length (A)","contrast (A^-2)","incoh. bkg (cm^-1)"}
44        Edit smear_parameters_Hcyl,smear_coef_Hcyl
45       
46        // output smeared intensity wave, dimensions are identical to experimental QSIG values
47        // make extra copy of experimental q-values for easy plotting
48        Duplicate/O $gQvals smeared_Hcyl,smeared_qvals
49        SetScale d,0,0,"1/cm",smeared_Hcyl     
50
51        smeared_Hcyl := SmearedHollowCylinderForm(smear_coef_Hcyl,$gQvals)
52        Display smeared_Hcyl vs smeared_qvals
53        ModifyGraph log=1,marker=29,msize=2,mode=4
54        Label bottom "q (\\S-1\\M)"
55        Label left "Intensity (cm\\S-1\\M)"
56        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
57End
58
59///////////////////////////////////////////////////////////////
60// unsmeared model calculation
61///////////////////////////
62Function HollowCylinderForm(w,x) : FitFunc
63        Wave w
64        Variable x
65
66//The input variables are (and output)
67        //[0] scale
68        //[1] cylinder CORE RADIUS (A)
69        //[2] cylinder shell radius (A)
70        //[3] total cylinder LENGTH (A)
71        //[4] contrast (A^-2)
72        //[5] background (cm^-1)
73        Variable scale,length,delrho,bkg,rcore,rshell,contrast
74        scale = w[0]
75        rcore = w[1]
76        rshell = w[2]
77        length = w[3]
78        contrast = w[4]
79        bkg = w[5]
80//
81// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
82//
83
84// local variables
85        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
86        Variable answer
87        String weightStr,zStr
88       
89        weightStr = "gauss76wt"
90        zStr = "gauss76z"
91
92       
93//      if wt,z waves don't exist, create them
94// 20 Gauss points is not enough for cylinder calculation
95       
96        if (WaveExists($weightStr) == 0) // wave reference is not valid,
97                Make/D/N=76 $weightStr,$zStr
98                Wave w76 = $weightStr
99                Wave z76 = $zStr                // wave references to pass
100                Make76GaussPoints(w76,z76)     
101        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
102        else
103                if(exists(weightStr) > 1)
104                         Abort "wave name is already in use"    // execute if condition is false
105                endif
106                Wave w76 = $weightStr
107                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
108        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
109        endif
110
111
112// set up the integration
113        // end points and weights
114        nord = 76
115        va = 0
116        vb = 1
117      halfheight = length/2.0
118
119// evaluate at Gauss points
120        // remember to index from 0,size-1
121
122        qq = x          //current x point is the q-value for evaluation
123      summ = 0.0                // initialize integral
124      ii=0
125      do
126                // Using 76 Gauss points
127                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
128                yyy = w76[ii] * Hollowcyl(qq, rcore, rshell, length, zi)
129                summ += yyy
130
131                ii+=1
132        while (ii<nord)                         // end of loop over quadrature points
133//   
134// calculate value of integral to return
135
136      answer = (vb-va)/2.0*summ
137     
138// multiply by the contrast
139        answer *= contrast*contrast
140
141//normalize by cylinder volume
142//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
143        vcyl=Pi*(rshell^2-rcore^2)*length
144        answer *= vcyl
145//convert to [cm-1]
146        answer *= 1.0e8
147//Scale
148        answer *= scale
149// add in the background
150        answer += bkg
151
152        Return (answer)
153       
154End             //End of function HollowCylinderForm()
155
156///////////////////////////////////////////////////////////////
157Function Hollowcyl(qq,r2,r1,h,theta)
158        Variable qq,r2,r1,h,theta
159       
160// qq is the q-value for the calculation (1/A)
161// r2 is the core radius of the cylinder (A)
162//r1 is the shell raduis
163// rho(n) are the respective SLD's
164// h is the total-LENGTH of the cylinder = L (A)
165// theta is the dummy variable for the integration (x in Feigin's notation)
166
167   //Local variables
168        Variable gamma,besarg1,besarg2,lam1,lam2,t2,retval,psi,sinarg
169       
170        gamma = r2/r1
171        besarg1 = qq*r1*sqrt(1-theta^2)
172        besarg2 = qq*r2*sqrt(1-theta^2)
173        lam1 = 2*bessJ(1,besarg1)/besarg1
174        lam2 = 2*bessJ(1,besarg2)/besarg2
175        psi = 1/(1-gamma^2)*(lam1 -  gamma^2*lam2)              //SRK 10/19/00
176       
177        sinarg = qq*h*theta/2
178        t2 = sin(sinarg)/sinarg
179       
180        retval = psi*psi*t2*t2
181       
182    return retval
183   
184End     //Function Hollowcyl()
185
186
187// this is all there is to the smeared calculation!
188Function SmearedHollowCylinderForm(w,x) :FitFunc
189        Wave w
190        Variable x
191       
192        Variable ans
193        SVAR sq = gSig_Q
194        SVAR qb = gQ_bar
195        SVAR sh = gShadow
196        SVAR gQ = gQVals
197       
198        //the name of your unsmeared model is the first argument
199        ans = Smear_Model_20(HollowCylinderForm,$sq,$qb,$sh,$gQ,w,x)
200
201        return(ans)
202End
Note: See TracBrowser for help on using the repository browser.