#! /usr/bin/env python
# -*- coding: iso-8859-15 -*-
"""
Created on Mon May 25 13:51:41 2015

@author: kervella
"""

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

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

import sys
import datetime
import numpy as np
from . import gravi_visual_class


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

    from reportlab.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
    
    
    #==============================================================================
    # Global parameters
    #==============================================================================
    
    sobjname =  get_key_withdefault(oifits.header,'HIERARCH ESO INS SOBJ NAME','')
    objname  =  get_key_withdefault(oifits.header,'HIERARCH ESO FT ROBJ NAME','')
    
    triplet_list = [str(oifits.t3_sc_staindex[0,:]),str(oifits.t3_sc_staindex[1,:]),\
                    str(oifits.t3_sc_staindex[2,:]),str(oifits.t3_sc_staindex[3,:])]

    minwave = np.nanmin(oifits.wave_sc) # minimum wavelength for FT plots
    maxwave = np.nanmax(oifits.wave_sc) # maximum wavelength for FT plots
    minwaveSC = float(minwaveSC) # different from FT if specified by user
    maxwaveSC = float(maxwaveSC)

    print("Wavelength limits for SC plots (microns):", minwaveSC, maxwaveSC)
    
    grandmaxvis = 1.4 # maximum of Y axis of VISAMP / VIS2 plots


    #==============================================================================
    # Selection of SC wavelength range for plots
    #==============================================================================

    waveok = np.squeeze(np.where((oifits.wave_sc>=minwaveSC) * (oifits.wave_sc<=maxwaveSC)))

    oifits.wave_sc = oifits.wave_sc[waveok]
    oifits.nwave_sc = len(waveok)
    oifits.spfreq_sc = oifits.spfreq_sc[:,waveok]
    oifits.t3_spfreq_sc = oifits.t3_spfreq_sc[:,waveok]
    
    if oifits.polarsplit == False:
        oifits.oi_flux_sc = oifits.oi_flux_sc[:,:,waveok]
        oifits.oi_vis2_sc_vis2data = oifits.oi_vis2_sc_vis2data[:,waveok]
        oifits.oi_vis2_sc_vis2err = oifits.oi_vis2_sc_vis2err[:,waveok]
        oifits.t3phi_sc = oifits.t3phi_sc[:,waveok]
        oifits.t3phierr_sc = oifits.t3phierr_sc[:,waveok]
        oifits.oi_vis_sc_visamp=oifits.oi_vis_sc_visamp[:,waveok]
        oifits.oi_vis_sc_visamperr=oifits.oi_vis_sc_visamperr[:,waveok]
        oifits.oi_vis_sc_visphi=oifits.oi_vis_sc_visphi[:,waveok]
        oifits.oi_vis_sc_visphierr=oifits.oi_vis_sc_visphierr[:,waveok]
    else:
        oifits.oi_flux_sc_s = oifits.oi_flux_sc_s[:,:,waveok]
        oifits.oi_flux_sc_p = oifits.oi_flux_sc_p[:,:,waveok]
        oifits.oi_vis2_sc_vis2data_s = oifits.oi_vis2_sc_vis2data_s[:,waveok]
        oifits.oi_vis2_sc_vis2err_s = oifits.oi_vis2_sc_vis2err_s[:,waveok]
        oifits.oi_vis2_sc_vis2data_p = oifits.oi_vis2_sc_vis2data_p[:,waveok]
        oifits.oi_vis2_sc_vis2err_p = oifits.oi_vis2_sc_vis2err_p[:,waveok]
        oifits.t3phi_sc_s = oifits.t3phi_sc_s[:,waveok]
        oifits.t3phierr_sc_s = oifits.t3phierr_sc_s[:,waveok]
        oifits.t3phi_sc_p = oifits.t3phi_sc_p[:,waveok]
        oifits.t3phierr_sc_p = oifits.t3phierr_sc_p[:,waveok]
        oifits.oi_vis_sc_visamp_s=oifits.oi_vis_sc_visamp_s[:,waveok]
        oifits.oi_vis_sc_visamp_p=oifits.oi_vis_sc_visamp_p[:,waveok]
        oifits.oi_vis_sc_visamperr_s=oifits.oi_vis_sc_visamperr_s[:,waveok]
        oifits.oi_vis_sc_visamperr_p=oifits.oi_vis_sc_visamperr_p[:,waveok]
        oifits.oi_vis_sc_visphi_s=oifits.oi_vis_sc_visphi_s[:,waveok]
        oifits.oi_vis_sc_visphi_p=oifits.oi_vis_sc_visphi_p[:,waveok]
        oifits.oi_vis_sc_visphierr_s=oifits.oi_vis_sc_visphierr_s[:,waveok]
        oifits.oi_vis_sc_visphierr_p=oifits.oi_vis_sc_visphierr_p[:,waveok]

    #==============================================================================
    # Resampling factor to avoid too large PDF files
    #==============================================================================
    
    resamp = 3 if oifits.nwave_sc > 500 else 1

    #==============================================================================
    # Matplotlib experimental plots
    #==============================================================================

    def add_subplot_axes(ax,rect,axisbg='w'):
        fig = plt.gcf()

        inax_position  = ax.transAxes.transform(rect[0:2])
        transFigure = fig.transFigure.inverted()
        infig_position = transFigure.transform(inax_position)    
        x = infig_position[0]
        y = infig_position[1]

        box = ax.get_position()
        width = box.width
        height = box.height
        width *= rect[2]
        height *= rect[3]
    
        subax = fig.add_axes([x,y,width,height],axisbg=axisbg)
        x_labelsize = subax.get_xticklabels()[0].get_size()
        y_labelsize = subax.get_yticklabels()[0].get_size()
        x_labelsize *= rect[2]**0.5
        y_labelsize *= rect[3]**0.5
        subax.xaxis.set_tick_params(labelsize=x_labelsize)
        subax.yaxis.set_tick_params(labelsize=y_labelsize)
        return subax

    def plot_Vis2_t3phi_UV_SC(oifits):
        ucoord = oifits.ucoord_sc
        vcoord = oifits.vcoord_sc
        hfont = {'fontname':'Helvetica'}
        resamp = 3 if oifits.nwave_sc > 500 else 1
