1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950 |
- # this module will allow us to call scattnlay
- import numpy as np
- #import PyMieScatt as ps
- #import matplotlib.pyplot as plt
- from scattnlay import scattnlay, fieldnlay
- from scipy import interpolate
- #make a materials dictionary
- matsdict = {
- 1: './materials/gold.dat',
- 2: './materials/silicon.dat',
- 3: './materials/silica.dat',
- 4: './materials/tio2.dat',
- 5: './materials/silver.dat'
- }
- def get_nk(datafile, wavelengths):
- """Reads the given file and returns the n+ik complex at
- the given wavelength after suitable interpolation
- :datafile: TODO
- :wavelength: TODO
- :returns: TODO
- """
- rawdisp = np.loadtxt(datafile)
- f_r = interpolate.interp1d(rawdisp[:,0], rawdisp[:,1])
- f_i = interpolate.interp1d(rawdisp[:,0], rawdisp[:,2])
- return f_r(wavelengths) + 1j*f_i(wavelengths)
- # make a function that will return the x and m matrices of
- # any given multilayered spherical particle
- # these can then be passed on to scattnlay
- def make_xm(sizes, mats, lams):
- x = np.ones((len(lams), len(sizes)), dtype = np.float64)
- m = x - x + 0 + 0*1j
- sizes = np.cumsum(sizes)
- for s, size in enumerate(sizes):
- x[:,s] = np.pi*size/lams
- m[:,s] = get_nk(matsdict[mats[s]], lams)
- #theta = np.linspace(0.0, 180, num_pts, dtype = np.float64)*np.pi/180.0
- return x, m
- def calc_spectrum(sizes, mats, lams):
- x, m = make_xm(sizes, mats, lams)
- terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
- return Qext
|