source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/PolyCoreShellCylinder_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: 11.1 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","parameters_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 (A\\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","smear_parameters_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.