#        l0 = 0 # first wavelength channel to display (to skip bad channels)
        
        imgdata = io.StringIO()
        plt.close('all')
        figwidth = 6
        figheight = 9
        plt.figure(figsize=(figwidth,figheight))
        colorlist = ['red', 'orange', 'blue', 'green', 'cyan', 'purple']
        gs = plt.GridSpec(2, 1, height_ratios=[1,1])
        plt.subplots_adjust(hspace=0)
        ax0 = plt.subplot(gs[0])
        for base in range(0,nbase):
            if oifits.polarsplit == False:
                ax0.errorbar(oifits.spfreq_sc[base,::resamp], oifits.oi_vis2_sc_vis2data[base,::resamp],\
                        yerr=oifits.oi_vis2_sc_vis2err[base,::resamp],\
                        fmt='o', color=colorlist[base], markersize=3, markeredgecolor='gray',\
                        ecolor=colorlist[base], capsize=1, capthick=1)
                ax0.set_xlim([0., 1.3*np.max(oifits.spfreq_sc)])
                ax0.set_ylim([0., 1.3*np.nanmin([1.,np.max(oifits.oi_vis2_sc_vis2data)])])
                #ax0.set_ylim([0., 0.7])
            else:
                ax0.errorbar(oifits.spfreq_sc[base,::resamp], oifits.oi_vis2_sc_vis2data_s[base,::resamp],\
                        yerr=oifits.oi_vis2_sc_vis2err_s[base,::resamp],\
                        fmt='o', color=colorlist[base], markersize=3, markeredgecolor='gray',\
                        ecolor=colorlist[base], capsize=1, capthick=1)
                ax0.errorbar(oifits.spfreq_sc[base,::resamp], oifits.oi_vis2_sc_vis2data_p[base,::resamp],\
                        yerr=oifits.oi_vis2_sc_vis2err_p[base,::resamp],\
                        fmt='o', color=colorlist[base], markersize=3, markeredgecolor='gray',\
                        ecolor=colorlist[base], capsize=1, capthick=1)
                ax0.set_xlim([0., 1.3*np.max(oifits.spfreq_sc)])
                ax0.set_ylim([0., 1.3*np.nanmin([1.,np.min(np.nanmax([oifits.oi_vis2_sc_vis2data_s,oifits.oi_vis2_sc_vis2data_p]))])])
        ax0.xaxis.label.set_fontsize(10)
        ax0.yaxis.label.set_fontsize(10)
    #        plotx = np.arange(0,1.5E8,0.01*np.mean(spfreqok))
    #        ax0.plot(plotx/cycarcsec,UDvis2(fit_ud[0],plotx), color='r', linewidth=1)
    #        ax0.plot(plotx/cycarcsec,UDvis2(fit_ud[0]+errtot,plotx), color='r', linestyle='--', linewidth=0.5)
    #        ax0.plot(plotx/cycarcsec,UDvis2(fit_ud[0]-errtot,plotx), color='r', linestyle='--', linewidth=0.5)
        #ax0.set_xticklabels([])
        ax0.set_ylabel('Squared visibility', **hfont)
        ax0.set_xlabel('Spatial frequency (1/arcsec)', **hfont)
        ax0.tick_params(axis='both',labelsize=10)
        ax0.grid(True)
        #ax0.text(15, 1.3, titletext, rotation='horizontal', color='k', size=14, **hfont)
        #ax0.text(15, 1.25, r"$\theta_{\mathrm{UD}}$ = %(val).3f $\pm$ %(err).3f  $\pm$ %(errsys).3f mas"%{"val":fit_ud[0], "err":fit_ud[1], "errsys":errsys}, rotation='horizontal', color='k', size=16, **hfont)
        rect = [0.77,0.82,0.25,0.32]
        ax0b = add_subplot_axes(ax0, rect)
        maxbase = 1.3*np.max(oifits.projbase_sc)
        ax0b.set_xlim([-maxbase, maxbase])
        ax0b.set_ylim([-maxbase, maxbase])
        for base in range(0,nbase):
            ax0b.plot(ucoord[base],vcoord[base], linestyle='None', marker='o', markersize=4, color=colorlist[base])
            ax0b.plot(-ucoord[base],-vcoord[base], linestyle='None', marker='o', markersize=4, color=colorlist[base])
        ax0b.tick_params(axis='both',labelsize=6)
    
        ax1 = plt.subplot(gs[1])
        for triplet in range(0,4):
            if oifits.polarsplit == False:
                ax1.errorbar(oifits.t3_spfreq_sc[triplet,::resamp], oifits.t3phi_sc[triplet,::resamp]*180./np.pi,\
                yerr=oifits.t3phierr_sc[triplet,::resamp]*180./np.pi,\
                fmt='o', color='red', markersize=3, markeredgecolor='gray', ecolor='red', capsize=1, capthick=1)
            else:
                ax1.errorbar(oifits.t3_spfreq_sc[triplet,::resamp], oifits.t3phi_sc_s[triplet,::resamp]*180./np.pi,\
                yerr=oifits.t3phierr_sc_s[triplet,::resamp]*180./np.pi,\
                fmt='o', color='red', markersize=3, markeredgecolor='gray', ecolor='red', capsize=1, capthick=1)

                ax1.errorbar(oifits.t3_spfreq_sc[triplet,::resamp], oifits.t3phi_sc_p[triplet,::resamp]*180./np.pi,\
                yerr=oifits.t3phierr_sc_p[triplet,::resamp]*180./np.pi,\
                fmt='o', color='orange', markersize=3, markeredgecolor='gray', ecolor='orange', capsize=1, capthick=1)
                
        ax1.set_ylabel('Closure phase (degrees)', **hfont)
        ax1.set_xlabel('Spatial frequency (1/arcsec)', **hfont)
        ax1.xaxis.label.set_fontsize(10)
        ax1.yaxis.label.set_fontsize(10)
        ax1.tick_params(axis='both',labelsize=10)
        if oifits.polarsplit == False:
            t3phi_max = np.nanpercentile(np.abs(oifits.t3phi_sc),99)*180./np.pi
            #t3phi_max = 60.
        else:
            t3phi_max = np.nanpercentile(np.abs([oifits.t3phi_sc_p,oifits.t3phi_sc_p]),99)*180./np.pi
            #t3phi_max = 60.
            
        ax1.set_xlim([0., 1.3*np.max(oifits.spfreq_sc)])
        ax1.set_ylim([-t3phi_max-10, t3phi_max+10])
        ax1.grid(True)
        #ax1.text(15, 1.3, titletext, rotation='horizontal', color='k', size=20, **hfont)
        plt.tight_layout()
        
        return imgdata
    
    def plot_Vis2_t3phi_UV_FT(oifits):
        ucoord = oifits.ucoord_ft
        vcoord = oifits.vcoord_ft
        hfont = {'fontname':'Helvetica'}
    
        imgdata = io.StringIO()
        plt.close('all')
        figwidth = 6
        figheight = 9
        plt.figure(figsize=(figwidth,figheight))
        colorlist = ['red', 'orange', 'blue', 'green', 'cyan', 'purple']
        gs = plt.GridSpec(2, 1, height_ratios=[1,1])
        plt.subplots_adjust(hspace=0)
        ax0 = plt.subplot(gs[0])
        for base in range(0,nbase):
            if oifits.polarsplit == False:
                ax0.errorbar(oifits.spfreq_ft[base,:], oifits.oi_vis2_ft_vis2data[base,:],\
                        yerr=oifits.oi_vis2_ft_vis2err[base,:],\
                        fmt='o', color=colorlist[base], markersize=3, linestyle='-', markeredgecolor='gray',\
                        ecolor=colorlist[base], capsize=1, capthick=1)
                ax0.set_xlim([0., 1.3*np.max(oifits.spfreq_sc)]) # same scale as SC
                ax0.set_ylim([0., 1.3*np.nanmin([1.,np.max(oifits.oi_vis2_ft_vis2data)])])
            else:
                ax0.errorbar(oifits.spfreq_ft[base,:], oifits.oi_vis2_ft_vis2data_s[base,:],\
                        yerr=oifits.oi_vis2_ft_vis2err_s[base,:],\
                        fmt='o', color=colorlist[base], markersize=3, linestyle='-', markeredgecolor='gray',\
                        ecolor=colorlist[base], capsize=1, capthick=1)
                ax0.errorbar(oifits.spfreq_ft[base,:], oifits.oi_vis2_ft_vis2data_p[base,:],\
                        yerr=oifits.oi_vis2_ft_vis2err_p[base,:],\
                        fmt='o', color=colorlist[base], markersize=3, linestyle='-', markeredgecolor='gray',\
                        ecolor=colorlist[base], capsize=1, capthick=1)
                ax0.set_xlim([0., 1.3*np.max(oifits.spfreq_sc)])
                ax0.set_ylim([0., 1.3*np.nanmin([1.,np.min(np.nanmax([oifits.oi_vis2_ft_vis2data_s,oifits.oi_vis2_ft_vis2data_p]))])])
        ax0.xaxis.label.set_fontsize(10)
        ax0.yaxis.label.set_fontsize(10)
        ax0.set_ylabel('Squared visibility', **hfont)
        ax0.set_xlabel('Spatial frequency (1/arcsec)', **hfont)
        ax0.tick_params(axis='both',labelsize=10)
        ax0.grid(True)
        #ax0.text(15, 1.3, titletext, rotation='horizontal', color='k', size=14, **hfont)
        #ax0.text(15, 1.25, r"$\theta_{\mathrm{UD}}$ = %(val).3f $\pm$ %(err).3f  $\pm$ %(errsys).3f mas"%{"val":fit_ud[0], "err":fit_ud[1], "errsys":errsys}, rotation='horizontal', color='k', size=16, **hfont)

        rect = [0.77,0.82,0.25,0.32]
        ax0b = add_subplot_axes(ax0, rect)
        maxbase = 1.3*np.max(oifits.projbase_ft)
        ax0b.set_xlim([-maxbase, maxbase])
        ax0b.set_ylim([-maxbase, maxbase])
        for base in range(0,nbase):
            ax0b.plot(ucoord[base],vcoord[base], linestyle='None', marker='o', markersize=4, color=colorlist[base])
            ax0b.plot(-ucoord[base],-vcoord[base], linestyle='None', marker='o', markersize=4, color=colorlist[base])
        ax0b.tick_params(axis='both',labelsize=6)
    
        ax1 = plt.subplot(gs[1])
        for triplet in range(0,4):
            if oifits.polarsplit == False:
                ax1.errorbar(oifits.t3_spfreq_ft[triplet,:], oifits.t3phi_ft[triplet,:]*180./np.pi,\
                yerr=oifits.t3phierr_ft[triplet,:]*180./np.pi,\
                fmt='o', color='red', markersize=3, linestyle='-', markeredgecolor='gray', ecolor='red', capsize=1, capthick=1)
            else:
                ax1.errorbar(oifits.t3_spfreq_ft[triplet,:], oifits.t3phi_ft_s[triplet,:]*180./np.pi,\
                yerr=oifits.t3phierr_ft_s[triplet,:]*180./np.pi,\
                fmt='o', color='red', markersize=3, linestyle='-', markeredgecolor='gray', ecolor='red', capsize=1, capthick=1)

                ax1.errorbar(oifits.t3_spfreq_ft[triplet,:], oifits.t3phi_ft_p[triplet,:]*180./np.pi,\
                yerr=oifits.t3phierr_ft_p[triplet,:]*180./np.pi,\
                fmt='o', color='orange', markersize=3, linestyle='-', markeredgecolor='gray', ecolor='orange', capsize=1, capthick=1)
                
        ax1.set_ylabel('Closure phase (degrees)', **hfont)
        ax1.set_xlabel('Spatial frequency (1/arcsec)', **hfont)
        ax1.xaxis.label.set_fontsize(10)
        ax1.yaxis.label.set_fontsize(10)
        ax1.tick_params(axis='both',labelsize=10)
        if oifits.polarsplit == False:
            t3phi_max = np.nanpercentile(np.abs(oifits.t3phi_ft),99)*180./np.pi
        else:
            t3phi_max = np.nanpercentile(np.abs([oifits.t3phi_ft_p,oifits.t3phi_ft_p]),99)*180./np.pi
            
        ax1.set_xlim([0., 1.3*np.max(oifits.spfreq_sc)])
        ax1.set_ylim([-t3phi_max-10, t3phi_max+10])
        ax1.grid(True)
        #ax1.text(15, 1.3, titletext, rotation='horizontal', color='k', size=20, **hfont)
        plt.tight_layout()
        
        return imgdata



    #==============================================================================
    # Start Story
    #==============================================================================

    Story = [Spacer(1,1*mm)]
    plotSummary (Story, filename, oifits, onTarget=True)

    #==============================================================================
    # Overview
    #==============================================================================
    
    plotTitle(Story,"Overview of processing results")
    plotSubtitle(Story,"(u,v) plane in meters, closure phase in degrees, spatial frequencies in megalambda.")
        
    # Plot of the (u,v) plane
    hsize = 6*cm
    vsize = 6*cm
    yminval = -150
    ymaxval = +150
    ticky = 50
    d = []
    data=[]
    for base in range(0,nbase):
        uvcoord = tuple([oifits.ucoord_sc[base], oifits.vcoord_sc[base]]),\
            tuple([-1.0*oifits.ucoord_sc[base], -1.0*oifits.vcoord_sc[base]])
        data.append( uvcoord )
    a = graphscatteraxes( data,yminval,ymaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,ticky,ticky,5,'(u,v) meters')
    d.append(a)

    # Plot of the SC squared visibility vs. spatial frequency
    spfreq_sc = np.zeros((nbase,oifits.nwave_sc))
    for base in range(0,nbase):
        spfreq_sc[base,:] = (oifits.projbase_sc[base]/(oifits.wave_sc/1.e6))/1.e6 # Spatial frequency in mega lambda

    xminval = int(np.min(spfreq_sc))-5
    xmaxval = int(np.max(spfreq_sc))+5
    yminval = 0.
    ymaxval = grandmaxvis
    tickx = (xmaxval-xminval)/10
    ticky = (ymaxval-yminval)/10.

    data=[]
    for base in range(0,nbase):
        if oifits.polarsplit == True:
            mean_vis2 = np.mean([oifits.oi_vis2_sc_vis2data_s[base,::resamp],oifits.oi_vis2_sc_vis2data_p[base,::resamp]],axis=0)
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(mean_vis2,yminval,ymaxval)) ))
        else:
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis2_sc_vis2data[base,::resamp],yminval,ymaxval)) ))
        a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,2,'SC VIS2')
    d.append(a)

    # Plot of the FT squared visibility vs. spatial frequency
    spfreq_ft = np.zeros((nbase,oifits.nwave_ft))
    for base in range(0,nbase):
        spfreq_ft[base,:] = (oifits.projbase_ft[base]/(oifits.wave_ft/1.e6))/1.e6 # Spatial frequency in mega lambda

    xminval = int(np.min(spfreq_sc))-5
    xmaxval = int(np.max(spfreq_sc))+5
    yminval = 0.
    ymaxval = grandmaxvis
    tickx = (xmaxval-xminval)/10
    ticky = (ymaxval-yminval)/10.

    data=[]
    for base in range(0,nbase):
        if oifits.polarsplit == True:
            mean_vis2 = np.mean([oifits.oi_vis2_ft_vis2data_s[base,:],oifits.oi_vis2_ft_vis2data_p[base,:]],axis=0)
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(mean_vis2,yminval,ymaxval)) ))
        else:
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis2_ft_vis2data[base,:],yminval,ymaxval)) ))
        a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VIS2')
    d.append(a)


    # Plot of the SC and FT closure phases  vs. spatial frequency
    nt3 = 4
    spfreq_sc_t3phi = np.zeros((nt3,oifits.nwave_sc))
    spfreq_ft_t3phi = np.zeros((nt3,oifits.nwave_ft))
    for triangle in range(0,nt3):
        spfreq_sc_t3phi[triangle,:] = (oifits.t3_baseavg[triangle]/(oifits.wave_sc/1.e6))/1.e6 # Average spatial frequency of triangle in mega lambda
        spfreq_ft_t3phi[triangle,:] = (oifits.t3_baseavg[triangle]/(oifits.wave_ft/1.e6))/1.e6 
        
    xminval = int(np.min(spfreq_sc_t3phi))-5
    xmaxval = int(np.max(spfreq_sc_t3phi))+5
    if oifits.polarsplit == True:
        ymaxval = int(np.abs(np.max([np.max(oifits.t3phi_sc_s),np.max(oifits.t3phi_sc_p),
                                        np.max(oifits.t3phi_ft_s),np.max(oifits.t3phi_ft_p)])*180./np.pi))+5
    else:
        ymaxval = int(np.abs(np.max([np.max(oifits.t3phi_sc),np.max(oifits.t3phi_ft)]))*180./np.pi)+5
    yminval = -ymaxval
    tickx = (xmaxval-xminval)/10
    ticky = (ymaxval-yminval)/10.
    
    data=[]
    for triangle in range(0,nt3):
        if oifits.polarsplit == True:
            mean_t3phi_sc = np.mean([oifits.t3phi_sc_s[triangle,::resamp]*180./np.pi,oifits.t3phi_sc_p[triangle,::resamp]*180./np.pi],axis=0)
            data.append( tuple(zip(spfreq_sc_t3phi[triangle,::resamp], clipdata(mean_t3phi_sc,yminval,ymaxval)) ))
            mean_t3phi_ft = np.mean([oifits.t3phi_ft_s[triangle,:]*180./np.pi,oifits.t3phi_ft_p[triangle,:]*180./np.pi],axis=0)
            data.append( tuple(zip(spfreq_ft_t3phi[triangle,:], clipdata(mean_t3phi_ft,yminval,ymaxval)) ))
        else:
            data.append( tuple(zip(spfreq_sc_t3phi[triangle,::resamp], clipdata(oifits.t3phi_sc[triangle,::resamp]*180./np.pi,yminval,ymaxval)) ))
            data.append( tuple(zip(spfreq_ft_t3phi[triangle,:], clipdata(oifits.t3phi_ft[triangle,:]*180./np.pi,yminval,ymaxval)) ))
        a = graphscatteraxes(data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,2,'SC and FT T3PHI')
    d.append(a)

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

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

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

    #==============================================================================
    # SC
    #==============================================================================

    plotTitle(Story,"SC VIS2 and T3PHI of "+sobjname)
    Story.append(Spacer(1,5*mm))
    
    imgdata = plot_Vis2_t3phi_UV_SC(oifits)
    widthfig = 15*cm
    heightfig= 22*cm
    plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
    pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
    plt.close('all')
    
    Story.append(pi1)
    Story.append(PageBreak())
    
    #==============================================================================
    # FT
    #==============================================================================
    
    plotTitle(Story,"FT VIS2 and T3PHI of "+objname)
    Story.append(Spacer(1,5*mm))
    
    imgdata = plot_Vis2_t3phi_UV_FT(oifits)
    widthfig = 15*cm
    heightfig= 22*cm
    plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
    pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
    plt.close('all')
    
    Story.append(pi1)
    Story.append(PageBreak())
    
    #==============================================================================
    # SC Flux signals vs wavelength
    #==============================================================================

    plotTitle(Story,"Average spectrum for SC vs. wavelength (OI_FLUX)")
    plotSubtitle(Story,"In photo-e-/DIT vs. wavelength (microns). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.")

    if oifits.polarsplit == False:
        nmeasure = oifits.oi_flux_sc_time.shape[0]
        avgflux = np.zeros((4,oifits.nwave_sc),dtype='d')
        for tel in range(0,4):
            avgflux[tel,:]=np.nanmean(oifits.oi_flux_sc[:,tel,:],axis=0) # average flux over time
    
        yminval = 0
        ymaxval = np.nanpercentile(avgflux,98)*1.1
        tickx = (maxwaveSC-minwaveSC)/10.
        ticky = (ymaxval-yminval)/10. # in photo-e-
        hsize = 16*cm
        vsize = 9*cm
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(oifits.wave_sc[::resamp], clipdata(avgflux[tel,::resamp],yminval,ymaxval)) ))
        a = graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2e',hsize,vsize,tickx,ticky,'Combined-polarization')
        Story.append(a)
        
    if oifits.polarsplit == True:
        nmeasure = oifits.oi_flux_sc_time.shape[0]
        avgfluxs = np.zeros((4,oifits.nwave_sc),dtype='d')
        avgfluxp = np.zeros((4,oifits.nwave_sc),dtype='d')
        for tel in range(0,4):
            avgfluxs[tel,:]=np.nanmean(oifits.oi_flux_sc_s[:,tel,:],axis=0) # average flux over time
            avgfluxp[tel,:]=np.nanmean(oifits.oi_flux_sc_p[:,tel,:],axis=0) # average flux over time
    
        yminval = 0
        ymaxval = np.nanpercentile(np.array([avgfluxs,avgfluxp]),98)*1.1
        tickx = (maxwaveSC-minwaveSC)/10.
        ticky = (ymaxval-yminval)/10 # in photo-e-
        hsize = 16*cm
        vsize = 4.5*cm
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(oifits.wave_sc[::resamp], clipdata(avgfluxs[tel,::resamp],yminval,ymaxval)) ))
        a = graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2e',hsize,vsize,tickx,ticky,'S-polarization')
        Story.append(a)
        
        Story.append(Spacer(1,7*mm))
    
        yminval = 0
        ymaxval = np.nanpercentile(np.array([avgfluxs,avgfluxp]),98)*1.1
        ticky = (ymaxval-yminval)/10 # in photo-e-
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(oifits.wave_sc[::resamp], clipdata(avgfluxp[tel,::resamp],yminval,ymaxval)) ))
        a = graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'P-polarization')
        Story.append(a)

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

    #==============================================================================
    # FT Flux signals vs wavelength
    #==============================================================================
    
    plotTitle(Story,"Average spectrum for FT vs. wavelength (OI_FLUX)")
    plotSubtitle(Story,"In photo-e-/DIT vs. wavelength (microns). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.")

    if oifits.polarsplit == False:
        #nmeasure = oifits.oi_flux_ft_time.shape[0]
        avgflux = np.zeros((4,oifits.nwave_ft),dtype='d')
        for tel in range(0,4):
            avgflux[tel,:]=np.nanmean(oifits.oi_flux_ft[:,tel,:],axis=0) # average flux over time
    
        hsize = 16*cm
        vsize = 9*cm
        yminval = 0
        ymaxval = np.max(avgflux)*1.2
        tickx = (maxwave-minwave)/10.
        ticky = (ymaxval-yminval)/10 # in photo-e-
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(oifits.wave_ft, clipdata(avgflux[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes(data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2e',hsize,vsize,tickx,ticky,'Combined-polarization')
        Story.append(a)

    if oifits.polarsplit == True:
        #nmeasure = oifits.oi_flux_ft_time_s.shape[0]
        avgfluxs = np.zeros((4,oifits.nwave_ft),dtype='d')
        avgfluxp = np.zeros((4,oifits.nwave_ft),dtype='d')
        for tel in range(0,4):
            avgfluxs[tel,:]=np.nanmean(oifits.oi_flux_ft_s[:,tel,:],axis=0) # average flux over time
            avgfluxp[tel,:]=np.nanmean(oifits.oi_flux_ft_p[:,tel,:],axis=0) # average flux over time
    
        hsize = 16*cm
        vsize = 4.5*cm
        yminval = 0
        ymaxval = np.max(avgfluxs)*1.2
        tickx = (maxwave-minwave)/10.
        ticky = (ymaxval-yminval)/10 # in photo-e-
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(oifits.wave_ft, clipdata(avgfluxs[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes(data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2e',hsize,vsize,tickx,ticky,'S-polarization')
        Story.append(a)
        
        Story.append(Spacer(1,7*mm))
    
        yminval = 0
        ymaxval = np.max(avgfluxp)*1.2
        ticky = (ymaxval-yminval)/10 # in photo-e-
        data = []
        for tel in range(0,4):
            data.append( tuple(zip(oifits.wave_ft, clipdata(avgfluxp[tel,:],yminval,ymaxval)) ))
        a = graphoutaxes(data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2e',hsize,vsize,tickx,ticky,'P-polarization')
        Story.append(a)

    Story.append(PageBreak())
    
    #==============================================================================
    # SC VISAMP vs spatial frequency
    #==============================================================================

    hsize = 16*cm
    vsize = 9*cm
    yminval = 0.0
    ymaxval = grandmaxvis
    ticky = 0.2

    plotTitle(Story,"SC visibility amplitude vs. spatial frequency (VISAMP SC)")
    plotSubtitle(Story,"Baselines: 12 = red, 13 = orange, 14 = blue, 23 = green, 24 = cyan, 34 = purple. X axis = megalambda")
    xminval = int(np.min(spfreq_sc))-5
    xmaxval = int(np.max(spfreq_sc))+5
    tickx = (xmaxval-xminval)/10
    yminval = 0.
    if 'LOW' in oifits.header['HIERARCH ESO INS SPEC RES']:
        ticksize = 4
    else:
        ticksize = 2

    if oifits.polarsplit==True:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis_sc_visamp_s,oifits.oi_vis_sc_visamp_p],98) + 0.1,grandmaxvis])
    else:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis_sc_visamp],98) + 0.1,grandmaxvis])


    if oifits.polarsplit == True:
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis_sc_visamp_s[base,::resamp],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,ticksize,'SC VISAMP S')
        plotSubtitle(Story,"Polarization S")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis_sc_visamp_p[base,::resamp],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,ticksize,'SC VISAMP P')
        plotSubtitle(Story,"Polarization P")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        Story.append(PageBreak())
        
    else: # combined polarization
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis_sc_visamp[base,::resamp],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i', yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,ticksize,'SC VISAMP')
        Story.append(a)
        Story.append(Spacer(1,7*mm))


    #==============================================================================
    # SC VIS2DATA vs spatial frequency
    #==============================================================================

    hsize = 16*cm
    vsize = 9*cm

    plotTitle(Story,"SC squared visibility vs. spatial frequency (VIS2DATA SC)")
    plotSubtitle(Story,"Baselines: 12 = red, 13 = orange, 14 = blue, 23 = green, 24 = cyan, 34 = purple. X axis = megalambda")
    
    xminval = int(np.min(spfreq_sc))-5
    xmaxval = int(np.max(spfreq_sc))+5
    tickx = (xmaxval-xminval)/10
    yminval = 0.

    if oifits.polarsplit==True:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis2_sc_vis2data_s,oifits.oi_vis2_sc_vis2data_p],98) + 0.1,grandmaxvis])
    else:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis2_sc_vis2data],98) + 0.1,grandmaxvis])

    if oifits.polarsplit == True:
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis2_sc_vis2data_s[base,::resamp],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,ticksize,'SC VIS2DATA S')
        plotSubtitle(Story,"Polarization S")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis2_sc_vis2data_p[base,::resamp],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,ticksize,'SC VIS2DATA P')
        plotSubtitle(Story,"Polarization P")
        Story.append(a)

        Story.append(PageBreak())
        
    else: # combined polarization
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_sc[base,::resamp], clipdata(oifits.oi_vis2_sc_vis2data[base,::resamp],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,ticksize,'SC VIS2DATA')
        Story.append(a)
        
        Story.append(PageBreak())

    #==============================================================================
    # FT VISAMP vs spatial frequency
    #==============================================================================

    hsize = 16*cm
    vsize = 9*cm

    plotTitle(Story,"FT visibility amplitude vs. spatial frequency (VISAMP FT)")
    plotSubtitle(Story,"Baselines: 12 = red, 13 = orange, 14 = blue, 23 = green, 24 = cyan, 34 = purple. X axis = megalambda")

    if oifits.polarsplit==True:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis_ft_visamp_s,oifits.oi_vis_ft_visamp_p],98) + 0.1,grandmaxvis])
    else:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis_ft_visamp],98) + 0.1,grandmaxvis])
        
    if oifits.polarsplit == True:
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis_ft_visamp_s[base,:],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VISAMP S')
        plotSubtitle(Story,"Polarization S")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis_ft_visamp_p[base,:],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VISAMP P')
        plotSubtitle(Story,"Polarization P")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        Story.append(PageBreak())
        
    else: # combined polarization
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis_ft_visamp[base,:],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VISAMP')
        Story.append(a)
        Story.append(Spacer(1,7*mm))
        
    #==============================================================================
    # FT VIS2DATA vs spatial frequency
    #==============================================================================

    hsize = 16*cm
    vsize = 9*cm

    plotTitle(Story,"FT squared visibility vs. spatial frequency (VIS2DATA FT)")
    plotSubtitle(Story,"Baselines: 12 = red, 13 = orange, 14 = blue, 23 = green, 24 = cyan, 34 = purple. X axis = megalambda")
