#
#--- module polar (atelier obspm, 18-21 dec 2006)
#
#--- written by fpaletou@ast.obs-mip.fr (oct 2006)
#
#------modified DEC 22, 2006 (mirror function added)
#
#------modified FEB 19, 2008 for numpy/scipy
#
#------modified 2018 (frederic.paletou@univ-tlse3.fr)
#
#      python3+ compliant now...

######################################
# import necessary external packages #
######################################

import numpy as np
import matplotlib.pyplot as plt

##################################################################
# returns the mueller matrix of a retarder of (phase) - DEGREES #
##################################################################
def matphase(phase):

    # converts DEGREES (input) to radians
    phase=phase*np.pi/180.

    mph =[ [1., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., np.cos(phase), np.sin(phase)],
           [0., 0.,-np.sin(phase), np.cos(phase)] ]

    return mph

##################################################
# returns a rotation matrix (2*theta) - DEGREES #
##################################################
def matrot(theta):

    # converts DEGREES (input) to radians
    theta=theta*np.pi/180.

    mrot  =[ [1., 0., 0., 0.],
             [0., np.cos(2.*theta), np.sin(2.*theta), 0.],
             [0.,-np.sin(2.*theta), np.cos(2.*theta), 0.],
             [0., 0., 0., 1.] ]

    return mrot

##############################################################################
# returns the Muller matrix of a Al inclined mirror in the 0.47-2 mics range #
##############################################################################
def mirror(inc, londe):

    # LONDE is wavelength in nanometers!!

    # converts DEGREES (input) to radians
    inc=inc*np.pi/180.

    # complex refractive index for Aluminium [from LIDE, Handbook]

    n=np.array([2.273, 1.770, 1.444, 1.264, 1.212, 1.201, 1.260, 1.468, 2.237,
       2.745, 2.625, 2.143, 1.741, 1.488, 1.304, 1.018, 0.826, 0.695])

    k=np.array([21.403, 18.328, 15.955, 14.021, 12.464, 11.181, 10.01, 8.949,
       8.212, 8.309, 8.597, 8.573, 8.205, 7.821, 7.479, 6.846, 6.283, 5.800])

    ee=np.array([0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
       1.9, 2., 2.2, 2.4, 2.6])

    # wavelength l is in nanometers
    l=1.2398/ee

    # interpolate LONDE in tabulated l-range WARNING WITH BOUNDARIES
    i=0
    while (l[i] > londe):
#        print(l[i], londe)
        i=i+1

    # simple linear interpolation on real/im parts of complex indices
    nn=0.5*(n[i]+n[i-1])
    kk=0.5*(k[i]+k[i-1])

    p=nn*nn - kk*kk - np.sin(inc)*np.sin(inc)
    q=2.*nn*kk
    
    rp=np.sqrt( p+np.sqrt(p*p+q*q))/np.sqrt(2.) 
    rm=np.sqrt(-p+np.sqrt(p*p+q*q))/np.sqrt(2.) 
    
    s=np.sin(inc)*np.tan(inc)
    
    rho2=(np.sqrt(p*p+q*q)+s*s-2.*s*rp)/(np.sqrt(p*p+q*q)+s*s+2.*s*rp)
    rho=np.sqrt(rho2)

    tdelta = 2.*s*rm/(np.sqrt(p*p+q*q)-s*s)
    delta=np.arctan(tdelta)

    mmir  =np.array([ [1.+rho2, 1.-rho2, 0., 0.],
             [1.-rho2, 1.+rho2, 0., 0.],
             [0., 0., -2.*rho*np.cos(delta), -2.*rho*np.sin(delta)],
             [0., 0.,  2.*rho*np.sin(delta), -2.*rho*np.cos(delta)]])

    mmir=0.5*mmir
 
    return mmir

