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

Last change on this file since 663 was 663, checked in by ilavsky, 4 years ago

Add User Form factors for Irena

  • Property svn:eol-style set to native
File size: 13.5 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 ForFactor 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// par 2 is shell thickness in A, and it is the same thickness everywhere on teh 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 implicitelyu multipled 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 creect 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//            yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
176          Endif
177//
178          if(npro ==0) //then   // oblate 
179                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
180//            yyy = w76[ii]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
181          Endif
182        Endif           //nfn = 1
183        //
184        if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
185          if(npro ==1) //then   //prolate
186             zi = ( z76[ii]*(vb-va) + vb + va )/2.0
187//            yyy = w76[ii]*gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
188          //printf "yyy = %g\r",yyy
189          Endif
190//
191          if(npro ==0) //then   //oblate
192                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
193                yyy = w76[ii]*IR1T_gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
194          Endif
195        Endif           //nfn <>1
196       
197        summ = yyy + summ               // get running total of integral
198        ii+=1
199        while (ii<nord)                         // end of loop over quadrature points
200//   
201// calculate value of integral to return
202
203      answer = (vb-va)/2.0*summ
204     
205      // normalize by particle volume
206      oblatevol = 4*Pi/3*trmaj*trmaj*trmin
207      answer /= oblatevol
208     
209      //convert answer [A-1] to [cm-1]
210      answer *= 1.0e8 
211      //scale
212      answer *= scale
213      // //then add background
214      answer += bkg
215
216        Return (answer)
217End
218//
219//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
220//                       BY (53) & (58-59) IN CHEN AND
221//                       KOTLARCHYK REFERENCE
222//
223//       <OBLATE ELLIPSOID>
224
225static Function IR1T_gfn4(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
226        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
227        // local variables
228        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43
229        Variable siq,sit
230       
231        PI43=4.0/3.0*PI
232        aa = crmaj
233        bb = crmin
234        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx))
235        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx))
236        uq = sqrt(u2)*qq
237        ut= sqrt(ut2)*qq
238        vc = PI43*aa*aa*bb
239        vt = PI43*trmaj*trmaj*trmin
240        if(uq == 0)
241                siq = 1/3
242        else
243                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq
244        endif
245        if(ut == 0)
246                sit = 1/3
247        else
248                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut
249        endif
250       
251        gfnc = 3.0*siq*vc*delpc
252        gfnt = 3.0*sit*vt*delps
253        tgfn = gfnc+gfnt
254        gfn4 = tgfn*tgfn
255       
256        return gfn4
257       
258End             // function gfn4 for oblate ellipsoids
259
260///////////////////////////////////////////////////////////////
261// unsmeared model calculation
262///////////////////////////
263static Function IR1T_fProlateForm(w,x) : FitFunc
264        Wave w
265        Variable x
266
267//The input variables are (and output)
268        //[0] scale
269        //[1] crmaj, major radius of core       [A]
270        //[2] crmin, minor radius of core
271        //[3] trmaj, overall major radius
272        //[4] trmin, overall minor radius
273        //[5] sld core, [A^-2]
274        //[6] sld shell,
275        //[7] sld solvent
276        //[8] bkg [cm-1]
277        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg,sldc,slds,sld
278        scale = w[0]
279        crmaj = w[1]
280        crmin = w[2]
281        trmaj = w[3]
282        trmin = w[4]
283        sldc = w[5]
284        slds = w[6]
285        sld = w[7]
286        bkg = w[8]
287
288        delpc = sldc - slds                     //core - shell
289        delps = slds - sld                      //shell - solvent
290// local variables
291        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,prolatevol
292        String weightStr,zStr
293       
294        weightStr = "gauss76wt"
295        zStr = "gauss76z"
296
297//      if wt,z waves don't exist, create them
298
299        if (WaveExists($weightStr) == 0) // wave reference is not valid,
300                Make/D/N=76 $weightStr,$zStr
301                Wave w76 = $weightStr
302                Wave z76 = $zStr                // wave references to pass
303                IR1T_Make76GaussPoints(w76,z76)
304        else
305                if(exists(weightStr) > 1)
306                         Abort "wave name is already in use"    // execute if condition is false
307                endif
308                Wave w76 = $weightStr
309                Wave z76 = $zStr
310        endif
311
312// set up the integration
313        // end points and weights
314        nord = 76
315        nfn = 2         //only <f^2> is calculated
316        npro = 1        // PROLATE ELLIPSOIDS
317        va =0
318        vb =1
319//move this zi(i) evaluation inside other nord loop, since I don't have an array
320//      i=0
321//      do
322//       zi[i] = ( z76[i]*(vb-va) + vb + va )/2.0
323 //       i +=1
324 //     while (i<nord)
325//
326// evaluate at Gauss points
327        // remember to index from 0,size-1
328
329        qq = x          //current x point is the q-value for evaluation
330        summ = 0.0
331        ii=0
332        do
333                //printf "top of nord loop, i = %g\r",i
334                if(nfn ==1) //then              // "f1" required for beta factor
335                        if(npro ==1) //then     // prolate
336                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0 
337//           yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
338                        Endif
339//
340                        if(npro ==0) //then     // oblate 
341                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
342//            yyy = w76[i]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
343                        Endif
344                Endif           //nfn = 1
345          //
346                if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
347                        if(npro ==1) //then     //prolate
348                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
349                                yyy = w76[ii]*IR1T_gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
350                                //printf "yyy = %g\r",yyy
351                        Endif
352//
353                        if(npro ==0) //then     //oblate
354                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
355//              yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
356                        Endif
357                Endif           //nfn <>1
358         
359                summ = yyy + summ               // get running total of integral
360                ii+=1
361        while (ii<nord)                         // end of loop over quadrature points
362//   
363// calculate value of integral to return
364
365        answer = (vb-va)/2.0*summ
366       
367        //normailze by particle volume
368        prolatevol = 4*Pi/3*trmaj*trmin*trmin
369        answer /= prolatevol
370       
371        // rescale from 1/A to 1/cm
372        answer *= 1.0e8
373        //scale (arb)
374        answer *= scale
375        ////then add in background
376        answer += bkg
377
378        Return (answer)
379End     //prolate form factor
380
381//
382//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
383//                       BY (53) AND (56,57) IN CHEN AND
384//                       KOTLARCHYK REFERENCE
385//
386//     <PROLATE ELLIPSOIDS>
387//
388static Function IR1T_gfn2(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
389        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
390        // local variables
391        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn
392        Variable siq,sit
393
394        PI43=4.0/3.0*PI
395        aa = crmaj
396        bb = crmin
397        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx))
398        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx))
399        uq = sqrt(u2)*qq
400        ut= sqrt(ut2)*qq
401        vc = PI43*aa*bb*bb
402        vt = PI43*trmaj*trmin*trmin
403       
404        if(uq == 0.0)
405                siq = 1.0/3.0
406   else
407                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq
408        endif
409        if(ut == 0.0)
410                sit = 1.0/3.0
411        else
412                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut
413        endif
414       
415        gfnc = 3.0*siq*vc*delpc
416        gfnt = 3.0*sit*vt*delps
417        gfn = gfnc+gfnt
418        gfn2 = gfn*gfn
419       
420        return gfn2
421End             //function gfn2 for prolate ellipsoids
422
423
Note: See TracBrowser for help on using the repository browser.