#    xminval = int(np.min(spfreq_ft))-5
#    xmaxval = int(np.max(spfreq_ft))+5
#    tickx = (xmaxval-xminval)/10
#    yminval = 0.

    if oifits.polarsplit==True:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis2_ft_vis2data_s,oifits.oi_vis2_ft_vis2data_p],98) + 0.1,grandmaxvis])
    else:
        ymaxval = np.nanmin([np.nanpercentile([oifits.oi_vis2_ft_vis2data],98) + 0.1,grandmaxvis])

    if oifits.polarsplit == True:
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis2_ft_vis2data_s[base,:],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VIS2DATA S')
        plotSubtitle(Story,"Polarization S")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis2_ft_vis2data_p[base,:],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VIS2DATA P')
        plotSubtitle(Story,"Polarization P")
        Story.append(a)
        Story.append(Spacer(1,7*mm))

        Story.append(PageBreak())
        
    else: # combined polarization
        data=[]
        for base in range(0,nbase):
            data.append( tuple(zip(spfreq_ft[base,:], clipdata(oifits.oi_vis2_ft_vis2data[base,:],yminval,ymaxval)) ))
            a = graphscatteraxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,4,'FT VIS2DATA')
        Story.append(a)
        Story.append(PageBreak())

    #==============================================================================
    # Time averaged spectrum of SC VISAMP
    #==============================================================================
    hsize = 16*cm
    vsize = 6.5*cm
    yminval = 0.0
    ymaxval = grandmaxvis
    ticky = 0.2
    tickx = 0.1

    plotTitle(Story,"Time averaged SC visibility amp. vs. wavelength (VISAMP SC)")
    plotSubtitle(Story,"Visibility vs. wavelength (microns).")
    
    d=()
    for baseline in range(0,3):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visamp[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visamp_s[baseline::6,::resamp], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visamp_p[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        a = graphoutaxes( data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'SC VISAMP '+base_list[baseline])
        #a = graphoutaxes( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline])
        Story.append(a)
        Story.append(Spacer(1,7*mm))
    Story.append(PageBreak())


    plotTitle(Story,"Time averaged SC visibility amp. vs. wavelength (VISAMP SC)")
    plotSubtitle(Story,"Visibility vs. wavelength (microns).")

    for baseline in range(3,6):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visamp[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visamp_s[baseline::6,::resamp], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visamp_p[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        a = graphoutaxes( data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'SC VISAMP '+base_list[baseline])
        #a = graphoutaxes( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline])
        Story.append(a)
        Story.append(Spacer(1,7*mm))

    Story.append(PageBreak())

    #==============================================================================
    # Time averaged spectrum of FT VISAMP
    #==============================================================================

    hsize = 8*cm
    vsize = 6.5*cm
    yminval = 0.0
    ymaxval = grandmaxvis
    ticky = 0.2

    plotTitle(Story,"FT visibility amp. vs. wavelength (VISAMP SC)")
    plotSubtitle(Story,"Visibility vs. wavelength (microns).")

    d=()
    for baseline in range(0,nbase):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis_ft_visamp[baseline::6,:], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis_ft_visamp_s[baseline::6,:], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis_ft_visamp_p[baseline::6,:], axis=0),yminval,ymaxval))))
        if(baseline == 4 or baseline == 2 or baseline == 0):
            d = d + (graphoutaxes( data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]),)
        else:
            d = d + (graphoutnoaxis( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]),)

    Story.append(Table( [[d[0],d[1]],[d[2],d[3]],[d[4],d[5]]], colWidths=hsize+1*mm, rowHeights=vsize+7*mm))
    Story.append(Spacer(1,2*mm))
    Story.append(PageBreak())

    #==============================================================================
    # Time averaged spectrum of SC VIS2
    #==============================================================================

    plotTitle(Story,"SC squared visibility vs. wavelength (VIS2DATA SC)")
    plotSubtitle(Story,"Squared visibility vs. wavelength (microns).")
    
    hsize = 16*cm
    vsize = 6.5*cm
    yminval = 0.0
    ymaxval = grandmaxvis
    ticky = 0.2
    tickx = 0.1

    d=()
    for baseline in range(0,3):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis2_sc_vis2data[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis2_sc_vis2data_s[baseline::6,::resamp], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis2_sc_vis2data_p[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        a = graphoutaxes( data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'SC VIS2DATA '+base_list[baseline])
        #a = graphoutaxes( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline])
        Story.append(a)
        Story.append(Spacer(1,7*mm))
    Story.append(PageBreak())


    plotTitle(Story,"SC squared visibility vs. wavelength (VIS2DATA SC)")
    plotSubtitle(Story,"Squared visibility vs. wavelength (microns).")

    for baseline in range(3,6):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis2_sc_vis2data[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis2_sc_vis2data_s[baseline::6,::resamp], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis2_sc_vis2data_p[baseline::6,::resamp], axis=0),yminval,ymaxval))))
        a = graphoutaxes( data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'SC VIS2DATA '+base_list[baseline])
        #a = graphoutaxes( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline])
        Story.append(a)
        Story.append(Spacer(1,7*mm))

    Story.append(PageBreak())



    #==============================================================================
    # Time averaged spectrum of FT VIS2
    #==============================================================================

    plotTitle(Story,"FT squared visibility sq. vs. wavelength (VIS2 FT)")
    plotSubtitle(Story,"Squared visibility vs. wavelength (microns).")

    hsize = 8*cm
    vsize = 6.5*cm
    yminval = 0.0
    ymaxval = grandmaxvis
    ticky = 0.2
    
    d=()
    for baseline in range(0,nbase):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis2_ft_vis2data[baseline::6,:], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis2_ft_vis2data_s[baseline::6,:], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis2_ft_vis2data_p[baseline::6,:], axis=0),yminval,ymaxval))))
        if(baseline == 4 or baseline == 2 or baseline == 0):
            d = d + (graphoutaxes(data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]),)
        else:
            d = d + (graphoutnoaxis( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]),)

    Story.append(Table( [[d[0],d[1]],[d[2],d[3]],[d[4],d[5]]], colWidths=hsize+1*mm, rowHeights=vsize+7*mm))
    Story.append(Spacer(1,2*mm))
    Story.append(PageBreak())
    

    #==============================================================================
    # SC T3PHI vs mean spatial frequency
    #==============================================================================

    hsize = 16*cm
    vsize = 9*cm

    plotTitle(Story,"SC closure phase vs. spatial frequency (T3PHI SC)")
    plotSubtitle(Story,"Closure phase in degrees vs. mean spatial frequency over triangle (megalambda).")

    xminval = int(np.min(spfreq_sc_t3phi))-5
    xmaxval = int(np.max(spfreq_sc_t3phi))+5
    if oifits.polarsplit == True:
        ymaxval = int(np.abs(np.percentile([oifits.t3phi_sc_s*180./np.pi,oifits.t3phi_sc_p*180./np.pi],98)))+10
    else:
        ymaxval = int(np.abs(np.percentile(oifits.t3phi_sc*180./np.pi,98)))+10
    yminval = -ymaxval
    tickx = (xmaxval-xminval)/10
    ticky = (ymaxval-yminval)/10.
    
    if oifits.polarsplit == True:
        data=[]
        for triangle in range(0,nt3):
            data.append( tuple(zip(spfreq_sc_t3phi[triangle,::resamp], clipdata(oifits.t3phi_sc_s[triangle,::resamp]*180./np.pi,yminval,ymaxval)) ))
        a = graphoutaxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,'SC T3PHI S')
        plotSubtitle(Story,"Polarization S.")
        Story.append(a)
        
        Story.append(Spacer(1,7*mm))

        data=[]
        for triangle in range(0,nt3):
            data.append( tuple(zip(spfreq_sc_t3phi[triangle,::resamp], clipdata(oifits.t3phi_sc_p[triangle,::resamp]*180./np.pi,yminval,ymaxval)) ))
        a = graphoutaxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,'SC T3PHI P')
        plotSubtitle(Story,"Polarization P.")
        Story.append(a)

        Story.append(PageBreak())

    else:
        data=[]
        for triangle in range(0,nt3):
            data.append( tuple(zip(spfreq_sc_t3phi[triangle,::resamp], clipdata(oifits.t3phi_sc[triangle,::resamp]*180./np.pi,yminval,ymaxval)) ))
        a = graphoutaxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,'SC T3PHI')
        Story.append(a)

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

    #==============================================================================
    # FT T3PHI vs mean spatial frequency
    #==============================================================================

    plotTitle(Story,"FT closure phase vs. spatial frequency (T3PHI FT)")
    plotSubtitle(Story,"Closure phase in degrees vs. mean spatial frequency over triangle (megalambda).")

    # xminval = int(np.min(spfreq_ft_t3phi))-5
    # xmaxval = int(np.max(spfreq_ft_t3phi))+5
    # if oifits.polarsplit == True:
    #     ymaxval = int(np.abs(np.percentile([oifits.t3phi_ft_s*180./np.pi,oifits.t3phi_ft_p*180./np.pi],98)))+5
    # else:
    #     ymaxval = int(np.abs(np.percentile(oifits.t3phi_ft*180./np.pi,98)))+5
    # yminval = -ymaxval
    # tickx = (xmaxval-xminval)/10
    # ticky = (ymaxval-yminval)/10.
    
    if oifits.polarsplit == True:
        data=[]
        for triangle in range(0,nt3):
            data.append( tuple(zip(spfreq_ft_t3phi[triangle,:], clipdata(oifits.t3phi_ft_s[triangle,:]*180./np.pi,yminval,ymaxval)) ))
        a = graphoutaxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,'FT T3PHI S')
        plotSubtitle(Story,"Polarization S.")
        Story.append(a)
        
        Story.append(Spacer(1,7*mm))

        data=[]
        for triangle in range(0,nt3):
            data.append( tuple(zip(spfreq_ft_t3phi[triangle,:], clipdata(oifits.t3phi_ft_p[triangle,:]*180./np.pi,yminval,ymaxval)) ))
        a = graphoutaxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,'FT T3PHI P')
        plotSubtitle(Story,"Polarization P.")
        Story.append(a)

        Story.append(PageBreak())

    else:
        data=[]
        for triangle in range(0,nt3):
            data.append( tuple(zip(spfreq_ft_t3phi[triangle,:], clipdata(oifits.t3phi_ft[triangle,:]*180./np.pi,yminval,ymaxval)) ))
        a = graphoutaxes( data,xminval,xmaxval,'%i',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,'FT T3PHI')
        Story.append(a)

        Story.append(PageBreak())
        
    #==============================================================================
    # T3PHI spectrum of SC
    #==============================================================================

    hsize = 16*cm
    vsize = 5*cm
    tickx = 0.1
    
    plotTitle(Story,"SC closure phase vs. wavelength (T3PHI SC)")
    plotSubtitle(Story,"Closure phase in degrees vs. wavelength (microns).")
     
    d=()
    for triplet in range(0,4):
        data=[]
        if oifits.polarsplit == False:
            # same scale as SC
            #yminval = -np.percentile(np.abs(oifits.t3phi_sc),99)*180./np.pi-2.
            #ymaxval = +np.percentile(np.abs(oifits.t3phi_sc),99)*180./np.pi+2.
            #ticky = (ymaxval-yminval)/10.
            data.append(tuple(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.t3phi_sc[triplet::4,::resamp],axis=0)*180./np.pi,yminval,ymaxval))))
        else:
            #yminval = -np.percentile([np.abs(oifits.t3phi_sc_s),np.abs(oifits.t3phi_sc_p)],99)*180./np.pi-2.
            #ymaxval = +np.percentile([np.abs(oifits.t3phi_sc_s),np.abs(oifits.t3phi_sc_p)],99)*180./np.pi+2.
            #ticky = (ymaxval-yminval)/10.
            data.append(tuple(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.t3phi_sc_s[triplet::4,::resamp],axis=0)*180./np.pi,yminval,ymaxval))))
            data.append(tuple(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.t3phi_sc_p[triplet::4,::resamp],axis=0)*180./np.pi,yminval,ymaxval))))
        a = graphoutaxes( data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,triplet_list[triplet])
