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