try:
   import pyfits
except:
   from astropy.io import fits as pyfits
import numpy as np
import matplotlib.pyplot as plt
import os

#file = "phase_sc_ft_normal.fits"
#file = "phase_sc_ft_1er.fits"
file = "phase_sc_ft_200um_1period.fits"
# file = "phase_sc_ft.fits"

hdulist = pyfits.open(file)

tel = [[0,1],[0,2],[0,3],[1,2],[1,3],[2,3]]


def smooth(y,n):
   w = np.ones(n)
   return np.convolve (w/w.sum(),y,mode='same')

def wrap(y):
   return np.angle (np.exp (1.j * y))

lbd_met = 1.908
lbd_sc  = 2.24
lbd_ft  = 2.2
tpi = np.pi * 2

for base in range(6):
    opd_sc  = hdulist['PHASE_SC'].data['PHASE'][base,:] / tpi * lbd_sc
    time_sc = hdulist['PHASE_SC'].data['TIME'][base,:] * 1e-6

    opd_ft  = hdulist['PHASE_FT'].data['PHASE'][base,:] / tpi * lbd_ft
    time_ft = hdulist['PHASE_FT'].data['TIME'][base,:] * 1e-6

    opl1 = hdulist['PHASE_MET'].data['OPL_FC_SC%i'%(tel[base][0]+1,)]
    opl2 = hdulist['PHASE_MET'].data['OPL_FC_SC%i'%(tel[base][1]+1,)]
    opd_met_sc = (opl1 - opl2) * 1e6
    opl1 = hdulist['PHASE_MET'].data['OPL_FC_FT%i'%(tel[base][0]+1,)]
    opl2 = hdulist['PHASE_MET'].data['OPL_FC_FT%i'%(tel[base][1]+1,)]
    opd_met_ft = (opl1 - opl2) * 1e6
    opl1 = hdulist['PHASE_MET'].data['OPD_FC%i'%(tel[base][0]+1,)]
    opl2 = hdulist['PHASE_MET'].data['OPD_FC%i'%(tel[base][1]+1,)]
    opd_met = (opl1 - opl2) * 1e6
    time_met = hdulist['PHASE_MET'].data['TIME'] * 1e-6

    # plt.plot (time_sc, opd_sc)
    # plt.plot (time_ft, opd_ft)
    # plt.plot (time_met, opd_met)

    opd_met = np.interp (time_sc, time_met, smooth(opd_met,10))
    opd_ft  = np.interp (time_sc, time_ft, smooth(opd_ft,10))
    opd_met_ft = np.interp (time_sc, time_met, smooth(opd_met_ft,10))
    opd_met_sc = np.interp (time_sc, time_met, smooth(opd_met_sc,10))
    
    #plt.plot (time_sc, opd_sc, 'r-')
    #plt.plot (time_sc, opd_met_sc, 'r--')
    #plt.plot (time_sc, opd_ft, 'b-')
    #plt.plot (time_sc, opd_met_ft, 'b--')
    #plt.plot (time_sc, opd_met, 'g-')
    #
    #plt.plot (time_sc, opd_sc - opd_met - opd_ft, 'r.')
    #plt.show()

    #plt.plot (time_sc, wrap (opd_sc/lbd_sc*tpi))
    #plt.plot (time_sc, opd_ft/lbd_ft*tpi)
    #plt.plot (time_sc, opd_met/lbd_met*tpi)

    # opd_met = np.interp (time_sc, time_met, opd_met)
    # opd_ft  = np.interp (time_sc, time_ft, opd_ft)

    # Residual phases
    phi_a = np.angle (np.exp (tpi*1.j * (opd_sc - opd_ft - opd_met)/lbd_sc))
    # phi_b = np.angle (np.exp (tpi*1.j * (opd_sc - opd_met_sc)/lbd_sc))
    plt.plot (time_sc, phi_a, "ob-")
    # plt.plot (time_sc, phi_b, "g")
    plt.show()

   #  plt.plot (tphase * 1e6, phi_unwrap)

hdulist.close()