#! /usr/bin/env python3
# -*- coding: iso-8859-15 -*-
"""
Created on Mon May 25 13:51:41 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
from . import bootstrap as boot
import numpy as np
import pdb
import warnings

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_astroreduced_report(astroreduced,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.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 io
    from io import StringIO, BytesIO
    from svglib.svglib import svg2rlg
    
    #==============================================================================
    # Start Story
    #==============================================================================
    Story = [Spacer(1,1*mm)]
    plotSummary (Story, filename, astroreduced, onTarget=True)
    
    if 'HIERARCH ESO INS MLC WAVELENG' in astroreduced.header:
        Lambda_met=astroreduced.header['HIERARCH ESO INS MLC WAVELENG']/1000. # Metrology wavelength in microns (1.908 microns)
    else: Lambda_met = 1.908287
    
    gravitybasenames = ['12','13','14','23','24','34']
    nframeplot = astroreduced.nframe_sc
    
    #==============================================================================
    # Photometry
    #==============================================================================
    style = styles["Heading2"]
    p = Paragraph("Average photom. flux (photo-e-/SC +/- std)", style)
    Story.append(p)
    
    Story.append(Spacer(1,2*mm))

    if hasattr(astroreduced,'oi_flux_sc') & (astroreduced.polarsplit == False):
        avgfluxsc = []
        avgfluxft = []
        for tel in range(0,4):
            avgfluxsc.append('%(val).2e ' % {"val":np.nanmean(astroreduced.oi_flux_sc[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmean(np.nanstd(astroreduced.oi_flux_sc[:,tel,:],axis=0))})
            # average flux over wavelength
        data = avgfluxsc
        data = np.vstack((['Tel1 '+astroreduced.stations[0],'Tel2 '+astroreduced.stations[1],\
                           'Tel3 '+astroreduced.stations[2],'Tel4 '+astroreduced.stations[3]],data))
        data = np.hstack(([["Telescope"],["SC"]],data))
        t=Table(data.tolist(),colWidths=3.5*cm)
        t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black),
                           ('BACKGROUND',(0,0),(0,0),colors.pink),
                           ('BOX', (0,0), (-1,-1), 2, colors.black),]))
        Story.append(t)
    
    if hasattr(astroreduced, 'oi_flux_sc_s') & (astroreduced.polarsplit == True):
        avgfluxs = []
        avgfluxp = []
        for tel in range(0,4):
            avgfluxs.append('%(val).2e ' % {"val":np.nanmean(astroreduced.oi_flux_sc_s[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmean(np.nanstd(astroreduced.oi_flux_sc_s[:,tel,:],axis=0))})
            avgfluxp.append('%(val).2e ' % {"val":np.nanmean(astroreduced.oi_flux_sc_p[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmean(np.nanstd(astroreduced.oi_flux_sc_p[:,tel,:],axis=0))})
            # average flux over wavelength
        data = avgfluxs
        data = np.vstack((['Tel1','Tel2','Tel3','Tel4'],data))
        data = np.vstack((data,avgfluxp))
        data = np.hstack(([["SC"],["Polar. S"],["Polar. P"]],data))
        t=Table(data.tolist(),colWidths=3.5*cm)
        t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black),
                           ('BACKGROUND',(0,0),(0,0),colors.pink),
                           ('BOX', (0,0), (-1,-1), 2, colors.black),]))
        Story.append(t)

    #==============================================================================
    # Grand average photometric spectrum of SC
    #==============================================================================

    plotTitle(Story,"SC Average spectrum of all outputs (OI_FLUX)")
    plotSubtitle(Story,"In photo-e- vs. wavelength (microns).")

    if hasattr(astroreduced,'oi_flux_sc') & (astroreduced.polarsplit == False):
        nmeasure = astroreduced.time_sc.shape[0]
        avgflux = np.zeros((astroreduced.nwave_sc),dtype='d')
        avgflux = np.nanmean(np.nanmean(astroreduced.oi_flux_sc,axis=0), axis=0) # average flux over time, then over telescopes
    
        yminval = 0 # np.nanpercentile(avgflux,1)-0.1*np.abs(np.nanpercentile(avgflux,99))
        ymaxval = 1.3*np.nanpercentile(avgflux,99)
        ticky = (ymaxval-yminval)/10 # in photo-e-
        hsize = 16*cm
        vsize = 4*cm
        data = [ tuple(zip(astroreduced.wave_sc, clipdata(avgflux[:],yminval,ymaxval)) ) ]
        a = graphoutaxes( data, min(astroreduced.wave_sc),max(astroreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined polarization')
        Story.append(a)
        
    if hasattr(astroreduced,'oi_flux_sc_s') & (astroreduced.polarsplit == True):
        nmeasure = astroreduced.time_sc.shape[0]
        avgfluxs = np.zeros((astroreduced.nwave_sc),dtype='d')
        avgfluxp = np.zeros((astroreduced.nwave_sc),dtype='d')
        avgfluxs =np.nanmean(np.nanmean(astroreduced.oi_flux_sc_s,axis=0), axis=0) # average flux over time, then over telescopes
        avgfluxp =np.nanmean(np.nanmean(astroreduced.oi_flux_sc_p,axis=0), axis=0) # average flux over time, then over telescopes
    
        yminval = 0 # np.nanpercentile(avgfluxs,1)-0.1*np.abs(np.nanpercentile(avgfluxs,99))
        ymaxval = 1.3*np.nanpercentile([avgfluxs,avgfluxp],99)
        ticky = (ymaxval-yminval)/5 # in photo-e-
        hsize = 16*cm
        vsize = 4*cm
        avgflux = np.mean([avgfluxs,avgfluxp],axis=0)
        data = [tuple(zip(astroreduced.wave_sc, clipdata(avgflux[:],yminval,ymaxval)) )]
        a = graphoutaxes( data, min(astroreduced.wave_sc),max(astroreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Average of S and P')
 
        Story.append(a)
        Story.append(Spacer(1,7*mm))

    #==============================================================================
    # SC rejection flag
    #==============================================================================

    plotTitle(Story,"SC Rejection flag of frames")
    plotSubtitle(Story,"Flag on any of the six baselines and all polarizations vs. time (seconds).")
    
    if astroreduced.polarsplit == True:
        hsize = 16*cm
        vsize = 4.0*cm
        global_flag = np.max(astroreduced.oi_vis_reject_flag_s,axis=1)
        ymaxval = np.max(global_flag)
    else:
        hsize = 16*cm
        vsize = 4.0*cm
        global_flag = np.max(astroreduced.oi_vis_reject_flag, axis=1)
        ymaxval = np.max(global_flag)
        
    yminval = 0
    ticky = 1
    data = [tuple(zip(astroreduced.time_sc, global_flag) )]
    
    if ymaxval <= 4:        
        a = graphoutaxes( data, min(astroreduced.time_sc),max(astroreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
        a.add(String(0.03*hsize, 0.22*vsize, '0=Good', fontSize=9))
        a.add(String(0.03*hsize, 0.42*vsize, '1=low FT SNR', fontSize=9))
        a.add(String(0.03*hsize, 0.62*vsize, '2=low FT Vfactor', fontSize=9))
        a.add(String(0.03*hsize, 0.82*vsize, '3=low FT Vfactor and low FT SNR', fontSize=9))
    else:
        a = graphoutaxes( data, min(astroreduced.time_sc),max(astroreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'')
    
    Story.append(a)

    Story.append(PageBreak())

        
    #==============================================================================
    # Summary of reduction parameters
    #==============================================================================

    plotReductionSummary(Story, astroreduced)
    Story.append(PageBreak())


    #==============================================================================
    # Pupil position ACQ_CAM OI_VIS_ACQ UVW OPD
    #==============================================================================
    if astroreduced.acq_in is True:
        
        plotTitle(Story,"Pupil UVW shift from ACQ_CAM in 2D")
        plotSubtitle(Story,"UV positions in meters")

        fig,ax = plt.subplots(nrows=2, ncols=2, figsize=(5,5), sharex='all', sharey='all')
        
        # adapted from code by Guillaume Bourdarot
        for k in range(4):
            # Plot the pupil
            pupil_U = astroreduced.oi_acq_pup_u[k::ntel] # in m
            pupil_V = astroreduced.oi_acq_pup_v[k::ntel]
            #pupil_W = astroreduced.oi_acq_pup_w[k::ntel]
            
            # Number of telescopes with valid pupils
            pupil_N = astroreduced.oi_acq_pup_n[k::ntel]
            pupOK = np.where(pupil_N>0)

            cplot = plt.cm.gnuplot(np.linspace(0.1,0.9,len(pupil_U[pupOK])))
            
            ax[k//2,k%2].scatter(pupil_U[pupOK],pupil_V[pupOK],s=6,c=cplot)
            
        # Plot the lateral and longitudinal tolerance
        if 'A' in astroreduced.telnames[0]: #if ATs
            pupil_tol = 0.04 / (8.2/1.8) # m
            xlim = np.max([1.2*np.nanpercentile(np.abs(astroreduced.oi_acq_pup_u),99),
                           1.2*np.nanpercentile(np.abs(astroreduced.oi_acq_pup_v),99),
                           0.05]) # 0.05 # m
            wlim = 5 # m
            opdlim = np.max([1.1E6*np.nanmax(np.abs(astroreduced.oi_acq_opd_pup)),1.00]) # 1.00 # micron
        else: # if UTs
            pupil_tol = 0.04 # m
            xlim = np.max([1.2*np.nanpercentile(np.abs(astroreduced.oi_acq_pup_u),99),
                           1.2*np.nanpercentile(np.abs(astroreduced.oi_acq_pup_v),99),
                           0.05]) # * (8.2/1.8) # m
            wlim = 5 # m
            opdlim = np.max([1.1E6*np.nanmax(np.abs(astroreduced.oi_acq_opd_pup)),1.00]) # 1.00 # micron

        # print(astroreduced.oi_acq_pup_u)                
        # print("MAX  = ",np.nanmax(astroreduced.oi_acq_pup_u))
        # print("XLIM = ",xlim)
        # print(np.abs(astroreduced.oi_acq_opd_pup))
        # print("OPDLIM = ",opdlim)

        ax[0,0].set_xlim(-xlim,xlim)
        ax[0,0].set_ylim(-xlim,xlim)

        for k in range(4):
            # add the requirements on the pupil for astrometry : 4cm for UTs
            circ = plt.Circle((0,0),radius=pupil_tol,fc='None',ec='k',linestyle='--')
            ax[k//2,k%2].add_patch(circ)
            ax[k//2,0].set_ylabel('PUPIL_V [m]',fontsize=8)
            ax[k//2,k%2].set_aspect('equal', 'box')
            ax[1,k%2].set_xlabel('PUPIL_U [m]',fontsize=8)
            ax[k//2,k%2].tick_params(labelsize=8)
            ax[k//2,k%2].scatter([],[],c='none',edgecolor='none',label='%s'%(astroreduced.telnames[k]))
            ax[k//2,k%2].legend(handlelength=0, handletextpad=0, fancybox=True,fontsize=10)
            
        fig.tight_layout()

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

        Story.append(PageBreak())

        
        plotTitle(Story,"Pupil UVW and OPD shift from ACQ_CAM versus time")
        plotSubtitle(Story,"UVW in [meter], and OPD in [micron] from vertical "+\
                    astroreduced.telnames[0] +" = red, "+\
                    astroreduced.telnames[1] +" = orange, "+\
                    astroreduced.telnames[2] +" = blue, "+\
                    astroreduced.telnames[3] +" = green.")

        time = astroreduced.oi_acq_time[0::ntel]
        xmin = np.nanmin (time) - 1
        xmax = np.nanmax (time) + 1

        data = []
        ymin = -xlim
        ymax = xlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_pup_u[tel::ntel] == astroreduced.oi_acq_pup_u[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids],clipdata(astroreduced.oi_acq_pup_u[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None, ymin, ymax, None,
                              16*cm,5*cm,None,None,'Pupil_U [m]')
            Story.append(a)
        Story.append(Spacer(1,7*mm))

        data = []
        ymin = -xlim
        ymax = xlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_pup_v[tel::ntel] == astroreduced.oi_acq_pup_v[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids],clipdata(astroreduced.oi_acq_pup_v[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None, ymin, ymax, None,
                              16*cm,5*cm,None,None,'Pupil_V [m]')
            Story.append(a)
        Story.append(Spacer(1,7*mm))

        data = []
        ymin = - wlim
        ymax = wlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_pup_w[tel::ntel] == astroreduced.oi_acq_pup_w[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids],clipdata(astroreduced.oi_acq_pup_w[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None,ymin,ymax, None,
                              16*cm,3.4*cm,None,None,'Pupil_W [m]')
            Story.append(a)
        Story.append(Spacer(1,14*mm))

        data = []
        ymin = -opdlim
        ymax = opdlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_opd_pup[tel::ntel] == astroreduced.oi_acq_opd_pup[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids], clipdata(1E6 * astroreduced.oi_acq_opd_pup[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None, ymin, ymax, None,
                              16*cm,5*cm,None,None,'Pupil_OPD [micron]')
            Story.append(a)
                                
        Story.append(PageBreak())
                
    #==============================================================================
    # ACQ_CAM OI_VIS_ACQ XYZR
    #==============================================================================
    if astroreduced.acq_in is True:
        plotTitle(Story,"Pupil XYZ and R from ACQ_CAM versus time")
        plotSubtitle(Story,"XYZ in [pix] from reference, and R in [deg] from vertical "+\
                    astroreduced.telnames[0] +" = red, "+\
                    astroreduced.telnames[1] +" = orange, "+\
                    astroreduced.telnames[2] +" = blue, "+\
                    astroreduced.telnames[3] +" = green.")

        time = astroreduced.oi_acq_time[0::ntel]
        xmin = np.nanmin (time) - 1
        xmax = np.nanmax (time) + 1
        pixlim = np.nanmax([1.1*np.nanpercentile(np.abs(astroreduced.oi_acq_pup_x),95),
                            1.1*np.nanpercentile(np.abs(astroreduced.oi_acq_pup_y),95),
                            3]) # plot limit in pixels
        
        data = []
        ymin = np.nanmedian(astroreduced.oi_acq_pup_x) - pixlim
        ymax = ymin + 2*pixlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_pup_x[tel::ntel] == astroreduced.oi_acq_pup_x[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids],clipdata(astroreduced.oi_acq_pup_x[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None, ymin, ymax, None,
                              16*cm,5*cm,None,None,'Pupil_X [pix]')
            Story.append(a)
        Story.append(Spacer(1,7*mm))

        data = []
        ymin = np.nanmedian(astroreduced.oi_acq_pup_y) - pixlim
        ymax = ymin + 2*pixlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_pup_y[tel::ntel] == astroreduced.oi_acq_pup_y[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids],clipdata(astroreduced.oi_acq_pup_y[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None, ymin, ymax, None,
                              16*cm,5*cm,None,None,'Pupil_Y [pix]')
            Story.append(a)
        Story.append(Spacer(1,7*mm))

        data = []
        ymin = np.nanmedian (astroreduced.oi_acq_pup_z) - pixlim
        ymax = ymin + 2*pixlim
        for tel in range(0,ntel):
            ids = (astroreduced.oi_acq_pup_z[tel::ntel] == astroreduced.oi_acq_pup_z[tel::ntel]) * \
                  (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
            if (np.sum(ids) == 0):
                continue
            data.append(list(zip(time[ids],clipdata(astroreduced.oi_acq_pup_z[tel::ntel][ids],ymin,ymax))))
        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None,ymin,ymax, None,
                              16*cm,3.4*cm,None,None,'Pupil_Z [pix]')
            Story.append(a)
        Story.append(Spacer(1,14*mm))


        # Indices where any of the telescope pupil diodes are visible
        diodes_on_ids = np.where(astroreduced.oi_acq_pup_n[0::ntel]+
                                 astroreduced.oi_acq_pup_n[1::ntel]+
                                 astroreduced.oi_acq_pup_n[2::ntel]+
                                 astroreduced.oi_acq_pup_n[3::ntel] > 0)[0]
        # First index with diodes on
        first_diode_on = diodes_on_ids[0]*ntel
        # print("First diode on = %i"%(first_diode_on))
        
        times_ok = time[first_diode_on::]
        acq_pup_N_on = astroreduced.oi_acq_pup_n[first_diode_on::]
        
        data = []
        for tel in range(0,ntel):
            # Step 2*ntel to skip the frames with diodes off
            data.append(list(zip(times_ok[tel::2*ntel], acq_pup_N_on[tel::2*ntel])))

        if (len(data)>0):
            a = graphoutaxes(data, xmin, xmax, None, -1, 17, None,
                              16*cm,5*cm, None, 2.0,'Pupil_Nspot')
            Story.append(a)

        # data = []
        # for tel in range(0,ntel):
        #     ids = astroreduced.oi_acq_pup_r[tel::ntel] == astroreduced.oi_acq_pup_r[tel::ntel] * \
        #           (astroreduced.oi_acq_pup_n[tel::ntel] > 0)
        #     if (np.sum(ids) == 0):
        #         continue
        #     data.append(list(zip(time[ids],astroreduced.oi_acq_pup_r[tel::ntel][ids])))
        # if (len(data)>0):
        #     a = graphoutaxes(data, xmin, xmax, None,np.nanmin (astroreduced.oi_acq_pup_r-1),
        #                       np.nanmax (astroreduced.oi_acq_pup_r+1), None,
        #                       16*cm,3.4*cm,None,None,'Pupil_R [deg]')
        #     Story.append(a)
        # Story.append(Spacer(1,7*mm))
                                
        Story.append(PageBreak())
        
    #==============================================================================
    # Average photometric spectra
    #==============================================================================

    plotTitle(Story,"SC Average spectrum (OI_FLUX)")
    plotSubtitle(Story,"In photo-e- vs. wavelength (microns). "+\
                astroreduced.stations[0] +" = red, "+\
                astroreduced.stations[1] +" = orange, "+\
                astroreduced.stations[2] +" = blue, "+\
                astroreduced.stations[3] +" = green.")

    if hasattr(astroreduced,'oi_flux_sc') & (astroreduced.polarsplit == False):
        nmeasure = astroreduced.time_sc.shape[0]
        avgflux = np.zeros((4,astroreduced.nwave_sc),dtype='d')
        for tel in range(0,4):
            avgflux[tel,:]=np.nanmean(astroreduced.oi_flux_sc[:,tel,:],axis=0) # average flux over time
    
        yminval = 0 # np.nanpercentile(avgflux,1)-0.1*np.abs(np.nanpercentile(avgflux,99))
        ymaxval = 1.1*np.nanpercentile(avgflux,99)
        ticky = (ymaxval-yminval)/10 # in photo-e-
        hsize = 16*cm
        vsize = 12*cm
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(astroreduced.wave_sc, clipdata(avgflux[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(astroreduced.wave_sc),max(astroreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined-polarization')
        Story.append(a)
        
    if hasattr(astroreduced,'oi_flux_sc_s') & (astroreduced.polarsplit == True):
        nmeasure = astroreduced.time_sc.shape[0]
        avgfluxs = np.zeros((4,astroreduced.nwave_sc),dtype='d')
        avgfluxp = np.zeros((4,astroreduced.nwave_sc),dtype='d')
        for tel in range(0,4):
            avgfluxs[tel,:]=np.nanmean(astroreduced.oi_flux_sc_s[:,tel,:],axis=0) # average flux over time
            avgfluxp[tel,:]=np.nanmean(astroreduced.oi_flux_sc_p[:,tel,:],axis=0) # average flux over time
    
        yminval = 0 # np.nanpercentile(avgfluxs,1)-0.1*np.abs(np.nanpercentile(avgfluxs,99))
        ymaxval = 1.1*np.nanpercentile([avgfluxs,avgfluxp],99)
        ticky = (ymaxval-yminval)/10 # in photo-e-
        hsize = 16*cm
        vsize = 9*cm
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(astroreduced.wave_sc, clipdata(avgfluxs[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(astroreduced.wave_sc),max(astroreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'S-polarization')
        Story.append(a)
        
        Story.append(Spacer(1,7*mm))
    
        ticky = (ymaxval-yminval)/10 # in photo-e-
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(astroreduced.wave_sc, clipdata(avgfluxp[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(astroreduced.wave_sc),max(astroreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'P-polarization')
        Story.append(a)

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

    Story.append(PageBreak())

    #==============================================================================
    # Image (waterfall) display of SC OIFLUX for the four telescopes
    #==============================================================================
    plotTitle(Story,"SC flux waterfall per telescope (OI_FLUX SC)")
    plotSubtitle(Story,"In photo-e-/DIT, wavelength channel number in horizontal axis, number of frame in sequence in vertical axis. The average over the sequence has been subtracted and is given above each image.")

    if astroreduced.polarsplit == False:
        meanflux = np.zeros((4))
        meansub = np.zeros((astroreduced.nframe_sc,4,astroreduced.nwave_sc))        
        for tel in range(0,4):
            meanflux[tel] = np.nanmean(astroreduced.oi_flux_sc[:,tel,:])
            meansub[:,tel,:] = astroreduced.oi_flux_sc[:,tel,:] - meanflux[tel]
            
        valmin = np.nanpercentile(meansub,1)
        valmax = np.nanpercentile(meansub,99)

        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7))
        im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[0].set_title('SC OIFLUX Telescope '+astroreduced.stations[0]+' - Comb - Mean = %.1e' % meanflux[0], fontsize=8)
        axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[1].set_title('SC OIFLUX Telescope '+astroreduced.stations[1]+' - Comb - Mean = %.1e' % meanflux[1], fontsize=8)
        axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[2].set_title('SC OIFLUX Telescope '+astroreduced.stations[2]+' - Comb - Mean = %.1e' % meanflux[2], fontsize=8)
        axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[3].set_title('SC OIFLUX Telescope '+astroreduced.stations[3]+' - Comb - Mean = %.1e' % meanflux[3], fontsize=8)
        fig.subplots_adjust(hspace=.3)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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

    if astroreduced.polarsplit == True:
        meanflux = np.zeros((4))
        meansub = np.zeros((astroreduced.nframe_sc,4,astroreduced.nwave_sc))        
        for tel in range(0,4):
            meanflux[tel] = np.nanmean(astroreduced.oi_flux_sc_s[:,tel,:])
            meansub[:,tel,:] = astroreduced.oi_flux_sc_s[:,tel,:] - meanflux[tel]
            
        valmin = np.nanpercentile(meansub,1)
        valmax = np.nanpercentile(meansub,99)

        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7))
        im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[0].set_title('SC OIFLUX Telescope '+astroreduced.stations[0]+' - S - Mean = %.1e' % meanflux[0], fontsize=8)
        axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[1].set_title('SC OIFLUX Telescope '+astroreduced.stations[1]+' - S - Mean = %.1e' % meanflux[1], fontsize=8)
        axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[2].set_title('SC OIFLUX Telescope '+astroreduced.stations[2]+' - S - Mean = %.1e' % meanflux[2], fontsize=8)
        axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[3].set_title('SC OIFLUX Telescope '+astroreduced.stations[3]+' - S - Mean = %.1e' % meanflux[3], fontsize=8)
        fig.subplots_adjust(hspace=.3)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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

        plotTitle(Story,"SC flux waterfall per telescope (OI_FLUX SC)")
        plotSubtitle(Story,"In photo-e-/DIT, wavelength channel number in horizontal axis, number of frame in sequence in vertical axis. The average over the sequence has been subtracted and is given above each image.")

        meansub = np.zeros((astroreduced.nframe_sc,4,astroreduced.nwave_sc))        
        meanflux = np.zeros((4))
        for tel in range(0,4):
            meanflux[tel] = np.nanmean(astroreduced.oi_flux_sc_p[:,tel,:])
            meansub[:,tel,:] = astroreduced.oi_flux_sc_p[:,tel,:] - meanflux[tel]
            
        valmin = np.nanpercentile(meansub,1)
        valmax = np.nanpercentile(meansub,99)

        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7))
        im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[0].set_title('SC OIFLUX Telescope '+astroreduced.stations[0]+' - P - Mean = %.1e' % meanflux[0], fontsize=8)
        axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[1].set_title('SC OIFLUX Telescope '+astroreduced.stations[1]+' - P - Mean = %.1e' % meanflux[1], fontsize=8)
        axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[2].set_title('SC OIFLUX Telescope '+astroreduced.stations[2]+' - P - Mean = %.1e' % meanflux[2], fontsize=8)
        axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[3].set_title('SC OIFLUX Telescope '+astroreduced.stations[3]+' - P - Mean = %.1e' % meanflux[3], fontsize=8)
        fig.subplots_adjust(hspace=.3)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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

    #==============================================================================
    # Flux signals
    #==============================================================================

    plotTitle(Story,"SC flux vs. time (OI_FLUX wavelength averaged)")
    plotSubtitle(Story,"In photo-e-/DIT vs. seconds. "+\
                astroreduced.stations[0] +" = red, "+\
                astroreduced.stations[1] +" = orange, "+\
                astroreduced.stations[2] +" = blue, "+\
                astroreduced.stations[3] +" = green.")
                
    if hasattr(astroreduced,'oi_flux_sc') & (astroreduced.polarsplit == False):
        nmeasure = astroreduced.time_sc.shape[0]
        avgflux = np.zeros((4,nmeasure),dtype='d')
        for tel in range(0,4):
            avgflux[tel,:]=np.nanmean(astroreduced.oi_flux_sc[:,tel,:],axis=1) # average flux over wavelength
    
        yminval = 0
        ymaxval = np.nanmax(avgflux)+0.1*np.abs(np.nanmax(avgflux))
        ticky = (ymaxval-yminval)/10 # in photo-e-
        hsize = 16*cm
        vsize = 10*cm
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(astroreduced.time_sc, clipdata(avgflux[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(astroreduced.time_sc),max(astroreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined-polarization')
        Story.append(a)
        
        Story.append(PageBreak())

    if hasattr(astroreduced,'oi_flux_sc_s') & (astroreduced.polarsplit == True):
        nmeasure = astroreduced.time_sc.shape[0]
        avgfluxs = np.zeros((4,nmeasure),dtype='d')
        avgfluxp = np.zeros((4,nmeasure),dtype='d')
        for tel in range(0,4):
            avgfluxs[tel,:]=np.nanmean(astroreduced.oi_flux_sc_s[:,tel,:],axis=1) # average flux over wavelength
            avgfluxp[tel,:]=np.nanmean(astroreduced.oi_flux_sc_p[:,tel,:],axis=1) # average flux over wavelength
    
        yminval = 0
        ymaxval = np.nanmax(avgfluxs)+0.1*np.abs(np.nanmax(avgfluxs))
        ticky = (ymaxval-yminval)/10 # in photo-e-
        hsize = 16*cm
        vsize = 10*cm
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(astroreduced.time_sc, clipdata(avgfluxs[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(astroreduced.time_sc),max(astroreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'S-polarization')
        Story.append(a)
        data = []

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

        for tel in range(0,4):
            data.append( tuple(zip(astroreduced.time_sc, clipdata(avgfluxp[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes( data, min(astroreduced.time_sc),max(astroreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'P-polarization')
        Story.append(a)
 
        Story.append(PageBreak())


    #==============================================================================
    # SC REJECTION_FLAG
    #==============================================================================
    plotTitle(Story,"SC REJECTION_FLAG per baseline")

    if astroreduced.polarsplit == False:
        ymaxval = np.max([4,np.max(astroreduced.oi_vis_reject_flag)])
    else:
        ymaxval = np.max([4,np.max(astroreduced.oi_vis_reject_flag_s)])

    if ymaxval <= 4:
        plotSubtitle(Story,"0= accepted frame, 1= Low FT tracking, 2= Low FT vFactor, 3= both")
    else:
        Story.append(Spacer(1,7*mm))
        
    resamp = 1 # take one value every resamp value for the plot
    yminval = -1
    ticky = 1 # visibility tick
    xval = astroreduced.time_sc[::resamp]
    xminval = np.min(xval)
    xmaxval = np.max(xval)

    d=()
    hsize=16*cm
    vsize=3*cm

    for baseline in range(0,nbase):
        data=[]
        if astroreduced.polarsplit == False:
            data.append (list(zip(xval, astroreduced.oi_vis_reject_flag[::resamp,baseline])))
        else:
            data.append (list(zip(xval, astroreduced.oi_vis_reject_flag_s[::resamp,baseline])))
            data.append (list(zip(xval, astroreduced.oi_vis_reject_flag_p[::resamp,baseline]+0.05)))
        
        a = graphoutaxes (data, xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+astroreduced.basenames[baseline])
        #a.add(String(0.05*hsize, 0.22*vsize, '0=Good', fontSize=9))
        #a.add(String(0.05*hsize, 0.42*vsize, '1=low FT tracking', fontSize=9))
        #a.add(String(0.05*hsize, 0.62*vsize, '2=low FT Vfactor', fontSize=9))
        #a.add(String(0.05*hsize, 0.82*vsize, '3=low FT Vfactor and low FT tracking', fontSize=9))
        
        Story.append(a)
        Story.append(Spacer(1,7*mm))
 
    Story.append(PageBreak())
    
    #==============================================================================
    # Plotting of the SC group delay signals
    #==============================================================================
    plotTitle(Story,"SC group delay vs. time")
    plotSubtitle(Story,"In OPD (microns) vs. time (seconds). Red = S, Orange = P.")

    if astroreduced.gdelay_in == True:

        timescale_sc = astroreduced.time_sc
      
        hsize=16*cm
        vsize=3*cm
        for baseline in range(0,nbase):
            if (astroreduced.polarsplit == False):
                tmp1 = astroreduced.gdelay_sc[baseline,:]
                yminval = np.nanmin(tmp1)
                ymaxval = np.nanmax(tmp1)
                data = [tuple(zip(timescale_sc[:], clipdata(tmp1,yminval,ymaxval)))]
            else:
                tmp1 = astroreduced.gdelay_sc_s[baseline,:]
                tmp2 = astroreduced.gdelay_sc_p[baseline,:]
                yminval = np.nanmin(tmp1)
                ymaxval = np.nanmax(tmp1)
                data = [tuple(zip(timescale_sc[:], clipdata(tmp1,yminval,ymaxval)))]
                data.append(tuple(zip(timescale_sc[:], clipdata(tmp2,yminval,ymaxval))))

            ticky = (ymaxval-yminval)/5
            yminval -= ticky
            ymaxval += ticky
            a = graphoutaxes( data, np.nanmin(timescale_sc[:]), np.nanmax(timescale_sc[:]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+astroreduced.basenames[baseline])
            Story.append(a)
            Story.append(Spacer(1,7*mm))
      
    Story.append(PageBreak())

    #==============================================================================
    # Plotting of the SC coherent flux SNR signals from PIPELINE
    #==============================================================================
    plotTitle(Story,"SC coherent flux SNR vs. time (OI_VIS/SNR)")
    plotSubtitle(Story,"SNR vs. time (seconds). Red = S, Orange = P.")
            
    if astroreduced.gdelay_in == True:
        
        timescale_sc = astroreduced.time_sc
        d=()
        data=[]
        hsize=16*cm
        vsize=3*cm

        for baseline in range(0,nbase):
            data = []
            yminval = 0
            if (astroreduced.polarsplit == False):
                ymaxval = np.nanpercentile(astroreduced.snr_sc[baseline,0:nframeplot],98)
                data.append(tuple(zip(timescale_sc[0:nframeplot], clipdata(astroreduced.snr_sc[baseline,0:nframeplot],yminval,ymaxval))))
            else:
                ymaxval = np.nanpercentile(astroreduced.snr_sc_s[baseline,0:nframeplot],98)
                data.append(tuple(zip(timescale_sc[0:nframeplot], clipdata(astroreduced.snr_sc_s[baseline,0:nframeplot],yminval,ymaxval))))
                data.append(tuple(zip(timescale_sc[0:nframeplot], clipdata(astroreduced.snr_sc_p[baseline,0:nframeplot],yminval,ymaxval))))
               
            ticky = (ymaxval-yminval)/5
            yminval -= ticky
            ymaxval += ticky
            
            a = graphoutaxes( data, np.nanmin(timescale_sc[0:nframeplot]), np.nanmax(timescale_sc[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+astroreduced.basenames[baseline])
            Story.append(a)
            Story.append(Spacer(1,7*mm))
   
    Story.append(PageBreak())

    #==============================================================================
    # Average V factor (quadratic quantity) computed for SC vs. wavelength
    #==============================================================================

    if hasattr(astroreduced,'oi_vis_vfactor_ft') or hasattr(astroreduced,'oi_vis_vfactor_ft_s'):
        plotTitle(Story,"V_FACTOR and V_FACTOR_FT vs. wavelength (averaged over frames)")
        
        if (astroreduced.polarsplit == False):
            plotSubtitle(Story,"Wavelength in microns, red = SC, orange = FT.")
        else:
            plotSubtitle(Story,"Wavelength in microns, red = SC S, orange = SC P, blue = FT S, green = FT P.")
            
        hsize = 6*cm
        vsize = 6*cm
        yminval = 0.0
        ymaxval = 1.2
        xminval = np.min(astroreduced.wave_sc)
        xmaxval = np.max(astroreduced.wave_sc)
        ticky = 0.25
        tickx = 0.25
        
        d = []
        for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength]
            data=[]
            if (astroreduced.polarsplit == False):
                data.append( tuple(zip(astroreduced.wave_sc, np.mean(astroreduced.oi_vis_vfactor[:,baseline,:],axis=0)) ))
                data.append( tuple(zip(astroreduced.wave_ft, np.mean(astroreduced.oi_vis_vfactor_ft[:,baseline,:],axis=0)) ))
            else:
                data.append( tuple(zip(astroreduced.wave_sc, np.mean(astroreduced.oi_vis_vfactor_s[:,baseline,:],axis=0)) ))
                data.append( tuple(zip(astroreduced.wave_sc, np.mean(astroreduced.oi_vis_vfactor_p[:,baseline,:],axis=0)) ))
                data.append( tuple(zip(astroreduced.wave_ft, np.mean(astroreduced.oi_vis_vfactor_ft_s[:,baseline,:],axis=0)) ))
                data.append( tuple(zip(astroreduced.wave_ft, np.mean(astroreduced.oi_vis_vfactor_ft_p[:,baseline,:],axis=0)) ))

            a = graphoutaxes( data,xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+astroreduced.basenames[baseline])
            d.append(a)
     
        plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ]
        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())

  
    #==============================================================================
    # Image (waterfall) display of model V_FACTOR (of SC)
    #==============================================================================
    plotTitle(Story,"V_FACTOR_SC waterfall vs. wavelength and time")
    plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of SC frame in vertical axis. Mean spectrum over time subtracted.")

    if astroreduced.polarsplit == False:

        meanv = np.zeros((nbase,astroreduced.nwave_sc))
        for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength]
            meanv[baseline,:] = np.nanmean(astroreduced.oi_vis_vfactor[:,baseline,:],axis=0)

        meansub = np.zeros((nbase,astroreduced.nframe_sc,astroreduced.nwave_sc))        
        for baseline in range(0,nbase):
            for frame in range(0,astroreduced.nframe_sc):
                meansub[baseline,frame,:] = astroreduced.oi_vis_vfactor[frame,baseline,:] - meanv[baseline,:]
           
        valmin = np.nanpercentile(meansub,1)-0.05
        valmax = np.nanpercentile(meansub,99)+0.05
        #print(valmin, valmax)

        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=nbase, ncols=1, sharex=True, sharey=True, figsize=(5,8))
        for baseline in range(0,nbase):
            im = axarr[baseline].imshow(meansub[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[baseline].set_title('Baseline '+astroreduced.basenames[baseline], fontsize=7)
        fig.subplots_adjust(hspace=.2)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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


    if astroreduced.polarsplit == True:

        meanvs = np.zeros((nbase,astroreduced.nwave_sc))
        meanvp = np.zeros((nbase,astroreduced.nwave_sc))

        for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength]
            meanvs[baseline,:] = np.nanmean(astroreduced.oi_vis_vfactor_s[:,baseline,:],axis=0)
            meanvp[baseline,:] = np.nanmean(astroreduced.oi_vis_vfactor_p[:,baseline,:],axis=0)

        meansubs = np.zeros((nbase,astroreduced.nframe_sc,astroreduced.nwave_sc))        
        meansubp = np.zeros((nbase,astroreduced.nframe_sc,astroreduced.nwave_sc))        
        for baseline in range(0,nbase):
            for frame in range(0,astroreduced.nframe_sc):
                meansubs[baseline,frame,:] = astroreduced.oi_vis_vfactor_s[frame,baseline,:] - meanvs[baseline,:]
                meansubp[baseline,frame,:] = astroreduced.oi_vis_vfactor_p[frame,baseline,:] - meanvp[baseline,:]
           
        valmin = np.nanpercentile([meansubs,meansubp],1)-0.05
        valmax = np.nanpercentile([meansubs,meansubp],99)+0.05

        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=nbase, ncols=2, sharex=True, sharey=True, figsize=(5,8))
        for baseline in range(0,nbase):
            im = axarr[baseline,0].imshow(meansubs[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[baseline,0].set_title('Baseline '+astroreduced.basenames[baseline]+'-S', fontsize=7)
            im = axarr[baseline,1].imshow(meansubp[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[baseline,1].set_title('Baseline '+astroreduced.basenames[baseline]+'-P', fontsize=7)
        fig.subplots_adjust(hspace=.2)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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

    #==============================================================================
    # V_FACTOR versus time
    #==============================================================================

    if hasattr(astroreduced,'oi_vis_vfactor') or hasattr(astroreduced,'oi_vis_vfactor_s'):
        plotTitle(Story,"V_FACTOR for SC versus time")
        plotSubtitle(Story,"The V_FACTOR is averaged over the band")
    
        timescale_sc = astroreduced.time_sc;
        d=()
        data=[]
        hsize=16*cm
        vsize=3*cm
        yminval = 0
    
        for baseline in range(0,nbase):
            data = []
            if (astroreduced.polarsplit == True):
                d1 = np.nanmean(astroreduced.oi_vis_vfactor_s[:,baseline,:],axis=0);
                d2 = np.nanmean(astroreduced.oi_vis_vfactor_p[:,baseline,:],axis=0);
                ymaxval = np.maximum (d1.max(),d2.max());
                data.append(tuple(zip(timescale_sc, clipdata(d1,yminval,ymaxval))));
                data.append(tuple(zip(timescale_sc, clipdata(d2,yminval,ymaxval))));
            else:
                d1 = np.nanmean(astroreduced.oi_vis_vfactor[:,baseline,:],axis=0);
                ymaxval = d1.max();
                data.append(tuple(zip(timescale_sc, clipdata(d1,yminval,ymaxval))))

            ticky = (ymaxval-yminval) / 5;
            
            a = graphoutaxes(data, np.nanmin(timescale_sc), np.nanmax(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+astroreduced.basenames[baseline])
            Story.append(a)
            Story.append(Spacer(1,7*mm))
        
        Story.append(PageBreak())

    #==============================================================================
    # P_FACTOR versus time
    #==============================================================================

    if hasattr(astroreduced,'oi_vis_pfactor') or hasattr(astroreduced,'oi_vis_pfactor_s'):
        plotTitle(Story,"P_FACTOR for SC versus time")
        plotSubtitle(Story,"The P_FACTOR is computed broad-band")
    
        timescale_sc = astroreduced.time_sc
        d=()
        data=[]
        hsize=16*cm
        vsize=3*cm
        yminval = -0.1
        ymaxval =  1.1
        ticky = 0.2
    
        for base in range(0,nbase):
            data = []
            if (astroreduced.polarsplit == True):
                data.append(tuple(zip(timescale_sc, clipdata(astroreduced.oi_vis_pfactor_s[base::nbase],yminval,ymaxval))))
                data.append(tuple(zip(timescale_sc, clipdata(astroreduced.oi_vis_pfactor_p[base::nbase],yminval,ymaxval))))
            else:
                data.append(tuple(zip(timescale_sc, clipdata(astroreduced.oi_vis_pfactor[base::nbase],yminval,ymaxval))))
               
            a = graphoutaxes( data, np.nanmin(timescale_sc), np.nanmax(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+astroreduced.basenames[base])
            Story.append(a)
            Story.append(Spacer(1,7*mm))
        
        Story.append(PageBreak())
    

    #==============================================================================
    # V factor (quadratic quantity) computed for SC versus VIS2
    #==============================================================================
    if True:
        if hasattr(astroreduced,'oi_vis_vfactor_ft') or hasattr(astroreduced,'oi_vis_vfactor_ft_s'):
            plotTitle(Story,"SC VIS2 versus V_FACTOR (OI_VIS2, wavelength averaged)")
            plotSubtitle(Story,"VIS2 (vertical) versus V_FACTOR (horizontal), Red = S or combined.")
        
            hsize = 6*cm
            vsize = 6*cm
            ymaxval = 1.2
            yminval = 0.0
            ticky = 0.25
            tickx = 0.25
            
            d = []
            for baseline in range(0,nbase):
                data=[]
                if (astroreduced.polarsplit == False):
                    data.append( tuple(zip(clipdata(np.mean(astroreduced.oi_vis_vfactor[:,baseline,:],axis=1),yminval,ymaxval), \
                        clipdata(np.mean(astroreduced.oi_vis2_sc[baseline,:,:],axis=1),yminval,ymaxval) )))
                else:
                    data.append( tuple(zip(clipdata(np.mean(astroreduced.oi_vis_vfactor_s[:,baseline,:],axis=1),yminval,ymaxval), \
                        clipdata(np.mean(astroreduced.oi_vis2_sc_s[baseline,:,:],axis=1),yminval,ymaxval) )))
                a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+astroreduced.basenames[baseline])
                d.append(a)
         
            plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ]
            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())


    #==============================================================================
    # VIS2 * F1 * F2 vs. F1 * F2 of SC
    #==============================================================================
    plotTitle(Story,"SC  VIS2F1F2 vs. F1*F2  (OI_VIS, wavelength averaged)")
    plotSubtitle(Story,"In photo-e-^2/SC DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.")

    resamp = 1
    hsize = 6*cm
    vsize = 6*cm
    d = []

    for baseline in range(0,nbase):

        if (astroreduced.polarsplit == True):
            f1f2 = np.mean(astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],:]*astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],:],axis=1)
            vis2 = np.mean(astroreduced.oi_vis2_sc_s[baseline,::resamp,:],axis=1)
        else:
            f1f2 = np.mean(astroreduced.oi_flux_sc[::resamp,tel_list[baseline,0],:]*astroreduced.oi_flux_sc[::resamp,tel_list[baseline,1],:],axis=1)
            vis2 = np.mean(astroreduced.oi_vis2_sc[baseline,::resamp,:],axis=1)

        ymaxval = 1.05*np.max(f1f2)
        yminval = -0.05*ymaxval
        ticky = (ymaxval-yminval)/3.
        tickx = ticky
        
        data=[]
        data.append( tuple(zip(f1f2, clipdata(vis2*f1f2,yminval,ymaxval)) ))
        a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+astroreduced.basenames[baseline])
        d.append(a)

    plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ]
    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())

    #==============================================================================
    # VIS2 * F1 * F2 vs. F1 * F2 of SC
    #==============================================================================
    plotTitle(Story,"SC  VIS2F1F2 vs. F1*F2*VFACTOR  (OI_VIS, wavelength averaged)")
    plotSubtitle(Story,"In photo-e-^2/SC DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.")

    resamp = 1
    hsize = 6*cm
    vsize = 6*cm
    d = []

    for baseline in range(0,nbase):

        if (astroreduced.polarsplit == True):
            f1f2 = np.mean(astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],:]*astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],:]*astroreduced.oi_vis_vfactor_s[::resamp,baseline,:],axis=1)
            vis2 = np.mean(astroreduced.oi_vis2_sc_s[baseline,::resamp,:]*astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],:]*astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],:],axis=1)
        else:
            f1f2 = np.mean(astroreduced.oi_flux_sc[::resamp,tel_list[baseline,0],:]*astroreduced.oi_flux_sc[::resamp,tel_list[baseline,1],:]*astroreduced.oi_vis_vfactor[::resamp,baseline,:],axis=1)
            vis2 = np.mean(astroreduced.oi_vis2_sc[baseline,::resamp,:]*astroreduced.oi_flux_sc[::resamp,tel_list[baseline,0],:]*astroreduced.oi_flux_sc[::resamp,tel_list[baseline,1],:],axis=1)

        ymaxval = 1.05*np.max(f1f2)
        yminval = -0.05*ymaxval
        ticky = (ymaxval-yminval)/3.
        tickx = ticky
        
        data=[]
        data.append( tuple(zip(f1f2, clipdata(vis2,yminval,ymaxval)) ))
        a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+astroreduced.basenames[baseline])
        d.append(a)

    plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ]
    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())
    
    #==============================================================================
    # VIS2 * F1 * F2 vs. F1 * F2 of SC
    #==============================================================================
    plotTitle(Story,"SC  VIS2_F1F2 vs. F1*F2  (OI_VIS, 2 channels)")
    plotSubtitle(Story,"In photo-e-^2/FT DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.")

    resamp = 1
    hsize = 6*cm
    vsize = 6*cm
    d = []
    st1 = int(astroreduced.nwave_sc / 4)
    sp1  = int(astroreduced.nwave_sc / 4 + 1)
    st2 = int(astroreduced.nwave_sc / 4 * 3)
    sp2  = int(astroreduced.nwave_sc / 4 * 3 + 1)

    for baseline in range(0,nbase):
        data=[]

        if (astroreduced.polarsplit == True):
            f1f2a = np.mean(astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],st1:sp1]*astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],st1:sp1],axis=1)
            vis2a = np.mean(astroreduced.oi_vis2_sc_s[baseline,::resamp,st1:sp1],axis=1)
            f1f2b = np.mean(astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],st2:sp2]*astroreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],st2:sp2],axis=1)
            vis2b = np.mean(astroreduced.oi_vis2_sc_s[baseline,::resamp,st2:sp2],axis=1)
        else:
            f1f2a = np.mean(astroreduced.oi_flux_sc[::resamp,tel_list[baseline,0],st1:sp1]*astroreduced.oi_flux_sc[::resamp,tel_list[baseline,1],st1:sp1],axis=1)
            vis2a = np.mean(astroreduced.oi_vis2_sc[baseline,::resamp,st1:sp1],axis=1)
            f1f2b = np.mean(astroreduced.oi_flux_sc[::resamp,tel_list[baseline,0],st2:sp2]*astroreduced.oi_flux_sc[::resamp,tel_list[baseline,1],st2:sp2],axis=1)
            vis2b = np.mean(astroreduced.oi_vis2_sc[baseline,::resamp,st2:sp2],axis=1)

        ymaxval = 1.05*np.max(f1f2a)
        yminval = -0.05*ymaxval
        ticky = (ymaxval-yminval)/3.
        tickx = ticky
        
        data.append( tuple(zip(f1f2a, clipdata(vis2a*f1f2a,yminval,ymaxval)) ))
        data.append( tuple(zip(f1f2b, clipdata(vis2b*f1f2b,yminval,ymaxval)) ))
        
        a = graphscatteraxes (data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+astroreduced.basenames[baseline])
        d.append(a)

    plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ]
    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())
          

    #==============================================================================
    # Image (waterfall) display of SC VISAMP
    #==============================================================================
    plotTitle(Story,"SC visibility amplitude waterfall (VISAMP)")
    plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.")

    if hasattr(astroreduced,'oi_vis_sc_visamp') & (astroreduced.polarsplit == False):
        valmin = np.max([0,np.nanpercentile(astroreduced.oi_vis_sc_visamp,1)])
        valmax = np.min([1.3,np.nanpercentile(astroreduced.oi_vis_sc_visamp,99)])
        
        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7))
        im = axarr[0].imshow(astroreduced.oi_vis_sc_visamp[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[0].set_title('SC VISAMP Baseline '+astroreduced.basenames[0]+' - Comb', fontsize=8)
        axarr[1].imshow(astroreduced.oi_vis_sc_visamp[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[1].set_title('SC VISAMP Baseline '+astroreduced.basenames[1]+' - Comb', fontsize=8)
        axarr[2].imshow(astroreduced.oi_vis_sc_visamp[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[2].set_title('SC VISAMP Baseline '+astroreduced.basenames[2]+' - Comb', fontsize=8)
        axarr[3].imshow(astroreduced.oi_vis_sc_visamp[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[3].set_title('SC VISAMP Baseline '+astroreduced.basenames[3]+' - Comb', fontsize=8)
        axarr[4].imshow(astroreduced.oi_vis_sc_visamp[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[4].set_title('SC VISAMP Baseline '+astroreduced.basenames[4]+' - Comb', fontsize=8)
        axarr[5].imshow(astroreduced.oi_vis_sc_visamp[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[5].set_title('SC VISAMP Baseline '+astroreduced.basenames[5]+' - Comb', fontsize=8)
        fig.subplots_adjust(hspace=.3)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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


    if hasattr(astroreduced,'oi_vis_sc_visamp_s') & (astroreduced.polarsplit == True):
        valmin = np.max([0,np.nanpercentile(astroreduced.oi_vis_sc_visamp_s,1)])
        valmax = np.min([1.3,np.nanpercentile(astroreduced.oi_vis_sc_visamp_s,99)])
        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7))
        im = axarr[0].imshow(astroreduced.oi_vis_sc_visamp_s[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[0].set_title('SC VISAMP Baseline '+astroreduced.basenames[0]+' - S', fontsize=8)
        axarr[1].imshow(astroreduced.oi_vis_sc_visamp_s[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[1].set_title('SC VISAMP Baseline '+astroreduced.basenames[1]+' - S', fontsize=8)
        axarr[2].imshow(astroreduced.oi_vis_sc_visamp_s[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[2].set_title('SC VISAMP Baseline '+astroreduced.basenames[2]+' - S', fontsize=8)
        axarr[3].imshow(astroreduced.oi_vis_sc_visamp_s[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[3].set_title('SC VISAMP Baseline '+astroreduced.basenames[3]+' - S', fontsize=8)
        axarr[4].imshow(astroreduced.oi_vis_sc_visamp_s[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[4].set_title('SC VISAMP Baseline '+astroreduced.basenames[4]+' - S', fontsize=8)
        axarr[5].imshow(astroreduced.oi_vis_sc_visamp_s[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[5].set_title('SC VISAMP Baseline '+astroreduced.basenames[5]+' - S', fontsize=8)
        fig.subplots_adjust(hspace=.3)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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

        plotTitle(Story,"SC visibility amplitude waterfall (VISAMP)")
        plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.")


        plt.close('all')
        plt.rc('xtick', labelsize='x-small')
        plt.rc('ytick', labelsize='x-small')
        fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7))
        im = axarr[0].imshow(astroreduced.oi_vis_sc_visamp_p[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[0].set_title('SC VISAMP Baseline '+astroreduced.basenames[0]+' - P', fontsize=8)
        axarr[1].imshow(astroreduced.oi_vis_sc_visamp_p[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[1].set_title('SC VISAMP Baseline '+astroreduced.basenames[1]+' - P', fontsize=8)
        axarr[2].imshow(astroreduced.oi_vis_sc_visamp_p[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[2].set_title('SC VISAMP Baseline '+astroreduced.basenames[2]+' - P', fontsize=8)
        axarr[3].imshow(astroreduced.oi_vis_sc_visamp_p[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[3].set_title('SC VISAMP Baseline '+astroreduced.basenames[3]+' - P', fontsize=8)
        axarr[4].imshow(astroreduced.oi_vis_sc_visamp_p[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[4].set_title('SC VISAMP Baseline '+astroreduced.basenames[4]+' - P', fontsize=8)
        axarr[5].imshow(astroreduced.oi_vis_sc_visamp_p[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
        axarr[5].set_title('SC VISAMP Baseline '+astroreduced.basenames[5]+' - P', fontsize=8)
        fig.subplots_adjust(hspace=.3)
        # Color bar plot
        fig.subplots_adjust(right=0.8)
        cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar(im, cax=cbar_ax)

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

    #==============================================================================
    # Wavelength average SC VISAMP vs. group delay
    #==============================================================================
    plotTitle(Story,"SC visibility amp. vs. group delay (VISAMP wavelength averaged)")
    plotSubtitle(Story,"Visibility vs. group delay (microns). Red = S, Orange = P.  (VISAMP/GDELAY)")

    if astroreduced.gdelay_in == True:
    
        d=()
        data=[]
        resamp = 1 # take one value every resamp value for the plot
        hsize=16*cm
        vsize=3*cm

        if (astroreduced.polarsplit == False):
            visamp_wavg = np.zeros((nbase,astroreduced.oi_vis_sc_visamp[baseline::6*resamp,:].shape[0]),dtype='d')
            
            for baseline in range(0,nbase):
                visamp_wavg[baseline,:] = np.nanmean(astroreduced.oi_vis_sc_visamp[baseline::6*resamp,:], axis=1)
                print("Baseline %s = %s mean VISAMP = %.3f"%(gravitybasenames[baseline],astroreduced.basenames[baseline],np.nanmean(visamp_wavg[baseline,:])))
                
            mingd = np.nanmin(astroreduced.gdelay_sc)
            maxgd = np.nanmax(astroreduced.gdelay_sc)
            for baseline in range(0,nbase):
                data=[]
                yminval = np.nanpercentile(visamp_wavg[baseline,:],2)-0.1
                ymaxval = np.nanpercentile(visamp_wavg[baseline,:],98)+0.1
                ticky = (ymaxval - yminval)/5.
                tickx = (maxgd-mingd)/5.
                data.append(list(zip(astroreduced.gdelay_sc[baseline,::resamp], clipdata(visamp_wavg[baseline,:],yminval,ymaxval))))
                a = graphscatteraxes(data, mingd, maxgd,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,astroreduced.basenames[baseline])
                Story.append(a)
                Story.append(Spacer(1,7*mm))
        
        if (astroreduced.polarsplit == True):
            visamp_wavg_s = np.zeros((nbase,astroreduced.oi_vis_sc_visamp_s[baseline::6*resamp,:].shape[0]),dtype='d')
            visamp_wavg_p = np.zeros((nbase,astroreduced.oi_vis_sc_visamp_p[baseline::6*resamp,:].shape[0]),dtype='d')
            
            for baseline in range(0,nbase):
                visamp_wavg_s[baseline,:] = np.nanmean(astroreduced.oi_vis_sc_visamp_s[baseline::6*resamp,:], axis=1)
                visamp_wavg_p[baseline,:] = np.nanmean(astroreduced.oi_vis_sc_visamp_p[baseline::6*resamp,:], axis=1)
                print("Baseline %s = %s S polar mean VISAMP = %.3f   P polar mean VISAMP = %.3f"%\
                            (gravitybasenames[baseline],astroreduced.basenames[baseline],\
                            np.nanmean(visamp_wavg_s[baseline,:]),np.nanmean(visamp_wavg_p[baseline,:])))
            
            mingd = np.nanmin([astroreduced.gdelay_sc_s,astroreduced.gdelay_sc_p])
            maxgd = np.nanmax([astroreduced.gdelay_sc_s,astroreduced.gdelay_sc_p])
            for baseline in range(0,nbase):
                data=[]
                yminval = np.nanpercentile([visamp_wavg_s[baseline,:],visamp_wavg_p[baseline,:]],2)-0.1
                ymaxval = np.nanpercentile([visamp_wavg_s[baseline,:],visamp_wavg_p[baseline,:]],98)+0.1
                ticky = (ymaxval - yminval)/5.
                tickx = (maxgd-mingd)/5.
                data.append(list(zip(astroreduced.gdelay_sc_s[baseline,::resamp], clipdata(visamp_wavg_s[baseline,:],yminval,ymaxval))))
                data.append(list(zip(astroreduced.gdelay_sc_p[baseline,::resamp], clipdata(visamp_wavg_p[baseline,:],yminval,ymaxval))))
                a = graphscatteraxes(data,mingd,maxgd,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,astroreduced.basenames[baseline])
                Story.append(a)
                Story.append(Spacer(1,7*mm))
    
    Story.append(PageBreak())

    #==============================================================================
    # Temporal quality of the PHASE_REF arg(VISDATA_SC * PHASE_REF)
    #==============================================================================
    if hasattr(astroreduced,'oi_vis_sc_exp_phase') or hasattr(astroreduced,'oi_vis_sc_exp_phase_s'):

        minw = int(0.4*astroreduced.nwave_sc)
        maxw = int(0.6*astroreduced.nwave_sc)
      
        plotTitle(Story,"arg(VISDATA_SC * PHASE_REF) vs. time")
        plotSubtitle(Story,"If I am not here but you want me, fix me in the python tools!")

        resamp = 1 # take one value every resamp value for the plot
        if hasattr(astroreduced,'oi_vis_sc_phi_met') & hasattr(astroreduced,'oi_vis_sc_visphi') & (astroreduced.polarsplit == False):
            plotSubtitle(Story,"If phase is not flat => DITSHIFT. In rad vs. time (seconds). Mean over wavelengths.")
            hsize=5*cm
            vsize=10*cm
            opd_sc = np.zeros((nbase,astroreduced.time_sc.shape[0]),dtype='d')
            timescale_sc = np.zeros((nbase,astroreduced.time_sc.shape[0]),dtype='d')
   
            for baseline in range(0,nbase):
                phi_met = astroreduced.oi_vis_sc_phi_met[baseline::nbase]
                phi_met = np.outer(phi_met - np.mean(phi_met),1.908287/astroreduced.wave_sc)
                # In pipeline: exp_phase = -phi_ft + phi_met      or     exp_phase = -phi_ft   if --use-met=FALSE
                f1 = (astroreduced.oi_vis_sc_visdata[baseline::nbase,:] * np.exp(1.J*(astroreduced.oi_vis_sc_exp_phase[baseline::nbase,:])))
                opd_sc[baseline,:] = np.angle(np.nanmean(f1[:,minw:maxw],axis=1))# mean over wavelengths
            timescale_sc = astroreduced.time_sc
            
            yminval = np.nanmin(opd_sc)-2
            ymaxval = np.nanmax(opd_sc)+2
            ticky = np.pi/2.
          
            d=()
            data=[]
            for baseline in range(0,nbase):
                data.append(list(zip(timescale_sc[0:nframeplot:resamp], clipdata(opd_sc[baseline,0:nframeplot:resamp],yminval,ymaxval))))
                if(baseline == 3):
                    a = graphoutaxes( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,astroreduced.basenames[baseline]),
                else:
                    a = graphoutnoaxis( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+astroreduced.basenames[baseline]),
                data = []
                d = d + a
            plotmatrix = [list(d[0:int(nbase/2)]),list(d[int(nbase/2):nbase])]
            t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize)
            Story.append(t)
            Story.append(Spacer(1,2*mm))
   
        if hasattr(astroreduced,'oi_vis_sc_visphi_s') & (astroreduced.polarsplit == True):
            plotSubtitle(Story,"If phase is not flat => DITSHIFT. In rad vs. time (s), mean over wavelengths.  Red = S, orange = P.")
    
            hsize=5*cm
            vsize=10*cm
            opd_sc_s = np.zeros((nbase,astroreduced.time_sc.shape[0]),dtype='d')
            opd_sc_p = np.zeros((nbase,astroreduced.time_sc.shape[0]),dtype='d')
            timescale_sc = np.zeros((nbase,astroreduced.time_sc.shape[0]),dtype='d')
            for baseline in range(0,nbase):
                f1 = (astroreduced.oi_vis_sc_visdata_s[baseline::nbase,:] * np.exp(1.J*(astroreduced.oi_vis_sc_exp_phase_s[baseline::nbase,:])))
                opd_sc_s[baseline,:] = np.angle(np.nanmean(f1[:,minw:maxw],axis=1))
                f1 = (astroreduced.oi_vis_sc_visdata_p[baseline::nbase,:] * np.exp(1.J*(astroreduced.oi_vis_sc_exp_phase_p[baseline::nbase,:])))
                opd_sc_p[baseline,:] = np.angle(np.nanmean(f1[:,minw:maxw],axis=1))
    
            timescale_sc = astroreduced.time_sc
            
            yminval = np.nanmin(opd_sc_s)-2
            ymaxval = np.nanmax(opd_sc_s)+2
            ticky = np.pi/2.
            
            d=()
            data=[]
            for baseline in range(0,nbase):
                data.append(list(zip(timescale_sc[0:nframeplot:resamp], clipdata(opd_sc_s[baseline,0:nframeplot:resamp],yminval,ymaxval))))
                data.append(list(zip(timescale_sc[0:nframeplot:resamp], clipdata(opd_sc_p[baseline,0:nframeplot:resamp],yminval,ymaxval))))
                #print data
                if(baseline == 3):
                    a = graphoutaxes( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,astroreduced.basenames[baseline]),
                else:
                    a = graphoutnoaxis( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+astroreduced.basenames[baseline]),
                data = []
                d = d + a
            plotmatrix = [list(d[0:int(nbase/2)]),list(d[int(nbase/2):nbase])]
            t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize)
            Story.append(t)
            Story.append(Spacer(1,2*mm))
    
        Story.append(PageBreak())

    #==============================================================================
    # Differential phase between S and P polarizations (birefringence test) vs. wavelength
    #==============================================================================

    if astroreduced.polarsplit == True:
        plotTitle(Story,"SC phasediff, difference between S and P polar. vs. wavelength")
        plotSubtitle(Story,"In phase (deg) vs. wavelength (microns).")
    
        d=()
        hsize=16*cm
        vsize=3*cm
        yminval=-180
        ymaxval=+180
        ticky = 180
    
        for baseline in range(0,nbase):
            data=[]
            tmp = np.angle(np.mean(astroreduced.oi_vis_sc_visdata_s[baseline::6,:] * np.conj(astroreduced.oi_vis_sc_visdata_p[baseline::6,:]), axis=0))*180./np.pi;
            data.append(list(zip(astroreduced.wave_sc, tmp)))
            a = graphoutaxes( data, min(astroreduced.wave_sc), max(astroreduced.wave_sc),None,yminval,ymaxval,None,
                              hsize,vsize,None,ticky,'SC diff. phase S-P Baseline '+astroreduced.basenames[baseline])
            Story.append(a)
            Story.append(Spacer(1,7*mm))

        Story.append(PageBreak())


    #==============================================================================
    # Metrology TAC algorithm VAMP
    #==============================================================================

    if hasattr(astroreduced,'met_vamp_tel_sc'):
        for tel in range(0,4):
            plotTitle(Story,"OI_VIS_MET - VAMP")
            plotSubtitle(Story,"TAC algo VAMP (Volts) for Tel %s versus time (s)"%(astroreduced.stations[tel],))
    
            hsize=16*cm
            vsize=1.8*cm
            resamp = 100
            xval = astroreduced.oi_vis_met_time[::resamp]
            xminval = min(astroreduced.oi_vis_met_time)
            xmaxval = max(astroreduced.oi_vis_met_time)
            yminval = 0.
            ymaxval = 0.15
            ticky = (ymaxval - yminval)/2.
            
            d = []
            for diode in range(0,4):
                tmp = astroreduced.met_vamp_tel_sc[::resamp,tel*4+diode]
                d.append( graphoutnoxaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))], xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s SC diode %i'%(astroreduced.stations[tel],diode+1)) )
                tmp = astroreduced.met_vamp_tel_ft[::resamp,tel*4+diode]
                d.append( graphoutnoxaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))], xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s FT diode %i'%(astroreduced.stations[tel],diode+1)) )
            yminval = 0.
            ymaxval = 1.
            ticky = (ymaxval - yminval)/2.

            # last diode is FC
            yminval = 0.
            ymaxval = 3.
            ticky = (ymaxval - yminval)/2.
            tmp = astroreduced.met_vamp_fc_sc[::resamp,tel]
            d.append( graphoutnoxaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))],xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s SC -- FC'%(astroreduced.stations[tel],)) )
            tmp = astroreduced.met_vamp_fc_ft[::resamp,tel]
            d.append( graphoutaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))],xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Tel %s FT -- FC'%(astroreduced.stations[tel],)) )
            
            plotmatrix = [ [d[0]],[d[1]],[d[2]],[d[3]],[d[4]],[d[5]],[d[6]],[d[7]],[d[8]],[d[9]] ]
            t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize+4*mm)
            t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')]))
            t.hAlign = 'LEFT'
            Story.append(t)
            
            Story.append(PageBreak())

    #==============================================================================
    # Create PDF report from Story
    #==============================================================================
    print("Create the PDF")

    gravi_visual_class.TITLE = "GRAVITY astroreduced Quality Control Report"
    gravi_visual_class.PAGEINFO = "astroreduced file: "+filename+".fits"
    reportname = filename+"-ASTRORED.pdf"
    
    doc = SimpleDocTemplate(reportname)
    doc.build(Story, onFirstPage=myFirstPage, onLaterPages=myLaterPages)
    print((" "+reportname)) 

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

if __name__ == '__main__':
    filename=''
    # filename = '../GRAVITY.2015-11-16T00-06-33_0001'
    # filename = 'GRAVITY.2016-02-28T06-03-18_astroreduced'
    filename='../../GRAVITY.2016-03-25T00-08-57_astroreduced'

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

    astroreduced = gravi_visual_class.Astroreduced(filename+'.fits')

    produce_astroreduced_report(astroreduced,filename)