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

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

Modified all procedure files to add to the keyword=value strings to identify the function, coefficients, and suffix once the model is plotted (and the objects will exist)

a one-liner in the Plot and PlotSmeared? macros.

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