source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/PolyCoreBicelle_v40.ipf @ 682

Last change on this file since 682 was 633, checked in by srkline, 13 years ago

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.

File size: 10.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
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//  7th Jan 2007  DS  Core-shell to include varying SLD in rim and face: Model 2
18//  Jan 2009  ACH - changed to run in Igor 6.03
19//
20//  May 2009 AJJ - Renamed to PolyCoreBicelle_v40.ipf
21/////////////////////////////////////////////////////////////////
22
23Proc PlotPolyCoreBicelle(num,qmin,qmax)
24        Variable num=100,qmin=0.001,qmax=0.7
25        Prompt num "Enter number of data points for model: "
26        Prompt qmin "Enter minimum q-value (A^-1) for model: "
27        Prompt qmax "Enter maximum q-value (A^-1) for model: "
28       
29        make/o/d/n=(num) xwave_PCBicelle,ywave_PCBicelle
30        xwave_PCBicelle = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
31        make/o/d coef_PCBicelle = {0.01,150,0.10,10,20.,10.,4.0e-6, 1.0e-6,2.0e-6,4.0e-6,0.001}
32        make/o/t parameters_PCBicelle = {"scale","mean CORE radius (A)","radial polydispersity (sigma)","CORE length (A)","radial shell thickness (A)","face shell thickness (A)","SLD core (A^-2)","SLD face (A^-2)","SLD rim(A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
33        Edit/W=(410,44,757,306)  parameters_PCBicelle,coef_PCBicelle
34        ModifyTable width(parameters_PCBicelle)=162
35       
36        Variable/G root:g_PCBicelle
37        g_PCBicelle := PolyCoreBicelle(coef_PCBicelle,ywave_PCBicelle,xwave_PCBicelle)
38        Display ywave_PCBicelle vs xwave_PCBicelle
39        ModifyGraph log=1, marker=29,msize=2,mode=4
40        Label bottom "q (A\\S-1\\M)"
41        Label left "Intensity (cm\\S-1\\M)"
42        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
43       
44        AddModelToStrings("PolyCoreBicelle","coef_PCBicelle","parameters_PCBicelle","PCBicelle")
45End
46
47
48Proc PlotSmearedPolyCoreBicelle(str)
49        String str
50        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)                                                 
51        //no input parameters necessary, it MUST use the experimental q-values
52        // from the experimental data read in from an AVE/QSIG data file
53       
54        // if no gQvals wave, data must not have been loaded => abort
55        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
56                Abort
57        endif
58       
59        SetDataFolder $("root:"+str)
60       
61        // Setup parameter table for model function
62        make/o/d smear_coef_PCBicelle = {0.01,150,0.10,10,20.,10.,4.0e-6,1.0e-6,2.0e-6,4.0e-6,0.001}
63        make/o/t smear_parameters_PCBicelle = {"scale","mean CORE radius (A)","radial polydispersity (sigma)","CORE length (A)","radial shell thickness (A)","face shell thickness (A)","SLD core (A^-2)","SLD face (A^-2)","SLD rim(A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
64        Edit smear_parameters_PCBicelle,smear_coef_PCBicelle
65       
66        // output smeared intensity wave, dimensions are identical to experimental QSIG values
67        // make extra copy of experimental q-values for easy plotting
68       
69        Duplicate/O $(str+"_q") smeared_PCBicelle,smeared_qvals                         
70        SetScale d,0,0,"1/cm",smeared_PCBicelle                                                 
71                                       
72        Variable/G gs_PCBicelle=0
73        gs_PCBicelle := fSmearedPolyCoreBicelle(smear_coef_PCBicelle,smeared_PCBicelle,smeared_qvals)
74               
75        Display smeared_PCBicelle vs smeared_qvals                                                                     
76        ModifyGraph log=1,marker=29,msize=2,mode=4
77        Label bottom "q (A\\S-1\\M)"
78        Label left "Intensity (cm\\S-1\\M)"
79        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
80       
81        SetDataFolder root:
82        AddModelToStrings("SmearedPolyCoreBicelle","smear_coef_PCBicelle","smear_parameters_PCBicelle","PCBicelle")
83End
84///////////////////////////////////////////////////////////////////////////////
85// unsmeared model calculation: function integrates for a polydisperse radius.
86//  Relies on the following two functions to return the monodisperse form factor.
87///////////////////////////////////////////////////////////////////////////////
88Function PolyCoreBicelle(cw,yw,xw) : FitFunc
89        Wave cw,yw,xw
90       
91#if exists("PolyCoreBicelleX")
92        yw = PolyCoreBicelleX(cw,xw)
93#else
94        yw = fPolyCoreBicelle(cw,xw)
95#endif
96        return(0)
97End
98
99Function  fPolyCoreBicelle(w,x) : FitFunc
100        Wave w
101        Variable x
102
103//The input variables are (and output)
104        //[0] scale
105        //[1] cylinder CORE RADIUS (A)
106        //[2] radial polydispersity (sigma)
107        //[3] cylinder CORE LENGTH (A)
108        //[4] radial shell Thickness (A)
109        //[5] face shell Thickness (A)
110        //[6] core SLD (A^-2)
111        //[7] face SLD (A^-2)
112        //[8] rim SLD (A^-2)
113        //[9] solvent SLD (A^-2)
114        //[10] background (cm^-1)
115        Variable scale,length,sigma,bkg,radius,radthick,facthick,rhoc,rhoh, rhor, 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        rhoh = w[7]
127        rhor=w[8]
128        rhosolv = w[9]
129        bkg = w[10]
130       
131        weightStr = "gauss20wt"
132        zStr = "gauss20z"
133       
134//      if wt,z waves don't exist, create them
135       
136        if (WaveExists($weightStr) == 0)                // wave reference is not valid,
137                Make/D/N=20 $weightStr,$zStr
138                Wave w20 = $weightStr
139                Wave z20 = $zStr                                // wave references to pass
140                Make20GaussPoints(w20,z20)     
141        else
142                if(exists(weightStr) > 1)
143                         Abort "wave name is already in use"            // execute if condition is false
144                endif
145                Wave w20 = $weightStr
146                Wave z20 = $zStr                       
147        endif
148
149/////////////////////////////////////////////////////////////////////////
150// This integration loop is for the radial polydispersity.
151//  The loop uses values from cylintegration to average
152// the scattering over a radial size distribution.
153/////////////////////////////////////////////////////////////////////////
154
155        nord = 20
156        va = exp(ln(radius)-(4.*sigma))
157        if (va<0)
158                va=0                                    //to avoid numerical error when  va<0 (-ve r value)
159        endif
160        vb = exp(ln(radius)+(4.*sigma))
161
162// evaluate at Gauss points
163// remember to index from 0,size-1
164        qq = x         
165        summ = 0.0                      // initialize integral
166        Rsqrsumm = 0.0
167       
168        ii=0
169        do
170                // Using 20 Gauss points
171                rad = ( z20[ii]*(vb-va) + vb + va )/2.0                 //make distribution points
172               
173                AR=(1/(rad*sigma*sqrt(2*Pi)))*exp(-(0.5*((ln(radius/rad))/sigma)*((ln(radius/rad))/sigma)))
174
175                yyy = w20[ii] * AR * BicelleIntegration(qq,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length)
176                Rsqryyy= w20[ii] * AR * (rad+radthick)*(rad+radthick)           //SRK normalize to total dimensions
177               
178                summ += yyy
179                Rsqrsumm +=  Rsqryyy
180                ii+=1
181        while (ii<nord)         // end of loop over quadrature points
182
183// calculate value of integral to return
184        fc = (vb-va)/2.0*summ
185        Rsqr=(vb-va)/2.0*Rsqrsumm
186
187//NOTE that for absolute intensity scaling you need to multiply by the
188// number density of particles. This is the vol frac of core particles
189// divided by the core volume.
190        vcyl=Pi*Rsqr*(length+2*facthick)                //SRK normalize to total dimensions
191        fc /= vcyl
192
193//convert to [cm-1]
194        fc *= 1.0e8
195//Scale
196        fc *= scale                     //scale will be the volume fraction of core particles.
197// add in the  incoherent background
198        fc += bkg
199
200        Return (fc)
201End
202
203////////////////////////////////////////////////////////////////////////////
204//BicelleIntegration calculates the Form factor for the Bicelle core shell cylinder
205////////////////////////////////////////////////////////////////////////////
206Function  BicelleIntegration(qq,rad,radthick,facthick,rhoc, rhoh, rhor, rhosolv, length)
207        Variable  qq,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length
208        Variable  answer,halfheight
209        Variable nord,ii,va,vb,summ,yyy,zi
210        String weightStr,zStr
211       
212        weightStr = "gauss76wt"
213        zStr = "gauss76z"
214
215//      if wt,z waves don't exist, create them
216//     20 Gauss points is not enough for cylinder calculation
217       
218        if (WaveExists($weightStr) == 0)                                // wave reference is not valid,
219                Make/D/N=76 $weightStr,$zStr
220                Wave w76 = $weightStr
221                Wave z76 = $zStr                                                        // wave references to pass
222                Make76GaussPoints(w76,z76)     
223        else
224                if(exists(weightStr) > 1)
225                         Abort "wave name is already in use"   
226                endif
227                Wave w76 = $weightStr
228                Wave z76 = $zStr       
229        endif
230
231// set up the integration end points
232        nord = 76
233        va = 0
234        vb = Pi/2
235        halfheight = length/2.0
236
237// evaluate at Gauss points
238// remember to index from 0,size-1
239
240        summ = 0.0                              // initialize integral
241        ii=0
242        do
243                // Using 76 Gauss points
244                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
245                yyy = w76[ii] * BicelleKernel(qq, rad, radthick, facthick, rhoc,rhoh, rhor,rhosolv, halfheight, zi)
246                summ += yyy
247     ii+=1
248        while(ii<nord)          // end of loop over quadrature points
249
250// calculate value of integral to return
251  answer = (vb-va)/2.0*summ
252        Return (answer)
253
254End                             //End of function cylintegration
255
256////////////////////////////////////////////////////////////////////////
257// F(qq, rcore, thick, rhoc,rhos,rhosolv, length, zi)  This returns the
258// arguments used for the integration over theta for orientational average
259////////////////////////////////////////////////////////////////////////
260Function BicelleKernel(qq, rad, radthick, facthick, rhoc,rhoh, rhor,rhosolv, length, dum)
261        Variable qq, rad, radthick, facthick, rhoc,rhoh,rhor,rhosolv, length, dum
262       
263// qq is the q-value for the calculation (1/A)
264// radius is the core radius of the cylinder (A)
265//  radthick and facthick are the radial and face layer thicknesses
266// rho(n) are the respective SLD's
267// length is the *Half* CORE-LENGTH of the cylinder
268// dum is the dummy variable for the integration (theta)
269
270        Variable dr1,dr2, dr3, besarg1,besarg2, vol1,vol2, vol3,sinarg1,sinarg2,t1,t2,t3, retval                //Local variables
271        Variable si1,si2,be1,be2
272
273        dr1 = rhoc-rhoh
274        dr2 = rhor-rhosolv
275        dr3=  rhoh-rhor
276        vol1 = Pi*rad*rad*(2*length)
277        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2*length+2*facthick)
278        vol3= Pi*(rad)*(rad)*(2*length+2*facthick)
279        besarg1 = qq*rad*sin(dum)
280        besarg2 = qq*(rad+radthick)*sin(dum)
281        sinarg1 = qq*length*cos(dum)
282        sinarg2 = qq*(length+facthick)*cos(dum)
283       
284        if(besarg1 == 0)
285                be1 = 0.5
286        else
287                be1 = bessJ(1,besarg1)/besarg1
288        endif
289        if(besarg2 == 0)
290                be2 = 0.5
291        else
292                be2 = bessJ(1,besarg2)/besarg2
293        endif
294        if(sinarg1 == 0)
295                si1 = 1
296        else
297                si1 = sin(sinarg1)/sinarg1
298        endif
299        if(sinarg2 == 0)
300                si2 = 1
301        else
302                si2 = sin(sinarg2)/sinarg2
303        endif
304       
305        t1 = 2*vol1*dr1*si1*be1
306        t2 = 2*vol2*dr2*si2*be2
307        t3 = 2*vol3*dr3*si2*be1
308       
309        retval = ((t1+t2+t3)^2)*sin(dum)
310        return retval
311   
312End
313
314// this is all there is to the smeared calculation!
315Function fSmearedPolyCoreBicelle(coefW,yW,xW)
316        Wave coefW,yW,xW
317       
318        String str = getWavesDataFolder(yW,0)
319        String DF="root:"+str+":"
320       
321        WAVE resW = $(DF+str+"_res")
322       
323        STRUCT ResSmearAAOStruct fs
324        WAVE fs.coefW = coefW   
325        WAVE fs.yW = yW
326        WAVE fs.xW = xW
327        WAVE fs.resW = resW
328       
329        Variable err
330        err = SmearedPolyCoreBicelle(fs)
331         
332        return (0)
333End
334
335// this is all there is to the smeared calculation!
336Function SmearedPolyCoreBicelle(s) :FitFunc
337        Struct ResSmearAAOStruct &s
338
339//      the name of your unsmeared model (AAO) is the first argument
340        Smear_Model_20(PolyCoreBicelle,s.coefW,s.xW,s.yW,s.resW)
341
342        return(0)
343End
344       
Note: See TracBrowser for help on using the repository browser.