# -*- coding: utf-8 -*-
"""
Created on Wed Apr 29 19:07:38 2015

@author: slacour
"""

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

dir = "../DATA/data_2015-05-26_Julien/"
sof = "p2vm_data.sof"

text=open(dir+sof,'r').readlines()
files={}
flat=0
dataf=0
for line in text:
    if line.split()[1] == 'DARK_RAW':
        files[line.split()[1]]=dir+line.split()[0]
    if line.split()[1] == 'FLAT_RAW':
        flat+=1
        files[line.split()[1]+str(flat)]=dir+line.split()[0]
    if line.split()[1] == 'WAVE':
        files[line.split()[1]]=dir+line.split()[0]
    if line.split()[1] == 'DARK':
        files[line.split()[1]]=dir+line.split()[0]
    if line.split()[1] == 'FLAT':
        files[line.split()[1]]=dir+line.split()[0]
    if line.split()[1] == 'P2VM':
        files[line.split()[1]]=dir+line.split()[0]
    if line.split()[1] == 'BAD':
        files[line.split()[1]]=dir+line.split()[0]
    if line.split()[1] == 'DUAL_SCIENCE_RAW':
        dataf+=1
        files[line.split()[1]+str(dataf)]=dir+line.split()[0]
        

f1= files['FLAT_RAW1']
f2= files['FLAT_RAW2']
f3= files['FLAT_RAW3']
f4= files['FLAT_RAW4']

dark=files['DARK_RAW']

dark_sc = fits.getdata(dark,"IMAGING_DATA_SC")
dark_ft = fits.getdata(dark,"IMAGING_DATA_FT").field("PIX")


# FOR DARK CALCULATION 

median_dark_sc=median(dark_sc,axis=0)
median_dark_ft=median(dark_ft,axis=0)
var_dark_sc=mean(((dark_sc-mean(dark_sc,axis=0))**2),axis=0)
var_dark_ft=mean(((dark_ft-mean(dark_ft,axis=0))**2),axis=0)

# FOR DARK QC PARAMETERS

QC_med_sc=median(median_dark_sc)
QC_rms_sc=sqrt(median(var_dark_sc))
print("QC Median SC = "+str(QC_med_sc))
print("QC Median RMS SC = "+str(QC_rms_sc))

QC_mea_ft=mean(median_dark_ft)
QC_rms_ft=sqrt(mean(var_dark_ft))
print("QC Mean FT = "+str(QC_mea_ft))
print("QC Mean RMS FT = "+str(QC_rms_ft))

sc1 = fits.getdata(f1,"IMAGING_DATA_SC")-median_dark_sc
sc2 = fits.getdata(f2,"IMAGING_DATA_SC")-median_dark_sc
sc3 = fits.getdata(f3,"IMAGING_DATA_SC")-median_dark_sc
sc4 = fits.getdata(f4,"IMAGING_DATA_SC")-median_dark_sc

mean_sc1=sc1.mean(axis=0)
var_sc1=((sc1-mean_sc1)**2).mean(axis=0)
mean_sc2=sc2.mean(axis=0)
var_sc2=((sc2-mean_sc2)**2).mean(axis=0)
mean_sc3=sc3.mean(axis=0)
var_sc3=((sc3-mean_sc3)**2).mean(axis=0)
mean_sc4=sc4.mean(axis=0)
var_sc4=((sc4-mean_sc4)**2).mean(axis=0)

ft1 = fits.getdata(f1,"IMAGING_DATA_FT").field("PIX")-median_dark_ft
ft2 = fits.getdata(f2,"IMAGING_DATA_FT").field("PIX")-median_dark_ft
ft3 = fits.getdata(f3,"IMAGING_DATA_FT").field("PIX")-median_dark_ft
ft4 = fits.getdata(f4,"IMAGING_DATA_FT").field("PIX")-median_dark_ft

mean_ft1=ft1.mean(axis=0)
var_ft1=((ft1-mean_ft1)**2).mean(axis=0)
mean_ft2=ft2.mean(axis=0)
var_ft2=((ft2-mean_ft2)**2).mean(axis=0)
mean_ft3=ft3.mean(axis=0)
var_ft3=((ft3-mean_ft3)**2).mean(axis=0)
mean_ft4=ft4.mean(axis=0)
var_ft4=((ft4-mean_ft4)**2).mean(axis=0)


# FOR BAD PIXEL CALCULATION

ran_mean  = QC_med_sc + 10 * QC_rms_sc
ran_mean2 = QC_med_sc - 10 * QC_rms_sc  
ran_var   = QC_rms_sc**2 * 10
ran_var2  = QC_rms_sc**2 * 0.1
    
bad_sc= median_dark_sc>ran_mean
bad_sc|=median_dark_sc<ran_mean2

bad_sc|=var_dark_sc   >ran_var
bad_sc|=var_dark_sc  <ran_var2


print("QC BadPix SC = "+str(sum(bad_sc)))    

Mean_sc=array((mean_sc1[~bad_sc],mean_sc2[~bad_sc],mean_sc3[~bad_sc],mean_sc4[~bad_sc]))
Var_sc=array((var_sc1[~bad_sc],var_sc2[~bad_sc],var_sc3[~bad_sc],var_sc4[~bad_sc]))
Gain1=(mean(Var_sc*Mean_sc)-mean(Var_sc)*mean(Mean_sc))/(mean(Mean_sc**2)-mean(Mean_sc)*mean(Mean_sc))

Mean_ft=array((mean_ft1,mean_ft2,mean_ft3,mean_ft4))
Var_ft=array((var_ft1,var_ft2,var_ft3,var_ft4))
Gain2=(mean(Var_ft*Mean_ft)-mean(Var_ft)*mean(Mean_ft))/(mean(Mean_ft**2)-mean(Mean_ft)*mean(Mean_ft))


print("QC MeanGain FT = "+str(Gain2))
print("QC MeanGain SC = "+str(Gain1))

profiles = [fits.getdata(dir+"gravi_profile_map.fits","PROFILE_DATA").field(i) for i in range(48)]
profiles=array(profiles)[:,0,]