#############################
# l/4 polarimeter (azimuth) #
#############################
def l4mod(npts):

    az=np.zeros((npts), 'Float32')
    xq=np.zeros((npts), 'Float32')
    xu=np.zeros((npts), 'Float32')
    xv=np.zeros((npts), 'Float32')
    
    for i in np.arange(npts):
        az[i]=180.*i/(npts-1.)
        mpol=np.dot(matrot(-az[i]),np.dot(matphase(90.),matrot(az[i])))
        xq[i]=mpol[1,1]
        xu[i]=mpol[1,2]
        xv[i]=mpol[1,3]

    print("|-------------------------------------|")
    print("| USAGE : xx[0] = tableau des azimuts |")
    print("|         xx[1] = tableau des Xq      |")
    print("|         xx[2] = tableau des Xu      |")
    print("|         xx[3] = tableau des Xv      |")
    print("|                                     |")
    print("| GRAPHIQUES:                         |")
    print("|   plt.plot(xx[0],xx[3])             |")
    print("|   plt.show()                        |")
    print("|-------------------------------------|")

    return az, xq, xu, xv

#############################
# l/2 polarimeter (azimuth) #
#############################
def l2mod(npts):

    az=np.zeros((npts), 'Float32')
    xq=np.zeros((npts), 'Float32')
    xu=np.zeros((npts), 'Float32')
    xv=np.zeros((npts), 'Float32')
    
    for i in np.arange(npts):
        az[i]=180.*i/(npts-1.)
        mpol=np.dot(matrot(-az[i]),np.dot(matphase(180.),matrot(az[i])))
        xq[i]=mpol[1,1]
        xu[i]=mpol[1,2]
        xv[i]=mpol[1,3]

    print("|-------------------------------------|")
    print("| USAGE : xx[0] = tableau des azimuts |")
    print("|         xx[1] = tableau des Xq      |")
    print("|         xx[2] = tableau des Xu      |")
    print("|         xx[3] = tableau des Xv      |")
    print("|                                     |")
    print("| GRAPHIQUES:                         |")
    print("|   plt.plot(xx[0],xx[3])             |")
    print("|   plt.show()                        |")
    print("|-------------------------------------|")

    return az, xq, xu, xv

#############################################
# THeMIS polarimeter (theta1,theta2,retard) #
#############################################
#
# - CAUTION: angles should be in DEGREES
#
# - returns MPOL==(4x4) before analyser so that
#   elements of 2nd row, 3 last columns MPOL[1, 1:3] contains
#   modulation paramaters for Q, U and V
#
# written by fpaletou@ast.obs-mip.fr [oct 2006]
#
#
def th_pol(theta1,theta2,faz):

    m1=np.dot(matrot(-theta1),np.dot(matphase(faz),matrot(theta1)))
    m2=np.dot(matrot(-theta2),np.dot(matphase(faz),matrot(theta2)))
    m_pol_th=np.dot(m2,m1)

    print("----------------------------------------------------")
    print("  Matrice de Mueller du polarimetre de THeMIS pour")
    print("  ** theta_1 =",theta1,"** theta_2 =",theta2,"** retard =",faz) 
    print("----------------------------------------------------")

    print(np.array_str(m_pol_th,precision=3,suppress_small=1))

    return m_pol_th

#######################################################################
# computes theoretical retardance of 400-700 nm Fichou achromatic QWP #  
#######################################################################
#
# - CAUTION: wavelength should be in Angstroems!
#
# - returns PHASE in degrees
#
# written by fpaletou@ast.obs-mip.fr [jul 2006]
#
def phase_400700(wavelength):

    #- indices for MgF2 plate
    no_mgf2= 1.36957 + 35.821/(wavelength - 1492.5)
    ne_mgf2= 1.381   + 37.415/(wavelength - 1494.7)
    
    #-indices for quartz plate
    no_qrtz= 1.52654 + 77.324/(wavelength - 1521.02)
    ne_qrtz= 1.535   + 80.288/(wavelength - 1514.10)

    #- birefringence for each material
    mu_mgf2= ne_mgf2 - no_mgf2
    mu_qrtz= ne_qrtz - no_qrtz

    #- nominal thickness of elementary plates
    ee_mgf2=202.438 * 1.e4
    ee_qrtz=245.664 * 1.e4
    
    ph_mgf2=360. * mu_mgf2 * ee_mgf2 / wavelength
    ph_qrtz=360. * mu_qrtz * ee_qrtz / wavelength
    
    #- crossed elementary plates gives final phase
    phase=ph_mgf2 - ph_qrtz

    return phase