dielectricfun.py

# -*- coding: utf-8 -*-
"""
Created on Mon Jan 11 09:31:32 2021

library to access n, k database. still in a very early form. At some point it 
would be good to have a common database format. But I am not sure what is the 
way to implement something like that. Also there are python implementations 
that enable access to the refractiveindex.info database. Now also some 
functions for calculations/rotations of dielectric tensors are included.

@author: m.rabe
"""
import specage.unit_convert as uc
import numpy as np
import pkg_resources
import sympy as _sp

DB_PATH = pkg_resources.resource_filename('specage', 'dielectricDB/')

def readnkdbfile(file, xunit):
    '''
    needs to be implemented

    Parameters
    ----------
    file : TYPE
        DESCRIPTION.
    xunit : TYPE
        DESCRIPTION.

    Returns
    -------
    None.

    '''
    pass

def writenkdbfile(file, x, n, k, source = '', cmt=''):
    '''
    needs to be implemented

    Parameters
    ----------
    file : TYPE
        DESCRIPTION.
    x : TYPE
        DESCRIPTION.
    n : TYPE
        DESCRIPTION.
    k : TYPE
        DESCRIPTION.
    source : TYPE, optional
        DESCRIPTION. The default is ''.
    cmt : TYPE, optional
        DESCRIPTION. The default is ''.

    Returns
    -------
    None.

    '''
    
    pass

def geBurnett_n(wavelength):
    '''
    n for Ge according to :
        J. H. Burnett, S. G. Kaplan, E. Stover, A. Phenis. 
        Refractive index measurements of Ge, Proc. SPIE 9974, 99740X (2016). 
        https://doi.org/10.1117/12.2237978

    Returns
    -------
    (n, nan)

    '''
    n=(1
      +0.4886331/(1-1.393959/wavelength**2)
      +14.5142535/(1-0.1626427/wavelength**2)
      +0.0091224/(1-752.190/wavelength**2)
      )**.5
    return n
    
def siChandler_Horowitz_n(wavelength):
    '''not functional at the moment '''
    # TODO: program that if needed
    # n=(11.67316
    #    +1/(x**2)
    #    +0.004482633/(x**2-1.108205**2))**.5
    pass

    
def D2O_MaxChapados_nk(wvlofinterest):
    fulldata=np.loadtxt(DB_PATH + 'nk_H2O_D2O_MaxChapados2009.txt', 
                        skiprows=4)
    #the dataset is in wavenumbers and the first value is 0, thus:
    wvl = uc.wavenum2nm(fulldata[1:,0])*1E-3
    
    #for the interpolation x must increase, thus it should be double checked
    if ~np.all(np.diff(wvl)>0):
        wvl = np.flipud(wvl)
        nk = np.flipud(fulldata[1:,3:5])
                
    return (np.interp(wvlofinterest, wvl, nk[:,0]),
            np.interp(wvlofinterest, wvl, nk[:,1]))

def H2O_MaxChapados_nk(wvlofinterest):
    fulldata=np.loadtxt(DB_PATH + 'nk_H2O_D2O_MaxChapados2009.txt', 
                        skiprows=4)
    #the dataset is in wavenumbers and the first value is 0, thus:
    wvl = uc.wavenum2nm(fulldata[1:,0])*1E-3
    
    #for the interpolation x must increase, thus it should be double checked
    if ~np.all(np.diff(wvl)>0):
        wvl = np.flipud(wvl)
        nk = np.flipud(fulldata[1:,1:3])
                
    return (np.interp(wvlofinterest, wvl, nk[:,0]),
            np.interp(wvlofinterest, wvl, nk[:,1]))

    
    
    
def nkdatabase(material, wavelength):
    '''
    return vavelength dependent refractive index and absorption coefficient 
    (n, k) values retrieved from database.

    Parameters
    ----------
    material : (str)
        Name of material. 
    wavelength : (numpy.array)
        Wavelength [micrometer].

    Returns
    -------
    n : (numpy.array)
        refractive index
    k : (numpy.array) 
        absorption coefficient

    '''
    # definition of the data base. this might be improved (maybe put in separate file?):
    
    alldatadict = {'Ge_n':geBurnett_n,
                   'Si_n':siChandler_Horowitz_n,
                   'D2O_nk':D2O_MaxChapados_nk,
                   'H2O_nk':H2O_MaxChapados_nk
                   }
    
    return alldatadict[material](wavelength)


def LorentzOscillator(w, eps_inf, S, w_res, fwhm):
    '''Returns epsilon for lorentz oscillator model with several resonances
    ... usage comes here
    works for frequencies (and wavenumbers) of any unit
    '''
    
    epsilon = eps_inf
    
    for (Sj, w_res_j, fwhm_j) in zip(S, w_res, fwhm):
        epsilon += (Sj*(w_res_j**2))/((w_res_j**2)-(w**2)-1j*w*fwhm_j)
        
    return epsilon

