def get_files(d,sof):
    text=open(d+sof,'r').readlines()
    files={}
    flat=0
    dataf=0
    for line in text:
        if len(line.split()) == 2:
            if line.split()[1] == 'DARK_RAW':
                files[line.split()[1]]=d+line.split()[0]
            if line.split()[1] == 'FLAT_RAW':
                flat+=1
                files[line.split()[1]+str(flat)]=d+line.split()[0]
            if line.split()[1] == 'WAVE':
                files[line.split()[1]]=d+line.split()[0]
            if line.split()[1] == 'DARK':
                files[line.split()[1]]=d+line.split()[0]
            if line.split()[1] == 'FLAT':
                files[line.split()[1]]=d+line.split()[0]
            if line.split()[1] == 'P2VM':
                files[line.split()[1]]=d+line.split()[0]
            if line.split()[1] == 'BAD':
                files[line.split()[1]]=d+line.split()[0]
            if line.split()[1] == 'DUAL_SCIENCE_RAW':
                dataf+=1
                files[line.split()[1]+str(dataf)]=d+line.split()[0]
                files["DUAL_SCIENCE_PREPROC"+str(dataf)]=d+"gravi_preproc_"+line.split()[0]
                files["DUAL_SCIENCE_REDUCED"+str(dataf)]=d+"dual_reduced_"+line.split()[0]
    return files
        
        
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 29 19:07:38 2015

@author: slacour
"""
    
from numpy import *   
import cmath
from matplotlib.pyplot import *    
import pyfits as fits
from scipy.interpolate import interp1d


from numpy import *   
from matplotlib.pyplot import *    
import pyfits as fits

d = "../DATA/data_2015-05-26_Julien/"
sof = "science_dual_data.sof"

files=get_files(d,sof)

dir="/Users/slacour/DRS/DATA/data_2015-04-30_opd_stability/"

file=files["DUAL_SCIENCE_REDUCED1"]
#
#res="vis_dual_science.fits"
#Phi  = angle(fits.getdata(dir+res,10).field('VISDATA'))
#Phi2  = angle(fits.getdata(dir+res,18).field('VISDATA'))
Lambda_met=1.908e-6

M_matrix = -1*array([1.,-1.,0.0,0.0,-1.,0.0,1.,0.0,-1.,0.0,0.0,1.,0.0,-1.,1.,0.0,0.0,-1.,0.0,1.,0.0,0.0,-1.,1.]);
M_matrix = M_matrix.reshape((6,4))
    
MJD=fits.getheader(file,0)['MJD-OBS']
print((fits.getheader(file,0)['DATE-OBS']))

# Je lis la metrologie
met = dot(M_matrix,fits.getdata(file,8).field('PHI').T).T
metT= fits.getdata(file,8).field('TIME')

# Puis les longueurs d'onde
Wave_SC1  = fits.getdata(file,4).field('EFF_WAVE')
Wave_FT1  = fits.getdata(file,6).field('EFF_WAVE')

# Enfin les "time stamps"
SC1T = fits.getdata(file,9).field('TIME')[::6]
FT1T = fits.getdata(file,15).field('TIME')[::6]
# Puis les donnees SC (pola1)
SC1  = fits.getdata(file,9).field('VISDATA').reshape(len(SC1T),6,241)
SC2  = fits.getdata(file,12).field('VISDATA').reshape(len(SC1T),6,241)
# Puis les donnees FT (pola1)
FT1  = fits.getdata(file,15).field('VISDATA').reshape(len(FT1T),6,5)
FT2  = fits.getdata(file,18).field('VISDATA').reshape(len(FT1T),6,5)

# Je selectionne l'echantillonage temporel    
T=SC1T[(SC1T>FT1T.min())&(SC1T<FT1T.max())]
# Et l'echantillonage en longeur d'onde (mais ici, plutot faire un fit, c'est mieux.)
W1=Wave_SC1[(Wave_SC1>Wave_FT1.min())&(Wave_SC1<Wave_FT1.max())]

# Puis je re-echantillone SC (mais on ne devrait pas en avoir dans le pipeline)
S1W=SC1[(SC1T>FT1T.min())&(SC1T<FT1T.max())]
S1W=SC1[(SC1T>FT1T.min())&(SC1T<FT1T.max())][:,:,(Wave_SC1>Wave_FT1.min())&(Wave_SC1<Wave_FT1.max())]

DIT=diff(T).mean()

# Et FT (ca, on en a besoin)
F1=array([FT1[(FT1T>T[t])&(FT1T<T[t]+DIT),:,:].mean(axis=0) for t in range(len(T))])
F2=array([FT2[(FT1T>T[t])&(FT1T<T[t]+DIT),:,:].mean(axis=0) for t in range(len(T))])
#F1C=F1.mean(axis=0).conj()
#F1*=F1.mean(axis=0).conj()
#a=angle(F1)

F1W=array([[interp1d(Wave_FT1,F1[t,b])(W1) for b in range(6)] for t in range(len(T))])

# Et enfin la metrologie
M=array([met[(metT>T[t])&(metT<T[t]+DIT),:].mean(axis=0) for t in range(len(T))])

# La, les choses serieuses commencent. C'est un peu complique a expliquer. Mais il suffit de faire cela:

#OPD_Mean est l'OPD moyenne sur le temps qui correspond a la metrologie.
OPD_Mean=M.mean(axis=0)*Lambda_met

#OPD_Met est l'OPD qui correspond a la metrologie.
OPD_Met=M*Lambda_met-OPD_Mean

# VISDATA est "l'ancien" VISDATA
phi0_ref= angle(F1W) - OPD_Met[:,:,None]/W1
VISDATA = (S1W*exp(1j* phi0_ref)).mean(axis=0)  

GD=(VISDATA[:,1:]*conj(VISDATA[:,:-1])).sum(axis=1)

#OPD_GD est l'OPD qui correspond a la pente de phase.
OPD_GD=(angle(GD)/diff(1/W1).mean())[:,None]/W1

VISDATA2=(VISDATA*exp(-1j*OPD_GD))
VISDATA2_Mean=VISDATA2.mean(axis=1) 
# VISDATA2_Mean correspond a une moyenne generale sur la phase

VISDATA3=VISDATA2*exp(-1j* angle(VISDATA2_Mean)[:,None])

Vis_Phi=angle(VISDATA3)+angle(VISDATA2_Mean)[:,None]-OPD_Mean[:,None]/W1 + OPD_GD
Vis_Amp=abs(VISDATA)
Vis_Data=Vis_Amp*exp(1j*Vis_Phi)