1 | # Example of conversion of scattering cross section from SANS in absolute |
---|
2 | # units into SESANS using a Hankel transformation |
---|
3 | # everything is in units of metres except specified otherwise |
---|
4 | # Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013 |
---|
5 | |
---|
6 | from __future__ import division |
---|
7 | |
---|
8 | from pylab import * |
---|
9 | from scipy.special import jv as besselj |
---|
10 | |
---|
11 | # q-range parameters |
---|
12 | |
---|
13 | q = arange(0.0003, 1.0, 0.0003); # [nm^-1] range wide enough for Hankel transform |
---|
14 | dq=(q[1]-q[0])*1e9; # [m^-1] step size in q, needed for integration |
---|
15 | nq=len(q); |
---|
16 | Lambda=2e-10; # [m] wavelength |
---|
17 | # sample parameters |
---|
18 | phi=0.1; # volume fraction |
---|
19 | R=100; # [nm] radius particles |
---|
20 | DeltaRho=6e14; # [m^-2] |
---|
21 | V=4/3*pi*R**3 * 1e-27; # [m^3] |
---|
22 | th=0.002; # [m] thickness sample |
---|
23 | |
---|
24 | #2 PHASE SYSTEM |
---|
25 | st= 1.5*Lambda**2*DeltaRho**2*th*phi*(1-phi)*R*1e-9 # scattering power in sesans formalism |
---|
26 | |
---|
27 | # Form factor solid sphere |
---|
28 | qr=q*R; |
---|
29 | P=(3.*(sin(qr)-qr*cos(qr)) / qr**3)**2; |
---|
30 | # Structure factor dilute |
---|
31 | S=1.; |
---|
32 | #2 PHASE SYSTEM |
---|
33 | # scattered intensity [m^-1] in absolute units according to SANS |
---|
34 | I=phi*(1-phi)*V*(DeltaRho**2)*P*S; |
---|
35 | |
---|
36 | clf() |
---|
37 | subplot(211) # plot the SANS calculation |
---|
38 | plot(q,I,'k') |
---|
39 | loglog(q,I) |
---|
40 | xlim([0.01, 1]) |
---|
41 | ylim([1, 1e9]) |
---|
42 | xlabel(r'$Q [nm^{-1}]$') |
---|
43 | ylabel(r'$d\Sigma/d\Omega [m^{-1}]$') |
---|
44 | |
---|
45 | |
---|
46 | # Hankel transform to nice range for plot |
---|
47 | nz=61; |
---|
48 | zz=linspace(0,240,nz); # [nm], should be less than reciprocal from q |
---|
49 | G=zeros(nz); |
---|
50 | for i in range(len(zz)): |
---|
51 | integr=besselj(0,q*zz[i])*I*q; |
---|
52 | G[i]=sum(integr); |
---|
53 | G=G*dq*1e9*2*pi; # integr step, conver q into [m**-1] and 2 pi circle integr |
---|
54 | # plot(zz,G); |
---|
55 | stt= th*Lambda**2/4/pi/pi*G[0] # scattering power according to SANS formalism |
---|
56 | PP=exp(th*Lambda**2/4/pi/pi*(G-G[0])); |
---|
57 | |
---|
58 | subplot(212) |
---|
59 | plot(zz,PP,'k',label="Hankel transform") # Hankel transform 1D |
---|
60 | xlabel('spin-echo length [nm]') |
---|
61 | ylabel('polarisation normalised') |
---|
62 | hold(True) |
---|
63 | |
---|
64 | # Cosine transformation of 2D scattering patern |
---|
65 | if False: |
---|
66 | qy,qz = meshgrid(q,q) |
---|
67 | qr=R*sqrt(qy**2 + qz**2); # reuse variable names Hankel transform, but now 2D |
---|
68 | P=(3.*(sin(qr)-qr*cos(qr)) / qr**3)**2; |
---|
69 | # Structure factor dilute |
---|
70 | S=1.; |
---|
71 | # scattered intensity [m^-1] in absolute units according to SANS |
---|
72 | I=phi*V*(DeltaRho**2)*P*S; |
---|
73 | GG=zeros(nz); |
---|
74 | for i in range(len(zz)): |
---|
75 | integr=cos(qz*zz[i])*I; |
---|
76 | GG[i]=sum(sum(integr)); |
---|
77 | GG=4*GG* dq**2; # take integration step into account take 4 quadrants |
---|
78 | # plot(zz,GG); |
---|
79 | sstt= th*Lambda**2/4/pi/pi*GG[0] # scattering power according to SANS formalism |
---|
80 | PPP=exp(th*Lambda**2/4/pi/pi*(GG-GG[0])); |
---|
81 | |
---|
82 | plot(zz,PPP,label="cosine transform") # cosine transform 2D |
---|
83 | |
---|
84 | # For comparison calculation in SESANS formalism, which overlaps perfectly |
---|
85 | def gsphere(z,r): |
---|
86 | """ |
---|
87 | Calculate SESANS-correlation function for a solid sphere. |
---|
88 | |
---|
89 | Wim Bouwman after formulae Timofei Kruglov J.Appl.Cryst. 2003 article |
---|
90 | """ |
---|
91 | d = z/r |
---|
92 | g = zeros_like(z) |
---|
93 | g[d==0] = 1. |
---|
94 | low = ((d > 0) & (d < 2)) |
---|
95 | dlow = d[low] |
---|
96 | dlow2 = dlow**2 |
---|
97 | print dlow.shape, dlow2.shape |
---|
98 | g[low] = sqrt(1-dlow2/4.)*(1+dlow2/8.) + dlow2/2.*(1-dlow2/16.)*log(dlow/(2.+sqrt(4.-dlow2))) |
---|
99 | return g |
---|
100 | |
---|
101 | if True: |
---|
102 | plot(zz,exp(st*(gsphere(zz,R)-1)),'r', label="analytical") |
---|
103 | legend() |
---|
104 | show() |
---|