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)
#
#file=files["DUAL_SCIENCE_REDUCED1"]
#file_f=d+"vis_dual_science.fits"
#
#
#dir="/Users/slacour/GRAVITY/Data/2015Aout10/p2vm_raw/"
#file_p2vm="p2vm_map.2015-08-07T15-32-53.fits"

f1="GRAVITY.2015-08-07T15-12-31.fits"

file_raw=dir+"gravi_preproc_"+f1
file=dir+"single_reduced_"+f1
file_f=dir+"vis_dual_science.fits"
#
#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.905e-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
WS1  = fits.getdata(file,4).field('EFF_WAVE')
WF1  = fits.getdata(file,6).field('EFF_WAVE')
WS2  = fits.getdata(file,5).field('EFF_WAVE')
WF2  = fits.getdata(file,7).field('EFF_WAVE')

# Enfin les "time stamps"
SC1T = fits.getdata(file,9).field('TIME')[::6]
FT1T = fits.getdata(file,15).field('TIME')[::6]
SC2T = fits.getdata(file,12).field('TIME')[::6]
FT2T = fits.getdata(file,18).field('TIME')[::6]
# Puis les donnees SC (pola1 &2)
Nw=fits.getdata(file,9).field('VISDATA').shape[1]
S1  = fits.getdata(file,9).field('VISDATA').reshape(len(SC1T),6,Nw)
S2  = fits.getdata(file,12).field('VISDATA').reshape(len(SC1T),6,Nw)
# Puis les donnees FT (pola1 &2)
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())]

# Puis je re-echantillone SC (mais on ne devrait pas en avoir dans le pipeline)
#S1W=SC1[(SC1T>FT1T.min())&(SC1T<FT1T.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))])

F1W=array([[angle(F1[i,j].mean())+np.poly1d(np.polyfit(WF1, angle(F1[i,j]*conj(F1[i,j].mean())), 2))(WS1) for j in range(6)] for i in range(len(T))])
F2W=array([[angle(F2[i,j].mean())+np.poly1d(np.polyfit(WF2, angle(F2[i,j]*conj(F2[i,j].mean())), 2))(WS2) for j in range(6)] for i in range(len(T))])

#S1[:,0,:]=conj(S1[:,0])
#S2[:,0,:]=conj(S2[:,0])

#F1C=F1.mean(axis=0).conj()
#F1*=F1.mean(axis=0).conj()
#a=angle(F1)

# 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-M.mean(axis=0))*Lambda_met

# VISDATA est "l'ancien" VISDATA
phi1= - F1W + (OPD_Met[:,:,None])/WS1
phi2= - F2W + (OPD_Met[:,:,None])/WS2
VISDATA1 = (S1*exp(1j* phi1)).mean(axis=0)  
VISDATA2 = (S2*exp(1j* phi2)).mean(axis=0)  
VISDATA1_T = (S1*exp(1j* (-F1W+M[:,:,None]*Lambda_met/WS1))).mean(axis=0)  
VISDATA2_T = (S2*exp(1j* (-F2W+M[:,:,None]*Lambda_met/WS2))).mean(axis=0)  

GD1=(VISDATA1[:,1:]*conj(VISDATA1[:,:-1])).sum(axis=1)
GD2=(VISDATA2[:,1:]*conj(VISDATA2[:,:-1])).sum(axis=1)

#OPD_GD est l'OPD qui correspond a la pente de phase.
OPD_GD1=(angle(GD1)/diff(1/WS1).mean())[:,None]/WS1
OPD_GD2=(angle(GD2)/diff(1/WS2).mean())[:,None]/WS2

VISDATA1_2=(VISDATA1*exp(-1j*OPD_GD1))
VISDATA1_2_Mean=VISDATA1_2.mean(axis=1) 


VISDATA2_2=(VISDATA2*exp(-1j*OPD_GD2))
VISDATA2_2_Mean=VISDATA2_2.mean(axis=1) 

# VISDATA2_Mean correspond a une moyenne generale sur la phase

VISDATA1_3=VISDATA1_2*exp(-1j* angle(VISDATA1_2_Mean)[:,None])
VISDATA2_3=VISDATA2_2*exp(-1j* angle(VISDATA2_2_Mean)[:,None])

Vis_Phi1=angle(VISDATA1_3)+angle(VISDATA1_2_Mean)[:,None]+OPD_Mean[:,None]/WS1 + OPD_GD1
Vis_Phi2=angle(VISDATA2_3)+angle(VISDATA2_2_Mean)[:,None]+OPD_Mean[:,None]/WS2 + OPD_GD2

#Vis_Amp=abs(VISDATA)
#Vis_Data=Vis_Amp*exp(1j*Vis_Phi)



#Check:
    
VIS_PHI1  = fits.getdata(file_f,11).field('VISPHI')
VIS_PHI2  = fits.getdata(file_f,19).field('VISPHI')

# TEST
S1met=array([interp(metT,SC1T,unwrap(angle(S1)[:,i,150])) for i in range(6)]).T
F1met=array([interp(metT,FT1T,unwrap(angle(FT1)[:,i,2])) for i in range(6)]).T

S2met=array([interp(metT,SC2T,unwrap(angle(S2)[:,i,150])) for i in range(6)]).T
F2met=array([interp(metT,FT2T,unwrap(angle(FT2)[:,i,2])) for i in range(6)]).T

clf()
plot(metT,met[:,1])
plot(metT,S1met[:,1])
plot(metT,F1met[:,1])


print(M.mean(axis=0)*Lambda_met/WS1[100])