source: trunk/User form factors for Irena/CoreShellEllipsoidsForIrena.ipf

Last change on this file was 1035, checked in by ilavsky, 22 months ago

Minor changes in Coreshellellipsoid
  • Property svn:eol-style set to native
File size: 13.2 KB
Line 
1#pragma TextEncoding = "UTF-8"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3#include "IR1_Loader"
4
5
6
7
8////////////////////////////////////////////////
9//NOTE: this code is copy of code from OblateCoreShell.ipf from NIST Igor package, modified to use in Irena. Below is original notes from the NIST package:
10////////////////////////////////////////////////
11//
12// this function is for the form factor of an oblate ellipsoid with a core-shell structure
13//
14// 06 NOV 98 SRK
15//
16// 2018-06 modified for use in Irena, all glory goes to Steven Kline!!!!
17//note: modified for use in Irena by reducing flexibility and changing parameters description
18
19//*************************************************************************************************
20//*************************************************************************************************
21//USE in Irena :
22//In Modeling II select User form factor
23//In panel put in "Name of FormFactor function this string:    IR1T_EllipsoidalCoreShell
24//In Panel put in Name of volume FF function this string:     IR1T_EllipsoidalVolume
25//
26// Par1 is the aspect ratio which for ellipsoids are defiend as rotational objects with dimensions R x R x AR*R, note, AR=1 may fail.
27// par2 is shell thickness in A, and it is the same thickness everywhere on the ellipsoid.
28// par3, 4 and 5 are contrasts as this is core shell system and contrasts are part of the form factor.
29// par3, 4 and 5 are implicitely multiplied by 10^10cm^-2, so insert only a number. These are rhos not, delta-rho-square
30// In main panel set contrast = 1 !!!!!
31//*************************************************************************************************
32//*************************************************************************************************
33 
34////////////////////////////////////////////////
35//****   IR1T_OblateCoreShell ***** is simply the form factor, which normalizes to 1 at Q=0
36//Threadsafe
37
38Function IR1T_EllipsoidalCoreShell(Qval,radius, par1,par2,par3,par4,par5)
39        variable Qval, radius, par1,par2,par3,par4,par5                                                                                         
40        //par1  aspect ratio
41        //par2  shell thickness (same everywhere)
42        //par3, par4, par5              SLD of core, shell, and sovent
43       
44        Make/D/Free coef_oef = {1.,200,20,250,30,1e-6,2e-6,6.3e-6,0}
45        //make/o/t parameters_oef = {"scale","major core (A)","minor core (A)","major shell (A)","minor shell (A)","SLD core (A-2)","SLD shell (A-2)","SLD solvent (A^-2)","bkg (cm-1)"}
46        //set scale to 1, background to 0,
47        //
48        if(par1 <= 1 )                  //oblate shape
49                coef_oef[0] =1                                                                  //scale, set to 1, Irena uses its own scale
50                coef_oef[1] = radius                                                    //this is the main dimension of the shape
51                coef_oef[2] = radius*par1                                       //minor, smaller, diemnsion = secondary (AR*R) dimension of the particle
52                coef_oef[3] = radius+par2                                       //major shell = radius+ShellThickness
53                coef_oef[4] = radius*par1+par2                          //minor shell (A) = radius*AR+ShellThickness
54                coef_oef[5] = par3//    *1e-6                                                   //SLD core, their input units are A^-2 and values around 10-6 A^-2, for Irena input units are N*10^10 cm^-2,
55                coef_oef[6] = par4      //*1e-6                                         //SLD shell, this       *1e-6 converts it into their original units...
56                coef_oef[7] = par5      //*1e-6                                         //SLD solvent
57                coef_oef[8] = 0                                                                 //background, irena uses its own background
58        else            //prolate
59                coef_oef[0] =1                                                                  //scale, set to 1, Irena uses its own scale
60                coef_oef[1] = radius*par1                                       //this is the main dimension of the shape
61                coef_oef[2] = radius                                                    //minor, smaller, diemnsion = secondary (AR*R) dimension of the particle
62                coef_oef[3] = radius*par1+par2                                  //major shell = radius+ShellThickness
63                coef_oef[4] = radius+par2                                               //minor shell (A) = radius*AR+ShellThickness
64                coef_oef[5] = par3//    *1e-6                                                   //SLD core, their input units are A^-2 and values around 10-6 A^-2, for Irena input units are N*10^10 cm^-2,
65                coef_oef[6] = par4      //*1e-6                                         //SLD shell, this       *1e-6 converts it into their original units...
66                coef_oef[7] = par5      //*1e-6                                         //SLD solvent
67                coef_oef[8] = 0                                                                 //background, irena uses its own background
68        endif
69        //
70       
71        if(par1 <= 1 )                  //oblate shape
72#if exists("OblateFormX")
73   //  this returns F^2, we need F, so square root that
74   //           after scaling by volume which needs to be converted to cm.
75   //   IR1T_EllispodalVolume(radius, par1, par2,par3,par4,par5)/1e8
76   //   oblatevol = IR1T_OblateVolume(trmaj, AspectRatio)
77   //   answer /= oblatevol                     -- this is needs to be taken out, Irena does its own volumehandling here... 
78   //   also correct for their conversion to [A-1] to [cm-1]
79   //   answer *= 1.0e8 
80   //   not needed, set to 1 scale
81   //  answer *= scale
82   // not needed, set to 0 then add background
83   //   answer += bkg
84        return sqrt(OblateFormX(coef_oef,Qval)/(IR1T_EllipsoidalVolume(radius, par1, par2,par3,par4,par5)/1e8))/1e8
85#else
86        return sqrt(IR1T_fOblateForm(coef_oef,Qval)/(IR1T_EllipsoidalVolume(radius, par1, par2,par3,par4,par5)/1e8))/1e8
87#endif
88        else            //prolate shape
89#if exists("ProlateFormX")
90        return sqrt(ProlateFormX(coef_oef,Qval)/(IR1T_EllipsoidalVolume(radius, par1, par2,par3,par4,par5)/1e8))/1e8
91#else
92        return sqrt(IR1T_fProlateForm(coef_oef,Qval)/(IR1T_EllipsoidalVolume(radius, par1, par2,par3,par4,par5)/1e8))/1e8
93#endif 
94        endif
95end
96
97Threadsafe Function IR1T_EllipsoidalVolume(radius, AspectRatio, par2,par3,par4,par5)
98        variable radius, AspectRatio, par2,par3,par4,par5
99        return 4*Pi/3*radius*radius*radius*AspectRatio
100end
101
102
103///////////////////////////////////////////////////////////////
104// Oblate functions unsmeared model calculation
105///////////////////////////
106static Function IR1T_fOblateForm(w,x) : FitFunc
107        Wave w
108        Variable x
109
110//The input variables are (and output)
111        //[0] scale
112        //[1] crmaj, major radius of core       [A]
113        //[2] crmin, minor radius of core
114        //[3] trmaj, overall major radius
115        //[4] trmin, overall minor radius
116        //[5] sldc, core [A-2]
117        //[6] slds, shell
118        //[7] sld (solvent)
119        //[8] bkg, [cm-1]
120        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg,sldc,slds,sld
121        scale = w[0]
122        crmaj = w[1]
123        crmin = w[2]
124        trmaj = w[3]
125        trmin = w[4]
126        sldc = w[5]
127        slds = w[6]
128        sld = w[7]
129        bkg = w[8]
130
131        delpc = sldc - slds                     //core - shell
132        delps = slds - sld                      //shell - solvent
133       
134// local variables
135        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,oblatevol
136        String weightStr,zStr
137       
138        weightStr = "gauss76wt"
139        zStr = "gauss76z"
140
141       
142//      if wt,z waves don't exist, create them
143
144        if (WaveExists($weightStr) == 0) // wave reference is not valid,
145                Make/D/N=76 $weightStr,$zStr
146                Wave w76 = $weightStr
147                Wave z76 = $zStr                // wave references to pass
148                IR1T_Make76GaussPoints(w76,z76)
149        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
150        else
151                if(exists(weightStr) > 1)
152                         Abort "wave name is already in use"    // execute if condition is false
153                endif
154                Wave w76 = $weightStr
155                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
156        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
157        endif
158
159// set up the integration
160        // end points and weights
161        nord = 76
162        nfn = 2         //only <f^2> is calculated
163        npro = 0        // OBLATE ELLIPSOIDS
164        va =0
165        vb =1
166
167        qq = x          //current x point is the q-value for evaluation
168      summ = 0.0
169      ii=0
170      do
171                //printf "top of nord loop, i = %g\r",i
172        if(nfn ==1) //then              // "f1" required for beta factor
173          if(npro ==1) //then   // prolate
174                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
175          Endif
176          if(npro ==0) //then   // oblate 
177                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
178          Endif
179        Endif           //nfn = 1
180         if(nfn !=1) //then             //calculate"f2" = <f^2> = averaged form factor
181          if(npro ==1) //then   //prolate
182             zi = ( z76[ii]*(vb-va) + vb + va )/2.0
183          Endif
184          if(npro ==0) //then   //oblate
185                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
186                yyy = w76[ii]*IR1T_gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
187          Endif
188        Endif           //nfn <>1
189       
190        summ = yyy + summ               // get running total of integral
191        ii+=1
192        while (ii<nord)                         // end of loop over quadrature points
193        // calculate value of integral to return
194
195      answer = (vb-va)/2.0*summ
196     
197      // normalize by particle volume
198      oblatevol = 4*Pi/3*trmaj*trmaj*trmin
199      answer /= oblatevol
200     
201      //convert answer [A-1] to [cm-1]
202      answer *= 1.0e8 
203      //scale
204      answer *= scale
205      // //then add background
206      answer += bkg
207
208        Return (answer)
209End
210//
211//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
212//                       BY (53) & (58-59) IN CHEN AND
213//                       KOTLARCHYK REFERENCE
214//
215//       <OBLATE ELLIPSOID>
216
217static Function IR1T_gfn4(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
218        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
219        // local variables
220        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43
221        Variable siq,sit
222       
223        PI43=4.0/3.0*PI
224        aa = crmaj
225        bb = crmin
226        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx))
227        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx))
228        uq = sqrt(u2)*qq
229        ut= sqrt(ut2)*qq
230        vc = PI43*aa*aa*bb
231        vt = PI43*trmaj*trmaj*trmin
232        if(uq == 0)
233                siq = 1/3
234        else
235                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq
236        endif
237        if(ut == 0)
238                sit = 1/3
239        else
240                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut
241        endif
242       
243        gfnc = 3.0*siq*vc*delpc
244        gfnt = 3.0*sit*vt*delps
245        tgfn = gfnc+gfnt
246        gfn4 = tgfn*tgfn
247       
248        return gfn4
249       
250End             // function gfn4 for oblate ellipsoids
251
252///////////////////////////////////////////////////////////////
253// unsmeared model calculation
254///////////////////////////
255static Function IR1T_fProlateForm(w,x) : FitFunc
256        Wave w
257        Variable x
258
259//The input variables are (and output)
260        //[0] scale
261        //[1] crmaj, major radius of core       [A]
262        //[2] crmin, minor radius of core
263        //[3] trmaj, overall major radius
264        //[4] trmin, overall minor radius
265        //[5] sld core, [A^-2]
266        //[6] sld shell,
267        //[7] sld solvent
268        //[8] bkg [cm-1]
269        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg,sldc,slds,sld
270        scale = w[0]
271        crmaj = w[1]
272        crmin = w[2]
273        trmaj = w[3]
274        trmin = w[4]
275        sldc = w[5]
276        slds = w[6]
277        sld = w[7]
278        bkg = w[8]
279
280        delpc = sldc - slds                     //core - shell
281        delps = slds - sld                      //shell - solvent
282// local variables
283        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,prolatevol
284        String weightStr,zStr
285       
286        weightStr = "gauss76wt"
287        zStr = "gauss76z"
288
289//      if wt,z waves don't exist, create them
290
291        if (WaveExists($weightStr) == 0) // wave reference is not valid,
292                Make/D/N=76 $weightStr,$zStr
293                Wave w76 = $weightStr
294                Wave z76 = $zStr                // wave references to pass
295                IR1T_Make76GaussPoints(w76,z76)
296        else
297                if(exists(weightStr) > 1)
298                         Abort "wave name is already in use"    // execute if condition is false
299                endif
300                Wave w76 = $weightStr
301                Wave z76 = $zStr
302        endif
303
304// set up the integration
305        // end points and weights
306        nord = 76
307        nfn = 2         //only <f^2> is calculated
308        npro = 1        // PROLATE ELLIPSOIDS
309        va =0
310        vb =1
311//move this zi(i) evaluation inside other nord loop, since I don't have an array
312//      i=0
313//      do
314//       zi[i] = ( z76[i]*(vb-va) + vb + va )/2.0
315 //       i +=1
316 //     while (i<nord)
317//
318// evaluate at Gauss points
319        // remember to index from 0,size-1
320
321        qq = x          //current x point is the q-value for evaluation
322        summ = 0.0
323        ii=0
324        do
325                //printf "top of nord loop, i = %g\r",i
326                if(nfn ==1) //then              // "f1" required for beta factor
327                        if(npro ==1) //then     // prolate
328                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0 
329//           yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
330                        Endif
331//
332                        if(npro ==0) //then     // oblate 
333                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
334//            yyy = w76[i]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
335                        Endif
336                Endif           //nfn = 1
337          //
338                if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
339                        if(npro ==1) //then     //prolate
340                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
341                                yyy = w76[ii]*IR1T_gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
342                                //printf "yyy = %g\r",yyy
343                        Endif
344//
345                        if(npro ==0) //then     //oblate
346                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
347//              yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
348                        Endif
349                Endif           //nfn <>1
350         
351                summ = yyy + summ               // get running total of integral
352                ii+=1
353        while (ii<nord)                         // end of loop over quadrature points
354//   
355// calculate value of integral to return
356
357        answer = (vb-va)/2.0*summ
358       
359        //normailze by particle volume
360        prolatevol = 4*Pi/3*trmaj*trmin*trmin
361        answer /= prolatevol
362       
363        // rescale from 1/A to 1/cm
364        answer *= 1.0e8
365        //scale (arb)
366        answer *= scale
367        ////then add in background
368        answer += bkg
369
370        Return (answer)
371End     //prolate form factor
372
373//
374//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
375//                       BY (53) AND (56,57) IN CHEN AND
376//                       KOTLARCHYK REFERENCE
377//
378//     <PROLATE ELLIPSOIDS>
379//
380static Function IR1T_gfn2(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
381        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
382        // local variables
383        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn
384        Variable siq,sit
385
386        PI43=4.0/3.0*PI
387        aa = crmaj
388        bb = crmin
389        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx))
390        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx))
391        uq = sqrt(u2)*qq
392        ut= sqrt(ut2)*qq
393        vc = PI43*aa*bb*bb
394        vt = PI43*trmaj*trmin*trmin
395       
396        if(uq == 0.0)
397                siq = 1.0/3.0
398   else
399                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq
400        endif
401        if(ut == 0.0)
402                sit = 1.0/3.0
403        else
404                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut
405        endif
406       
407        gfnc = 3.0*siq*vc*delpc
408        gfnt = 3.0*sit*vt*delps
409        gfn = gfnc+gfnt
410        gfn2 = gfn*gfn
411       
412        return gfn2
413End             //function gfn2 for prolate ellipsoids
414
415
Note: See TracBrowser for help on using the repository browser.