def LorentzianAbsorbance(kmax, nu_0, hwhh, n_inf, nu):
    """
    Returns refractive KK conform index of refraction (n) and absorption 
    coefficient from an antisymmetric linear combination of the Lorentzian 
    functions presented by Huang and Urban (Appl. Spectrosc., 1992, 46, 1666–1672.).

    Parameters
    ----------
    kmax : np.array
        maximum (peak) k values at nu_0
    nu_0 : np.array
        wavenumbers of k maxima (peaks)
    hwhh : np.array
        Half widths at half maximum (peak widths)
    n_inf : numeric scalar
        refractive index far away from peaks (background n)
    nu : numpy.array
        wavenumbers. 

    Returns
    -------
    n : numpy.array
        refractive indeces at nu
    k : numpy.array
        absorption coefficients at nu

    """
    k = np.zeros(nu.shape)
    n = n_inf*np.ones(nu.shape)
    for (km, u_0, hw) in zip(kmax, nu_0, hwhh):
        k += ( ((km*hw**2)/((nu-u_0)**2 + hw**2))
              -((km*hw**2)/((nu+u_0)**2 + hw**2)) )
        n += ( - (((nu-u_0)*km*hw)/((nu-u_0)**2 + hw**2)) 
               + (((nu+u_0)*km*hw)/((nu+u_0)**2 + hw**2)) )
    
    return n, k

def euler_rotmat(Phi, Theta, Psi):
    """
    Returns the orthogonal rotation matrix from the Euler angles φ, θ, and ψ
    (phi, theta, and psi) for active rotation. As defined in Goldstein, Poole 
    and Safko, "Classical Mechanics" (3rd Edition), 2002, p.153, eq(4.47). 
    This is the zx'z' convention of definition of Euler angles. Positive angles
    mean counter-clockwise rotations watching along the positive rotation axis 
    towards the origin. 

    Parameters
    ----------
    Phi : numeric or symbolic
        Azimuth angle in radians.
    Theta : numeric or symbolic
        Tilt angle in radians.
    Psi : numeric or symbolic
        Twist angle in radians.

    Returns
    -------
    sympy matrix
        Rotation matrix for active rotation by Euler angles.

    """
    phi, theta, psi = _sp.symbols('φ θ ψ')
    s=_sp.sin
    c=_sp.cos
    R = _sp.Matrix([
                    [c(psi)*c(phi) - c(theta)*s(phi)*s(psi), #1st row
                     -s(psi)*c(phi) - c(theta)*s(phi)*c(psi),
                     s(theta)*s(phi)],
                    [c(psi)*s(phi) + c(theta)*c(phi)*s(psi), # 2nd row
                     -s(psi)*s(phi) + c(theta)*c(phi)*c(psi),
                     -s(theta)*c(phi)],
                    [s(theta)*s(psi), #3rd row
                     s(theta)*c(psi),
                     c(theta)]])
    return R.subs([(phi, Phi), (theta, Theta), (psi, Psi)])


def euler_rotmat_zyz(Phi, Theta, Psi):
    """
    Returns the orthogonal rotation matrix from the Euler angles φ, θ, and ψ
    (phi, theta, and psi) for active rotation. As defined in Wilson, Decius, 
    Cross „Molecular Vibrations“ , 1955 p. 285. 
    This is the zy'z' convention of definition of Euler angles. Positive angles
    mean counter-clockwise rotations watching along the positive rotation axis 
    towards the origin.

    Parameters
    ----------
    Phi : numeric or symbolic
        Azimuth angle in radians.
    Theta : numeric or symbolic
        Tilt angle in radians.
    Psi : numeric or symbolic
        Twist angle in radians.

    Returns
    -------
    sympy matrix
        Rotation matrix for active rotation by Euler angles.

    """
    phi, theta, psi = _sp.symbols('φ θ ψ')
    s=_sp.sin
    c=_sp.cos
    R = _sp.Matrix([
                    [c(theta)*c(phi)*c(psi)-s(phi)*s(psi), #1st row
                     c(theta)*s(phi)*c(psi)+c(phi)*s(psi),
                     -s(theta)*c(psi)],
                    [-c(theta)*c(phi)*s(psi)-s(phi)*c(psi), # 2nd row
                     -c(theta)*s(phi)*s(psi)+c(phi)*c(psi),
                     s(theta)*s(psi)],
                    [s(theta)*c(phi), #3rd row
                     s(theta)*s(phi),
                     c(theta)]])
    return R.transpose().subs([(phi, Phi), (theta, Theta), (psi, Psi)]) #we need the transpose since we need it for active rotation