#        a = graphoutaxes( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,triplet_list[triplet])
        Story.append(a)
        if triplet <=2: Story.append(Spacer(1,7*mm))

    Story.append(PageBreak())

    #==============================================================================
    # T3PHI spectrum of FT
    #==============================================================================

    hsize = 16*cm
    vsize = 5*cm

    plotTitle(Story,"FT closure phase vs. wavelength (T3PHI FT)")
    plotSubtitle(Story,"Closure phase in degrees vs. wavelength (microns).")
     
    d=()
    for triplet in range(0,4):
        data=[]
        if oifits.polarsplit == False:
            # same limits as SC
            # yminval = -np.percentile(np.abs(oifits.t3phi_ft),99)*180./np.pi-2.
            # ymaxval = +np.percentile(np.abs(oifits.t3phi_ft),99)*180./np.pi+2.
            # ticky = (ymaxval-yminval)/10.
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.t3phi_sc[triplet::4,:],axis=0)*180./np.pi,yminval,ymaxval))))
        else:
            # same limits as SC
            # yminval = -np.percentile([np.abs(oifits.t3phi_ft_s),np.abs(oifits.t3phi_ft_p)],99)*180./np.pi-2.
            # ymaxval = +np.percentile([np.abs(oifits.t3phi_ft_s),np.abs(oifits.t3phi_ft_p)],99)*180./np.pi+2.
            # ticky = (ymaxval-yminval)/10.
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.t3phi_ft_s[triplet::4,:],axis=0)*180./np.pi,yminval,ymaxval))))
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.t3phi_ft_p[triplet::4,:],axis=0)*180./np.pi,yminval,ymaxval))))
        a = graphscatteraxes( data,minwave,maxwave,'%.2f',yminval,ymaxval,'%i',hsize,vsize,tickx,ticky,6,triplet_list[triplet])
        #a = graphoutaxes( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,triplet_list[triplet])
        Story.append(a)
        if triplet <=2: Story.append(Spacer(1,7*mm))

    Story.append(PageBreak())


    #==============================================================================
    # VISPHI spectrum of SC
    #==============================================================================

    plotTitle(Story,"SC differential phase vs. wavelength (VISPHI SC)")
    plotSubtitle(Story,"Closure phase in degrees vs. wavelength (microns).")

    hsize = 8*cm
    vsize = 6.5*cm
    if oifits.polarsplit == False:
        ymaxval = np.nanpercentile(np.abs(oifits.oi_vis_sc_visphi)*180./np.pi,98)*2.
        yminval = -ymaxval
    else:
        ymaxval = np.nanpercentile(np.abs(oifits.oi_vis_sc_visphi_s)*180./np.pi,98)*2.
        yminval = -ymaxval

    ticky = (ymaxval-yminval)/10.

    d=()
    for baseline in range(0,nbase):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visphi[baseline::6,::resamp]*180./np.pi, axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visphi_s[baseline::6,::resamp]*180./np.pi, axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(np.nanmean(oifits.oi_vis_sc_visphi_p[baseline::6,::resamp]*180./np.pi, axis=0),yminval,ymaxval))))
        if(baseline == 4 or baseline == 2 or baseline == 0):
            d = d + (graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]),)
        else:
            d = d + (graphoutnoaxis(data,minwaveSC,maxwaveSC,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]),)

    Story.append(Table( [[d[0],d[1]],[d[2],d[3]],[d[4],d[5]]], colWidths=hsize+1*mm, rowHeights=vsize+7*mm))
    Story.append(PageBreak())


    #==============================================================================
    # VISPHI spectrum of FT
    #==============================================================================

    plotTitle(Story,"FT differential phase vs. wavelength (VISPHI SC)")
    plotSubtitle(Story,"Phase in degrees vs. wavelength (microns).")

    # Same limits as SC
