snlay.py 1.5 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950
  1. # this module will allow us to call scattnlay
  2. import numpy as np
  3. #import PyMieScatt as ps
  4. #import matplotlib.pyplot as plt
  5. from scattnlay import scattnlay, fieldnlay
  6. from scipy import interpolate
  7. #make a materials dictionary
  8. matsdict = {
  9. 1: './materials/gold.dat',
  10. 2: './materials/silicon.dat',
  11. 3: './materials/silica.dat',
  12. 4: './materials/tio2.dat',
  13. 5: './materials/silver.dat'
  14. }
  15. def get_nk(datafile, wavelengths):
  16. """Reads the given file and returns the n+ik complex at
  17. the given wavelength after suitable interpolation
  18. :datafile: TODO
  19. :wavelength: TODO
  20. :returns: TODO
  21. """
  22. rawdisp = np.loadtxt(datafile)
  23. f_r = interpolate.interp1d(rawdisp[:,0], rawdisp[:,1])
  24. f_i = interpolate.interp1d(rawdisp[:,0], rawdisp[:,2])
  25. return f_r(wavelengths) + 1j*f_i(wavelengths)
  26. # make a function that will return the x and m matrices of
  27. # any given multilayered spherical particle
  28. # these can then be passed on to scattnlay
  29. def make_xm(sizes, mats, lams):
  30. x = np.ones((len(lams), len(sizes)), dtype = np.float64)
  31. m = x - x + 0 + 0*1j
  32. sizes = np.cumsum(sizes)
  33. for s, size in enumerate(sizes):
  34. x[:,s] = np.pi*size/lams
  35. m[:,s] = get_nk(matsdict[mats[s]], lams)
  36. #theta = np.linspace(0.0, 180, num_pts, dtype = np.float64)*np.pi/180.0
  37. return x, m
  38. def calc_spectrum(sizes, mats, lams):
  39. x, m = make_xm(sizes, mats, lams)
  40. terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
  41. return Qext