HoneyComb_3_2x2_compare.py

a.glavic, 13 Apr 2015 23:47

Download (5.17 KB)

 
1
#-*- coding: utf8 -*-
2
"""
3
Spheres on two hexagonal close packed layers
4
"""
5
import numpy
6
import matplotlib
7
import pylab
8
from numpy import sin, cos, arange, pi, linspace, array, ones, abs, meshgrid, savetxt
9
from bornagain import *
10

    
11
alpha_i=0.35
12
lambda_i=5.5
13

    
14
phi_f_min, phi_f_max=-1.9, 1.75
15
tth_min, tth_max=-0.08, 4.4
16

    
17
lattice_rotation=0.
18

    
19
def get_sample(double_cell):
20
    """
21
    Build and return the sample representing spheres on two hexagonal close packed layers
22
    """
23
    film_thickness=13.*nanometer
24
    lattice_length=20.0*nanometer
25
    radius=lattice_length*0.6
26
    cylinder_ff=FormFactorCylinder(radius, film_thickness)
27
    surface_fraction=pi*radius**2/(3.*lattice_length**2*sin(60.*degree))
28

    
29
    m_air=HomogeneousMaterial("Air", 0.0, 0.0)
30
    m_substrate=HomogeneousMaterial("Silicon", 2.07261e-06, 0.)
31
    m_particle=HomogeneousMaterial("PermalloyHole", (-1.-surface_fraction)*9.09538e-06, 0.)
32
    m_layer=HomogeneousMaterial("PermalloyLayer", (1.-surface_fraction)*9.09538e-06, 0.)
33

    
34
    cylinder=Particle(m_particle, cylinder_ff)
35

    
36
    origin=kvector_t(0.0, 0.0, 0.0)
37
    if double_cell:
38
      basis_vec1=kvector_t(cos(lattice_rotation*degree)*lattice_length*1.73,
39
                           sin(lattice_rotation*degree)*lattice_length*1.73, 0.)
40
      basis_vec2=kvector_t(cos((lattice_rotation+120.)*degree)*lattice_length*1.73,
41
                           sin((lattice_rotation+120.)*degree)*lattice_length*1.73, 0.)
42
      basis=LatticeBasis()
43
      basis.addParticle(cylinder, [origin, basis_vec1, basis_vec2, basis_vec1+basis_vec2])
44
      particle_layout=ParticleLayout()
45
      particle_layout.addParticle(basis)
46

    
47
      interference=InterferenceFunction2DLattice.createHexagonal(2.*lattice_length*1.73,
48
                                                                 lattice_rotation*degree)
49
    else:
50
      basis=LatticeBasis()
51
      basis.addParticle(cylinder, [origin])
52
      particle_layout=ParticleLayout()
53
      particle_layout.addParticle(basis)
54

    
55
      interference=InterferenceFunction2DLattice.createHexagonal(lattice_length*1.73,
56
                                                                 lattice_rotation*degree)
57
    pdf=FTDistribution2DCauchy(150*nanometer, 150*nanometer)
58
    interference.setProbabilityDistribution(pdf)
59

    
60
    particle_layout.addInterferenceFunction(interference)
61

    
62
    multi_layer=MultiLayer()
63
    air_layer=Layer(m_air)
64
    multi_layer.addLayer(air_layer)
65

    
66
    part_layer=Layer(m_layer, film_thickness)
67
    part_layer.addLayout(particle_layout)
68
    multi_layer.addLayer(part_layer)
69

    
70
    substrate_layer=Layer(m_substrate, 0)
71
    multi_layer.addLayer(substrate_layer)
72
    return multi_layer
73

    
74

    
75
def get_simulation():
76
    """
77
    Create and return GISAXS simulation with beam and detector defined
78
    """
79
    simulation=Simulation()
80
    simulation.setDetectorParameters(248, phi_f_min*degree, phi_f_max*degree, 296,
81
                                     (tth_min-alpha_i)*degree, (tth_max-alpha_i)*degree)
82
    simulation.setBeamParameters(lambda_i*angstrom, alpha_i*degree, 0.0*degree)
83

    
84
    return simulation
85

    
86
def run_simulation_powder():
87
    """
88
    Run simulation and plot results
89
    """
90
    global lattice_rotation
91

    
92
    sample=get_sample(False)
93
    simulation=get_simulation()
94
    simulation.setSample(sample)
95
    simulation.runSimulation()
96
    result=simulation.getIntensityData().getArray()
97
#    for lattice_rotation in arange(5., 180., 5.):
98
#      print lattice_rotation
99
#      sample=get_sample(False)
100
#      simulation=get_simulation()
101
#      simulation.setSample(sample)
102
#      simulation.runSimulation()
103
#      result+=simulation.getIntensityData().getArray()
104

    
105
    # showing the result
106
    im=pylab.imshow(numpy.rot90(result+1., 1), norm=matplotlib.colors.LogNorm(1e2, 1e8),
107
                 extent=[2*pi/lambda_i*phi_f_min*degree, 2*pi/lambda_i*phi_f_max*degree,
108
                         2*pi/lambda_i*tth_min*degree, 2*pi/lambda_i*tth_max*degree],
109
                    aspect='auto', cmap='gist_ncar')
110
    cb=pylab.colorbar(im)
111
    pylab.title('single particle cell')
112
    cb.set_label(r'Intensity (arb. u.)', fontsize=16)
113
    pylab.xlabel(r'$Q_y (\AA^{-1})$', fontsize=16)
114
    pylab.ylabel(r'$Q_z (\AA^{-1})$', fontsize=16)
115

    
116
    sample=get_sample(True)
117
    simulation=get_simulation()
118
    simulation.setSample(sample)
119
    simulation.runSimulation()
120
    result=simulation.getIntensityData().getArray()
121
#    for lattice_rotation in arange(5., 180., 5.):
122
#      print lattice_rotation
123
#      sample=get_sample(True)
124
#      simulation=get_simulation()
125
#      simulation.setSample(sample)
126
#      simulation.runSimulation()
127
#      result+=simulation.getIntensityData().getArray()
128

    
129
    # showing the result
130
    pylab.figure()
131
    im=pylab.imshow(numpy.rot90(result+1., 1), norm=matplotlib.colors.LogNorm(1e2, 1e8),
132
                 extent=[2*pi/lambda_i*phi_f_min*degree, 2*pi/lambda_i*phi_f_max*degree,
133
                         2*pi/lambda_i*tth_min*degree, 2*pi/lambda_i*tth_max*degree],
134
                    aspect='auto', cmap='gist_ncar')
135
    cb=pylab.colorbar(im)
136
    pylab.title('2x2 particle cell')
137
    cb.set_label(r'Intensity (arb. u.)', fontsize=16)
138
    pylab.xlabel(r'$Q_y (\AA^{-1})$', fontsize=16)
139
    pylab.ylabel(r'$Q_z (\AA^{-1})$', fontsize=16)
140

    
141
    pylab.show()
142

    
143

    
144
if __name__=='__main__':
145
    run_simulation_powder()
146

    
147