source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/PolyCoreBicelle_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: 10.7 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//  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
272        dr1 = rhoc-rhoh
273        dr2 = rhor-rhosolv
274        dr3=  rhoh-rhor
275        vol1 = Pi*rad*rad*(2*length)
276        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2*length+2*facthick)
277        vol3= Pi*(rad)*(rad)*(2*length+2*facthick)
278        besarg1 = qq*rad*sin(dum)
279        besarg2 = qq*(rad+radthick)*sin(dum)
280        sinarg1 = qq*length*cos(dum)
281        sinarg2 = qq*(length+facthick)*cos(dum)
282       
283        t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1
284        t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2
285        t3 = 2*vol3*dr3*sin(sinarg2)/sinarg2*bessJ(1,besarg1)/besarg1
286       
287        retval = ((t1+t2+t3)^2)*sin(dum)
288        return retval
289   
290End
291
292// this is all there is to the smeared calculation!
293Function fSmearedPolyCoreBicelle(coefW,yW,xW)
294        Wave coefW,yW,xW
295       
296        String str = getWavesDataFolder(yW,0)
297        String DF="root:"+str+":"
298       
299        WAVE resW = $(DF+str+"_res")
300       
301        STRUCT ResSmearAAOStruct fs
302        WAVE fs.coefW = coefW   
303        WAVE fs.yW = yW
304        WAVE fs.xW = xW
305        WAVE fs.resW = resW
306       
307        Variable err
308        err = SmearedPolyCoreBicelle(fs)
309         
310        return (0)
311End
312
313// this is all there is to the smeared calculation!
314Function SmearedPolyCoreBicelle(s) :FitFunc
315        Struct ResSmearAAOStruct &s
316
317//      the name of your unsmeared model (AAO) is the first argument
318        Smear_Model_20(PolyCoreBicelle,s.coefW,s.xW,s.yW,s.resW)
319
320        return(0)
321End
322       
Note: See TracBrowser for help on using the repository browser.