# 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