This tutorial describes how to add your own structor factor data to Monte Carlo calculations
Copyright 2016, Vinothan N. Manoharan, Victoria Hwang, Annie Stephenson
This file is part of the structural-color python package.
This package is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This package is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this package. If not, see http://www.gnu.org/licenses/.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import structcol as sc
import structcol.refractive_index as ri
from structcol import montecarlo as mc
from structcol import detector as det
from structcol import model
from structcol import structure
%matplotlib inline
set parameters
In [7]:
wavelengths = sc.Quantity(np.arange(400, 800, 20), 'nm') # wavelengths
radius = sc.Quantity('0.5 um') # particle radius
volume_fraction = sc.Quantity(0.5, '') # volume fraction of particles
n_particle = ri.n('fused silica', wavelengths)
n_matrix = ri.n('vacuum', wavelengths) # called from the refractive_index module. n_matrix is the
n_medium = ri.n('vacuum', wavelengths) # space within sample. n_medium is outside the sample.
# n_particle and n_matrix can have complex indices if absorption is desired
thickness = sc.Quantity('50 um') # thickness of the sample film
Construct the structure factor data
Here, we use discrete points from the percus-yevick approximation for structure factor, as an example. In practice, you will most likely use actual structure factor data imported from your own file
In [12]:
qd_data = np.arange(0,75, 0.1)
s_data = structure.factor_py(qd_data, volume_fraction.magnitude)
plot the structure factor data and interpolated function
In [9]:
qd = np.arange(0,70, 0.1)# works up to qd = 72
s = structure.factor_data(qd, s_data, qd_data)
plt.figure()
plt.plot(qd, s, label = 'interpolated')
plt.plot(qd_data, s_data,'.', label = 'data')
plt.legend()
plt.xlabel('qd')
plt.ylabel('structure factor')
Out[9]:
Calculate reflectance
In [24]:
reflectance=np.zeros(len(wavelengths))
for i in range(len(wavelengths)):
reflectance[i],_,_,_,_ = sc.model.reflection(n_particle[i], n_matrix[i], n_medium[i], wavelengths[i],
radius, volume_fraction,
thickness=thickness,
structure_type='data',
structure_s_data=s_data,
structure_qd_data=qd_data)
plot
In [25]:
plt.figure()
plt.plot(wavelengths, reflectance)
plt.ylim([0,0.1])
plt.ylabel('Reflectance')
plt.xlabel('wavelength (nm)')
Out[25]:
set parameters
In [2]:
ntrajectories = 500 # number of trajectories
nevents = 500 # number of scattering events in each trajectory
wavelengths = sc.Quantity(np.arange(400, 800, 20), 'nm') # wavelengths
radius = sc.Quantity('0.5 um') # particle radius
volume_fraction = sc.Quantity(0.5, '') # volume fraction of particles
n_particle = ri.n('fused silica', wavelengths)
n_matrix = ri.n('vacuum', wavelengths) # called from the refractive_index module. n_matrix is the
n_medium = ri.n('vacuum', wavelengths) # space within sample. n_medium is outside the sample.
# n_particle and n_matrix can have complex indices if absorption is desired
boundary = 'film' # geometry of sample, can be 'film' or 'sphere', see below for tutorial
# on sphere case
thickness = sc.Quantity('50 um') # thickness of the sample film
Construct the structure factor data
Here, we use discrete points from the percus-yevick approximation for structure factor, as an example. In practice, you will most likely use actual structure factor data imported from your own file
In [3]:
qd_data = np.arange(0,75, 0.1)
s_data = structure.factor_py(qd_data, volume_fraction.magnitude)
plot the structure factor data and interpolated function
In [4]:
qd = np.arange(0,70, 0.1)# works up to qd = 72
s = structure.factor_data(qd, s_data, qd_data)
plt.figure()
plt.plot(qd, s, label = 'interpolated')
plt.plot(qd_data, s_data,'.', label = 'data')
plt.legend()
plt.xlabel('qd')
plt.ylabel('structure factor')
Out[4]:
Calculate reflectance
In [5]:
reflectance = np.zeros(wavelengths.size)
for i in range(wavelengths.size):
# calculate n_sample
n_sample = ri.n_eff(n_particle[i], n_matrix[i], volume_fraction)
# Calculate the phase function and scattering and absorption coefficients from the single scattering model
p, mu_scat, mu_abs = mc.calc_scat(radius, n_particle[i], n_sample, volume_fraction, wavelengths[i],
structure_type = 'data',
structure_s_data = s_data,
structure_qd_data = qd_data)
# Initialize the trajectories
r0, k0, W0 = mc.initialize(nevents, ntrajectories, n_medium[i], n_sample, boundary)
r0 = sc.Quantity(r0, 'um')
k0 = sc.Quantity(k0, '')
W0 = sc.Quantity(W0, '')
# Generate a matrix of all the randomly sampled angles first
sintheta, costheta, sinphi, cosphi, _, _ = mc.sample_angles(nevents, ntrajectories, p)
# Create step size distribution
step = mc.sample_step(nevents, ntrajectories, mu_scat)
# Create trajectories object
trajectories = mc.Trajectory(r0, k0, W0)
# Run photons
trajectories.absorb(mu_abs, step)
trajectories.scatter(sintheta, costheta, sinphi, cosphi)
trajectories.move(step)
reflectance[i], transmittance = det.calc_refl_trans(trajectories, thickness, n_medium[i], n_sample, boundary)
plot
In [6]:
plt.figure()
plt.plot(wavelengths, reflectance)
plt.ylim([0,1])
plt.ylabel('Reflectance')
plt.xlabel('wavelength (nm)')
Out[6]:
In [ ]: