source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models_v3.00/CoreShellCylinder.ipf @ 381

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

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

File size: 6.8 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 core/shell scattering length density profile
8//
9// the core dimensions are given and a constant shell thickness is added to the radius and to dach and of the length
10// this way, the scattering amplitude is simply the difference between two cylinders of different dimensions
11//
12// 06 NOV 98 SRK
13////////////////////////////////////////////////
14
15Proc PlotCoreShellCylinderForm(num,qmin,qmax)
16        Variable num=128,qmin=0.001,qmax=0.7
17        Prompt num "Enter number of data points for model: "
18        Prompt qmin "Enter minimum q-value (A^-1) for model: "
19        Prompt qmax "Enter maximum q-value (A^-1) for model: "
20       
21        make/o/d/n=(num) xwave_cscyl,ywave_cscyl
22        xwave_cscyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
23        make/o/d coef_cscyl = {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01}
24        make/o/t parameters_cscyl = {"scale","core radius (A)","shell THICKNESS (A)","CORE length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
25        Edit parameters_cscyl,coef_cscyl
26        ywave_cscyl := CoreShellCylinderForm(coef_cscyl,xwave_cscyl)
27        Display ywave_cscyl vs xwave_cscyl
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)
32End
33///////////////////////////////////////////////////////////
34
35Proc PlotSmearedCSCylinderForm()       
36        //no input parameters necessary, it MUST use the experimental q-values
37        // from the experimental data read in from an AVE/QSIG data file
38       
39        // if no gQvals wave, data must not have been loaded => abort
40        if(ResolutionWavesMissing())
41                Abort
42        endif
43       
44        // Setup parameter table for model function
45        make/o/d smear_coef_cscyl =  {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01}
46        make/o/t smear_parameters_cscyl = {"scale","core radius (A)","shell radius (A)","length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
47        Edit smear_parameters_cscyl,smear_coef_cscyl
48       
49        // output smeared intensity wave, dimensions are identical to experimental QSIG values
50        // make extra copy of experimental q-values for easy plotting
51        Duplicate/O $gQvals smeared_cscyl,smeared_qvals
52        SetScale d,0,0,"1/cm",smeared_cscyl     
53
54        smeared_cscyl := SmearedCoreShellCylinderForm(smear_coef_cscyl,$gQvals)
55        Display smeared_cscyl vs smeared_qvals
56        ModifyGraph log=1,marker=29,msize=2,mode=4
57        Label bottom "q (A\\S-1\\M)"
58        Label left "Intensity (cm\\S-1\\M)"
59        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
60End
61
62///////////////////////////////////////////////////////////////
63// unsmeared model calculation
64///////////////////////////
65Function CoreShellCylinderForm(w,x) : FitFunc
66        Wave w
67        Variable x
68
69//The input variables are (and output)
70        //[0] scale
71        //[1] cylinder CORE RADIUS (A)
72        //[2] shell Thickness (A)
73        //[3]  cylinder CORE LENGTH (A)
74        //[4] core SLD (A^-2)
75        //[5] shell SLD (A^-2)
76        //[6] solvent SLD (A^-2)
77        //[7] background (cm^-1)
78        Variable scale,length,delrho,bkg,rcore,thick,rhoc,rhos,rhosolv
79        scale = w[0]
80        rcore = w[1]
81        thick = w[2]
82        length = w[3]
83        rhoc = w[4]
84        rhos = w[5]
85        rhosolv = w[6]
86        bkg = w[7]
87//
88// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
89//
90
91// local variables
92        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
93        Variable answer
94        String weightStr,zStr
95       
96        weightStr = "gauss76wt"
97        zStr = "gauss76z"
98
99       
100//      if wt,z waves don't exist, create them
101// 20 Gauss points is not enough for cylinder calculation
102       
103        if (WaveExists($weightStr) == 0) // wave reference is not valid,
104                Make/D/N=76 $weightStr,$zStr
105                Wave w76 = $weightStr
106                Wave z76 = $zStr                // wave references to pass
107                Make76GaussPoints(w76,z76)     
108        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
109        else
110                if(exists(weightStr) > 1)
111                         Abort "wave name is already in use"    // execute if condition is false
112                endif
113                Wave w76 = $weightStr
114                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
115        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
116        endif
117
118
119// set up the integration
120        // end points and weights
121        nord = 76
122        va = 0
123        vb = Pi/2
124      halfheight = length/2.0
125
126// evaluate at Gauss points
127        // remember to index from 0,size-1
128
129        qq = x          //current x point is the q-value for evaluation
130      summ = 0.0                // initialize integral
131      ii=0
132      do
133                // Using 76 Gauss points
134                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
135                yyy = w76[ii] * CoreShellcyl(qq, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi)
136                summ += yyy
137
138                ii+=1
139        while (ii<nord)                         // end of loop over quadrature points
140//   
141// calculate value of integral to return
142
143      answer = (vb-va)/2.0*summ
144     
145// contrast is now explicitly included in the core-shell calculation
146
147//normalize by cylinder volume
148//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
149//calculate TOTAL volume
150// length is the total core length
151        vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2*thick)
152        answer /= vcyl
153//convert to [cm-1]
154        answer *= 1.0e8
155//Scale
156        answer *= scale
157// add in the background
158        answer += bkg
159
160        Return (answer)
161       
162End             //End of function CoreShellCylinderForm()
163
164///////////////////////////////////////////////////////////////
165// F(qq, rcore, thick, rhoc,rhos,rhosolv, length, zi)
166//
167Function CoreShellcyl(qq, rcore, thick, rhoc,rhos,rhosolv, length, dum)
168        Variable qq, rcore, thick, rhoc,rhos,rhosolv, length, dum
169       
170// qq is the q-value for the calculation (1/A)
171// rcore is the core radius of the cylinder (A)
172//thick is the uniform thickness
173// rho(n) are the respective SLD's
174
175// length is the *Half* CORE-LENGTH of the cylinder = L (A)
176
177// dum is the dummy variable for the integration (x in Feigin's notation)
178
179   //Local variables
180        Variable dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval
181       
182        dr1 = rhoc-rhos
183        dr2 = rhos-rhosolv
184        vol1 = Pi*rcore*rcore*(2*length)
185        vol2 = Pi*(rcore+thick)*(rcore+thick)*(2*length+2*thick)
186       
187        besarg1 = qq*rcore*sin(dum)
188        besarg2 = qq*(rcore+thick)*sin(dum)
189        sinarg1 = qq*length*cos(dum)
190        sinarg2 = qq*(length+thick)*cos(dum)
191       
192        t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1
193        t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2
194       
195        retval = ((t1+t2)^2)*sin(dum)
196       
197    return retval
198   
199End     //Function CoreShellcyl()
200
201// this is all there is to the smeared calculation!
202Function SmearedCoreShellCylinderForm(w,x) :FitFunc
203        Wave w
204        Variable x
205       
206        Variable ans
207        SVAR sq = gSig_Q
208        SVAR qb = gQ_bar
209        SVAR sh = gShadow
210        SVAR gQ = gQVals
211       
212        //the name of your unsmeared model is the first argument
213        ans = Smear_Model_20(CoreShellCylinderForm,$sq,$qb,$sh,$gQ,w,x)
214
215        return(ans)
216End
217
Note: See TracBrowser for help on using the repository browser.