#    if oifits.polarsplit == False:
#        ymaxval = np.nanpercentile(np.abs(oifits.oi_vis_ft_visphi)*180./np.pi,98)*2.
#        yminval = -ymaxval
#    else:
#        ymaxval = np.nanpercentile(np.abs(oifits.oi_vis_ft_visphi_s)*180./np.pi,98)*2.
#        yminval = -ymaxval

    ticky = (ymaxval-yminval)/10.     
    d=()
    for baseline in range(0,nbase):
        data = []
        if oifits.polarsplit == False:
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis_ft_visphi[baseline::6,:], axis=0),yminval,ymaxval))))
        else:
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis_ft_visphi_s[baseline::6,:], axis=0),yminval,ymaxval))))
            data.append(list(zip(oifits.wave_ft, clipdata(np.nanmean(oifits.oi_vis_ft_visphi_p[baseline::6,:], axis=0),yminval,ymaxval))))
        if(baseline == 4 or baseline == 2 or baseline == 0):
            d = d + (graphoutaxes(data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]),)
        else:
            d = d + (graphoutnoaxis( data, minwave,maxwave,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]),)

    Story.append(Table( [[d[0],d[1]],[d[2],d[3]],[d[4],d[5]]], colWidths=hsize+1*mm, rowHeights=vsize+7*mm))
    Story.append(PageBreak())

    #==============================================================================
    # Polar diff VISDATA
    #==============================================================================
    hsize = 8*cm
    vsize = 6.5*cm
    yminval = -180
    ymaxval = +180
    ticky = 45.
 
    if oifits.polarsplit == True:
        plotTitle(Story,"Polar differential phase VISDATA_s*conj(VISDATA_p)")
        plotSubtitle(Story,"Phase in degrees vs. wavelength (microns). Red for SC and Orange for FT")
     
        for baseline in range(0,nbase):
            data = []
            phasor = np.angle(np.nanmean(oifits.oi_vis_sc_visdata_s[baseline::6,:] * np.conj(oifits.oi_vis_sc_visdata_p[baseline::6,:]),axis=0))*180./np.pi
            data.append(list(zip(oifits.wave_sc[::resamp],phasor[::resamp])))
            phasor = np.angle(np.nanmean(oifits.oi_vis_ft_visdata_s[baseline::6,:] * np.conj(oifits.oi_vis_ft_visdata_p[baseline::6,:]),axis=0))*180./np.pi
            data.append(list(zip(oifits.wave_ft,clipdata(phasor,yminval,ymaxval))))
            if(baseline == 4 or baseline == 2 or baseline == 0):
                d = d + (graphoutaxes(data,minwaveSC,maxwaveSC,'%i',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]),)
            else:
                d = d + (graphoutnoaxis( data,minwaveSC,maxwaveSC,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]),)
    
        Story.append(Table( [[d[0],d[1]],[d[2],d[3]],[d[4],d[5]]], colWidths=hsize+1*mm, rowHeights=vsize+7*mm))
    
        Story.append(PageBreak())


    #==============================================================================
    #  Time average SC-FT dOPD vs. wavelength (SC VISPHI)
    #==============================================================================

    #print oifits.polarsplit
    #if (oifits.single==True)
    if oifits.polarsplit == False:
 
        plotTitle(Story,"Time average SC-FT dOPD vs. wavelength (OI_VIS/VISPHI)")
        plotSubtitle(Story,"In microns dOPD vs. wavelength (microns).")
    
        Story.append(Spacer(1,2*mm))
        
        d=()
        data=[]
        Story.append(Spacer(1,5*mm))
        hsize=5*cm
        vsize=10*cm
        yminval = np.nanpercentile(oifits.oi_vis_sc_visphi*oifits.wave_sc/(2*np.pi),2)
        ymaxval = np.nanpercentile(oifits.oi_vis_sc_visphi*oifits.wave_sc/(2*np.pi),99)
        ticky = (ymaxval-yminval)/10 # visibility tick
        visphi = np.zeros((nbase,oifits.wave_sc.shape[0]),dtype='d')
        for baseline in range(0,nbase):
            visphi[baseline,:] = np.nanmean(oifits.oi_vis_sc_visphi[baseline::6,:], axis=0)*oifits.wave_sc/(2*np.pi)
            #visphi[baseline,:] = np.angle(np.nanmean(oifits.oi_vis_sc_visdata[baseline::6,:], axis=0))*180./np.pi
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(visphi[baseline,::resamp],yminval,ymaxval))))
            #print data
            if(baseline == 3):
                a = graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]),
            else:
                a = graphoutnoaxis( data,minwaveSC,maxwaveSC,yminval,ymaxval,hsize,vsize,ticky,base_list[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 oifits.polarsplit == True:
 
        plotTitle(Story,"Time average SC-FT dOPD vs. wavelength (OI_VIS/VISPHI)")
        plotSubtitle(Story,"In microns dOPD vs. wavelength (microns).")
    
        Story.append(Spacer(1,2*mm))
        
        d=()
        data=[]
        Story.append(Spacer(1,5*mm))

        hsize=5*cm
        vsize=5*cm
        yminval = np.nanpercentile(oifits.oi_vis_sc_visphi_s*oifits.wave_sc/(2*np.pi),2)
        ymaxval = np.nanpercentile(oifits.oi_vis_sc_visphi_s*oifits.wave_sc/(2*np.pi),99)
        ticky = (ymaxval-yminval)/10 # visibility tick
        visphi_s = np.zeros((nbase,oifits.wave_sc.shape[0]),dtype='d')
        for baseline in range(0,nbase):
            visphi_s[baseline,:] = np.nanmean(oifits.oi_vis_sc_visphi_s[baseline::6,:]*oifits.wave_sc/(2*np.pi), axis=0)
            #visphi_s[baseline,:] = np.angle(np.nanmean(oifits.oi_vis_sc_visdata_s[baseline::6,:], axis=0))*180./np.pi
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(visphi_s[baseline,::resamp],yminval,ymaxval))))
            a = graphoutnoaxis( data,minwaveSC,maxwaveSC,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]+'-S'),
            data = []
            d = d + a

        yminval = np.nanpercentile(oifits.oi_vis_sc_visphi_p*oifits.wave_sc/(2*np.pi),2)
        ymaxval = np.nanpercentile(oifits.oi_vis_sc_visphi_p*oifits.wave_sc/(2*np.pi),99)
        ticky = (ymaxval-yminval)/10 # visibility tick
        visphi_p = np.zeros((nbase,oifits.wave_sc.shape[0]),dtype='d')
        for baseline in range(0,nbase):
            visphi_p[baseline,:] = np.nanmean(oifits.oi_vis_sc_visphi_p[baseline::6,:]*oifits.wave_sc/(2*np.pi), axis=0)
            #visphi_p[baseline,:] = np.angle(np.nanmean(oifits.oi_vis_sc_visdata_p[baseline::6,:], axis=0))*180./np.pi
            data.append(list(zip(oifits.wave_sc[::resamp], clipdata(visphi_p[baseline,::resamp],yminval,ymaxval))))
            if(baseline == 3):
                a = graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,base_list[baseline]+'-P'),
            else:
                a = graphoutnoaxis( data,minwaveSC,maxwaveSC,yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]+'-P'),
            data = []
            d = d + a

        plotmatrix = [list(d[0:int(nbase/2)]),list(d[int(nbase/2):nbase]),list(d[nbase:nbase+int(nbase/2)]),list(d[nbase+int(nbase/2):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 oifits.polarsplit == True:
        plotTitle(Story,"Difference between S and P polar. vs. wavelength for FT")
        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(oifits.oi_vis_ft_visdata_s[baseline::6,:] * np.conj(oifits.oi_vis_ft_visdata_p[baseline::6,:]), axis=0))*180./np.pi;
            data.append(list(zip(oifits.wave_ft, tmp)))
            a = graphoutaxes(data,minwave,maxwave,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'FT diff. phase S-P Baseline '+base_list[baseline])
            Story.append(a)
            Story.append(Spacer(1,7*mm))

        Story.append(PageBreak())

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

    if oifits.polarsplit == True:
        plotTitle(Story,"Difference between S and P polar. vs. wavelength for SC")
        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(oifits.oi_vis_sc_visdata_s[baseline::6,:] * np.conj(oifits.oi_vis_sc_visdata_p[baseline::6,:]), axis=0))*180./np.pi;
            data.append(list(zip(oifits.wave_sc, tmp)))
            a = graphoutaxes(data,minwaveSC,maxwaveSC,'%.2f',yminval,ymaxval,'%.2f',hsize,vsize,tickx,ticky,'SC diff. phase S-P Baseline '+base_list[baseline])
            Story.append(a)
            if baseline != (nbase-1): Story.append(Spacer(1,7*mm));

    #==============================================================================
    #==============================================================================
    #==============================================================================
    if oifits.nrecords > 1: # plot time dependent quantities only if several observations are in the file
    #==============================================================================
    #==============================================================================
    #==============================================================================

        #==============================================================================
        # Flux signals vs time
        #==============================================================================
 
        plotTitle(Story,"Wavelength averaged flux values for SC vs. time (OI_FLUX)")
        plotSubtitle(Story,"In photo-e- vs. time (decimal hours). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.")
        
        Story.append(Spacer(1,5*mm))
    
        if oifits.polarsplit == False:
            avgflux = np.zeros((4,oifits.nrecords),dtype='d')
            for tel in range(0,4):
                avgflux[tel,:]=np.nanmean(oifits.oi_flux_sc[:,tel,:],axis=1) # average flux over lambda
        
            yminval = 0
            ymaxval = np.max(avgflux)
            ticky = (ymaxval-yminval)/10 # in photo-e-
            hsize = 16*cm
            vsize = 9*cm
            data = []
            for tel in range(0,4):
                data.append( tuple(zip(oifits.time, clipdata(avgflux[tel,:],yminval,ymaxval)) ))
            a = graphoutaxes( data, min(oifits.time),max(oifits.time),yminval,ymaxval,hsize,vsize,ticky,'Combined-polarization')
            Story.append(a)
            
        if oifits.polarsplit == True:
            avgfluxs = np.zeros((4,oifits.nrecords),dtype='d')
            avgfluxp = np.zeros((4,oifits.nrecords),dtype='d')
    
            for tel in range(0,4):
                avgfluxs[tel,:]=np.nanmean(oifits.oi_flux_sc_s[:,tel,:],axis=1) # average flux over lambda
                avgfluxp[tel,:]=np.nanmean(oifits.oi_flux_sc_p[:,tel,:],axis=1) # average flux over lambda
        
            yminval = 0
            ymaxval = np.max(np.array([avgfluxs,avgfluxp]))
            ticky = (ymaxval-yminval)/10 # in photo-e-
            hsize = 16*cm
            vsize = 4.5*cm
            data = []
            for tel in range(0,4):
                data.append( tuple(zip(oifits.time, clipdata(avgfluxs[tel,:],yminval,ymaxval)) ))
            a = graphoutaxes( data, min(oifits.time),max(oifits.time),yminval,ymaxval,hsize,vsize,ticky,'S-polarization')
            Story.append(a)
            
            Story.append(Spacer(1,7*mm))
        
            yminval = 0
            ymaxval = np.max(avgfluxp)
            ticky = (ymaxval-yminval)/10 # in photo-e-
            data = []
            for tel in range(0,4):
                data.append( tuple(zip(oifits.time, clipdata(avgfluxp[tel,:],yminval,ymaxval)) ))
            a = graphoutaxes( data, min(oifits.time),max(oifits.time),yminval,ymaxval,hsize,vsize,ticky,'P-polarization')
            Story.append(a)
    
        Story.append(Spacer(1,5*mm))
    
        plotTitle(Story,"Wavelength averaged flux values for FT vs. time (OI_FLUX)")
        plotSubtitle(Story,"In photo-e- vs. time (decimal hours). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.")
        
        Story.append(Spacer(1,5*mm))
    
        if oifits.polarsplit == False:
            nmeasure = oifits.oi_flux_ft_time.shape[0]
            avgflux = np.zeros((4,nmeasure),dtype='d')
            for tel in range(0,4):
                avgflux[tel,:]=np.nanmean(oifits.oi_flux_ft[:,tel,:],axis=1) # average flux over lambda
        
            hsize = 16*cm
            vsize = 9*cm
            yminval = 0
            ymaxval = np.max(avgflux)*1.2
            ticky = (ymaxval-yminval)/10 # in photo-e-
            data = []
            for tel in range(0,4):
                data.append( tuple(zip(oifits.time_ft, clipdata(avgflux[tel,:],yminval,ymaxval)) ))
            a = graphoutaxes( data, min(oifits.time_ft),max(oifits.time_ft),yminval,ymaxval,hsize,vsize,ticky,'Combined-polarization')
            Story.append(a)
    
        if oifits.polarsplit == True:
            nmeasure = oifits.oi_flux_ft_time.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(oifits.oi_flux_ft_s[:,tel,:],axis=1) # average flux over lambda
                avgfluxp[tel,:]=np.nanmean(oifits.oi_flux_ft_p[:,tel,:],axis=1) # average flux over lambda
        
            hsize = 16*cm
            vsize = 4.5*cm
            yminval = 0
            ymaxval = np.max(avgfluxs)*1.2
            ticky = (ymaxval-yminval)/10 # in photo-e-
            data = []
            for tel in range(0,4):
                data.append( tuple(zip(oifits.time, clipdata(avgfluxs[tel,:],yminval,ymaxval)) ))
            a = graphoutaxes( data, min(oifits.time),max(oifits.time),yminval,ymaxval,hsize,vsize,ticky,'S-polarization')
            Story.append(a)
            
            Story.append(Spacer(1,7*mm))
        
            yminval = 0
            ymaxval = np.max(avgfluxp)*1.2
            ticky = (ymaxval-yminval)/10 # in photo-e-
            data = []
            for tel in range(0,4):
                data.append( tuple(zip(oifits.time, clipdata(avgfluxp[tel,:],yminval,ymaxval)) ))
            a = graphoutaxes( data, min(oifits.time),max(oifits.time),yminval,ymaxval,hsize,vsize,ticky,'P-polarization')
            Story.append(a)
    
        Story.append(PageBreak())
    
        #==============================================================================
        # Wavelength median SC-FT dOPD vs. time (SC VISPHI)
        #==============================================================================
    
        if oifits.polarsplit == False:
            plotTitle(Story,"Wavelength median SC-FT dOPD vs. time (OI_VIS/VISPHI)")
            plotSubtitle(Story,"In radians vs. time (decimal hours).")
            
            d=()
            data=[]
            Story.append(Spacer(1,5*mm))
            #print oifits.polarsplit
            
            yminval = np.nanpercentile(oifits.oi_vis_sc_visphi*oifits.wave_sc/(2*np.pi),2)
            ymaxval = np.nanpercentile(oifits.oi_vis_sc_visphi*oifits.wave_sc/(2*np.pi),98)
            ticky = (ymaxval-yminval)/10 # visibility tick
            hsize=5*cm
            vsize=10*cm
            visphi = np.zeros((nbase,oifits.oi_vis_sc_time.shape[0]),dtype='d')
            for baseline in range(0,nbase):
                visphi[baseline,:] = np.nanmedian(oifits.oi_vis_sc_visphi[baseline::6,:]*oifits.wave_sc/(2*np.pi), axis=1)
                data.append(list(zip(oifits.oi_vis_sc_time, clipdata(visphi[baseline,:],yminval,ymaxval))))
                #print data
                if(baseline == 3):
                    a = graphoutaxes( data, min(oifits.oi_vis_sc_time), max(oifits.oi_vis_sc_time),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]),
                else:
                    a = graphoutnoaxis( data, min(oifits.oi_vis_sc_time), max(oifits.oi_vis_sc_time),yminval,ymaxval,hsize,vsize,ticky,base_list[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 oifits.polarsplit == True:
            plotTitle(Story,"Wavelength median SC-FT dOPD vs. time (OI_VIS/VISPHI)")
            plotSubtitle(Story,"In microns dOPD vs. time (decimal hours).")
        
            d=()
            data=[]
            Story.append(Spacer(1,5*mm))
            #print oifits.polarsplit
            
            hsize=5*cm
            vsize=5*cm
            visphi_s = np.zeros((nbase,oifits.oi_vis_sc_time.shape[0]),dtype='d')
            visphi_p = np.zeros((nbase,oifits.oi_vis_sc_time.shape[0]),dtype='d')
    
            yminval = np.nanpercentile(oifits.oi_vis_sc_visphi_s*oifits.wave_sc/(2*np.pi),2)
            ymaxval = np.nanpercentile(oifits.oi_vis_sc_visphi_s*oifits.wave_sc/(2*np.pi),98)
            ticky = (ymaxval-yminval)/10 # visibility tick
    
            for baseline in range(0,nbase):
                visphi_s[baseline,:] = np.nanmedian(oifits.oi_vis_sc_visphi_s[baseline::6,:]*oifits.wave_sc/(2*np.pi), axis=1)
                data.append(list(zip(oifits.oi_vis_sc_time, clipdata(visphi_s[baseline,:],yminval,ymaxval))))
                a = graphoutnoaxis( data, min(oifits.oi_vis_sc_time), max(oifits.oi_vis_sc_time),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]+'-S'),
                data = []
                d = d + a
    
            yminval = np.nanpercentile(oifits.oi_vis_sc_visphi_p*oifits.wave_sc/(2*np.pi),2)
            ymaxval = np.nanpercentile(oifits.oi_vis_sc_visphi_p*oifits.wave_sc/(2*np.pi),98)
            ticky = (ymaxval-yminval)/10 # visibility tick
            for baseline in range(0,nbase):
                visphi_p[baseline,:] = np.nanmedian(oifits.oi_vis_sc_visphi_p[baseline::6,:]*oifits.wave_sc/(2*np.pi), axis=1)
                data.append(list(zip(oifits.oi_vis_sc_time, clipdata(visphi_p[baseline,:],yminval,ymaxval))))
                #print data
                if(baseline == 3):
                    a = graphoutaxes( data, min(oifits.oi_vis_sc_time), max(oifits.oi_vis_sc_time),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]+'-P'),
                else:
                    a = graphoutnoaxis( data, min(oifits.oi_vis_sc_time), max(oifits.oi_vis_sc_time),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]+'-P'),
                data = []
                d = d + a
            plotmatrix = [list(d[0:int(nbase/2)]),list(d[int(nbase/2):nbase]),list(d[nbase:nbase+int(nbase/2)]),list(d[nbase+int(nbase/2):2*nbase])]
            t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize)
            Story.append(t)
            Story.append(Spacer(1,2*mm))
    
        Story.append(PageBreak())
    
    
        #==============================================================================
        # Image (waterfall) display of FT+SC-MET phase vs. time
        #==============================================================================
        
        plotTitle(Story,"SC-FT dOPD waterfall view (VISPHI SC)")
        plotSubtitle(Story,"In microns dOPD. Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.")
        
        if oifits.polarsplit == False:
            valmin = np.nanpercentile(oifits.oi_vis_sc_visphi*oifits.wave_sc/(2*np.pi),1)
            valmax = np.nanpercentile(oifits.oi_vis_sc_visphi*oifits.wave_sc/(2*np.pi),99)
            imgdata = io.StringIO()
            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(oifits.oi_vis_sc_visphi[0::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[0].set_title('SC-FT dOPD - Baseline 12 - Comb', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_sc_visphi[1::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[1].set_title('SC-FT dOPD - Baseline 13 - Comb', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_sc_visphi[2::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[2].set_title('SC-FT dOPD - Baseline 14 - Comb', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_sc_visphi[3::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[3].set_title('SC-FT dOPD - Baseline 23 - Comb', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_sc_visphi[4::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[4].set_title('SC-FT dOPD - Baseline 24 - Comb', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_sc_visphi[5::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[5].set_title('SC-FT dOPD - Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
    
        
        if oifits.polarsplit == True:
            imgdata = io.StringIO()
            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))
    
            valmin = np.nanpercentile(oifits.oi_vis_sc_visphi_s*oifits.wave_sc/(2*np.pi),1)
            valmax = np.nanpercentile(oifits.oi_vis_sc_visphi_s*oifits.wave_sc/(2*np.pi),99)
    
            im = axarr[0].imshow(oifits.oi_vis_sc_visphi_s[0::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[0].set_title('SC-FT dOPD - Baseline 12 - S', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_sc_visphi_s[1::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[1].set_title('SC-FT dOPD - Baseline 13 - S', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_sc_visphi_s[2::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[2].set_title('SC-FT dOPD - Baseline 14 - S', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_sc_visphi_s[3::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[3].set_title('SC-FT dOPD - Baseline 23 - S', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_sc_visphi_s[4::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[4].set_title('SC-FT dOPD - Baseline 24 - S', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_sc_visphi_s[5::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[5].set_title('SC-FT dOPD - Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            Story.append(PageBreak())
    
            imgdata = io.StringIO()
            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))
    
            valmin = np.nanpercentile(oifits.oi_vis_sc_visphi_p*oifits.wave_sc/(2*np.pi),1)
            valmax = np.nanpercentile(oifits.oi_vis_sc_visphi_p*oifits.wave_sc/(2*np.pi),99)
    
            im = axarr[0].imshow(oifits.oi_vis_sc_visphi_p[0::6,:]*oifits.wave_sc,vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[0].set_title('SC-FT dOPD - Baseline 12 - P', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_sc_visphi_p[1::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[1].set_title('SC-FT dOPD - Baseline 13 - P', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_sc_visphi_p[2::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[2].set_title('SC-FT dOPD - Baseline 14 - P', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_sc_visphi_p[3::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[3].set_title('SC-FT dOPD - Baseline 23 - P', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_sc_visphi_p[4::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[4].set_title('SC-FT dOPD - Baseline 24 - P', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_sc_visphi_p[5::6,:]*oifits.wave_sc/(2*np.pi),vmin=valmin,vmax=valmax,cmap='jet', interpolation='nearest', aspect='auto')
            axarr[5].set_title('SC-FT dOPD - Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            
        Story.append(PageBreak())
    
        #==============================================================================
        # Image (waterfall) display of SC VISAMP
        #==============================================================================
        
        plotTitle(Story,"SC visibility amplitude waterfall view (VISAMP SC)")
        plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.")
    
        #fig = plt.figure(figsize=(3,7),dpi=250)
        #valmin = np.nanpercentile(oifits.oi_vis_visamp,1)
        #valmax = np.nanpercentile(oifits.oi_vis_visamp,95)
        valmin = -0.2
        valmax = 2.
        #figaspect = heightfig/widthfig
        #plt.rc('font', family='serif')
    
        if oifits.polarsplit == False:
            imgdata = io.StringIO()
            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(oifits.oi_vis_sc_visamp[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('SC VISAMP Baseline 12 - Comb', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_sc_visamp[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('SC VISAMP Baseline 13 - Comb', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_sc_visamp[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('SC VISAMP Baseline 14 - Comb', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_sc_visamp[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('SC VISAMP Baseline 23 - Comb', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_sc_visamp[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('SC VISAMP Baseline 24 - Comb', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_sc_visamp[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('SC VISAMP Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
    
    
        if oifits.polarsplit == True:
            imgdata = io.StringIO()
            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(oifits.oi_vis_sc_visamp_s[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('SC VISAMP Baseline 12 - S', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_sc_visamp_s[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('SC VISAMP Baseline 13 - S', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_sc_visamp_s[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('SC VISAMP Baseline 14 - S', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_sc_visamp_s[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('SC VISAMP Baseline 23 - S', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_sc_visamp_s[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('SC VISAMP Baseline 24 - S', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_sc_visamp_s[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('SC VISAMP Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            Story.append(PageBreak())
    
            imgdata = io.StringIO()
            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(oifits.oi_vis_sc_visamp_p[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('SC VISAMP Baseline 12 - P', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_sc_visamp_p[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('SC VISAMP Baseline 13 - P', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_sc_visamp_p[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('SC VISAMP Baseline 14 - P', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_sc_visamp_p[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('SC VISAMP Baseline 23 - P', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_sc_visamp_p[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('SC VISAMP Baseline 24 - P', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_sc_visamp_p[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('SC VISAMP Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            
        Story.append(PageBreak())
        
    
    
    
    
    
        #==============================================================================
        # Image (waterfall) display of FT VISAMP
        #==============================================================================
    
        plotTitle(Story,"FT visibility amplitude waterfall view (VISAMP FT)")
        plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.")
    
        valmin = -0.4
        valmax = 2.
    
        if oifits.polarsplit == False:
            imgdata = io.StringIO()
            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(oifits.oi_vis_ft_visamp[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('FT VISAMP Baseline 12 - Comb', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_ft_visamp[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('FT VISAMP Baseline 13 - Comb', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_ft_visamp[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('FT VISAMP Baseline 14 - Comb', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_ft_visamp[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('FT VISAMP Baseline 23 - Comb', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_ft_visamp[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('FT VISAMP Baseline 24 - Comb', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_ft_visamp[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('FT VISAMP Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
    
        if oifits.polarsplit == True:
            imgdata = io.StringIO()
            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(oifits.oi_vis_ft_visamp_s[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('FT VISAMP Baseline 12 - S', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_ft_visamp_s[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('FT VISAMP Baseline 13 - S', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_ft_visamp_s[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('FT VISAMP Baseline 14 - S', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_ft_visamp_s[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('FT VISAMP Baseline 23 - S', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_ft_visamp_s[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('FT VISAMP Baseline 24 - S', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_ft_visamp_s[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('FT VISAMP Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            Story.append(PageBreak())
    
            imgdata = io.StringIO()
            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(oifits.oi_vis_ft_visamp_p[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('FT VISAMP Baseline 12 - P', fontsize=8)
            axarr[1].imshow(oifits.oi_vis_ft_visamp_p[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('FT VISAMP Baseline 13 - P', fontsize=8)
            axarr[2].imshow(oifits.oi_vis_ft_visamp_p[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('FT VISAMP Baseline 14 - P', fontsize=8)
            axarr[3].imshow(oifits.oi_vis_ft_visamp_p[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('FT VISAMP Baseline 23 - P', fontsize=8)
            axarr[4].imshow(oifits.oi_vis_ft_visamp_p[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('FT VISAMP Baseline 24 - P', fontsize=8)
            axarr[5].imshow(oifits.oi_vis_ft_visamp_p[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('FT VISAMP Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            
        Story.append(PageBreak())
    
    
        #==============================================================================
        # Image (waterfall) display of FT VIS2DATA
        #==============================================================================
    
        plotTitle(Story,"FT squared visibility amplitude waterfall view (VIS2DATA FT)")
        plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.")
    
        valmin = -0.4
        valmax = 2.
    
        if oifits.polarsplit == False:
            imgdata = io.StringIO()
            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(oifits.oi_vis2_ft_vis2data[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('FT VIS2DATA Baseline 12 - Comb', fontsize=8)
            axarr[1].imshow(oifits.oi_vis2_ft_vis2data[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('FT VIS2DATA Baseline 13 - Comb', fontsize=8)
            axarr[2].imshow(oifits.oi_vis2_ft_vis2data[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('FT VIS2DATA Baseline 14 - Comb', fontsize=8)
            axarr[3].imshow(oifits.oi_vis2_ft_vis2data[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('FT VIS2DATA Baseline 23 - Comb', fontsize=8)
            axarr[4].imshow(oifits.oi_vis2_ft_vis2data[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('FT VIS2DATA Baseline 24 - Comb', fontsize=8)
            axarr[5].imshow(oifits.oi_vis2_ft_vis2data[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('FT VIS2DATA Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
    
        if oifits.polarsplit == True:
            imgdata = io.StringIO()
            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(oifits.oi_vis2_ft_vis2data_s[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('FT VIS2DATA Baseline 12 - S', fontsize=8)
            axarr[1].imshow(oifits.oi_vis2_ft_vis2data_s[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('FT VIS2DATA Baseline 13 - S', fontsize=8)
            axarr[2].imshow(oifits.oi_vis2_ft_vis2data_s[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('FT VIS2DATA Baseline 14 - S', fontsize=8)
            axarr[3].imshow(oifits.oi_vis2_ft_vis2data_s[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('FT VIS2DATA Baseline 23 - S', fontsize=8)
            axarr[4].imshow(oifits.oi_vis2_ft_vis2data_s[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('FT VIS2DATA Baseline 24 - S', fontsize=8)
            axarr[5].imshow(oifits.oi_vis2_ft_vis2data_s[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('FT VIS2DATA Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            Story.append(PageBreak())
    
            imgdata = io.StringIO()
            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(oifits.oi_vis2_ft_vis2data_p[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[0].set_title('FT VIS2DATA Baseline 12 - P', fontsize=8)
            axarr[1].imshow(oifits.oi_vis2_ft_vis2data_p[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[1].set_title('FT VIS2DATA Baseline 13 - P', fontsize=8)
            axarr[2].imshow(oifits.oi_vis2_ft_vis2data_p[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[2].set_title('FT VIS2DATA Baseline 14 - P', fontsize=8)
            axarr[3].imshow(oifits.oi_vis2_ft_vis2data_p[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[3].set_title('FT VIS2DATA Baseline 23 - P', fontsize=8)
            axarr[4].imshow(oifits.oi_vis2_ft_vis2data_p[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[4].set_title('FT VIS2DATA Baseline 24 - P', fontsize=8)
            axarr[5].imshow(oifits.oi_vis2_ft_vis2data_p[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto')
            axarr[5].set_title('FT VIS2DATA Baseline 34 - 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)
            #fig.tight_layout()
            #plt.colorbar()
            widthfig = 16 * cm
            heightfig= 23 * cm
            plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight')
            pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound')
            plt.close('all')
            Story.append(pi1)
            
    Story.append(PageBreak())
    

    #==============================================================================
    # Tables of numerical values
    #==============================================================================

    plotTitle(Story,"Average photometric flux (photo-e-/s/sp.channel +/- std)")
    
    Story.append(Spacer(1,1*mm))

    if oifits.polarsplit == False:
        avgfluxsc = []
        avgfluxft = []
        for tel in range(0,4):
            avgfluxsc.append('%(val).2e ' % {"val":np.nanmean(oifits.oi_flux_sc[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmedian(oifits.oi_flux_err_sc[:,tel,:])})
            avgfluxft.append('%(val).2e ' % {"val":np.nanmean(oifits.oi_flux_ft[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmedian(oifits.oi_flux_err_ft[:,tel,:])})
            # average flux over wavelength
        data = avgfluxsc
        data = np.vstack((['Tel1','Tel2','Tel3','Tel4'],data))
        data = np.vstack((data,avgfluxft))
        data = np.hstack(([["Telescope"],["SC"],["FT"]],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 oifits.polarsplit == True:
        avgfluxs = []
        avgfluxp = []
        for tel in range(0,4):
            avgfluxs.append('%(val).2e ' % {"val":np.nanmean(oifits.oi_flux_sc_s[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmedian(oifits.oi_flux_err_sc_s[:,tel,:])})
            avgfluxp.append('%(val).2e ' % {"val":np.nanmean(oifits.oi_flux_sc_p[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmedian(oifits.oi_flux_err_sc_p[:,tel,:])})
            # 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)
    
        Story.append(Spacer(1,2*mm))
    
        avgfluxs = []
        avgfluxp = []
        for tel in range(0,4):
            avgfluxs.append('%(val).2e ' % {"val":np.nanmean(oifits.oi_flux_ft_s[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmedian(oifits.oi_flux_err_ft_s[:,tel,:])})
            avgfluxp.append('%(val).2e ' % {"val":np.nanmean(oifits.oi_flux_ft_p[:,tel,:])}+ '\u00B1' +' %(err).1e' % 
                {"err":np.nanmedian(oifits.oi_flux_err_ft_p[:,tel,:])})
        # average flux over wavelength
        data = avgfluxs
        data = np.vstack((['Tel1','Tel2','Tel3','Tel4'],data))
        data = np.vstack((data,avgfluxp))
        data = np.hstack(([["FT"],["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)

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

    plotTitle(Story,"Average visibility amplitude per baseline (vis +/- std)")
    
    Story.append(Spacer(1,1*mm))

    if oifits.polarsplit == False:
        avgvisampsc = []
        avgvisampft = []
        for base in range(0,nbase):
            visamp_sc = nanaverage(oifits.oi_vis_sc_visamp[base::6,:],oifits.oi_vis_sc_visamperr[base::6,:]**(-2))
            errvisamp_sc = np.nanmedian(oifits.oi_vis_sc_visamperr[base::6,:])
            visamp_ft = nanaverage(oifits.oi_vis_ft_visamp[base::6,:],oifits.oi_vis_ft_visamperr[base::6,:]**(-2))
            errvisamp_ft = np.nanmedian(oifits.oi_vis_ft_visamperr[base::6,:])
            if visamp_sc > 3 or visamp_sc <-3 or errvisamp_sc > 3 or errvisamp_sc<-3:
                avgvisampsc.append('bad')
            else:
                avgvisampsc.append('%(val).3f ' % {"val":visamp_sc}+ '\u00B1' +' %(err).3f' % {"err":errvisamp_sc})
            if visamp_ft > 3 or visamp_ft < -3 or errvisamp_ft > 3 or errvisamp_ft < -3:
                avgvisampft.append('bad')
            else:
                avgvisampft.append('%(val).3f ' % {"val":visamp_ft}+ '\u00B1' +' %(err).3f' % {"err":errvisamp_ft})
            # average flux over wavelength
        data = avgvisampsc
        data = np.vstack((base_list,data))
        data = np.vstack((data,avgvisampft))
        data = np.hstack(([["Baseline"],["SC"],["FT"]],data))
        t=Table(data.tolist(),colWidths=2.6*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 oifits.polarsplit == True:
        #visamp_s[baseline,:] = np.nanmean(oifits.oi_vis_sc_visamp_s[baseline::6,:], axis=0)

        avgvisamps = []
        avgvisampp = []

        for base in range(0,nbase):
            visamp_s = nanaverage(oifits.oi_vis_sc_visamp_s[base::6,:],oifits.oi_vis_sc_visamperr_s[base::6,:]**(-2))
            errvisamp_s = np.nanmedian(oifits.oi_vis_sc_visamperr_s[base::6,:])
            visamp_p = nanaverage(oifits.oi_vis_sc_visamp_p[base::6,:],oifits.oi_vis_sc_visamperr_p[base::6,:]**(-2))
            errvisamp_p = np.nanmedian(oifits.oi_vis_sc_visamperr_p[base::6,:])
            if visamp_s > 3 or visamp_s <-3 or errvisamp_s > 3 or errvisamp_s<-3:
                avgvisamps.append('bad')
            else:
                avgvisamps.append('%(val).3f ' % {"val":visamp_s}+ '\u00B1' +' %(err).3f' % {"err":errvisamp_s})
            if visamp_p > 3 or visamp_p < -3 or errvisamp_p > 3 or errvisamp_p < -3:
                avgvisampp.append('bad')
            else:
                avgvisampp.append('%(val).3f ' % {"val":visamp_p}+ '\u00B1' +' %(err).3f' % {"err":errvisamp_p})
            # average flux over wavelength
        data = avgvisamps
        data = np.vstack((base_list,data))
        data = np.vstack((data,avgvisampp))
        data = np.hstack(([["SC"],["Polar. S"],["Polar. P"]],data))
        t=Table(data.tolist(),colWidths=2.6*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)
    
        Story.append(Spacer(1,2*mm))
    
        avgvisamps = []
        avgvisampp = []
        for base in range(0,nbase):
            visamp_s = nanaverage(oifits.oi_vis_ft_visamp_s[base::6,:],oifits.oi_vis_ft_visamperr_s[base::6,:]**(-2))
            errvisamp_s = np.nanmedian(oifits.oi_vis_ft_visamperr_s[base::6,:])
            visamp_p = nanaverage(oifits.oi_vis_ft_visamp_p[base::6,:],oifits.oi_vis_ft_visamperr_p[base::6,:]**(-2))
            errvisamp_p = np.nanmedian(oifits.oi_vis_ft_visamperr_p[base::6,:])
            if visamp_s > 3 or visamp_s <-3 or errvisamp_s > 3 or errvisamp_s<-3:
                avgvisamps.append('bad')
            else:
                avgvisamps.append('%(val).3f ' % {"val":visamp_s}+ '\u00B1' +' %(err).3f' % {"err":errvisamp_s})
                
            if visamp_p > 3 or visamp_p < -3 or errvisamp_p > 3 or errvisamp_p < -3:
                avgvisampp.append('bad')
            else:               
                avgvisampp.append('%(val).3f ' % {"val":visamp_p}+ '\u00B1' +' %(err).3f' % {"err":errvisamp_p})
            # average flux over wavelength
        data = avgvisamps
        data = np.vstack((base_list,data))
        data = np.vstack((data,avgvisampp))
        data = np.hstack(([["FT"],["Polar. S"],["Polar. P"]],data))
        t=Table(data.tolist(),colWidths=2.6*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)

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

    plotTitle(Story,"Average squared visibility per baseline (vis^2 +/- std)")
    
    Story.append(Spacer(1,1*mm))

    if oifits.polarsplit == False:
        avgvis2sc = []
        avgvis2ft = []
        for base in range(0,nbase):
            vis2_sc = nanaverage(oifits.oi_vis2_sc_vis2data[base::6,:],oifits.oi_vis2_sc_vis2err[base::6,:]**(-2))
            errvis2_sc = np.nanmedian(oifits.oi_vis2_sc_vis2err[base::6,:])
            vis2_ft = nanaverage(oifits.oi_vis2_ft_vis2data[base::6,:],oifits.oi_vis2_ft_vis2err[base::6,:]**(-2))
            errvis2_ft = np.nanmedian(oifits.oi_vis2_ft_vis2err[base::6,:])
            if vis2_sc > 3 or vis2_sc <-3 or errvis2_sc > 3 or errvis2_sc<-3:
                avgvis2sc.append('bad')
            else:
                avgvis2sc.append('%(val).3f ' % {"val":vis2_sc}+ '\u00B1' +' %(err).3f' % {"err":errvis2_sc})
            if vis2_ft > 3 or vis2_ft < -3 or errvis2_ft > 3 or errvis2_ft < -3:
                avgvis2ft.append('bad')
            else:
                avgvis2ft.append('%(val).3f ' % {"val":vis2_ft}+ '\u00B1' +' %(err).3f' % {"err":errvis2_ft})
            # average flux over wavelength
        data = avgvis2sc
        data = np.vstack((base_list,data))
        data = np.vstack((data,avgvis2ft))
        data = np.hstack(([["Baseline"],["SC"],["FT"]],data))
        t=Table(data.tolist(),colWidths=2.6*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 oifits.polarsplit == True:
        #vis2_s[baseline,:] = np.nanmean(oifits.oi_vis_sc_vis2_s[baseline::6,:], axis=0)

        avgvis2s = []
        avgvis2p = []

        for base in range(0,nbase):
            vis2_s = nanaverage(oifits.oi_vis2_sc_vis2data_s[base::6,:],oifits.oi_vis2_sc_vis2err_s[base::6,:]**(-2))
            errvis2_s = np.nanmedian(oifits.oi_vis2_sc_vis2err_s[base::6,:])
            vis2_p = nanaverage(oifits.oi_vis2_sc_vis2data_p[base::6,:],oifits.oi_vis2_sc_vis2err_p[base::6,:]**(-2))
            errvis2_p = np.nanmedian(oifits.oi_vis2_sc_vis2err_p[base::6,:])
            if vis2_s > 3 or vis2_s <-3 or errvis2_s > 3 or errvis2_s<-3:
                avgvis2s.append('bad')
            else:
                avgvis2s.append('%(val).3f ' % {"val":vis2_s}+ '\u00B1' +' %(err).3f' % {"err":errvis2_s})
            if vis2_p > 3 or vis2_p < -3 or errvis2_p > 3 or errvis2_p < -3:
                avgvis2p.append('bad')
            else:
                avgvis2p.append('%(val).3f ' % {"val":vis2_p}+ '\u00B1' +' %(err).3f' % {"err":errvis2_p})
            # average flux over wavelength
        data = avgvis2s
        data = np.vstack((base_list,data))
        data = np.vstack((data,avgvis2p))
        data = np.hstack(([["SC"],["Polar. S"],["Polar. P"]],data))
        t=Table(data.tolist(),colWidths=2.6*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)
    
        Story.append(Spacer(1,2*mm))
    
        avgvis2s = []
        avgvis2p = []
        for base in range(0,nbase):
            vis2_s = nanaverage(oifits.oi_vis2_ft_vis2data_s[base::6,:],oifits.oi_vis2_ft_vis2err_s[base::6,:]**(-2))
            errvis2_s = np.nanmedian(oifits.oi_vis2_ft_vis2err_s[base::6,:])
            vis2_p = nanaverage(oifits.oi_vis2_ft_vis2data_p[base::6,:],oifits.oi_vis2_ft_vis2err_p[base::6,:]**(-2))
            errvis2_p = np.nanmedian(oifits.oi_vis2_ft_vis2err_p[base::6,:])
            if vis2_s > 3 or vis2_s <-3 or errvis2_s > 3 or errvis2_s<-3:
                avgvis2s.append('bad')
            else:
                avgvis2s.append('%(val).3f ' % {"val":vis2_s}+ '\u00B1' +' %(err).3f' % {"err":errvis2_s})
                
            if vis2_p > 3 or vis2_p < -3 or errvis2_p > 3 or errvis2_p < -3:
                avgvis2p.append('bad')
            else:               
                avgvis2p.append('%(val).3f ' % {"val":vis2_p}+ '\u00B1' +' %(err).3f' % {"err":errvis2_p})
            # average flux over wavelength
        data = avgvis2s
        data = np.vstack((base_list,data))
        data = np.vstack((data,avgvis2p))
        data = np.hstack(([["FT"],["Polar. S"],["Polar. P"]],data))
        t=Table(data.tolist(),colWidths=2.6*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)

        Story.append(PageBreak())

    plotTitle(Story,"Average closure phase per triplet (t3phi +/- std), in degrees")
    
    Story.append(Spacer(1,1*mm))

    if oifits.polarsplit == False:
        avgt3sc = []
        avgt3ft = []
        for triplet in range(0,4):
            t3_sc = nanaverage(oifits.t3phi_sc[triplet,:],oifits.t3phierr_sc[triplet,:]**(-2))
            stdt3_sc = np.nanmedian(oifits.t3phierr_sc[triplet,:])
            t3_ft = nanaverage(oifits.t3phi_ft[triplet,:],oifits.t3phierr_ft[triplet,:]**(-2))
            stdt3_ft = np.nanmedian(oifits.t3phierr_ft[triplet,:])
            avgt3sc.append('%(val).3f ' % {"val":t3_sc}+ '\u00B1' +' %(err).3f' % {"err":stdt3_sc})
            avgt3ft.append('%(val).3f ' % {"val":t3_ft}+ '\u00B1' +' %(err).3f' % {"err":stdt3_ft})
            # average flux over wavelength
        data = avgt3sc
        data = np.vstack((triplet_list,data))
        data = np.vstack((data,avgt3ft))
        data = np.hstack(([["Triplet"],["SC (degrees)"],["FT (degrees)"]],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 oifits.polarsplit == True:
        avgt3scs = []
        avgt3scp = []
        for triplet in range(0,4):
            t3_scs = nanaverage(oifits.t3phi_sc_s[triplet,:],oifits.t3phierr_sc_s[triplet,:]**(-2))
            stdt3_scs = np.nanmedian(oifits.t3phierr_sc_s[triplet,:])
            avgt3scs.append('%(val).3f ' % {"val":t3_scs}+ '\u00B1' +' %(err).3f' % {"err":stdt3_scs})
            t3_scp = nanaverage(oifits.t3phi_sc_p[triplet,:],oifits.t3phierr_sc_p[triplet,:]**(-2))
            stdt3_scp = np.nanmedian(oifits.t3phierr_sc_p[triplet,:])
            avgt3scp.append('%(val).3f ' % {"val":t3_scp}+ '\u00B1' +' %(err).3f' % {"err":stdt3_scp})
            # average flux over wavelength
        data = avgt3scs
        data = np.vstack((data,avgt3scp))
        data = np.vstack((triplet_list,data))
        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)
        Story.append(Spacer(1,2*mm))
        avgt3fts = []
        avgt3ftp = []
        for triplet in range(0,4):
            t3_fts = nanaverage(oifits.t3phi_ft_s[triplet,:],oifits.t3phierr_ft_s[triplet,:]**(-2))
            stdt3_fts = np.nanmedian(oifits.t3phierr_ft_s[triplet,:])
            avgt3fts.append('%(val).3f ' % {"val":t3_fts}+ '\u00B1' +' %(err).3f' % {"err":stdt3_fts})
            t3_ftp = nanaverage(oifits.t3phi_ft_p[triplet,:],oifits.t3phierr_ft_p[triplet,:]**(-2))
            stdt3_ftp = np.nanmedian(oifits.t3phierr_ft_p[triplet,:])
            avgt3ftp.append('%(val).3f ' % {"val":t3_ftp}+ '\u00B1' +' %(err).3f' % {"err":stdt3_ftp})
            # average flux over wavelength
        data = avgt3fts
        data = np.vstack((data,avgt3ftp))
        data = np.vstack((triplet_list,data))
        data = np.hstack(([["FT"],["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)
    
    Story.append(Spacer(1,5*mm))


    plotTitle(Story,"Average closure amplitude per triplet (t3phi +/- std), in degrees")
    
    Story.append(Spacer(1,1*mm))

    if oifits.polarsplit == False:
        avgt3sc = []
        avgt3ft = []
        for triplet in range(0,4):
            t3_sc = nanaverage(oifits.t3amp_sc[triplet,:],oifits.t3amperr_sc[triplet,:]**(-2))
            stdt3_sc = np.nanmedian(oifits.t3amperr_sc[triplet,:])
            t3_ft = nanaverage(oifits.t3amp_ft[triplet,:],oifits.t3amperr_ft[triplet,:]**(-2))
            stdt3_ft = np.nanmedian(oifits.t3amperr_ft[triplet,:])
            avgt3sc.append('%(val).3f ' % {"val":t3_sc}+ '\u00B1' +' %(err).3f' % {"err":stdt3_sc})
            avgt3ft.append('%(val).3f ' % {"val":t3_ft}+ '\u00B1' +' %(err).3f' % {"err":stdt3_ft})
            # average flux over wavelength
        data = avgt3sc
        data = np.vstack((triplet_list,data))
        data = np.vstack((data,avgt3ft))
        data = np.hstack(([["Triplet"],["SC (degrees)"],["FT (degrees)"]],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 oifits.polarsplit == True:
        avgt3scs = []
        avgt3scp = []
        for triplet in range(0,4):
            t3_scs = nanaverage(oifits.t3amp_sc_s[triplet,:],oifits.t3amperr_sc_s[triplet,:]**(-2))
            stdt3_scs = np.nanmedian(oifits.t3amperr_sc_s[triplet,:])
            avgt3scs.append('%(val).3f ' % {"val":t3_scs}+ '\u00B1' +' %(err).3f' % {"err":stdt3_scs})
            t3_scp = nanaverage(oifits.t3amp_sc_p[triplet,:],oifits.t3amperr_sc_p[triplet,:]**(-2))
            stdt3_scp = np.nanmedian(oifits.t3amperr_sc_p[triplet,:])
            avgt3scp.append('%(val).3f ' % {"val":t3_scp}+ '\u00B1' +' %(err).3f' % {"err":stdt3_scp})
            # average flux over wavelength
        data = avgt3scs
        data = np.vstack((data,avgt3scp))
        data = np.vstack((triplet_list,data))
        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)
        Story.append(Spacer(1,2*mm))
        avgt3fts = []
        avgt3ftp = []
        for triplet in range(0,4):
            t3_fts = nanaverage(oifits.t3amp_ft_s[triplet,:],oifits.t3amperr_ft_s[triplet,:]**(-2))
            stdt3_fts = np.nanmedian(oifits.t3amperr_ft_s[triplet,:])
            avgt3fts.append('%(val).3f ' % {"val":t3_fts}+ '\u00B1' +' %(err).3f' % {"err":stdt3_fts})
            t3_ftp = nanaverage(oifits.t3amp_ft_p[triplet,:],oifits.t3amperr_ft_p[triplet,:]**(-2))
            stdt3_ftp = np.nanmedian(oifits.t3amperr_ft_p[triplet,:])
            avgt3ftp.append('%(val).3f ' % {"val":t3_ftp}+ '\u00B1' +' %(err).3f' % {"err":stdt3_ftp})
            # average flux over wavelength
        data = avgt3fts
        data = np.vstack((data,avgt3ftp))
        data = np.vstack((triplet_list,data))
        data = np.hstack(([["FT"],["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)

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

    print((" "+reportname)) 

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

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

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

    oifits = gravi_visual_class.Oifits(filename+'.fits')

    produce_oifits_report(oifits,filename)