source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/PolyCoreShellCylinder.ipf @ 131

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

Typo in dialog in every model...

File size: 10.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// This function is for the form factor of a right circular
6// cylinder with core/shell scattering length density profile.
7// Note that a different shell thickness is added on the edge of
8// the particle, compared to the face. 
9// Furthermore the scattering is convoluted by a log normal (or Schultz)
10// distribution that creates polydispersity for the radius of the
11// particle core.
12//
13// 30 Apr 2003 Andrew Nelson
14//
15// 17 MAY 2006 SRK - changed to normalize to total particle dimensions
16//                                              (core+shell)
17//
18// The Gaussian quadrature routines are based on those in the
19// current NIST macros.
20/////////////////////////////////////////////////////////////////
21
22Proc PlotPolyCoShCylinder(num,qmin,qmax)
23        Variable num=100,qmin=0.001,qmax=0.7
24        Prompt num "Enter number of data points for model: "
25        Prompt qmin "Enter minimum q-value (A^-1) for model: "
26        Prompt qmax "Enter maximum q-value (A^-1) for model: "
27       
28        make/o/d/n=(num) xwave_CSCpr,ywave_CSCpr
29        xwave_CSCpr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
30        make/o/d coef_CSCpr = {0.01,150,0.10,10,20.,10.,4.0e-6,1.0e-6,4.0e-6,0.001}
31        make/o/t parameters_CSCpr = {"scale","mean CORE radius (A)","radial polydispersity (sigma)","CORE length (A)","radial shell thickness (A)","face shell thickness (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
32        Edit/W=(410,44,757,306)  parameters_CSCpr,coef_CSCpr
33        ModifyTable width(parameters_CSCpr)=162
34       
35        Variable/G root:g_CSCpr
36        g_CSCpr := PolyCoShCylinder(coef_CSCpr,ywave_CSCpr,xwave_CSCpr)
37        Display ywave_CSCpr vs xwave_CSCpr
38        ModifyGraph log=1,marker=29,msize=2,mode=4
39        Label bottom "q (A\\S-1\\M)"
40        Label left "Intensity (cm\\S-1\\M)"
41        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
42End
43
44// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
45Proc PlotSmearedPolyCoShCylinder(str)                                                           
46        String str
47        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
48       
49        // if any of the resolution waves are missing => abort
50        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
51                Abort
52        endif
53       
54        SetDataFolder $("root:"+str)
55       
56        // Setup parameter table for model function
57        make/o/d smear_coef_CSCpr = {0.01,150,0.10,10,20.,10.,4.0e-6,1.0e-6,4.0e-6,0.001}
58        make/o/t smear_parameters_CSCpr = {"scale","mean CORE radius (A)","radial polydispersity (sigma)","CORE length (A)","radial shell thickness (A)","face shell thickness (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
59        Edit smear_parameters_CSCpr,smear_coef_CSCpr
60       
61        // output smeared intensity wave, dimensions are identical to experimental QSIG values
62        // make extra copy of experimental q-values for easy plotting
63       
64        Duplicate/O $(str+"_q") smeared_CSCpr,smeared_qvals                             
65        SetScale d,0,0,"1/cm",smeared_CSCpr                                                     
66                                       
67        Variable/G gs_CSCpr=0
68        gs_CSCpr := fSmearedPolyCoShCylinder(smear_coef_CSCpr,smeared_CSCpr,smeared_qvals)      //this wrapper fills the STRUCT
69       
70        Display smeared_CSCpr vs smeared_qvals                                                                 
71        ModifyGraph log=1,marker=29,msize=2,mode=4
72        Label bottom "q (\\S-1\\M)"
73        Label left "Intensity (cm\\S-1\\M)"
74        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
75       
76        SetDataFolder root:
77End
78       
79
80
81//AAO version, uses XOP if available
82// simply calls the original single point calculation with
83// a wave assignment (this will behave nicely if given point ranges)
84Function PolyCoShCylinder(cw,yw,xw) : FitFunc
85        Wave cw,yw,xw
86       
87#if exists("PolyCoShCylinderX")
88        yw = PolyCoShCylinderX(cw,xw)
89#else
90        yw = fPolyCoShCylinder(cw,xw)
91#endif
92        return(0)
93End
94
95///////////////////////////////////////////////////////////////////////////////
96// unsmeared model calculation: function integrates for a polydisperse radius.
97//  Relies on the following two functions to return the monodisperse form factor.
98///////////////////////////////////////////////////////////////////////////////
99
100Function fPolyCoShCylinder(w,x) : FitFunc
101        Wave w
102        Variable x
103
104//The input variables are (and output)
105        //[0] scale
106        //[1] cylinder CORE RADIUS (A)
107        //[2] radial polydispersity (sigma)
108        //[3] cylinder CORE LENGTH (A)
109        //[4] radial shell Thickness (A)
110        //[5] face shell Thickness (A)
111        //[6] core SLD (A^-2)
112        //[7] shell SLD (A^-2)
113        //[8] solvent SLD (A^-2)
114        //[9] background (cm^-1)
115        Variable scale,length,sigma,bkg,radius,radthick,facthick,rhoc,rhos,rhosolv
116        Variable fc, vcyl,qq
117        Variable nord,ii,va,vb,summ,yyy,rad,AR,lgAR,zed,Rsqr,lgRsqr,Rsqrsumm,Rsqryyy,tot
118        String weightStr,zStr
119        scale = w[0]
120        radius = w[1]
121        sigma = w[2]                            //sigma is the standard mean deviation
122        length = w[3]
123        radthick = w[4]
124        facthick= w[5]
125        rhoc = w[6]
126        rhos = w[7]
127        rhosolv = w[8]
128        bkg = w[9]
129       
130        weightStr = "gauss20wt"
131        zStr = "gauss20z"
132       
133//      if wt,z waves don't exist, create them
134       
135        if (WaveExists($weightStr) == 0)                // wave reference is not valid,
136                Make/D/N=20 $weightStr,$zStr
137                Wave w20 = $weightStr
138                Wave z20 = $zStr                                // wave references to pass
139                Make20GaussPoints(w20,z20)     
140        else
141                if(exists(weightStr) > 1)
142                         Abort "wave name is already in use"            // execute if condition is false
143                endif
144                Wave w20 = $weightStr
145                Wave z20 = $zStr                       
146        endif
147
148/////////////////////////////////////////////////////////////////////////
149// This integration loop is for the radial polydispersity.
150//  The loop uses values from cylintegration to average
151// the scattering over a radial size distribution.
152/////////////////////////////////////////////////////////////////////////
153
154        nord = 20
155        va = exp(ln(radius)-(4.*sigma))
156        if (va<0)
157                va=0                                    //to avoid numerical error when  va<0 (-ve r value)
158        endif
159        vb = exp(ln(radius)+(4.*sigma))
160
161//      zed = ((radius*radius)/(sigma*sigma))-1         // If you want to use a Schultz distribution instead
162
163// evaluate at Gauss points
164// remember to index from 0,size-1
165        qq = x         
166        summ = 0.0                      // initialize integral
167        Rsqrsumm = 0.0
168       
169        ii=0
170        do
171                // Using 20 Gauss points
172                rad = ( z20[ii]*(vb-va) + vb + va )/2.0                 //make distribution points
173
174//              lgAR = (zed*ln(rad))-((rad*(zed+1))/radius)-((zed+1)*ln(radius/(zed+1)))-gammln(zed+1)
175                //create Schultz distribution
176//              AR = exp(lgAR)                                  //invert Schultz to prevent overflow/underflow
177               
178//              AR=(1/(rad*sigma*sqrt(2*Pi)))*exp(-(0.5*((ln(radius/rad))/sigma)*((ln(radius/rad))/sigma)))
179                AR=(1/(rad*sigma*sqrt(2*Pi)))*exp(-(0.5*((ln(rad/radius))/sigma)*((ln(rad/radius))/sigma)))
180
181                yyy = w20[ii] * AR * cylintegration(qq,rad,radthick,facthick,rhoc,rhos,rhosolv,length)
182//              Rsqryyy= w20[ii] * AR * rad*rad         //A.Nelson, original does not include shell
183                Rsqryyy= w20[ii] * AR * (rad+radthick)*(rad+radthick)           //SRK normalize to total dimensions
184               
185                summ += yyy
186                Rsqrsumm +=  Rsqryyy
187                ii+=1
188        while (ii<nord)         // end of loop over quadrature points
189
190
191// calculate value of integral to return
192        fc = (vb-va)/2.0*summ
193        Rsqr=(vb-va)/2.0*Rsqrsumm
194
195//NOTE that for absolute intensity scaling you need to multiply by the
196// number density of particles. This is the vol frac of core particles
197// divided by the core volume.
198
199//      lgRsqr=2*ln(radius/(zed+1))+gammln(zed+3)-gammln(zed+1)
200//      Rsqr=exp(lgRsqr)                                                                                                                                       
201
202//      vcyl=Pi*Rsqr*length             //but you have to multiply by <R2> not <R>2.
203        vcyl=Pi*Rsqr*(length+2*facthick)                //SRK normalize to total dimensions
204        fc /= vcyl
205
206//convert to [cm-1]
207        fc *= 1.0e8
208//Scale
209        fc *= scale                     //scale will be the volume fraction of core particles.
210// add in the  incoherent background
211        fc += bkg
212
213        Return (fc)
214End
215
216////////////////////////////////////////////////////////////////////////////
217//Cylintegration calculates the Form factor for the monodisperse core shell
218////////////////////////////////////////////////////////////////////////////
219Function cylintegration(qq,rad,radthick,facthick,rhoc,rhos,rhosolv,length)
220        Variable  qq,rad,radthick,facthick,rhoc,rhos,rhosolv,length
221        Variable  answer,halfheight
222        Variable nord,ii,va,vb,summ,yyy,zi
223        String weightStr,zStr
224       
225        weightStr = "gauss76wt"
226        zStr = "gauss76z"
227
228//      if wt,z waves don't exist, create them
229//     20 Gauss points is not enough for cylinder calculation
230       
231        if (WaveExists($weightStr) == 0)                                // wave reference is not valid,
232                Make/D/N=76 $weightStr,$zStr
233                Wave w76 = $weightStr
234                Wave z76 = $zStr                                                        // wave references to pass
235                Make76GaussPoints(w76,z76)     
236        else
237                if(exists(weightStr) > 1)
238                         Abort "wave name is already in use"   
239                endif
240                Wave w76 = $weightStr
241                Wave z76 = $zStr       
242        endif
243
244// set up the integration end points
245        nord = 76
246        va = 0
247        vb = Pi/2
248        halfheight = length/2.0
249
250// evaluate at Gauss points
251// remember to index from 0,size-1
252
253        summ = 0.0                              // initialize integral
254        ii=0
255        do
256                // Using 76 Gauss points
257                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
258                yyy = w76[ii] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi)
259                summ += yyy
260     ii+=1
261        while(ii<nord)          // end of loop over quadrature points
262
263// calculate value of integral to return
264  answer = (vb-va)/2.0*summ
265        Return (answer)
266
267End                             //End of function cylintegration
268
269////////////////////////////////////////////////////////////////////////
270// F(qq, rcore, thick, rhoc,rhos,rhosolv, length, zi)  This returns the
271// arguments used for the integration over theta.
272////////////////////////////////////////////////////////////////////////
273Function CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, length, dum)
274        Variable qq, rad, radthick, facthick, rhoc,rhos,rhosolv, length, dum
275       
276// qq is the q-value for the calculation (1/A)
277// radius is the core radius of the cylinder (A)
278//  radthick and facthick are the radial and face layer thicknesses
279// rho(n) are the respective SLD's
280// length is the *Half* CORE-LENGTH of the cylinder
281// dum is the dummy variable for the integration (theta)
282
283        Variable dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval         //Local variables
284
285        dr1 = rhoc-rhos
286        dr2 = rhos-rhosolv
287        vol1 = Pi*rad*rad*(2*length)
288        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2*length+2*facthick)
289       
290        besarg1 = qq*rad*sin(dum)
291        besarg2 = qq*(rad+radthick)*sin(dum)
292        sinarg1 = qq*length*cos(dum)
293        sinarg2 = qq*(length+facthick)*cos(dum)
294       
295        t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1
296        t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2
297       
298        retval = ((t1+t2)^2)*sin(dum)
299        return retval
300   
301End     //Function CScyl()
302
303//wrapper to calculate the smeared model as an AAO-Struct
304// fills the struct and calls the ususal function with the STRUCT parameter
305//
306// used only for the dependency, not for fitting
307//
308Function fSmearedPolyCoShCylinder(coefW,yW,xW)
309        Wave coefW,yW,xW
310       
311        String str = getWavesDataFolder(yW,0)
312        String DF="root:"+str+":"
313       
314        WAVE resW = $(DF+str+"_res")
315       
316        STRUCT ResSmearAAOStruct fs
317        WAVE fs.coefW = coefW   
318        WAVE fs.yW = yW
319        WAVE fs.xW = xW
320        WAVE fs.resW = resW
321       
322        Variable err
323        err = SmearedPolyCoShCylinder(fs)
324       
325        return (0)
326End
327
328// this is all there is to the smeared calculation!
329Function SmearedPolyCoShCylinder(s) :FitFunc
330        Struct ResSmearAAOStruct &s
331
332//      the name of your unsmeared model (AAO) is the first argument
333        s.yW = Smear_Model_20(PolyCoShCylinder,s.coefW,s.xW,s.resW)
334
335        return(0)
336End
337       
Note: See TracBrowser for help on using the repository browser.