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

@author: kervella
"""

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 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.graphics import renderPDF
    from reportlab.lib.units import inch, cm, mm
    from reportlab.lib import colors
    
    import io
    from io import StringIO, BytesIO
    from svglib.svglib import svg2rlg
    
    #==============================================================================
    # 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))

    plt.close('all')
    fig=plt.figure(figsize=(6,13))
    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((20,60), 200, 900, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV1',xy=(50, 45),  xycoords='data', fontsize=6, color='r')
    fig.gca().add_artist(text)

    rect = plt.Rectangle((270,60), 200, 900, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV2',xy=(300, 45),  xycoords='data', fontsize=6, color='r')
    fig.gca().add_artist(text)

    rect = plt.Rectangle((520,60), 200, 900, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV3',xy=(550, 45),  xycoords='data', fontsize=6, color='r')
    fig.gca().add_artist(text)

    rect = plt.Rectangle((770,60), 200, 900, color='r', linewidth=0.5, fill=False)
    fig.gca().add_artist(rect)
    text = plt.annotate('Beam GV4',xy=(800, 45),  xycoords='data', fontsize=6, color='r')
    fig.gca().add_artist(text)
    
    plt.colorbar(fraction=0.15,aspect=15,orientation='horizontal',shrink=0.7,pad=0.05)

    imgdata = BytesIO()
    fig.savefig(imgdata, format='svg', bbox_inches='tight')
    imgdata.seek(0)  # rewind the data
    drawing = svg2rlg(imgdata)
    Story.append(drawing)

    plt.close('all')    

    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

    plt.close('all')

    fig=plt.figure(figsize=(5,10))
    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)

    imgdata = BytesIO()
    fig.savefig(imgdata, format='svg', bbox_inches='tight')
    imgdata.seek(0)  # rewind the data
    drawing = svg2rlg(imgdata)
    Story.append(drawing)

    plt.close('all')    

    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')

    fig=plt.figure(figsize=(5,10))
    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)

    imgdata = BytesIO()
    fig.savefig(imgdata, format='svg', bbox_inches='tight')
    imgdata.seek(0)  # rewind the data
    drawing = svg2rlg(imgdata)
    Story.append(drawing)

    plt.close('all')    
    
    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())

    #==============================================================================
    # SNR in Metrology signals 
    #==============================================================================

    plotTitle(Story,"SNR telescope diodes (METROLOGY) Finge Tracker")
    plotSubtitle(Story,"In time in seconds vs. volts.")

    def averaging(x, N, median=False):
        if N == 1:
            return x
        if x.ndim == 2:
            res = np.zeros((x.shape[0], x.shape[1]//N))
            for idx in range(x.shape[0]):
                if median:
                    res[idx] = np.nanmedian(x[idx].reshape(-1, N), axis=1)
                else:
                    res[idx] = np.nanmean(x[idx].reshape(-1, N), axis=1)
            return res
        elif x.ndim == 1:
            ll = x.shape[0]
            xx = np.pad(x, (0, N - ll % N), constant_values=np.nan)
            if median:
                res = np.nanmedian(xx.reshape(-1, N), axis=1)
            else:
                res = np.nanmean(xx.reshape(-1, N), axis=1)
            return res

    V = raw.metrology_volt*1e3
    time = raw.metrology_time#*1e-6
    
    smoo = 100
    V = np.array([np.convolve(V[:, i], np.ones(smoo)/smoo, 'same')
                      for i in range(80)]).T
    
    VC = V[:, 1::2] + 1j * V[:, ::2]
    VTEL = VC[:, :-8]
    VTELFT = VTEL[:, :16].reshape(-1, 4, 4)
    VTELSC = VTEL[:, 16:].reshape(-1, 4, 4)

    hsize = 4*cm
    vsize = 4*cm
    yminval = 2
    ymaxval = 500

    d = []
    for tel in range(4):
        data = []
        data.append( tuple(zip(averaging(time, smoo)[::5], averaging(np.abs(VTELFT[:,tel].mean(1)), smoo)[::5])) )
        data.append( tuple(zip(averaging(time, smoo)[::5], averaging(np.abs(VTELSC[:,tel].mean(1)), smoo)[::5])) )
        _title = f'{raw.telnames[tel]}, Average'
        a = graphoutnoaxis(data, min(time),max(time),yminval,ymaxval,hsize,vsize,None,_title,log_y=True)
        d.append(a)
    plotmatrix = [ [d[0], d[1], d[2],d[3]] ]
    t = Table(plotmatrix,colWidths=hsize+2*mm,rowHeights=vsize+2*mm)
    t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')]))
    t.hAlign = 'LEFT'
    # Story.append(t)
    # Story.append(Spacer(1,3*mm))
    
    for dio in range(4):
        d = []
        for tel in range(4):
            data = []
            data.append( tuple(zip(averaging(time, smoo)[::5], averaging(np.abs(VTELFT[:,tel,dio]), smoo)[::5])) )
            # data.append( tuple(zip(averaging(time, smoo)[::5], averaging(np.abs(VTELSC[:,tel,dio]), smoo)[::5])) )
            _title = f'{raw.telnames[tel]}, Diode {dio+1}'
            if dio == 3 and tel == 0:
                a = graphoutaxes(data, min(time),max(time),None,yminval,ymaxval,None,hsize,vsize,None,None,_title,log_y=True)
            else:
                a = graphoutnoaxis(data, min(time),max(time),yminval,ymaxval,hsize,vsize,None,_title,log_y=True)
            d.append(a)
        plotmatrix = [ [d[0], d[1], d[2],d[3]] ]
        t = Table(plotmatrix,colWidths=hsize+2*mm,rowHeights=vsize+2*mm)
        t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')]))
        t.hAlign = 'LEFT'
        Story.append(t)        
        # Story.append(Spacer(1,3*mm))

    Story.append(PageBreak())
        

    plotTitle(Story,"SNR telescope diodes (METROLOGY) Science")
    plotSubtitle(Story,"In time in seconds vs. volts.")
    for dio in range(4):
        d = []
        for tel in range(4):
            data = []
            # data.append( tuple(zip(averaging(time, smoo)[::5], averaging(np.abs(VTELFT[:,tel,dio]), smoo)[::5]*0)) )
            data.append( tuple(zip(averaging(time, smoo)[::5], averaging(np.abs(VTELSC[:,tel,dio]), smoo)[::5])) )
            _title = f'{raw.telnames[tel]}, Diode {dio+1}'
            if dio == 3 and tel == 0:
                a = graphoutaxes(data, min(time),max(time),None,yminval,ymaxval,None,hsize,vsize,None,None,_title,log_y=True)
            else:
                a = graphoutnoaxis(data, min(time),max(time),yminval,ymaxval,hsize,vsize,None,_title,log_y=True)
            d.append(a)
        plotmatrix = [ [d[0], d[1], d[2],d[3]] ]
        t = Table(plotmatrix,colWidths=hsize+2*mm,rowHeights=vsize+2*mm)
        t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')]))
        t.hAlign = 'LEFT'
        Story.append(t)        
        # Story.append(Spacer(1,3*mm))

    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)

    #==============================================================================
    # Create PDF report from Story
    #==============================================================================
    print("Create the 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)) 

#==============================================================================
# MAIN PROGRAM    
#==============================================================================

if __name__ == '__main__':
    filename=''
    filename='../GRAVITY.2015-09-14T21-25-16'

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

    raw = gravi_visual_class.Rawdata(filename+'.fits')

    produce_raw_report(raw,filename)