# Tutorial for using structure factor data as the structure factor used in the structural-color package

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

``````

### For the single scattering model

set parameters

``````

In [7]:

wavelengths = sc.Quantity(np.arange(400, 800, 20), 'nm') # wavelengths
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]:

<matplotlib.text.Text at 0x119086f28>

``````

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],
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]:

<matplotlib.text.Text at 0x1195789e8>

``````

### For the Monte Carlo model

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
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]:

<matplotlib.text.Text at 0x118fab978>

``````

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]:

<matplotlib.text.Text at 0x119032f60>

``````
``````

In [ ]:

``````