source: sasmodels/sesansdemo.py @ a42b091

Last change on this file since a42b091 was c97724e, checked in by pkienzle, 10 years ago

add sesans support to bumps model

  • Property mode set to 100644
File size: 3.1 KB
Line 
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
6from __future__ import division
7
8from pylab import *
9from scipy.special import jv as besselj
10
11# q-range parameters
12
13q = arange(0.0003, 1.0, 0.0003);    # [nm^-1] range wide enough for  Hankel transform
14dq=(q[1]-q[0])*1e9;   # [m^-1] step size in q, needed for integration
15nq=len(q);
16Lambda=2e-10;   # [m] wavelength
17# sample parameters
18phi=0.1;   # volume fraction
19R=100;       # [nm] radius particles
20DeltaRho=6e14;  # [m^-2]
21V=4/3*pi*R**3 * 1e-27; # [m^3]
22th=0.002;    # [m] thickness sample
23
24#2 PHASE SYSTEM
25st= 1.5*Lambda**2*DeltaRho**2*th*phi*(1-phi)*R*1e-9  # scattering power in sesans formalism
26
27# Form factor solid sphere
28qr=q*R;
29P=(3.*(sin(qr)-qr*cos(qr)) / qr**3)**2;
30# Structure factor dilute
31S=1.;
32#2 PHASE SYSTEM
33# scattered intensity [m^-1] in absolute units according to SANS
34I=phi*(1-phi)*V*(DeltaRho**2)*P*S;
35
36clf()
37subplot(211)  # plot the SANS calculation
38plot(q,I,'k')
39loglog(q,I)
40xlim([0.01, 1])
41ylim([1, 1e9])
42xlabel(r'$Q [nm^{-1}]$')
43ylabel(r'$d\Sigma/d\Omega [m^{-1}]$')
44
45
46# Hankel transform to nice range for plot
47nz=61;
48zz=linspace(0,240,nz); # [nm], should be less than reciprocal from q
49G=zeros(nz);
50for i in range(len(zz)):
51    integr=besselj(0,q*zz[i])*I*q;
52    G[i]=sum(integr);
53G=G*dq*1e9*2*pi; # integr step, conver q into [m**-1] and 2 pi circle integr
54# plot(zz,G);
55stt= th*Lambda**2/4/pi/pi*G[0]  # scattering power according to SANS formalism
56PP=exp(th*Lambda**2/4/pi/pi*(G-G[0]));
57
58subplot(212)
59plot(zz,PP,'k',label="Hankel transform") # Hankel transform 1D
60xlabel('spin-echo length [nm]')
61ylabel('polarisation normalised')
62hold(True)
63
64# Cosine transformation of 2D scattering patern
65if 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
85def 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
101if True:
102    plot(zz,exp(st*(gsphere(zz,R)-1)),'r', label="analytical")
103legend()
104show()
Note: See TracBrowser for help on using the repository browser.