# -*- coding: iso-8859-15 -*-ls
#! /usr/bin/env python
"""
Created on Fri Jul 17 2015

@author: kervella
"""

# ATTENTION: necessary to install pdfrw library (sudo port install py27-pdfrw)

import matplotlib as mpl
mpl.use('Agg') # necessary to be able to generate the reports even without X server

import matplotlib.pyplot as plt
import sys
import datetime
import numpy as np
from . import gravi_visual_class

try:
    from scipy.signal import welch
    welchpresent = True
except ImportError as e:
    welchpresent = False
    pass # module doesn't exist, deal with it.

#==============================================================================
# Preparation of the report using reportlab
#==============================================================================
def produce_raw_report(raw,filename):
    # From gravi_visual_class
    from . import gravi_visual_class
    from .gravi_visual_class import PdfImage, myFirstPage, myLaterPages, styles
    from .gravi_visual_class import clipdata, transbarchartout, graphoutaxes, graphoutnoaxis
    from .gravi_visual_class import graphoutnoxaxes, graphscatteraxes, graphaxesplt
    from .gravi_visual_class import plotTitle, plotSubtitle
    from .gravi_visual_class import get_key_withdefault, create_array_from_list, baseline_phases
    from .gravi_visual_class import mean_angle, std_angle, clean_unwrap_phase, nanaverage
    from .gravi_visual_class import clean_gdelay_is, clean_gdelay_fft, clean_gdelay_full
    from .gravi_visual_class import base_list, tel_list, nbase, ntel
    from .gravi_visual_class import plotReductionSummary, plotSummary

    # From reportlab
    from reportlab.platypus import SimpleDocTemplate, Paragraph, Spacer, Table, TableStyle, Image, PageBreak
    from reportlab.graphics.shapes import String
    from reportlab.lib.units import inch, cm, mm
    from reportlab.lib import colors
    
    # import PIL
    import io
    from io import StringIO
    
    #==============================================================================
    # Global parameters
    #==============================================================================
    
    satlevel_sc = 28000 # Saturation level of SC in ADU
    pixwidth = 1500 # Number of elements for plots in horizontal direction

    ftfreq = 1./np.mean(np.diff(raw.time_ft))


    if 'HIERARCH ESO INS MLC WAVELENG' in raw.header:
        Lambda_met=raw.header['HIERARCH ESO INS MLC WAVELENG']/1000. # Metrology wavelength in microns (1.908 microns)
    else: Lambda_met = 1.908
        
    # Subtraction of the initial time to avoid full "MJD seconds" in graphs ALL SYNC TO FT TIME
    raw.metrology_time -=raw.time_ft[0]
    raw.fddl_time -= raw.time_ft[0]
    raw.time_sc -= raw.time_ft[0]
    raw.time_ft -= raw.time_ft[0]
    
    #==============================================================================
    # Start Story
    #==============================================================================
    
    Story = [Spacer(1,1*mm)]
    plotSummary (Story, filename, raw, onTarget=True)
    Story.append(Spacer(1,3*mm))

    #==============================================================================
    # Table with checks
    #==============================================================================
    
    # Saturation checks:
    saturation_flag = max(np.percentile(np.percentile(raw.data_sc,99,axis=0),99,axis=0)) > satlevel_sc
    saturation      = 'YES' if saturation_flag else 'No'
 
    # Data structure basic checks:
    quality_flag = False  if raw.datacheck == "# Good" else True
    quality      = "GOOD" if quality_flag else "BAD"

    # Colored table
    data = [['Saturation check:', saturation], ['File basic checks:', quality]]
    t=Table(data, colWidths=[5*cm, 11*cm])
    t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black),
                           ('BOX', (0,0), (-1,-1), 2, colors.black),]))
    t.setStyle(TableStyle([('TEXTCOLOR',(1,0),(1,0), colors.red if saturation_flag else colors.lightgreen)]))
    t.setStyle(TableStyle([('TEXTCOLOR',(1,1),(1,1), colors.lightgreen if quality_flag else colors.red)]))
    Story.append(t)
    Story.append(Spacer(1,3*mm))
    
    # Other warnings
    if (quality == "BAD"):
        Story.append(Spacer(1,5*mm))
        for a in raw.datacheck:
            p = Paragraph(a, styles["Normal"])
            Story.append(p)
            Story.append(Spacer(1,5*mm))

    Story.append(PageBreak())
    
    #==============================================================================
    # Image display of the acquisition camera average frame
    #==============================================================================
    
    plotTitle(Story, "Average of acquisition camera images")
    plotSubtitle(Story,"Color scale: square root of the image values in ADU.")

    # replace negative values by zero for plot
    raw.acq[raw.acq<0]=0
    minval = np.sqrt(np.percentile(np.nan_to_num(raw.acq),5))
    maxval = np.sqrt(np.percentile(np.nan_to_num(raw.acq),99.9))

    imgdata = io.StringIO()
    plt.close('all')
    fig=plt.figure(figsize=(3,6),dpi=300)
    plt.rc('xtick', labelsize=5)
    plt.rc('ytick', labelsize=5)
    figaspect=1
    
    if raw.acq.shape[0]==1: plotimg = np.sqrt(np.squeeze(raw.acq)) # if only one frame
    if raw.acq.shape[0]>1: plotimg = np.sqrt(np.mean(raw.acq,axis=0)) # several frames

    # figaspect = (1.0*raw.acq.shape[2]/raw.acq.shape[1])*2.0
    plt.imshow(plotimg,vmin=minval,vmax=maxval,cmap='cubehelix_r', interpolation='nearest', aspect=figaspect)
    #circ = plt.Circle((790,1400), radius=100, color='r', linewidth=0.5, fill=False)

    rect = plt.Rectangle((130,60), 300, 1450, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV1',xy=(200, 45),  xycoords='data', fontsize=3, color='r')
    fig.gca().add_artist(text)

    rect = plt.Rectangle((630,60), 300, 1450, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV2',xy=(700, 45),  xycoords='data', fontsize=3, color='r')
    fig.gca().add_artist(text)

    rect = plt.Rectangle((1070,60), 300, 1450, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV3',xy=(1140, 45),  xycoords='data', fontsize=3, color='r')
    fig.gca().add_artist(text)

    rect = plt.Rectangle((1520,60), 300, 1450, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV4',xy=(1590, 45),  xycoords='data', fontsize=3, color='r')
    fig.gca().add_artist(text)
    
    plt.colorbar(fraction=0.15,aspect=15,orientation='horizontal',shrink=0.7,pad=0.05)
    plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
    widthfig = 16 * cm
    heightfig= 22 * cm
    pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
    plt.close('all')    
    Story.append(pi1)

    Story.append(PageBreak())

    #==============================================================================
    # Fringe tracker photometry
    #==============================================================================

    plotTitle(Story, "Fringe Tracker mean photometric flux per output (ADU/s)")
    
    hsize = 19*cm
    vsize = 6*cm
    # Grand average transmission
    # raw.data_ft[time,output,wave_channel]
    # trans_avg_ft = np.mean(np.mean(raw.data_ft,axis=0).T,axis=1)
    trans_ft = np.mean(np.median(raw.data_ft,axis=0),axis=1)/raw.exptime_ft
    yminval=np.min([0,np.min(trans_ft)-0.1*np.abs(np.min(trans_ft))])
    ymaxval=np.max(trans_ft)+0.1*np.abs(np.min(trans_ft))
    ystep=(ymaxval-yminval)/10.
    if raw.polarsplit == True:
        transdata = [tuple(clipdata(trans_ft[::2],yminval,ymaxval)),tuple(clipdata(trans_ft[1::2],yminval,ymaxval))]
        #transdata = [tuple(clipdata(trans_avg_ft[:,0],yminval,ymaxval)),tuple(clipdata(trans_avg_ft[:,1],yminval,ymaxval))]
        labelbase = raw.regname_ft[::2].astype('|S4').tolist()
        bc = transbarchartout(transdata,labelbase,yminval,ymaxval,ystep,hsize,vsize,colors.aquamarine,colors.cornflower)
    else:
        transdata = [tuple(clipdata(trans_ft,yminval,ymaxval))]
        labelbase = raw.regname_ft.astype('|S4').tolist()
        bc = transbarchartout(transdata,labelbase,yminval,ymaxval,ystep,hsize,vsize,colors.aquamarine,colors.cornflower)
    Story.append(bc)
        
    #==============================================================================
    # FT median photometry vs. time
    #==============================================================================

    plotTitle(Story,"FT mean flux per pixel vs. time (ADU/s)")
    plotSubtitle(Story,"In ADU/s vs. time (seconds). The background is not subtracted.")
    Story.append(Spacer(1,2*mm))
    
    hsize = 16*cm
    vsize = 6*cm

    meanflux_ft = []
    for i in range(0,raw.nframe_ft):    
        meanflux_ft.append(np.median(raw.data_ft[i,:,:])/raw.exptime_ft)
    
    #resamp = 50
    resamp = int(np.max([int(len(raw.time_ft)/pixwidth),1]))

    yminval = np.min(meanflux_ft[::resamp])
    ymaxval = np.max(meanflux_ft[::resamp])
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append(tuple(zip(raw.time_ft[::resamp], clipdata(meanflux_ft[::resamp],yminval,ymaxval))))
    a = graphoutaxes(data, np.min(raw.time_ft[::resamp]), np.max(raw.time_ft[::resamp]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Mean FT flux ; Avg = %(val).3e' % {"val":np.mean(meanflux_ft)}+' ADU/s')
    Story.append(a)
    Story.append(Spacer(1,3*mm))

    #==============================================================================
    # Power spectral density of the average FT flux vs. time
    #==============================================================================

    if welchpresent == True:
        style = styles["Heading2"]
        p = Paragraph("Fringe Tracker mean flux per pixel periodogram", style)
        Story.append(p)
        style = styles["Normal"]
        p = Paragraph("In ADU RMS vs. Hz.", style)
        Story.append(p)
        Story.append(Spacer(1,2*mm))
    
        PSD=[1,2,3,4]    
        
        avgflux_ft = np.mean(np.mean(raw.data_ft,axis=1),axis=1)
        
        f, PSD = welch(avgflux_ft[:], fs=ftfreq, detrend='linear', nperseg=1024, scaling='spectrum')
        yminval = 0
        ymaxval = np.sqrt(np.max(PSD))
        hsize = 16*cm
        vsize = 6*cm
        ticky = (ymaxval-yminval)/10
        
        data = []
        data.append(tuple(zip(f, clipdata(np.sqrt(PSD), yminval,ymaxval))))
        a = graphoutaxes( data, min(f),max(f),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Average of all FT pixels')
        Story.append(a)
        Story.append(Spacer(1,5*mm))

    Story.append(PageBreak())

    if raw.gdelay_in == True:
        #==============================================================================
        # OPDC signals
        #==============================================================================
        
        style = styles["Heading2"]
        p = Paragraph("OPDC signals (OPDC)", style)
        Story.append(p)
    
        Story.append(Spacer(1,5*mm))
    
        resamp = 50 # take one value every resamp value for the plot
    
        hsize = 16*cm
        vsize = 5.5*cm
    
        piezo_time_plot = raw.opdc_time[::resamp]
        state_plot = raw.opdc_state[::resamp]
        
        style = styles["Normal"]
        p = Paragraph("OPDC state", style)
        Story.append(p)
        
        Story.append(Spacer(1,3*mm))
    
        data = []
        yminval = 0 
        ymaxval = 11
        ticky = 1
        data.append( tuple(zip(piezo_time_plot, clipdata(state_plot,yminval,ymaxval)) ))
        a = graphoutaxes( data, min(piezo_time_plot),max(piezo_time_plot),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
        a.add(String(0.05*hsize, 0.10*vsize, '1=IDLE', fontSize=9))
        a.add(String(0.05*hsize, 0.19*vsize, '2=GD tracking', fontSize=9))
        a.add(String(0.05*hsize, 0.28*vsize, '3=Phase tracking', fontSize=9))
        a.add(String(0.05*hsize, 0.37*vsize, '4=Search', fontSize=9))
        a.add(String(0.05*hsize, 0.46*vsize, '5=Calib VLTI DL (open)', fontSize=9))
        a.add(String(0.05*hsize, 0.55*vsize, '6=Calib Piezo (open)', fontSize=9))
        a.add(String(0.05*hsize, 0.92*vsize, '10=No flux', fontSize=9))
        Story.append(a)
        Story.append(Spacer(1,10*mm))
    
        style = styles["Normal"]
        p = Paragraph("Piezo offsets (Volts). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.", style)
        Story.append(p)
        
        Story.append(Spacer(1,3*mm))
    
        data = []
        yminval = np.nanmin(raw.opdc_piezo_offset)-1  
        ymaxval = np.nanmax(raw.opdc_piezo_offset)+1 
        ticky = (ymaxval-yminval)/10. # in microns
        for tel in range(0,ntel):
            data.append( tuple(zip(piezo_time_plot, clipdata(raw.opdc_piezo_offset[::resamp,tel],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(piezo_time_plot),max(piezo_time_plot),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
        Story.append(a)
        Story.append(Spacer(1,10*mm))
    
        style = styles["Normal"]
        p = Paragraph("FDDL offsets (microns). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.", style)
        Story.append(p)
        
        Story.append(Spacer(1,3*mm))
    
        data = []
        yminval = np.nanmin(raw.opdc_vltidl_offset)-1  
        ymaxval = np.nanmax(raw.opdc_vltidl_offset)+1 
        ticky = (ymaxval-yminval)/10. # in microns
        for tel in range(0,ntel):
            data.append( tuple(zip(piezo_time_plot, clipdata(raw.opdc_vltidl_offset[::resamp,tel],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(piezo_time_plot),max(piezo_time_plot),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
        Story.append(a)
    
        Story.append(PageBreak())
    

    #==============================================================================
    # Science combiner average vertical cut
    #==============================================================================

    style = styles["Heading2"]
    p = Paragraph("SC average frame vertical cut vs. line number (ADU/s)", style)
    Story.append(p)
    Story.append(Spacer(1,3*mm))
    
    if raw.header['MJD-OBS'] < 57241:
        if 'LOW' in raw.header['HIERARCH ESO INS SPEC RES']:
            minx = 950
            maxx = 1200 # TO BE CHECKED 
        if 'MEDIUM' in raw.header['HIERARCH ESO INS SPEC RES']:
            minx = 950
            maxx = 1200
        if 'HIGH' in raw.header['HIERARCH ESO INS SPEC RES']:
            minx = 0
            maxx = 2047 # TO BE CHECKED
    else:
        minx = 0
        maxx = raw.data_sc.shape[1]
        
    hsize = 16*cm
    vsize = 6*cm
    # Grand average transmission
    cut_avg_sc = np.mean(np.mean(raw.data_sc,axis=0)[:,minx:maxx],axis=1)/raw.exptime_sc

    pixscale = list(range(0,cut_avg_sc.shape[0]))
    yminval = np.min(cut_avg_sc)
    ymaxval = np.max(cut_avg_sc)
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append( tuple(zip(pixscale, clipdata(cut_avg_sc,yminval,ymaxval))))
    a = graphoutaxes(data, np.min(pixscale), np.max(pixscale),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
    Story.append(a)

    Story.append(Spacer(1,7*mm))

    #==============================================================================
    # Science combiner spectrum (95% percentile over all outputs)
    #==============================================================================

    style = styles["Heading2"]
    p = Paragraph("SC average spectrum (ADU/s)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In ADU/s vs. pixel number along dispersion (95% percentile of all outputs).", style)
    Story.append(p)
    Story.append(Spacer(1,3*mm))
    
    hsize = 16*cm
    vsize = 6*cm
    # Grand average spectrum. 95% percentile allows to exclude the hot badpix
    avg_sc_spectrum = np.percentile(np.mean(raw.data_sc,axis=0),95,axis=0)/raw.exptime_sc
    
    pixscale = list(range(0,avg_sc_spectrum.shape[0]))
    yminval = 0 #np.min(avg_sc_spectrum)
    ymaxval = np.max(avg_sc_spectrum)*1.3
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append( tuple(zip(pixscale, clipdata(avg_sc_spectrum,yminval,ymaxval))))
    a = graphoutaxes(data, np.min(pixscale), np.max(pixscale),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
    Story.append(a)

    Story.append(Spacer(1,7*mm))

    #==============================================================================
    # Science combiner total flux and background vs. time
    #==============================================================================

    style = styles["Heading2"]
    p = Paragraph("SC total flux vs. frame number (ADU/s)", style)
    Story.append(p)
    Story.append(Spacer(1,3*mm))
    
    # background extimation window located right of the spectra
    if raw.header['MJD-OBS'] < 57240: # for older files without windowing
        minx = 0
        maxx = 10
        if 'LOW' in raw.header['HIERARCH ESO INS SPEC RES']:
            minxbg = 5
            maxxbg = 100 # TO BE CHECKED 
        if 'MEDIUM' in raw.header['HIERARCH ESO INS SPEC RES']:
            minxbg = 5
            maxxbg = 100
        if 'HIGH' in raw.header['HIERARCH ESO INS SPEC RES']:
            minxbg = 5
            maxxbg = 50 # TO BE CHECKED
    else:
        minx = 1
        maxx = 10
        if 'LOW' in raw.header['HIERARCH ESO INS SPEC RES']:
            minxbg = raw.data_sc.shape[2]-10
            maxxbg = raw.data_sc.shape[2]-1 # TO BE CHECKED 
        if 'MEDIUM' in raw.header['HIERARCH ESO INS SPEC RES']:
            minxbg = raw.data_sc.shape[2]-10
            maxxbg = raw.data_sc.shape[2]-1
        if 'HIGH' in raw.header['HIERARCH ESO INS SPEC RES']:
            minxbg = raw.data_sc.shape[2]-10
            maxxbg = raw.data_sc.shape[2]-1 # TO BE CHECKED 
    
    hsize = 16*cm
    vsize = 6*cm
    flux_sc = np.zeros(raw.data_sc.shape[0])
    bglevel = np.zeros(raw.data_sc.shape[0])
    bgrms = np.zeros(raw.data_sc.shape[0])
    maxx = raw.data_sc.shape[2]
    
    images_sc = np.copy(raw.data_sc)

    # correction of outlier pixels    
    max_ok = np.percentile(raw.data_sc,99.5)
    min_ok = np.percentile(raw.data_sc,1)
    images_sc[np.where(images_sc>max_ok)]=0    
    images_sc[np.where(images_sc<min_ok)]=0    

    avgbgimage = np.mean(images_sc[:,:,minxbg:maxxbg],axis=0)
    
    # Grand total flux (bg corrected)
    for i in range(0,raw.nframe_sc):
        bglevel[i] = np.median(images_sc[i,:,minxbg:maxxbg])/raw.exptime_sc
        flux_sc[i] = np.sum(images_sc[i,:,:minxbg]/raw.exptime_sc-bglevel[i])
        bgrms[i]   = np.std(images_sc[i,:,minxbg:maxxbg]-avgbgimage)
    
    yminval = 0
    ymaxval = np.max(flux_sc)*1.5
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append(tuple(zip(list(range(0,raw.nframe_sc)), clipdata(flux_sc,yminval,ymaxval))))
    a = graphoutaxes(data, 0, raw.nframe_sc,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Total SC flux ; Mean = %(val).3e' % {"val":np.mean(flux_sc)}+' ADU/s')
    Story.append(a)

    Story.append(PageBreak())

    style = styles["Heading2"]
    p = Paragraph("SC background vs. frame number (ADU/s)", style)
    Story.append(p)
    Story.append(Spacer(1,3*mm))

    hsize = 16*cm
    vsize = 6*cm
    
    yminval = np.min([np.min(bglevel),0])
    ymaxval = np.max(bglevel)*1.5
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append(tuple(zip(list(range(0,raw.nframe_sc)), clipdata(bglevel,yminval,ymaxval))))
    a = graphoutaxes(data, 0, raw.nframe_sc,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Median SC background level ; Mean = %(val).0f' % {"val":np.mean(bglevel)}+' ADU/s')
    Story.append(a)

    Story.append(Spacer(1,7*mm))

    style = styles["Heading2"]
    p = Paragraph("SC background RMS vs. time (ADU)", style)
    Story.append(p)
    Story.append(Spacer(1,3*mm))
    
    yminval = 0
    ymaxval = np.max(bgrms)*1.5
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append(tuple(zip(list(range(0,raw.nframe_sc)), clipdata(bgrms,yminval,ymaxval))))
    a = graphoutaxes(data, 0, raw.nframe_sc,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'SC background time RMS ; Median = %(val).0f' % {"val":np.median(bgrms)}+' ADU')
    Story.append(a)

    Story.append(Spacer(1,7*mm))

    style = styles["Heading2"]
    p = Paragraph("SC detector saturation check (ADU)", style)
    Story.append(p)
    Story.append(Spacer(1,3*mm))
    
    hsize = 16*cm
    vsize = 6*cm
    # Maximum spectrum. 99% percentile allows to exclude the hot badpix
    max_sc_spectrum = np.percentile(np.percentile(raw.data_sc,99,axis=0),99,axis=0)
    
    pixscale = list(range(0,max_sc_spectrum.shape[0]))
    yminval = 0 #np.min(avg_sc_spectrum)
    ymaxval = np.max(max_sc_spectrum)*1.3
    ticky = (ymaxval-yminval)/5.

    data = []
    data.append( tuple(zip(pixscale, clipdata(max_sc_spectrum,yminval,ymaxval))))
    data.append( tuple(zip([0,len(max_sc_spectrum)],[satlevel_sc,satlevel_sc])))
    a = graphoutaxes(data, np.min(pixscale), np.max(pixscale),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'2.8e4 ADU Maximum tolerated')
    Story.append(a)

    Story.append(PageBreak())

    #==============================================================================
    # Image display of the SC average frame
    #==============================================================================

    style = styles["Heading2"]
    p = Paragraph("Science Combiner average detector frame (ADU/s/pix)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("Color scale in ADU/s.", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))
    
    minval = np.percentile(raw.data_sc,0.1)/raw.exptime_sc
    maxval = np.percentile(raw.data_sc,99.9)/raw.exptime_sc

    imgdata = io.StringIO()
    plt.close('all')
    plt.figure(figsize=(3,6),dpi=300)
    plt.rc('xtick', labelsize=5)
    plt.rc('ytick', labelsize=5)
    #figaspect=1
    figaspect = (0.9*raw.data_sc.shape[2]/raw.data_sc.shape[1])*2.0
    plt.imshow(np.mean(raw.data_sc,axis=0)/raw.exptime_sc,vmin=minval,vmax=maxval,cmap='cubehelix_r', interpolation='nearest', aspect=figaspect)
    plt.colorbar(fraction=0.03,aspect=15,orientation='vertical',shrink=0.7,pad=0.05)
    plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
    widthfig = 15 * cm
    heightfig= 23 * cm
    pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
    plt.close('all')    
    Story.append(pi1)

    Story.append(PageBreak())

    #==============================================================================
    # Image display of the FT average frame
    #==============================================================================
    
    style = styles["Heading2"]
    p = Paragraph("Fringe Tracker average detector frame (ADU/s/pix)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("Color scale in ADU/s.", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))
    
    minval = np.percentile(np.mean(raw.data_ft,axis=0),0.1)/raw.exptime_ft
    maxval = np.percentile(np.mean(raw.data_ft,axis=0),99.9)/raw.exptime_ft

    imgdata = io.StringIO()
    plt.close('all')
    plt.figure(figsize=(3,6),dpi=250)
    plt.rc('xtick', labelsize=5)
    plt.rc('ytick', labelsize=5)
    #figaspect=1
    figaspect = (0.9*raw.data_ft.shape[2]/raw.data_ft.shape[1])*2.0
    plt.imshow(np.mean(raw.data_ft,axis=0)/raw.exptime_ft,vmin=minval,vmax=maxval,cmap='cubehelix', interpolation='nearest', aspect=figaspect)
    plt.colorbar(fraction=0.03,aspect=15,orientation='vertical',shrink=0.7,pad=0.05)
    plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
    widthfig = 15 * cm
    heightfig= 23 * cm
    pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
    plt.close('all')    
    Story.append(pi1)    
    
    Story.append(PageBreak())

    #==============================================================================
    # Metrology plots
    #==============================================================================

    resamp = int(np.max([int(len(raw.metrology_volt)/pixwidth),1]))*2

    # The last 16 signals are from the fiber couplers
    labeltable = ['FT GV1', 'FT GV2', 'FT GV3', 'FT GV4',
                  'SC GV1', 'SC GV2', 'SC GV3', 'SC GV4']

    #==============================================================================
    # Fiber coupler and telescope metrology diode signals in circles
    #==============================================================================

    ymaxval = +np.percentile(np.abs(raw.metrology_volt[::resamp,0:64]),99.5)
    yminval = -ymaxval
    ticky = (ymaxval-yminval)/5.
    yminvalfc = -np.max(np.abs(raw.metrology_volt[::resamp,64::]))
    ymaxvalfc = +np.max(np.abs(raw.metrology_volt[::resamp,64::]))
    tickyfc = (ymaxvalfc-yminvalfc)/5.
    
    for tel in range(0,4):
        plotTitle(Story,"Telescope and fiber coupler diodes GV"+str(tel+1)+" = "+raw.telnames[tel]+" (METROLOGY)")
        plotSubtitle(Story,"In volts vs. volts. X=cosinus, Y=sinus. FT = Red, SC = Orange.")
    
        hsize = 6*cm
        vsize = 6*cm
    
        j=1 # index for label text
        d=[]
        # Telescope diodes first: order: 32 FT telescope diodes: T1D1X, T1D1Y, T1D2X, ...,T4D4Y 
        for i in range(0+tel*4,4+tel*4):
            data = []
            data.append( tuple(zip(clipdata(raw.metrology_volt[::resamp,2*i+1],yminval,ymaxval),\
                clipdata(raw.metrology_volt[::resamp,2*i],yminval,ymaxval))))
            data.append( tuple(zip(clipdata(raw.metrology_volt[::resamp,2*i+33],yminval,ymaxval),\
                clipdata(raw.metrology_volt[::resamp,2*i+32],yminval,ymaxval))))
            a = graphscatteraxes(data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,ticky,ticky,2,\
                ' Tel. '+raw.telnames[tel]+' Diode '+str(j))
            d.append(a)
            j=j+1
    
        # FC diodes
        data = []
        data.append( tuple(zip(clipdata(raw.metrology_volt[::resamp,65+tel*2],yminvalfc,ymaxvalfc),\
            clipdata(raw.metrology_volt[::resamp,64+tel*2],yminvalfc,ymaxvalfc))))
        data.append( tuple(zip(clipdata(raw.metrology_volt[::resamp,73+tel*2],yminvalfc,ymaxvalfc),\
            clipdata(raw.metrology_volt[::resamp,72+tel*2],yminvalfc,ymaxvalfc))))
        a = graphscatteraxes(data,yminvalfc,ymaxvalfc,None,yminvalfc,ymaxvalfc,None,hsize,vsize,tickyfc,tickyfc,2,\
            ' Beam GV'+str(tel+1)+' Fiber Coupler Diode')
        d.append(a)
    
        plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4]] ]
        t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm)
        t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')]))
        t.hAlign = 'LEFT'
        Story.append(t)
          
        Story.append(PageBreak())
        
    #==============================================================================
    # Fiber coupler unwrapped metrology OPL 
    #==============================================================================
            
    style = styles["Heading2"]
    p = Paragraph("Unwrapped FT OPL from Fiber Coupler (METROLOGY)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In microns vs. time (seconds).", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))

    hsize = 16*cm
    vsize = 4.7*cm
    j=0 # index for labeltable
    for i in range(64,72,2):
        data = []
        phase = np.arctan2(raw.metrology_volt[:,i], raw.metrology_volt[:,i+1])
        unwrapped_opl = np.unwrap(phase)*Lambda_met/(2*np.pi)
        yminval = np.min(unwrapped_opl)
        ymaxval = np.max(unwrapped_opl)
        ticky = (ymaxval-yminval)/5.
        data.append( tuple(zip(raw.metrology_time[::resamp], clipdata(unwrapped_opl[::resamp],yminval,ymaxval))))
        a = graphoutaxes(data, min(raw.metrology_time),max(raw.metrology_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Unwrapped '+labeltable[j])
        Story.append(a)
        Story.append(Spacer(1,7*mm))
        j=j+1

    Story.append(PageBreak())

    style = styles["Heading2"]
    p = Paragraph("Unwrapped SC OPL from Fiber Coupler (METROLOGY)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In microns vs. time (seconds).", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))

    for i in range(72,80,2):
        data = []
        phase = np.arctan2(raw.metrology_volt[:,i], raw.metrology_volt[:,i+1])
        unwrapped_opl = np.unwrap(phase)*Lambda_met/(2*np.pi)
        yminval = np.min(unwrapped_opl)
        ymaxval = np.max(unwrapped_opl)
        ticky = (ymaxval-yminval)/5.
        data.append( tuple(zip(raw.metrology_time[::resamp], clipdata(unwrapped_opl[::resamp],yminval,ymaxval))))
        a = graphoutaxes(data, min(raw.metrology_time),max(raw.metrology_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Unwrapped '+labeltable[j])
        Story.append(a)
        Story.append(Spacer(1,7*mm))
        j=j+1

    Story.append(PageBreak())

    #==============================================================================
    # Telescope diodes unwrapped metrology OPL 
    #==============================================================================
            
    for tel in range(0,4):
        style = styles["Heading2"]
        p = Paragraph("Unwrapped FT OPL from Telescope %s diodes (METROLOGY)" % str(tel+1), style)
        Story.append(p)
        style = styles["Normal"]
        p = Paragraph("In microns vs. time (seconds).", style)
        Story.append(p)
        Story.append(Spacer(1,2*mm))
        
        for i in range(8*tel,8*tel+8,2):
            data = []
            phase = np.arctan2(raw.metrology_volt[:,i], raw.metrology_volt[:,i+1])
            #phase = clean_unwrap_phase(phase,4,2*np.pi) # cleanup of the telescope diode phases from "soft jumps"
            unwrapped_opl = np.unwrap(phase)*Lambda_met/(2*np.pi)
            yminval = np.min(unwrapped_opl)
            ymaxval = np.max(unwrapped_opl)
            ticky = (ymaxval-yminval)/5.
            data.append( tuple(zip(raw.metrology_time[::resamp], clipdata(unwrapped_opl[::resamp],yminval,ymaxval))))
            a = graphoutaxes(data, min(raw.metrology_time),max(raw.metrology_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
            Story.append(a)
            Story.append(Spacer(1,7*mm))
            
        Story.append(PageBreak())
        
        style = styles["Heading2"]
        p = Paragraph("Unwrapped SC OPL from Telescope %s diodes (METROLOGY)" % str(tel+1), style)
        Story.append(p)
        style = styles["Normal"]
        p = Paragraph("In microns vs. time (seconds).", style)
        Story.append(p)
        Story.append(Spacer(1,2*mm))

        for i in range(8*tel+32,8*tel+40,2):
            data = []
            phase = np.arctan2(raw.metrology_volt[:,i], raw.metrology_volt[:,i+1])
            #phase = clean_unwrap_phase(phase,4,2*np.pi) # cleanup of the telescope diode phases from "soft jumps"
            unwrapped_opl = np.unwrap(phase)*Lambda_met/(2*np.pi)
            yminval = np.min(unwrapped_opl)
            ymaxval = np.max(unwrapped_opl)
            ticky = (ymaxval-yminval)/5.
            data.append( tuple(zip(raw.metrology_time[::resamp], clipdata(unwrapped_opl[::resamp],yminval,ymaxval))))
            a = graphoutaxes(data, min(raw.metrology_time),max(raw.metrology_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
            Story.append(a)
            Story.append(Spacer(1,7*mm))
            
        Story.append(PageBreak())
        
    #==============================================================================
    # Fiber coupler metrology diode signals 
    #==============================================================================

    style = styles["Heading2"]
    p = Paragraph("Fiber coupler diodes (METROLOGY) - Full range", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In volts vs. time (seconds). Colors: sinus=red, cosinus=orange.", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))

    resamp = int(np.max([int(len(raw.metrology_time)/pixwidth),1]))

    ncurves = 2 # Number of curves per subplot
    # take one value every resamp value for the plot, to have pixwidth points on page width
    yminval = np.min(raw.metrology_volt)
    ymaxval = np.max(raw.metrology_volt)
    ticky = (ymaxval-yminval)/5.
    hsize = 16*cm
    vsize = 2.7*cm

    j=0 # index for labeltable
    for i in range(int(64/ncurves),int(80/ncurves)):
        data = []
        for channel in range(i*ncurves,(i+1)*ncurves):
            data.append( tuple(zip(raw.metrology_time[::resamp], clipdata(raw.metrology_volt[::resamp,channel],yminval,ymaxval))))
        if ((i>0)*((i+1)%8==0) or (i==int(80/ncurves))):
            a = graphoutaxes(data, min(raw.metrology_time),max(raw.metrology_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Fiber coupler '+labeltable[j])
            Story.append(a)
            Story.append(PageBreak())
        else:
            a = graphoutnoaxis(data, min(raw.metrology_time),max(raw.metrology_time),yminval,ymaxval,hsize,vsize,ticky,'Fiber coupler '+labeltable[j])
            Story.append(a)
        j=j+1
        

    maxtime = 10 # maximum time for display in seconds
    style = styles["Heading2"]
    p = Paragraph("Fiber coupler diodes (METROLOGY) - First "+str(maxtime)+" seconds", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In volts vs. time (seconds). Colors: sinus=red, cosinus=orange.", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))

    # Enlargement of the first 10 seconds without resampling
    maxindex = np.max(np.where(raw.metrology_time<maxtime))
    j=0 # index for labeltable
    for i in range(int(64/ncurves),int(80/ncurves)):
        data = []
        for channel in range(i*ncurves,(i+1)*ncurves):
            data.append( tuple(zip(raw.metrology_time[0:maxindex], clipdata(raw.metrology_volt[0:maxindex,channel],yminval,ymaxval))))
        if ((i>0)*((i+1)%8==0) or (i==int(80/ncurves))):
            a = graphoutaxes(data, min(raw.metrology_time),maxtime,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Fiber coupler '+labeltable[j])
            Story.append(a)
            Story.append(PageBreak())
        else:
            a = graphoutnoaxis(data, min(raw.metrology_time),maxtime,yminval,ymaxval,hsize,vsize,ticky,'Fiber coupler '+labeltable[j])
            Story.append(a)
        j=j+1

    #==============================================================================
    # FDDL signals
    #==============================================================================

    plotTitle(Story,"SC FDDL voltage (FDDL)")
    plotSubtitle(Story,"In volts vs. time (seconds).")
    Story.append(Spacer(1,2*mm))
    
    resamp = int(np.max([int(raw.fddl_sc_pos.shape[0]/pixwidth),1]))

    for tel in range(0,4):
        yminval = np.min(raw.fddl_sc_pos[:,tel])-0.1*(np.max(raw.fddl_sc_pos[:,tel])-np.min(raw.fddl_sc_pos[:,tel]))
        ymaxval = np.max(raw.fddl_sc_pos[:,tel])+0.1*(np.max(raw.fddl_sc_pos[:,tel])-np.min(raw.fddl_sc_pos[:,tel]))
        ticky = (ymaxval-yminval)/10. # in volts
        hsize = 16*cm
        vsize = 4.5*cm
        data = []
        data.append( tuple(zip(raw.fddl_time[::resamp], clipdata(raw.fddl_sc_pos[::resamp,tel],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(raw.fddl_time),max(raw.fddl_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'SC FDDL Telescope '+str(tel+1))
        Story.append(a)        
        Story.append(Spacer(1,7*mm))

    Story.append(PageBreak())

    style = styles["Heading2"]
    p = Paragraph("FT FDDL voltage (FDDL)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In volts vs. time (seconds).", style)
    Story.append(p)
    
    Story.append(Spacer(1,2*mm))
    for tel in range(0,4):
        yminval = np.min(raw.fddl_ft_pos[:,tel])-0.1*(np.max(raw.fddl_ft_pos[:,tel])-np.min(raw.fddl_ft_pos[:,tel]))
        ymaxval = np.max(raw.fddl_ft_pos[:,tel])+0.1*(np.max(raw.fddl_ft_pos[:,tel])-np.min(raw.fddl_ft_pos[:,tel]))
        ticky = (ymaxval-yminval)/10. # in volts
        hsize = 16*cm
        vsize = 4.5*cm
        data = []
        data.append( tuple(zip(raw.fddl_time[::resamp], clipdata(raw.fddl_ft_pos[::resamp,tel],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(raw.fddl_time),max(raw.fddl_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'FT FDDL Telescope '+str(tel+1))
        Story.append(a)
        Story.append(Spacer(1,7*mm))

    Story.append(PageBreak())

    style = styles["Heading2"]
    p = Paragraph("Metrology laser power (METROLOGY)", style)
    Story.append(p)
    style = styles["Normal"]
    p = Paragraph("In millivolts vs. time (seconds).", style)
    Story.append(p)
    Story.append(Spacer(1,2*mm))

    resamp = int(np.max([int(raw.metrology_powerlaser.shape[0]/pixwidth),1]))
    # take one value every resamp value for the plot, to have pixwidth points on page width
    yminval = np.min(raw.metrology_powerlaser)
    ymaxval = np.max(raw.metrology_powerlaser)
    ticky = (ymaxval-yminval)/10.
    hsize = 16*cm
    vsize = 10*cm
    #self.oi_vis_met_time = oi_vis_met.field('TIME')
    #self.oi_vis_met_phase = oi_vis_met.field('PHI')
    data = []
    data.append( tuple(zip(raw.metrology_time[::resamp], clipdata(raw.metrology_powerlaser[::resamp],yminval,ymaxval))))
    a = graphoutaxes(data, min(raw.metrology_time),max(raw.metrology_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
    Story.append(a)
    Story.append(PageBreak())

    #==============================================================================
    # Create PDF report from Story
    #==============================================================================
    print("Create regular PDF")
    
    gravi_visual_class.TITLE    = "GRAVITY RAW Quality Control Report"
    gravi_visual_class.PAGEINFO = "RAW file: "+filename+".fits"
    reportname = filename+"-RAW.pdf"
    
    doc = SimpleDocTemplate(reportname)
    doc.build(Story, onFirstPage=myFirstPage, onLaterPages=myLaterPages)

    print((" "+reportname)) 
	
	#==============================================================================
    # MetroPy
    #==============================================================================
	
    import os
    import subprocess
	# adjust this path variable to where the MetroPy_SVN.sh script is located
    metropy_path = os.path.join(os.sep, 'tera', '5', 'fhenningsen', 'python_tools', 'gravi_visual', 'MetroPy')
	
    metropy = input("\nCreate MetroPy checkplots? [y/n]: ")
    while True:
	
		if metropy == "y":
			print("\nCreate MetroPy PDF")
			subprocess.call([metropy_path + os.sep + 'MetroPy_SVN.sh', str(filename) + '.fits'])
			break
			
		elif metropy == "n":
			print("\nNo MetroPy PDF will be created")
			break
			
		else:
			metropy = input("\tInvalid input. Please state 'y' or 'n': ")
			
		continue
		
#==============================================================================
# MAIN PROGRAM    
#==============================================================================

if __name__ == '__main__':
    filename=''

    if len(sys.argv) == 2 :
        filename=sys.argv[1]
        
    if filename == '' :
        filename=input(" Enter RAW file name (without .fits) : ")

    # produce raw report
    raw = gravi_visual_class.Rawdata(filename+'.fits')
    produce_raw_report(raw,filename)