#! /usr/bin/env python # -*- coding: iso-8859-15 -*- """ Created on Mon May 25 13:51:41 2015 @author: kervella """ # ATTENTION: it is 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 from . import bootstrap as boot import numpy as np import pdb import warnings try: from scipy.signal import welch welchpresent = True except ImportError as e: welchpresent = False pass # module doesn't exist, deal with it. #============================================================================== # Preparation of the report using reportlab #============================================================================== def produce_p2vmreduced_report(p2vmreduced,filename): # From gravi_visual_class from . import gravi_visual_class from .gravi_visual_class import PdfImage, myFirstPage, myLaterPages, styles from .gravi_visual_class import clipdata, transbarchartout, graphoutaxes, graphoutnoaxis from .gravi_visual_class import graphoutnoxaxes, graphscatteraxes, graphaxesplt from .gravi_visual_class import plotTitle, plotSubtitle from .gravi_visual_class import get_key_withdefault, create_array_from_list, baseline_phases from .gravi_visual_class import mean_angle, std_angle, clean_unwrap_phase, nanaverage from .gravi_visual_class import clean_gdelay_is, clean_gdelay_fft, clean_gdelay_full from .gravi_visual_class import base_list, tel_list, nbase, ntel from .gravi_visual_class import plotReductionSummary, plotSummary from reportlab.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 #============================================================================== # Start Story #============================================================================== Story = [Spacer(1,1*mm)] plotSummary (Story, filename, p2vmreduced, onTarget=True) if 'HIERARCH ESO INS MLC WAVELENG' in p2vmreduced.header: Lambda_met=p2vmreduced.header['HIERARCH ESO INS MLC WAVELENG']/1000. # Metrology wavelength in microns (1.908 microns) else: Lambda_met = 1.908287 gravitybasenames = ['12','13','14','23','24','34'] nframeplot = p2vmreduced.nframe_ft #============================================================================== # Photometry #============================================================================== style = styles["Heading2"] p = Paragraph("Average photometric flux (photo-e-/DIT/sp.channel +/- std)", style) Story.append(p) Story.append(Spacer(1,2*mm)) if hasattr(p2vmreduced,'oi_flux_sc') & (p2vmreduced.polarsplit == False): avgfluxsc = [] avgfluxft = [] for tel in range(0,4): avgfluxsc.append('%(val).2e ' % {"val":np.nanmean(p2vmreduced.oi_flux_sc[:,tel,:])}+ '\u00B1' +' %(err).1e' % {"err":np.nanmean(np.nanstd(p2vmreduced.oi_flux_sc[:,tel,:],axis=0))}) avgfluxft.append('%(val).2e ' % {"val":np.nanmean(p2vmreduced.oi_flux_ft[:,tel,:])}+ '\u00B1' +' %(err).1e' % {"err":np.nanmean(np.nanstd(p2vmreduced.oi_flux_ft[:,tel,:],axis=0))}) # average flux over wavelength data = avgfluxsc data = np.vstack((['Tel1 '+p2vmreduced.stations[0],'Tel2 '+p2vmreduced.stations[1],\ 'Tel3 '+p2vmreduced.stations[2],'Tel4 '+p2vmreduced.stations[3]],data)) # print data # pdb.set_trace() 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 hasattr(p2vmreduced, 'oi_flux_sc_s') & (p2vmreduced.polarsplit == True): avgfluxs = [] avgfluxp = [] for tel in range(0,4): avgfluxs.append('%(val).2e ' % {"val":np.nanmean(p2vmreduced.oi_flux_sc_s[:,tel,:])}+ '\u00B1' +' %(err).1e' % {"err":np.nanmean(np.nanstd(p2vmreduced.oi_flux_sc_s[:,tel,:],axis=0))}) avgfluxp.append('%(val).2e ' % {"val":np.nanmean(p2vmreduced.oi_flux_sc_p[:,tel,:])}+ '\u00B1' +' %(err).1e' % {"err":np.nanmean(np.nanstd(p2vmreduced.oi_flux_sc_p[:,tel,:],axis=0))}) # average flux over wavelength data = avgfluxs data = np.vstack((['Tel1','Tel2','Tel3','Tel4'],data)) data = np.vstack((data,avgfluxp)) data = np.hstack(([["SC"],["Polar. S"],["Polar. P"]],data)) t=Table(data.tolist(),colWidths=3.5*cm) t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black), ('BACKGROUND',(0,0),(0,0),colors.pink), ('BOX', (0,0), (-1,-1), 2, colors.black),])) Story.append(t) Story.append(Spacer(1,7*mm)) avgfluxs = [] avgfluxp = [] for tel in range(0,4): avgfluxs.append('%(val).2e ' % {"val":np.nanmean(p2vmreduced.oi_flux_ft_s[:,tel,:])}+ '\u00B1' +' %(err).1e' % {"err":np.nanmean(np.nanstd(p2vmreduced.oi_flux_ft_s[:,tel,:],axis=0))}) avgfluxp.append('%(val).2e ' % {"val":np.nanmean(p2vmreduced.oi_flux_ft_p[:,tel,:])}+ '\u00B1' +' %(err).1e' % {"err":np.nanmean(np.nanstd(p2vmreduced.oi_flux_ft_p[:,tel,:],axis=0))}) # average flux over wavelength data = avgfluxs data = np.vstack((['Tel1','Tel2','Tel3','Tel4'],data)) data = np.vstack((data,avgfluxp)) data = np.hstack(([["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) #============================================================================== # Grand average photometric spectrum of SC #============================================================================== plotTitle(Story,"SC Average spectrum of all outputs (OI_FLUX)") plotSubtitle(Story,"In photo-e- vs. wavelength (microns).") if hasattr(p2vmreduced,'oi_flux_sc') & (p2vmreduced.polarsplit == False): nmeasure = p2vmreduced.time_sc.shape[0] avgflux = np.zeros((p2vmreduced.nwave_sc),dtype='d') avgflux = np.nanmean(np.nanmean(p2vmreduced.oi_flux_sc,axis=0), axis=0) # average flux over time, then over telescopes yminval = 0 # np.nanpercentile(avgflux,1)-0.1*np.abs(np.nanpercentile(avgflux,99)) ymaxval = 1.3*np.nanpercentile(avgflux,99) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 4*cm data = [ tuple(zip(p2vmreduced.wave_sc, clipdata(avgflux[:],yminval,ymaxval)) ) ] a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined polarization') Story.append(a) if hasattr(p2vmreduced,'oi_flux_sc_s') & (p2vmreduced.polarsplit == True): nmeasure = p2vmreduced.time_sc.shape[0] avgfluxs = np.zeros((p2vmreduced.nwave_sc),dtype='d') avgfluxp = np.zeros((p2vmreduced.nwave_sc),dtype='d') avgfluxs =np.nanmean(np.nanmean(p2vmreduced.oi_flux_sc_s,axis=0), axis=0) # average flux over time, then over telescopes avgfluxp =np.nanmean(np.nanmean(p2vmreduced.oi_flux_sc_p,axis=0), axis=0) # average flux over time, then over telescopes yminval = 0 # np.nanpercentile(avgfluxs,1)-0.1*np.abs(np.nanpercentile(avgfluxs,99)) ymaxval = 1.3*np.nanpercentile([avgfluxs,avgfluxp],99) ticky = (ymaxval-yminval)/5 # in photo-e- hsize = 16*cm vsize = 4*cm avgflux = np.mean([avgfluxs,avgfluxp],axis=0) data = [tuple(zip(p2vmreduced.wave_sc, clipdata(avgflux[:],yminval,ymaxval)) )] a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Average of S and P') Story.append(a) Story.append(Spacer(1,5*mm)) #============================================================================== # SC rejection flag #============================================================================== plotTitle(Story,"SC Rejection flag of frames") plotSubtitle(Story,"Flag on any of the six baselines and all polarizations vs. time (seconds).") yminval = -1 ymaxval = 4 ticky = 1 if p2vmreduced.polarsplit == True: hsize = 16*cm vsize = 2.0*cm global_flag = np.max(p2vmreduced.oi_vis_reject_flag_s,axis=1) else: hsize = 16*cm vsize = 4.0*cm global_flag = np.max(p2vmreduced.oi_vis_reject_flag,axis=1) data = [tuple(zip(p2vmreduced.time_sc, global_flag) )] a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'') a.add(String(0.03*hsize, 0.22*vsize, '0=Good', fontSize=9)) a.add(String(0.03*hsize, 0.42*vsize, '1=low FT SNR', fontSize=9)) a.add(String(0.03*hsize, 0.62*vsize, '2=low FT Vfactor', fontSize=9)) a.add(String(0.03*hsize, 0.82*vsize, '3=low FT Vfactor and low FT SNR', fontSize=9)) Story.append(a) Story.append(PageBreak()) #============================================================================== # Summary of reduction parameters #============================================================================== plotReductionSummary(Story, p2vmreduced) Story.append(PageBreak()) #============================================================================== # Average photometric spectra #============================================================================== plotTitle(Story,"SC Average spectrum (OI_FLUX)") plotSubtitle(Story,"In photo-e- vs. wavelength (microns). "+\ p2vmreduced.stations[0] +" = red, "+\ p2vmreduced.stations[1] +" = orange, "+\ p2vmreduced.stations[2] +" = blue, "+\ p2vmreduced.stations[3] +" = green.") if hasattr(p2vmreduced,'oi_flux_sc') & (p2vmreduced.polarsplit == False): nmeasure = p2vmreduced.time_sc.shape[0] avgflux = np.zeros((4,p2vmreduced.nwave_sc),dtype='d') for tel in range(0,4): avgflux[tel,:]=np.nanmean(p2vmreduced.oi_flux_sc[:,tel,:],axis=0) # average flux over time yminval = 0 # np.nanpercentile(avgflux,1)-0.1*np.abs(np.nanpercentile(avgflux,99)) ymaxval = 1.1*np.nanpercentile(avgflux,99) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 9*cm data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.wave_sc, clipdata(avgflux[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined-polarization') Story.append(a) if hasattr(p2vmreduced,'oi_flux_sc_s') & (p2vmreduced.polarsplit == True): nmeasure = p2vmreduced.time_sc.shape[0] avgfluxs = np.zeros((4,p2vmreduced.nwave_sc),dtype='d') avgfluxp = np.zeros((4,p2vmreduced.nwave_sc),dtype='d') for tel in range(0,4): avgfluxs[tel,:]=np.nanmean(p2vmreduced.oi_flux_sc_s[:,tel,:],axis=0) # average flux over time avgfluxp[tel,:]=np.nanmean(p2vmreduced.oi_flux_sc_p[:,tel,:],axis=0) # average flux over time yminval = 0 # np.nanpercentile(avgfluxs,1)-0.1*np.abs(np.nanpercentile(avgfluxs,99)) ymaxval = 1.1*np.nanpercentile([avgfluxs,avgfluxp],99) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 4.5*cm data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.wave_sc, clipdata(avgfluxs[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'S-polarization') Story.append(a) Story.append(Spacer(1,7*mm)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.wave_sc, clipdata(avgfluxp[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'P-polarization') Story.append(a) Story.append(Spacer(1,5*mm)) plotTitle(Story,"FT Average spectrum (OI_FLUX)") plotSubtitle(Story,"In photo-e- vs. wavelength (microns). "+\ p2vmreduced.stations[0] +" = red, "+\ p2vmreduced.stations[1] +" = orange, "+\ p2vmreduced.stations[2] +" = blue, "+\ p2vmreduced.stations[3] +" = green.") if (p2vmreduced.polarsplit == False): nmeasure = p2vmreduced.time_ft.shape[0] avgflux = np.zeros((4,p2vmreduced.nwave_ft),dtype='d') for tel in range(0,4): avgflux[tel,:]=np.nanmean(p2vmreduced.oi_flux_ft[:,tel,:],axis=0) # average flux over time hsize = 16*cm vsize = 9*cm yminval = 0 # np.nanmin(avgflux)-0.1*np.abs(np.nanmax(avgflux)) ymaxval = 1.1*np.nanmax(avgflux) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.wave_ft, clipdata(avgflux[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'SCombined-polarization') Story.append(a) if (p2vmreduced.polarsplit == True): nmeasure = p2vmreduced.time_ft.shape[0] avgfluxs = np.zeros((4,p2vmreduced.nwave_ft),dtype='d') avgfluxp = np.zeros((4,p2vmreduced.nwave_ft),dtype='d') for tel in range(0,4): avgfluxs[tel,:]=np.nanmean(p2vmreduced.oi_flux_ft_s[:,tel,:],axis=0) # average flux over time avgfluxp[tel,:]=np.nanmean(p2vmreduced.oi_flux_ft_p[:,tel,:],axis=0) # average flux over time hsize = 16*cm vsize = 4.5*cm yminval = 0 # np.nanmin(avgfluxs)-0.1*np.abs(np.nanmax(avgfluxs)) ymaxval = 1.1*np.nanmax([avgfluxs,avgfluxp]) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.wave_ft, clipdata(avgfluxs[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'S-polarization') Story.append(a) Story.append(Spacer(1,7*mm)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.wave_ft, clipdata(avgfluxp[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.wave_sc),max(p2vmreduced.wave_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'P-polarization') Story.append(a) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of SC OIFLUX for the four telescopes #============================================================================== plotTitle(Story,"SC flux waterfall per telescope (OI_FLUX SC)") plotSubtitle(Story,"In photo-e-/DIT, wavelength channel number in horizontal axis, number of frame in sequence in vertical axis. The average over the sequence has been subtracted and is given above each image.") if p2vmreduced.polarsplit == False: meanflux = np.zeros((4)) meansub = np.zeros((p2vmreduced.nframe_sc,4,p2vmreduced.nwave_sc)) for tel in range(0,4): meanflux[tel] = np.nanmean(p2vmreduced.oi_flux_sc[:,tel,:]) meansub[:,tel,:] = p2vmreduced.oi_flux_sc[:,tel,:] - meanflux[tel] valmin = np.nanpercentile(meansub,1) valmax = np.nanpercentile(meansub,99) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[0]+' - Comb - Mean = %.1e' % meanflux[0], fontsize=8) axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[1]+' - Comb - Mean = %.1e' % meanflux[1], fontsize=8) axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[2]+' - Comb - Mean = %.1e' % meanflux[2], fontsize=8) axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[3]+' - Comb - Mean = %.1e' % meanflux[3], fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) if p2vmreduced.polarsplit == True: meanflux = np.zeros((4)) meansub = np.zeros((p2vmreduced.nframe_sc,4,p2vmreduced.nwave_sc)) for tel in range(0,4): meanflux[tel] = np.nanmean(p2vmreduced.oi_flux_sc_s[:,tel,:]) meansub[:,tel,:] = p2vmreduced.oi_flux_sc_s[:,tel,:] - meanflux[tel] valmin = np.nanpercentile(meansub,1) valmax = np.nanpercentile(meansub,99) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[0]+' - S - Mean = %.1e' % meanflux[0], fontsize=8) axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[1]+' - S - Mean = %.1e' % meanflux[1], fontsize=8) axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[2]+' - S - Mean = %.1e' % meanflux[2], fontsize=8) axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[3]+' - S - Mean = %.1e' % meanflux[3], fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) plotTitle(Story,"SC flux waterfall per telescope (OI_FLUX SC)") plotSubtitle(Story,"In photo-e-/DIT, wavelength channel number in horizontal axis, number of frame in sequence in vertical axis. The average over the sequence has been subtracted and is given above each image.") meansub = np.zeros((p2vmreduced.nframe_sc,4,p2vmreduced.nwave_sc)) meanflux = np.zeros((4)) for tel in range(0,4): meanflux[tel] = np.nanmean(p2vmreduced.oi_flux_sc_p[:,tel,:]) meansub[:,tel,:] = p2vmreduced.oi_flux_sc_p[:,tel,:] - meanflux[tel] valmin = np.nanpercentile(meansub,1) valmax = np.nanpercentile(meansub,99) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[0]+' - P - Mean = %.1e' % meanflux[0], fontsize=8) axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[1]+' - P - Mean = %.1e' % meanflux[1], fontsize=8) axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[2]+' - P - Mean = %.1e' % meanflux[2], fontsize=8) axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC OIFLUX Telescope '+p2vmreduced.stations[3]+' - P - Mean = %.1e' % meanflux[3], fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, 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 OIFLUX for the four telescopes #============================================================================== plotTitle(Story,"FT flux waterfall per telescope (OIFLUX FT)") plotSubtitle(Story,"In photo-e-/DIT, mean value subtracted. Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.") if p2vmreduced.polarsplit == False: meansub = np.zeros((p2vmreduced.nframe_ft,4,p2vmreduced.nwave_ft)) meanflux = np.zeros((4)) for tel in range(0,4): meanflux[tel] = np.nanmean(p2vmreduced.oi_flux_ft[:,tel,:]) meansub[:,tel,:] = p2vmreduced.oi_flux_ft[:,tel,:] - meanflux[tel] valmin = np.nanpercentile(meansub,1) valmax = np.nanpercentile(meansub,99) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(meansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('FT Tel. '+p2vmreduced.stations[0]+' - Comb - Mean = %.1e' % meanflux[0], fontsize=8) axarr[1].imshow(meansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('FT Tel. '+p2vmreduced.stations[1]+' - Comb - Mean = %.1e' % meanflux[1], fontsize=8) axarr[2].imshow(meansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('FT Tel. '+p2vmreduced.stations[2]+' - Comb - Mean = %.1e' % meanflux[2], fontsize=8) axarr[3].imshow(meansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('FT Tel. '+p2vmreduced.stations[3]+' - Comb - Mean = %.1e' % meanflux[3], fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) if p2vmreduced.polarsplit == True: meansubs = np.zeros((p2vmreduced.nframe_ft,4,p2vmreduced.nwave_ft)) meansubp = np.zeros((p2vmreduced.nframe_ft,4,p2vmreduced.nwave_ft)) meanfluxs = np.zeros((4)) meanfluxp = np.zeros((4)) for tel in range(0,4): meanfluxs[tel] = np.nanmean(p2vmreduced.oi_flux_ft_s[:,tel,:]) meansubs[:,tel,:] = p2vmreduced.oi_flux_ft_s[:,tel,:] - meanfluxs[tel] meanfluxp[tel] = np.nanmean(p2vmreduced.oi_flux_ft_p[:,tel,:]) meansubp[:,tel,:] = p2vmreduced.oi_flux_ft_p[:,tel,:] - meanfluxp[tel] valmin = np.nanpercentile(meansubs,1) valmax = np.nanpercentile(meansubs,99) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=4, ncols=2, sharex=True, sharey=True, figsize=(5,8)) for tel in range(0,4): im = axarr[tel,0].imshow(meansubs[:,tel,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel,0].set_title('FT Tel. '+p2vmreduced.stations[tel]+'-S - Mean = %.1e' % meanfluxs[tel], fontsize=7) im = axarr[tel,1].imshow(meansubp[:,tel,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel,1].set_title('FT Tel. '+p2vmreduced.stations[tel]+'-P - Mean = %.1e' % meanfluxp[tel], fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Flux signals #============================================================================== plotTitle(Story,"SC flux vs. time (OI_FLUX wavelength averaged)") plotSubtitle(Story,"In photo-e-/DIT vs. seconds. "+\ p2vmreduced.stations[0] +" = red, "+\ p2vmreduced.stations[1] +" = orange, "+\ p2vmreduced.stations[2] +" = blue, "+\ p2vmreduced.stations[3] +" = green.") if hasattr(p2vmreduced,'oi_flux_sc') & (p2vmreduced.polarsplit == False): nmeasure = p2vmreduced.time_sc.shape[0] avgflux = np.zeros((4,nmeasure),dtype='d') for tel in range(0,4): avgflux[tel,:]=np.nanmean(p2vmreduced.oi_flux_sc[:,tel,:],axis=1) # average flux over wavelength yminval = 0 ymaxval = np.nanmax(avgflux)+0.1*np.abs(np.nanmax(avgflux)) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 10*cm data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.time_sc, clipdata(avgflux[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined-polarization') Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(p2vmreduced,'oi_flux_sc_s') & (p2vmreduced.polarsplit == True): nmeasure = p2vmreduced.time_sc.shape[0] avgfluxs = np.zeros((4,nmeasure),dtype='d') avgfluxp = np.zeros((4,nmeasure),dtype='d') for tel in range(0,4): avgfluxs[tel,:]=np.nanmean(p2vmreduced.oi_flux_sc_s[:,tel,:],axis=1) # average flux over wavelength avgfluxp[tel,:]=np.nanmean(p2vmreduced.oi_flux_sc_p[:,tel,:],axis=1) # average flux over wavelength yminval = 0 ymaxval = np.nanmax(avgfluxs)+0.1*np.abs(np.nanmax(avgfluxs)) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 10*cm data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.time_sc, clipdata(avgfluxs[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'S-polarization') Story.append(a) data = [] Story.append(Spacer(1,7*mm)) for tel in range(0,4): data.append( tuple(zip(p2vmreduced.time_sc, clipdata(avgfluxp[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'P-polarization') Story.append(a) Story.append(PageBreak()) #============================================================================== plotTitle(Story,"FT flux vs. time (OI_FLUX wavelength averaged)") plotSubtitle(Story,"In photo-e-/DIT vs. seconds. "+\ p2vmreduced.stations[0] +" = red, "+\ p2vmreduced.stations[1] +" = orange, "+\ p2vmreduced.stations[2] +" = blue, "+\ p2vmreduced.stations[3] +" = green.") resamp = nframeplot // 500 # take one value every resamp value for the plot boxsmooth_ft = 100 if hasattr(p2vmreduced,'oi_flux_ft') & (p2vmreduced.polarsplit == False): nmeasure = p2vmreduced.time_ft.shape[0] avgflux = np.zeros((4,nmeasure),dtype='d') avgfluxsmooth = np.zeros((4,nmeasure),dtype='d') for tel in range(0,4): avgflux[tel,:]=np.nanmean(p2vmreduced.oi_flux_ft[:,tel,:],axis=1) # average flux over wavelength avgfluxsmooth[tel,:]=np.convolve(avgflux[tel,:], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') yminval = np.nanmin(avgflux)-0.1*np.abs(np.nanmax(avgflux)) ymaxval = np.nanmax(avgflux)+0.1*np.abs(np.nanmax(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(p2vmreduced.time_ft[::resamp], clipdata(avgflux[tel,::resamp],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_ft),max(p2vmreduced.time_ft),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Combined-polarization') Story.append(a) if hasattr(p2vmreduced,'oi_flux_ft_s') & (p2vmreduced.polarsplit == True): nmeasure = p2vmreduced.time_ft.shape[0] avgfluxs = np.zeros((4,nmeasure),dtype='d') avgfluxp = np.zeros((4,nmeasure),dtype='d') avgflux = np.zeros((4,nmeasure),dtype='d') # average flux vector of both polarizations for the power spectrum computation avgfluxsmooths = np.zeros((4,nmeasure),dtype='d') avgfluxsmoothp = np.zeros((4,nmeasure),dtype='d') for tel in range(0,4): avgfluxs[tel,:]= np.nanmean(p2vmreduced.oi_flux_ft_s[:,tel,:],axis=1) # average flux over wavelength avgfluxp[tel,:]= np.nanmean(p2vmreduced.oi_flux_ft_p[:,tel,:],axis=1) # average flux over wavelength avgflux[tel,:] = np.nanmean([avgfluxs[tel,:],avgfluxp[tel,:]],axis=0) # average of both polarizations for the power spectrum computation avgfluxsmooths[tel,:]=np.convolve(avgfluxs[tel,:], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') avgfluxsmoothp[tel,:]=np.convolve(avgfluxp[tel,:], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') # np.savetxt('1_Avgflux-FT-per-tel.txt', avgflux.T, fmt='%+1.4e') # np.savetxt('0_Time-FT.txt', p2vmreduced.time_ft, fmt='%+1.4e') hsize = 16*cm vsize = 9*cm yminval = np.nanmin(avgfluxs)-0.1*np.abs(np.nanmax(avgfluxs)) ymaxval = np.nanmax(avgfluxs)+0.1*np.abs(np.nanmax(avgfluxs)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.time_ft[::resamp], clipdata(avgfluxsmooths[tel,::resamp],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_ft),max(p2vmreduced.time_ft),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'S-polarization') Story.append(a) Story.append(Spacer(1,7*mm)) yminval = np.nanmin(avgfluxp)-0.1*np.abs(np.nanmax(avgfluxp)) ymaxval = np.nanmax(avgfluxp)+0.1*np.abs(np.nanmax(avgfluxp)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,4): data.append( tuple(zip(p2vmreduced.time_ft[::resamp], clipdata(avgfluxsmoothp[tel,::resamp],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_ft),max(p2vmreduced.time_ft),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'P-polarization') Story.append(a) Story.append(PageBreak()) #============================================================================== # Correlation of the flux SC and FT #============================================================================== if hasattr(p2vmreduced,'totalflux_sc') or hasattr(p2vmreduced,'totalflux_sc_s'): robj_mag = get_key_withdefault(p2vmreduced.header,'HIERARCH ESO FT ROBJ MAG',-10.) sobj_mag = get_key_withdefault(p2vmreduced.header,'HIERARCH ESO INS SOBJ MAG',robj_mag) plotTitle(Story,"Total flux on SC vs FT (OI_FLUX)") plotSubtitle(Story,"transmission, FT horizontal (mag=%.1f), SC vertical (mag=%.1f)."%(robj_mag,sobj_mag)) # Total flux is in [e/sc_dit] dt_sc = p2vmreduced.time_sc[1]-p2vmreduced.time_sc[0] # time between each frame of the sc norm_ft = 1.71173e9 * 10**(-robj_mag / 2.5) * (np.pi * 0.9**2) * dt_sc norm_sc = 1.71173e9 * 10**(-sobj_mag / 2.5) * (np.pi * 0.9**2) * dt_sc hsize = 6*cm vsize = 6*cm if (p2vmreduced.polarsplit == False): yminval = 0. ymaxval = 1.2*np.nanmax([np.nanmax(p2vmreduced.totalflux_sc)/norm_sc,np.nanmax(p2vmreduced.totalflux_ft)/norm_ft]) ticky = (ymaxval-yminval)/5. d = [] for tel in range(0,4): data=[] data.append( tuple(zip(p2vmreduced.totalflux_ft[:,tel]/norm_ft, p2vmreduced.totalflux_sc[:,tel]/norm_sc) )) a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,ticky,ticky,2,'Combined Tel '+p2vmreduced.stations[tel]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]] ] 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) if (p2vmreduced.polarsplit == True): yminval = 0 ymaxval = 1.2*np.nanmax([np.nanmax(p2vmreduced.totalflux_sc_s)/norm_sc,np.nanmax(p2vmreduced.totalflux_ft_s)/norm_ft]) ticky = (ymaxval-yminval)/5. d = [] for tel in range(0,4): data=[] data.append( tuple(zip(p2vmreduced.totalflux_ft_s[:,tel]/norm_ft, p2vmreduced.totalflux_sc_s[:,tel]/norm_sc) )) a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,ticky,ticky,2,'S-polar Tel '+p2vmreduced.stations[tel]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]] ] 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()) #============================================================================== # Power spectral density of the FT flux signal vs. time #============================================================================== plotTitle(Story,"FT flux periodogram (OI_FLUX FT)") plotSubtitle(Story,"In flux^2 RMS vs. Hz., power spectrum in red, reverse cumulated in orange.") if welchpresent == True: if (p2vmreduced.polarsplit == False): PSD=[1,2,3,4] revcum=[1,2,3,4] for tel in range(0,4): f, PSD[tel] = welch(avgflux[tel,:]/avgflux[tel,:].max(), fs=(1./np.mean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') revcum[tel] = (PSD[tel]-np.mean(PSD[tel][-100:]))[::-1].cumsum()[::-1] yminval = 0 ymaxval = 1e-3 hsize = 16*cm vsize = 4.5*cm for tel in range(0,4): data = [] #ymaxval = np.nanmean(np.sqrt(PSD[tel]))*10 ticky = (ymaxval-yminval)/5 data.append(tuple(zip(f, clipdata(PSD[tel], yminval,ymaxval)))) data.append(tuple(zip(f, clipdata(revcum[tel], yminval, ymaxval)))) #a = graphoutaxes( data, min(f),max(f),yminval,ymaxval,hsize,vsize,ticky,'Telescope '+str(tel+1)) a = graphoutaxes( data, np.nanmin(f),np.nanmax(f),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Telescope '+p2vmreduced.stations[tel]) Story.append(a) Story.append(Spacer(1,7*mm)) if (p2vmreduced.polarsplit == True): PSD_s=[1,2,3,4] PSD_p=[1,2,3,4] revcum_s=[1,2,3,4] revcum_p=[1,2,3,4] for tel in range(0,4): f, PSD_s[tel] = welch(avgfluxs[tel,:]/avgfluxs[tel,:].max(), fs=(1./np.mean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') f, PSD_p[tel] = welch(avgfluxp[tel,:]/avgfluxp[tel,:].max(), fs=(1./np.mean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') revcum_s[tel] = (PSD_s[tel]-np.mean(PSD_s[tel][-100:]))[::-1].cumsum()[::-1] revcum_p[tel] = (PSD_p[tel]-np.mean(PSD_p[tel][-100:]))[::-1].cumsum()[::-1] yminval = 0 ymaxval = 1e-3 hsize = 16*cm vsize = 4.5*cm for tel in range(0,4): data = [] #ymaxval = np.nanmean(np.sqrt(PSD_s[tel]))*10 ticky = (ymaxval-yminval)/5 data.append(tuple(zip(f, clipdata((PSD_s[tel]+PSD_p[tel])/2.0, yminval,ymaxval)))) data.append(tuple(zip(f, clipdata((revcum_s[tel]+revcum_p[tel])/2.0, yminval, ymaxval)))) a = graphoutaxes( data, np.nanmin(f),np.nanmax(f),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Telescope '+p2vmreduced.stations[tel]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Power spectral density of the FT unwrapped phase vs. time #============================================================================== plotTitle(Story,"FT OPD periodogram - all frames (OI_VIS FT)") plotSubtitle(Story,"In um^2 vs. Hz., power spectrum in red, reverse cumulated in orange.") if welchpresent == True: if (p2vmreduced.polarsplit == True): visphi_s = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') visphi_p = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') for baseline in range(0,nbase): demod_s = np.nanmean(p2vmreduced.oi_vis_ft_visdata_s[baseline::nbase,:], axis=1) * np.exp(-1j * p2vmreduced.oi_vis_ft_target_phase_s[:,baseline]) demod_p = np.nanmean(p2vmreduced.oi_vis_ft_visdata_p[baseline::nbase,:], axis=1) * np.exp(-1j * p2vmreduced.oi_vis_ft_target_phase_p[:,baseline]) visphi_s[baseline,:] = np.angle(demod_s * np.conj(np.nanmean(demod_s))) * p2vmreduced.wave_ft.mean() / (2.*np.pi) visphi_p[baseline,:] = np.angle(demod_p * np.conj(np.nanmean(demod_p))) * p2vmreduced.wave_ft.mean() / (2.*np.pi) PSD_s=[1,2,3,4,5,6] PSD_p=[1,2,3,4,5,6] revcum_s = [1,2,3,4,5,6] revcum_p = [1,2,3,4,5,6] for baseline in range(0,nbase): f_s, PSD_s[baseline] = welch(visphi_s[baseline,:], fs=(1./np.mean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') f_p, PSD_p[baseline] = welch(visphi_p[baseline,:], fs=(1./np.mean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') PSD_s[baseline] -= np.mean(PSD_s[baseline][-100:]) PSD_p[baseline] -= np.mean(PSD_p[baseline][-100:]) revcum_s[baseline] = PSD_s[baseline][::-1].cumsum()[::-1] revcum_p[baseline] = PSD_p[baseline][::-1].cumsum()[::-1] yminval = 0. ymaxval = 0.1 ticky = (ymaxval-yminval)/5 hsize = 16*cm vsize = 6.5*cm for baseline in range(0,nbase): data = [] data.append(tuple(zip(f_s, clipdata((PSD_s[baseline]+PSD_p[baseline])/2.0, yminval,ymaxval)))) data.append(tuple(zip(f_s, clipdata((revcum_s[baseline]+revcum_p[baseline])/2.0, yminval, ymaxval)))) a = graphoutaxes( data, np.nanmin(f_s),np.nanmax(f_s),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ 'Baseline {0} ({1:.0f} nm rms)'.format(p2vmreduced.basenames[baseline], 1e3*visphi_s[baseline,:].std())) Story.append(a) Story.append(Spacer(1,7*mm)) if (p2vmreduced.polarsplit == False): visphi = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') for baseline in range(0,nbase): demod = np.nanmean(p2vmreduced.oi_vis_ft_visdata[baseline::nbase,:], axis=1) * np.exp(-1j * p2vmreduced.oi_vis_ft_target_phase[:,baseline]) visphi[baseline,:] = np.angle(demod * np.conj(np.nanmean(demod))) * p2vmreduced.wave_ft.mean() / (2.*np.pi) PSD=[1,2,3,4,5,6] revcum = [1,2,3,4,5,6] for baseline in range(0,nbase): f, PSD[baseline] = welch(visphi[baseline,:], fs=(1./np.mean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') PSD[baseline] -= np.mean(PSD[baseline][-100:]) revcum[baseline] = PSD[baseline][::-1].cumsum()[::-1] yminval = 0. ymaxval = 0.1 ticky = (ymaxval-yminval)/5 hsize = 16*cm vsize = 6*cm for baseline in range(0,nbase): data = [] data.append(tuple(zip(f, clipdata(PSD[baseline], yminval,ymaxval)))) data.append(tuple(zip(f, clipdata(revcum[baseline], yminval, ymaxval)))) a = graphoutaxes( data, np.nanmin(f),np.nanmax(f),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ 'Baseline {0} ({1:.0f} nm rms)'.format(p2vmreduced.basenames[baseline], 1e3*visphi[baseline,:].std())) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Power spectral density of the FT unwrapped phase vs. time #============================================================================== plotTitle(Story,"FT OPD periodogram - only good FT frames (OI_VIS FT)") plotSubtitle(Story,"In um^2 vs. Hz., power spectrum in red, reverse cumulated in orange.") if welchpresent == True: if (p2vmreduced.polarsplit == True): visphi_s = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') visphi_p = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') visphi_std = np.zeros((nbase),dtype='d') for baseline in range(0,nbase): demod_s = np.nanmean(p2vmreduced.oi_vis_ft_visdata_s[baseline::nbase,:], axis=1)\ * np.exp(-1j * p2vmreduced.oi_vis_ft_target_phase_s[:,baseline]) demod_p = np.nanmean(p2vmreduced.oi_vis_ft_visdata_p[baseline::nbase,:], axis=1)\ * np.exp(-1j * p2vmreduced.oi_vis_ft_target_phase_p[:,baseline]) visphi_s[baseline,:] = np.angle(demod_s * np.conj(np.nanmean(demod_s)))\ * p2vmreduced.wave_ft.mean() / (2.*np.pi) visphi_p[baseline,:] = np.angle(demod_p * np.conj(np.nanmean(demod_p)))\ * p2vmreduced.wave_ft.mean() / (2.*np.pi) goodframes = np.squeeze(np.where((p2vmreduced.oi_vis_ft_reject_flag_s[:,baseline]==0) * (p2vmreduced.oi_vis_ft_reject_flag_p[:,baseline]==0))) if len(goodframes) > 100: PSD_s=[1,2,3,4,5,6] PSD_p=[1,2,3,4,5,6] revcum_s = [1,2,3,4,5,6] revcum_p = [1,2,3,4,5,6] for baseline in range(0,nbase): visphi_s_ok = visphi_s[baseline,goodframes] visphi_p_ok = visphi_p[baseline,goodframes] visphi_std[baseline] = visphi_s_ok.std() # the std is the same for both polar f_s, PSD_s[baseline] = welch(visphi_s_ok, fs=(1./np.nanmean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') f_p, PSD_p[baseline] = welch(visphi_p_ok, fs=(1./np.nanmean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') PSD_s[baseline] -= np.nanmean(PSD_s[baseline][-100:]) PSD_p[baseline] -= np.nanmean(PSD_p[baseline][-100:]) revcum_s[baseline] = PSD_s[baseline][::-1].cumsum()[::-1] revcum_p[baseline] = PSD_p[baseline][::-1].cumsum()[::-1] yminval = 0. ymaxval = 0.1 ticky = (ymaxval-yminval)/5 hsize = 16*cm vsize = 6.5*cm for baseline in range(0,nbase): data = [] data.append(tuple(zip(f_s, clipdata((PSD_s[baseline]+PSD_p[baseline])/2.0, yminval,ymaxval)))) data.append(tuple(zip(f_s, clipdata((revcum_s[baseline]+revcum_p[baseline])/2.0, yminval, ymaxval)))) a = graphoutaxes( data, np.nanmin(f_s),np.nanmax(f_s),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ 'Baseline {0} ({1:.0f} nm rms)'.format(p2vmreduced.basenames[baseline], 1e3*visphi_std[baseline])) Story.append(a) Story.append(Spacer(1,7*mm)) else: plotSubtitle(Story,"Not enough good FT frames to compute PSD only for good frames.") if (p2vmreduced.polarsplit == False): visphi = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') visphi_std = np.zeros((nbase),dtype='d') for baseline in range(0,nbase): demod = np.nanmean(p2vmreduced.oi_vis_ft_visdata[baseline::nbase,:], axis=1) * np.exp(-1j * p2vmreduced.oi_vis_ft_target_phase[:,baseline]) visphi[baseline,:] = np.angle(demod * np.conj(np.nanmean(demod))) * p2vmreduced.wave_ft.mean() / (2.*np.pi) goodframes = np.squeeze(np.where(p2vmreduced.oi_vis_ft_reject_flag[:,baseline]==0)) if len(goodframes) > 100: PSD=[1,2,3,4,5,6] revcum = [1,2,3,4,5,6] for baseline in range(0,nbase): visphi_ok = visphi[baseline,goodframes] visphi_std[baseline] = visphi_ok.std() f, PSD[baseline] = welch(visphi_ok, fs=(1./np.nanmean(np.diff(p2vmreduced.time_ft))), detrend='linear', nperseg=1024, scaling='spectrum') PSD[baseline] -= np.nanmean(PSD[baseline][-100:]) revcum[baseline] = PSD[baseline][::-1].cumsum()[::-1] yminval = 0. ymaxval = 0.1 ticky = (ymaxval-yminval)/5 hsize = 16*cm vsize = 3*cm for baseline in range(0,nbase): data = [] data.append(tuple(zip(f, clipdata(PSD[baseline], yminval,ymaxval)))) data.append(tuple(zip(f, clipdata(revcum[baseline], yminval, ymaxval)))) a = graphoutaxes( data, np.nanmin(f),np.nanmax(f),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ 'Baseline {0} ({1:.0f} nm rms)'.format(p2vmreduced.basenames[baseline], 1e3*visphi_std[baseline])) Story.append(a) Story.append(Spacer(1,7*mm)) else: plotSubtitle(Story,"Not enough good FT frames to compute PSD only for good frames.") Story.append(PageBreak()) #============================================================================== # OPDC signals #============================================================================== plotTitle(Story,"OPDC global signals (OPDC)") if p2vmreduced.gdelay_in == True: resamp = int(p2vmreduced.opdc_time.shape[0]/1000) # take one value every resamp value for the plot hsize = 16*cm vsize = 5.5*cm piezo_time_plot = p2vmreduced.opdc_time[::resamp] state_plot = p2vmreduced.opdc_state[::resamp] fraction_pd = 100.*(np.sum(p2vmreduced.opdc_state==3)/(1.0*p2vmreduced.opdc_state.shape[0])) fraction_gd = 100.*(np.sum(p2vmreduced.opdc_state==2)/(1.0*p2vmreduced.opdc_state.shape[0])) style = styles["Normal"] plotSubtitle(Story,"OPDC state : Phase tracking = %3.1f %% / Group delay tracking = %3.1f %% / Total tracking = %3.1f %%." %(fraction_pd,fraction_gd,fraction_pd+fraction_gd)) data = [] yminval = 0 ymaxval = 11 ticky = 1 data.append( tuple(zip(piezo_time_plot, clipdata(state_plot,yminval,ymaxval)) )) a = graphoutaxes( data, min(piezo_time_plot),max(piezo_time_plot),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'') a.add(String(0.05*hsize, 0.10*vsize, '1=IDLE', fontSize=9)) a.add(String(0.05*hsize, 0.19*vsize, '2=GD_TRACKING (at least one tel)', fontSize=9)) a.add(String(0.05*hsize, 0.28*vsize, '3=PHASE_TRACKING', fontSize=9)) a.add(String(0.05*hsize, 0.37*vsize, '4=SEARCHING (at least one tel)', fontSize=9)) a.add(String(0.05*hsize, 0.46*vsize, '5=SCAN_SEARCH (with VLTI DL)', fontSize=9)) a.add(String(0.05*hsize, 0.55*vsize, '6=CALIB_FT (calibrate piezo)', fontSize=9)) a.add(String(0.05*hsize, 0.64*vsize, '7=CALIB_SNR', fontSize=9)) a.add(String(0.05*hsize, 0.73*vsize, '10=SKY (no flux)', fontSize=9)) a.add(String(0.05*hsize, 0.82*vsize, '11=PSF_SCAN (fiber-injection map)', fontSize=9)) Story.append(a) Story.append(Spacer(1,10*mm)) plotSubtitle(Story,"Piezo offsets (Volts). "+\ p2vmreduced.stations[0] +" = red (gv beam 1), "+\ p2vmreduced.stations[1] +" = orange (gv beam 2), "+\ p2vmreduced.stations[2] +" = blue (gv beam 3), "+\ p2vmreduced.stations[3] +" = green (gv beam 4).") data = [] yminval = np.nanmin(p2vmreduced.opdc_piezo_offset)-1 ymaxval = np.nanmax(p2vmreduced.opdc_piezo_offset)+1 ticky = (ymaxval-yminval)/10. # in microns for tel in range(0,ntel): data.append( tuple(zip(piezo_time_plot, clipdata(p2vmreduced.opdc_piezo_offset[::resamp,tel],yminval,ymaxval)) )) a = graphoutaxes( data, np.nanmin(piezo_time_plot),np.nanmax(piezo_time_plot),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'') Story.append(a) Story.append(Spacer(1,7*mm)) value = p2vmreduced.opdc_vltidl_offset[::resamp,:] * 1e3 #print(value.shape) means = np.nanmean(value,0) #print(means.shape) value -= means[None,:] plotSubtitle(Story,"VLTI DL offsets (mm). means = %.2f %.2f %.2f %.2f [mm]"%(means[0],means[1],means[2],means[3])); data = [] yminval = np.nanmin(value)-0.05 ymaxval = np.nanmax(value)+0.05 ticky = (ymaxval-yminval)/10. # in mm for tel in range(0,ntel): data.append( tuple(zip(piezo_time_plot, clipdata(value[:,tel],yminval,ymaxval)) )) a = graphoutaxes( data, np.nanmin(piezo_time_plot),np.nanmax(piezo_time_plot),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'') Story.append(a) Story.append(PageBreak()) #============================================================================== # FT OPDC STATE #============================================================================== plotTitle(Story,"FT OPDC STATE per baseline") plotSubtitle(Story,"computed from OPDC BASE_STATE") resamp = nframeplot // 500 # take one value every resamp value for the plot yminval = -2 ymaxval = 6 ticky = 1 # visibility tick xval = p2vmreduced.time_ft[::resamp] xminval = np.min(xval) xmaxval = np.max(xval) d=() hsize=16*cm vsize=3*cm for baseline in range(0,nbase): data=[] if p2vmreduced.polarsplit == False: data.append( list(zip(xval, clipdata(p2vmreduced.oi_vis_ft_state[::resamp,baseline],yminval,ymaxval))) ) else: data.append( list(zip(xval, clipdata(p2vmreduced.oi_vis_ft_state_s[::resamp,baseline],yminval,ymaxval))) ) data.append( list(zip(xval, clipdata(p2vmreduced.oi_vis_ft_state_p[::resamp,baseline]+0.05,yminval,ymaxval))) ) a = graphoutaxes( data, xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # FT REJECTION_FLAG #============================================================================== plotTitle(Story,"FT REJECTION_FLAG per baseline") plotSubtitle(Story,"0= accepted frame, 1= Low SNR, 2= Low STATE, 3= both") resamp = nframeplot // 500 # take one value every resamp value for the plot yminval = -1 ymaxval = 4 ticky = 1 # visibility tick xval = p2vmreduced.time_ft[::resamp] xminval = np.min(xval) xmaxval = np.max(xval) d=() hsize=16*cm vsize=3*cm for baseline in range(0,nbase): data=[] if p2vmreduced.polarsplit == False: data.append (list(zip(xval, p2vmreduced.oi_vis_ft_reject_flag[::resamp,baseline]))) else: data.append (list(zip(xval, p2vmreduced.oi_vis_ft_reject_flag_s[::resamp,baseline]))) data.append (list(zip(xval, p2vmreduced.oi_vis_ft_reject_flag_p[::resamp,baseline]+0.05))) a = graphoutaxes (data, xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the FT group delay signals from PIPELINE #============================================================================== plotTitle(Story, "FT group delay vs. time (GDELAY/GDELAY_BOOT)") plotSubtitle(Story, "OPD (microns) vs. time (seconds). Red= GD, Orange = GDB.") if p2vmreduced.gdelay_in == True: resamp = nframeplot // 500 # take one value every resamp value for the plot timescale_ft = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') timescale_ft = p2vmreduced.time_ft hsize=16*cm vsize=3*cm for baseline in range(0,nbase): data = [] if (p2vmreduced.polarsplit == False): yminval = np.nanpercentile(p2vmreduced.gdelay_ft[baseline,0:nframeplot],1) ymaxval = np.nanpercentile(p2vmreduced.gdelay_ft[baseline,0:nframeplot],99) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.gdelay_ft[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.gdelay_boot_ft[baseline,0:nframeplot:resamp],yminval,ymaxval)))) else: yminval = np.nanpercentile([p2vmreduced.gdelay_ft_s[baseline,0:nframeplot],p2vmreduced.gdelay_ft_p[baseline,0:nframeplot]],1)-2 ymaxval = np.nanpercentile([p2vmreduced.gdelay_ft_s[baseline,0:nframeplot],p2vmreduced.gdelay_ft_p[baseline,0:nframeplot]],99)+2 data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.gdelay_ft_s[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.gdelay_ft_p[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.gdelay_boot_ft_s[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.gdelay_boot_ft_p[baseline,0:nframeplot:resamp],yminval,ymaxval)))) ticky = (ymaxval-yminval)/5 a = graphoutaxes( data, np.nanmin(timescale_ft[0:nframeplot]),np.nanmax(timescale_ft[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the FT phase signals #============================================================================== plotTitle(Story, "FT phase vs. time (phase of wavelength averaged OI_VIS)") plotSubtitle(Story, "Phase (deg) vs. time (seconds). red=S and orange=P") if p2vmreduced.gdelay_in == True: resamp = nframeplot // 1000 # take one value every resamp value for the plot timescale_ft = p2vmreduced.time_ft[0::resamp] xminval = np.min(timescale_ft) xmaxval = np.max(timescale_ft) yminval=-180 ymaxval=+180 ticky = (ymaxval-yminval)/2 hsize=16*cm vsize=3*cm boxsmooth_ft = 2 for baseline in range(0,nbase): data = [] if (p2vmreduced.polarsplit == False): data.append(tuple(zip(timescale_ft, np.angle(np.nanmean(p2vmreduced.oi_vis_ft_visdata[baseline::nbase*resamp,:], axis=1))*180./np.pi))) else: data.append(tuple(zip(timescale_ft, np.angle(np.nanmean(p2vmreduced.oi_vis_ft_visdata_s[baseline::nbase*resamp,:], axis=1))*180./np.pi))) data.append(tuple(zip(timescale_ft, np.angle(np.nanmean(p2vmreduced.oi_vis_ft_visdata_p[baseline::nbase*resamp,:], axis=1))*180./np.pi))) a = graphoutaxes( data,xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the demodulated FT phase signals #============================================================================== plotTitle(Story, "Demodulated FT phase vs. time (VISPHI - TARGET_PHASE)") plotSubtitle(Story, "Phase (deg) vs. time (seconds). red=S and orange=P, mean subtracted") if p2vmreduced.gdelay_in == True: resamp = nframeplot // 1000 # take one value every resamp value for the plot timescale_ft = p2vmreduced.time_ft[0::resamp] xminval = np.min(timescale_ft) xmaxval = np.max(timescale_ft) yminval=-180 ymaxval=+180 ticky = (ymaxval-yminval)/2 hsize=16*cm vsize=3*cm boxsmooth_ft = 2 rmsjitter = np.zeros(6) for baseline in range(0,nbase): data = [] if (p2vmreduced.polarsplit == False): tmp = np.nanmean(p2vmreduced.oi_vis_ft_visdata[baseline::nbase*resamp,:], axis=1) * np.exp(-1.J * p2vmreduced.oi_vis_ft_target_phase[::resamp,baseline]) demod = np.angle(tmp * np.conj(np.nanmean(tmp)))*180./np.pi rmsjitter[baseline] = np.std(demod)*np.mean(p2vmreduced.wave_ft)/360. data.append(tuple(zip(timescale_ft, demod))) else: tmp = np.nanmean(p2vmreduced.oi_vis_ft_visdata_s[baseline::nbase*resamp,:], axis=1) * np.exp(-1.J * p2vmreduced.oi_vis_ft_target_phase_s[::resamp,baseline]) demod_s = np.angle(tmp * np.conj(np.nanmean(tmp)))*180./np.pi data.append(tuple(zip(timescale_ft, demod_s))) tmp = np.nanmean(p2vmreduced.oi_vis_ft_visdata_p[baseline::nbase*resamp,:], axis=1) * np.exp(-1.J * p2vmreduced.oi_vis_ft_target_phase_p[::resamp,baseline]) demod_p = np.angle(tmp * np.conj(np.nanmean(tmp)))*180./np.pi data.append(tuple(zip(timescale_ft, demod_p))) rmsjitter[baseline] = np.std([demod_p,demod_s])*np.mean(p2vmreduced.wave_ft)/360. a = graphoutaxes(data,xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+' RMS jitter = %(jitter).3f mic'%{"jitter":rmsjitter[baseline]}) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the FT coherent flux SNR signals #============================================================================== plotTitle(Story, "FT coherent flux SNR vs. time (OI_VIS/SNR)") if p2vmreduced.gdelay_in == True: resamp = nframeplot // 500 # take one value every resamp value for the plot timescale_ft = p2vmreduced.time_ft d=() hsize=16*cm vsize=3*cm if (p2vmreduced.polarsplit == False): plotSubtitle(Story,"SNR coh. flux vs. time (seconds). Red = real-time, Orange = bootstr.") for baseline in range(0,nbase): data = [] yminval = 0 ymaxval = np.nanmax([p2vmreduced.snr_ft[baseline,0:nframeplot],p2vmreduced.snr_boot_ft[baseline,0:nframeplot]]) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.snr_ft[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.snr_boot_ft[baseline,0:nframeplot:resamp],yminval,ymaxval)))) ticky = (ymaxval-yminval)/5 a = graphoutaxes(data,np.nanmin(timescale_ft[0:nframeplot]),np.nanmax(timescale_ft[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if (p2vmreduced.polarsplit == True): plotSubtitle(Story,"SNR vs. time (seconds). Red = S, Orange = P, Blue = bootstr. S, Green = bootstr. P.") for baseline in range(0,nbase): data = [] yminval = 0 ymaxval = np.nanmax([p2vmreduced.snr_ft_s[baseline,0:nframeplot],\ p2vmreduced.snr_ft_p[baseline,0:nframeplot],\ p2vmreduced.snr_boot_ft_s[baseline,0:nframeplot],\ p2vmreduced.snr_boot_ft_p[baseline,0:nframeplot]]) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.snr_ft_s[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.snr_ft_p[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.snr_boot_ft_s[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], clipdata(p2vmreduced.snr_boot_ft_p[baseline,0:nframeplot:resamp],yminval,ymaxval)))) ticky = (ymaxval-yminval)/5 a = graphoutaxes( data, np.nanmin(timescale_ft[0:nframeplot]), np.nanmax(timescale_ft[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the Expected coherent flux signals for FT #============================================================================== plotTitle(Story,"FT normalisation vs. time (OI_FLUX, wavelength averaged)") plotSubtitle(Story,"sqrt(F1*F2) vs. time (seconds). Red = S, Orange = P") if p2vmreduced.gdelay_in == True: resamp = nframeplot // 500 # take one value every resamp value for the plot hsize=16*cm vsize=3*cm timescale_ft = p2vmreduced.time_ft for baseline in range(0,nbase): data = [] yminval = 0 if (p2vmreduced.polarsplit == False): f1f2 = np.sqrt( np.maximum( np.sum(p2vmreduced.oi_flux_ft[:,tel_list[baseline,0],:],axis=1) * np.sum(p2vmreduced.oi_flux_ft[:,tel_list[baseline,1],:],axis=1), 0.0 ) ) ymaxval = np.nanmax(f1f2) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], f1f2[0:nframeplot:resamp]))) else: f1f2_s = np.sqrt( np.maximum( np.sum(p2vmreduced.oi_flux_ft_s[:,tel_list[baseline,0],:],axis=1)*np.sum(p2vmreduced.oi_flux_ft_s[:,tel_list[baseline,1],:],axis=1), 0.0 )) f1f2_p = np.sqrt( np.maximum( np.sum(p2vmreduced.oi_flux_ft_p[:,tel_list[baseline,0],:],axis=1)*np.sum(p2vmreduced.oi_flux_ft_p[:,tel_list[baseline,1],:],axis=1), 0.0 )) ymaxval = np.nanmax(f1f2_s) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], f1f2_s[0:nframeplot:resamp]))) data.append(tuple(zip(timescale_ft[0:nframeplot:resamp], f1f2_p[0:nframeplot:resamp]))) ticky = (ymaxval-yminval)/5 a = graphoutaxes( data, np.nanmin(timescale_ft[0:nframeplot:resamp]), np.nanmax(timescale_ft[0:nframeplot:resamp]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # SC REJECTION_FLAG #============================================================================== plotTitle(Story,"SC REJECTION_FLAG per baseline") plotSubtitle(Story,"0= accepted frame, 1= Low FT tracking, 2= Low FT vFactor, 3= both") resamp = 1 # take one value every resamp value for the plot yminval = -1 ymaxval = 4 ticky = 1 # visibility tick xval = p2vmreduced.time_sc[::resamp] xminval = np.min(xval) xmaxval = np.max(xval) d=() hsize=16*cm vsize=3*cm for baseline in range(0,nbase): data=[] if p2vmreduced.polarsplit == False: data.append (list(zip(xval, p2vmreduced.oi_vis_reject_flag[::resamp,baseline]))) else: data.append (list(zip(xval, p2vmreduced.oi_vis_reject_flag_s[::resamp,baseline]))) data.append (list(zip(xval, p2vmreduced.oi_vis_reject_flag_p[::resamp,baseline]+0.05))) a = graphoutaxes (data, xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) #a.add(String(0.05*hsize, 0.22*vsize, '0=Good', fontSize=9)) #a.add(String(0.05*hsize, 0.42*vsize, '1=low FT tracking', fontSize=9)) #a.add(String(0.05*hsize, 0.62*vsize, '2=low FT Vfactor', fontSize=9)) #a.add(String(0.05*hsize, 0.82*vsize, '3=low FT Vfactor and low FT tracking', fontSize=9)) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the SC group delay signals #============================================================================== plotTitle(Story,"SC group delay vs. time") plotSubtitle(Story,"In OPD (microns) vs. time (seconds). Red = S, Orange = P.") if p2vmreduced.gdelay_in == True: timescale_sc = p2vmreduced.time_sc hsize=16*cm vsize=3*cm for baseline in range(0,nbase): if (p2vmreduced.polarsplit == False): tmp1 = p2vmreduced.gdelay_sc[baseline,:] yminval = np.nanmin(tmp1) ymaxval = np.nanmax(tmp1) data = [tuple(zip(timescale_sc[:], clipdata(tmp1,yminval,ymaxval)))] else: tmp1 = p2vmreduced.gdelay_sc_s[baseline,:] tmp2 = p2vmreduced.gdelay_sc_p[baseline,:] yminval = np.nanmin(tmp1) ymaxval = np.nanmax(tmp1) data = [tuple(zip(timescale_sc[:], clipdata(tmp1,yminval,ymaxval)))] data.append(tuple(zip(timescale_sc[:], clipdata(tmp2,yminval,ymaxval)))) ticky = (ymaxval-yminval)/5 yminval -= ticky ymaxval += ticky a = graphoutaxes( data, np.nanmin(timescale_sc[:]), np.nanmax(timescale_sc[:]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Plotting of the SC coherent flux SNR signals from PIPELINE #============================================================================== plotTitle(Story,"SC coherent flux SNR vs. time (OI_VIS/SNR)") plotSubtitle(Story,"SNR vs. time (seconds). Red = S, Orange = P.") if p2vmreduced.gdelay_in == True: timescale_sc = p2vmreduced.time_sc d=() data=[] hsize=16*cm vsize=3*cm for baseline in range(0,nbase): data = [] yminval = 0 if (p2vmreduced.polarsplit == False): ymaxval = np.nanpercentile(p2vmreduced.snr_sc[baseline,0:nframeplot],98) data.append(tuple(zip(timescale_sc[0:nframeplot], clipdata(p2vmreduced.snr_sc[baseline,0:nframeplot],yminval,ymaxval)))) else: ymaxval = np.nanpercentile(p2vmreduced.snr_sc_s[baseline,0:nframeplot],98) data.append(tuple(zip(timescale_sc[0:nframeplot], clipdata(p2vmreduced.snr_sc_s[baseline,0:nframeplot],yminval,ymaxval)))) data.append(tuple(zip(timescale_sc[0:nframeplot], clipdata(p2vmreduced.snr_sc_p[baseline,0:nframeplot],yminval,ymaxval)))) ticky = (ymaxval-yminval)/5 yminval -= ticky ymaxval += ticky a = graphoutaxes( data, np.nanmin(timescale_sc[0:nframeplot]), np.nanmax(timescale_sc[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Average V factor (quadratic quantity) computed for SC vs. wavelength #============================================================================== if hasattr(p2vmreduced,'oi_vis_vfactor_ft') or hasattr(p2vmreduced,'oi_vis_vfactor_ft_s'): plotTitle(Story,"V_FACTOR and V_FACTOR_FT vs. wavelength (averaged over frames)") if (p2vmreduced.polarsplit == False): plotSubtitle(Story,"Wavelength in microns, red = SC, orange = FT.") else: plotSubtitle(Story,"Wavelength in microns, red = SC S, orange = SC P, blue = FT S, green = FT P.") hsize = 6*cm vsize = 6*cm yminval = 0.0 ymaxval = 1.2 xminval = np.min(p2vmreduced.wave_sc) xmaxval = np.max(p2vmreduced.wave_sc) ticky = 0.25 tickx = 0.25 d = [] for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength] data=[] if (p2vmreduced.polarsplit == False): data.append( tuple(zip(p2vmreduced.wave_sc, np.mean(p2vmreduced.oi_vis_vfactor[:,baseline,:],axis=0)) )) data.append( tuple(zip(p2vmreduced.wave_ft, np.mean(p2vmreduced.oi_vis_vfactor_ft[:,baseline,:],axis=0)) )) else: data.append( tuple(zip(p2vmreduced.wave_sc, np.mean(p2vmreduced.oi_vis_vfactor_s[:,baseline,:],axis=0)) )) data.append( tuple(zip(p2vmreduced.wave_sc, np.mean(p2vmreduced.oi_vis_vfactor_p[:,baseline,:],axis=0)) )) data.append( tuple(zip(p2vmreduced.wave_ft, np.mean(p2vmreduced.oi_vis_vfactor_ft_s[:,baseline,:],axis=0)) )) data.append( tuple(zip(p2vmreduced.wave_ft, np.mean(p2vmreduced.oi_vis_vfactor_ft_p[:,baseline,:],axis=0)) )) a = graphoutaxes( data,xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ] t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of V_FACTOR_FT #============================================================================== plotTitle(Story,"V_FACTOR_FT waterfall vs. wavelength and time") plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of SC frame in vertical axis. Mean spectrum over time subtracted.") if p2vmreduced.polarsplit == False: meanv = np.zeros((nbase,p2vmreduced.nwave_ft)) for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength] meanv[baseline,:] = np.nanmean(p2vmreduced.oi_vis_vfactor_ft[:,baseline,:],axis=0) meansub = np.zeros((nbase,p2vmreduced.nframe_sc,p2vmreduced.nwave_ft)) for baseline in range(0,nbase): for frame in range(0,p2vmreduced.nframe_sc): meansub[baseline,frame,:] = p2vmreduced.oi_vis_vfactor_ft[frame,baseline,:] - meanv[baseline,:] valmin = np.nanpercentile(meansub,1)-0.05 valmax = np.nanpercentile(meansub,99)+0.05 imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=nbase, ncols=1, sharex=True, sharey=True, figsize=(5,8)) for baseline in range(0,nbase): im = axarr[baseline].imshow(meansub[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline].set_title('Baseline '+p2vmreduced.basenames[baseline], fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if p2vmreduced.polarsplit == True: meanvs = np.zeros((nbase,p2vmreduced.nwave_ft)) meanvp = np.zeros((nbase,p2vmreduced.nwave_ft)) for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength] meanvs[baseline,:] = np.nanmean(p2vmreduced.oi_vis_vfactor_ft_s[:,baseline,:],axis=0) meanvp[baseline,:] = np.nanmean(p2vmreduced.oi_vis_vfactor_ft_p[:,baseline,:],axis=0) meansubs = np.zeros((nbase,p2vmreduced.nframe_sc,p2vmreduced.nwave_ft)) meansubp = np.zeros((nbase,p2vmreduced.nframe_sc,p2vmreduced.nwave_ft)) for baseline in range(0,nbase): for frame in range(0,p2vmreduced.nframe_sc): meansubs[baseline,frame,:] = p2vmreduced.oi_vis_vfactor_ft_s[frame,baseline,:] - meanvs[baseline,:] meansubp[baseline,frame,:] = p2vmreduced.oi_vis_vfactor_ft_p[frame,baseline,:] - meanvp[baseline,:] valmin = np.nanpercentile([meansubs,meansubp],1)-0.05 valmax = np.nanpercentile([meansubs,meansubp],99)+0.05 imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=nbase, ncols=2, sharex=True, sharey=True, figsize=(5,8)) for baseline in range(0,nbase): im = axarr[baseline,0].imshow(meansubs[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,0].set_title('Baseline '+p2vmreduced.basenames[baseline]+'-S', fontsize=7) im = axarr[baseline,1].imshow(meansubp[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,1].set_title('Baseline '+p2vmreduced.basenames[baseline]+'-P', fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, 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 model V_FACTOR (of SC) #============================================================================== plotTitle(Story,"V_FACTOR_SC waterfall vs. wavelength and time") plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of SC frame in vertical axis. Mean spectrum over time subtracted.") if p2vmreduced.polarsplit == False: meanv = np.zeros((nbase,p2vmreduced.nwave_sc)) for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength] meanv[baseline,:] = np.nanmean(p2vmreduced.oi_vis_vfactor[:,baseline,:],axis=0) meansub = np.zeros((nbase,p2vmreduced.nframe_sc,p2vmreduced.nwave_sc)) for baseline in range(0,nbase): for frame in range(0,p2vmreduced.nframe_sc): meansub[baseline,frame,:] = p2vmreduced.oi_vis_vfactor[frame,baseline,:] - meanv[baseline,:] valmin = np.nanpercentile(meansub,1)-0.05 valmax = np.nanpercentile(meansub,99)+0.05 #print(valmin, valmax) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=nbase, ncols=1, sharex=True, sharey=True, figsize=(5,8)) for baseline in range(0,nbase): im = axarr[baseline].imshow(meansub[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline].set_title('Baseline '+p2vmreduced.basenames[baseline], fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if p2vmreduced.polarsplit == True: meanvs = np.zeros((nbase,p2vmreduced.nwave_sc)) meanvp = np.zeros((nbase,p2vmreduced.nwave_sc)) for baseline in range(0,nbase): # V_FACTOR[frame,baseline,wavelength] meanvs[baseline,:] = np.nanmean(p2vmreduced.oi_vis_vfactor_s[:,baseline,:],axis=0) meanvp[baseline,:] = np.nanmean(p2vmreduced.oi_vis_vfactor_p[:,baseline,:],axis=0) meansubs = np.zeros((nbase,p2vmreduced.nframe_sc,p2vmreduced.nwave_sc)) meansubp = np.zeros((nbase,p2vmreduced.nframe_sc,p2vmreduced.nwave_sc)) for baseline in range(0,nbase): for frame in range(0,p2vmreduced.nframe_sc): meansubs[baseline,frame,:] = p2vmreduced.oi_vis_vfactor_s[frame,baseline,:] - meanvs[baseline,:] meansubp[baseline,frame,:] = p2vmreduced.oi_vis_vfactor_p[frame,baseline,:] - meanvp[baseline,:] valmin = np.nanpercentile([meansubs,meansubp],1)-0.05 valmax = np.nanpercentile([meansubs,meansubp],99)+0.05 #print(valmin, valmax) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=nbase, ncols=2, sharex=True, sharey=True, figsize=(5,8)) for baseline in range(0,nbase): im = axarr[baseline,0].imshow(meansubs[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,0].set_title('Baseline '+p2vmreduced.basenames[baseline]+'-S', fontsize=7) im = axarr[baseline,1].imshow(meansubp[baseline,:,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,1].set_title('Baseline '+p2vmreduced.basenames[baseline]+'-P', fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # P_FACTOR versus time #============================================================================== if hasattr(p2vmreduced,'oi_vis_pfactor') or hasattr(p2vmreduced,'oi_vis_pfactor_s'): plotTitle(Story,"P_FACTOR for SC versus time") plotSubtitle(Story,"The P_FACTOR is computed broad-band") timescale_sc = p2vmreduced.time_sc d=() data=[] hsize=16*cm vsize=3*cm yminval = -0.1 ymaxval = 1.1 ticky = 0.2 for base in range(0,nbase): data = [] if (p2vmreduced.polarsplit == True): data.append(tuple(zip(timescale_sc, clipdata(p2vmreduced.oi_vis_pfactor_s[base::nbase],yminval,ymaxval)))) data.append(tuple(zip(timescale_sc, clipdata(p2vmreduced.oi_vis_pfactor_p[base::nbase],yminval,ymaxval)))) else: data.append(tuple(zip(timescale_sc, clipdata(p2vmreduced.oi_vis_pfactor[base::nbase],yminval,ymaxval)))) a = graphoutaxes( data, np.nanmin(timescale_sc), np.nanmax(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[base]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # V factor (quadratic quantity) computed for SC versus VIS2 #============================================================================== if True: if hasattr(p2vmreduced,'oi_vis_vfactor_ft') or hasattr(p2vmreduced,'oi_vis_vfactor_ft_s'): plotTitle(Story,"SC VIS2 versus V_FACTOR (OI_VIS2, wavelength averaged)") plotSubtitle(Story,"VIS2 (vertical) versus V_FACTOR (horizontal), Red = S or combined.") hsize = 6*cm vsize = 6*cm ymaxval = 1.2 yminval = 0.0 ticky = 0.25 tickx = 0.25 d = [] for baseline in range(0,nbase): data=[] if (p2vmreduced.polarsplit == False): data.append( tuple(zip(clipdata(np.mean(p2vmreduced.oi_vis_vfactor[:,baseline,:],axis=1),yminval,ymaxval), \ clipdata(np.mean(p2vmreduced.oi_vis2_sc[baseline,:,:],axis=1),yminval,ymaxval) ))) else: data.append( tuple(zip(clipdata(np.mean(p2vmreduced.oi_vis_vfactor_s[:,baseline,:],axis=1),yminval,ymaxval), \ clipdata(np.mean(p2vmreduced.oi_vis2_sc_s[baseline,:,:],axis=1),yminval,ymaxval) ))) a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+p2vmreduced.basenames[baseline]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ] t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # VIS2 * F1 * F2 vs. F1 * F2 of FT #============================================================================== plotTitle(Story,"FT VIS2_F1F2 vs. F1*F2 (OI_VIS, wavelength averaged)") plotSubtitle(Story,"In photo-e-^2/FT DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.") resamp = 100 hsize = 6*cm vsize = 6*cm d = [] for baseline in range(0,nbase): if (p2vmreduced.polarsplit == True): f1f2 = np.mean(p2vmreduced.oi_flux_ft_s[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_ft_s[::resamp,tel_list[baseline,1],:],axis=1) vis2 = np.mean(p2vmreduced.oi_vis2_ft_s[baseline,::resamp,:],axis=1) else: f1f2 = np.mean(p2vmreduced.oi_flux_ft[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_ft[::resamp,tel_list[baseline,1],:],axis=1) vis2 = np.mean(p2vmreduced.oi_vis2_ft[baseline,::resamp,:],axis=1) ymaxval = 1.05*np.max(f1f2) yminval = -0.05*ymaxval ticky = (ymaxval-yminval)/3. tickx = ticky data=[] data.append( tuple(zip(f1f2, clipdata(vis2*f1f2,yminval,ymaxval)) )) a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+p2vmreduced.basenames[baseline]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ] t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # VIS2 * F1 * F2 vs. F1 * F2 of SC #============================================================================== plotTitle(Story,"SC VIS2F1F2 vs. F1*F2 (OI_VIS, wavelength averaged)") plotSubtitle(Story,"In photo-e-^2/SC DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.") resamp = 1 hsize = 6*cm vsize = 6*cm d = [] for baseline in range(0,nbase): if (p2vmreduced.polarsplit == True): f1f2 = np.mean(p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],:],axis=1) vis2 = np.mean(p2vmreduced.oi_vis2_sc_s[baseline,::resamp,:],axis=1) else: f1f2 = np.mean(p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,1],:],axis=1) vis2 = np.mean(p2vmreduced.oi_vis2_sc[baseline,::resamp,:],axis=1) ymaxval = 1.05*np.max(f1f2) yminval = -0.05*ymaxval ticky = (ymaxval-yminval)/3. tickx = ticky data=[] data.append( tuple(zip(f1f2, clipdata(vis2*f1f2,yminval,ymaxval)) )) a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+p2vmreduced.basenames[baseline]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ] t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # VIS2 * F1 * F2 vs. F1 * F2 of SC #============================================================================== plotTitle(Story,"SC VIS2F1F2 vs. F1*F2*VFACTOR (OI_VIS, wavelength averaged)") plotSubtitle(Story,"In photo-e-^2/SC DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.") resamp = 1 hsize = 6*cm vsize = 6*cm d = [] for baseline in range(0,nbase): if (p2vmreduced.polarsplit == True): f1f2 = np.mean(p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],:]*p2vmreduced.oi_vis_vfactor_s[::resamp,baseline,:],axis=1) vis2 = np.mean(p2vmreduced.oi_vis2_sc_s[baseline,::resamp,:]*p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],:],axis=1) else: f1f2 = np.mean(p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,1],:]*p2vmreduced.oi_vis_vfactor[::resamp,baseline,:],axis=1) vis2 = np.mean(p2vmreduced.oi_vis2_sc[baseline,::resamp,:]*p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,0],:]*p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,1],:],axis=1) ymaxval = 1.05*np.max(f1f2) yminval = -0.05*ymaxval ticky = (ymaxval-yminval)/3. tickx = ticky data=[] data.append( tuple(zip(f1f2, clipdata(vis2,yminval,ymaxval)) )) a = graphscatteraxes( data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+p2vmreduced.basenames[baseline]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ] t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # VIS2 * F1 * F2 vs. F1 * F2 of SC #============================================================================== plotTitle(Story,"SC VIS2_F1F2 vs. F1*F2 (OI_VIS, 2 channels)") plotSubtitle(Story,"In photo-e-^2/FT DIT, F1*F2 in horizontal axis, VIS2_F1F2 in vertical axis.") resamp = 1 hsize = 6*cm vsize = 6*cm d = [] st1 = p2vmreduced.nwave_sc / 4 sp1 = p2vmreduced.nwave_sc / 4 + 1 st2 = p2vmreduced.nwave_sc / 4 * 3 sp2 = p2vmreduced.nwave_sc / 4 * 3 + 1 for baseline in range(0,nbase): data=[] if (p2vmreduced.polarsplit == True): f1f2a = np.mean(p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],st1:sp1]*p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],st1:sp1],axis=1) vis2a = np.mean(p2vmreduced.oi_vis2_sc_s[baseline,::resamp,st1:sp1],axis=1) f1f2b = np.mean(p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,0],st2:sp2]*p2vmreduced.oi_flux_sc_s[::resamp,tel_list[baseline,1],st2:sp2],axis=1) vis2b = np.mean(p2vmreduced.oi_vis2_sc_s[baseline,::resamp,st2:sp2],axis=1) else: f1f2a = np.mean(p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,0],st1:sp1]*p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,1],st1:sp1],axis=1) vis2a = np.mean(p2vmreduced.oi_vis2_sc[baseline,::resamp,st1:sp1],axis=1) f1f2b = np.mean(p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,0],st2:sp2]*p2vmreduced.oi_flux_sc[::resamp,tel_list[baseline,1],st2:sp2],axis=1) vis2b = np.mean(p2vmreduced.oi_vis2_sc[baseline,::resamp,st2:sp2],axis=1) ymaxval = 1.05*np.max(f1f2a) yminval = -0.05*ymaxval ticky = (ymaxval-yminval)/3. tickx = ticky data.append( tuple(zip(f1f2a, clipdata(vis2a*f1f2a,yminval,ymaxval)) )) data.append( tuple(zip(f1f2b, clipdata(vis2b*f1f2b,yminval,ymaxval)) )) a = graphscatteraxes (data,yminval,ymaxval,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,'Baseline '+p2vmreduced.basenames[baseline]) d.append(a) plotmatrix = [ [d[0],d[1]],[d[2],d[3]],[d[4],d[5]] ] t = Table(plotmatrix,colWidths=hsize+2*cm,rowHeights=vsize+1*cm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of FT VISAMP #============================================================================== plotTitle(Story,"FT visibility amplitude waterfall (VISAMP)") plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.") if hasattr(p2vmreduced,'oi_vis_ft_visamp') & (p2vmreduced.polarsplit == False): valmin = np.max([0,np.nanpercentile(p2vmreduced.oi_vis_ft_visamp,1)]) valmax = np.min([1.3,np.nanpercentile(p2vmreduced.oi_vis_ft_visamp,95)]) 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(p2vmreduced.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(p2vmreduced.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(p2vmreduced.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(p2vmreduced.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(p2vmreduced.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(p2vmreduced.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 = 20 * cm heightfig= 21 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) if hasattr(p2vmreduced,'oi_vis_ft_visamp_s') & (p2vmreduced.polarsplit == True): valmin = np.nanmax([0,np.nanpercentile(p2vmreduced.oi_vis_ft_visamp_s,1)]) valmax = np.nanmin([1.3,np.nanpercentile(p2vmreduced.oi_vis_ft_visamp_s,95)]) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=3, ncols=2, sharex=True, sharey=True, figsize=(5,8)) for baseline in range(0,3): im = axarr[baseline,0].imshow(p2vmreduced.oi_vis_ft_visamp_s[baseline::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,0].set_title('FT Vis. Amp. Baseline '+p2vmreduced.basenames[baseline]+'-S', fontsize=7) im = axarr[baseline,1].imshow(p2vmreduced.oi_vis_ft_visamp_p[baseline::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,1].set_title('FT Vis. Amp. Baseline '+p2vmreduced.basenames[baseline]+'-P', fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 21 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) plotTitle(Story,"FT visibility amplitude waterfall (VISAMP)") plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.") imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize='x-small') plt.rc('ytick', labelsize='x-small') fig, axarr = plt.subplots(nrows=3, ncols=2, sharex=True, sharey=True, figsize=(5,8)) for baseline in range(0,3): im = axarr[baseline,0].imshow(p2vmreduced.oi_vis_ft_visamp_s[baseline+3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,0].set_title('FT Vis. Amp. Baseline '+p2vmreduced.basenames[baseline+3]+'-S', fontsize=7) im = axarr[baseline,1].imshow(p2vmreduced.oi_vis_ft_visamp_p[baseline+3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,1].set_title('FT Vis. Amp. Baseline '+p2vmreduced.basenames[baseline+3]+'-P', fontsize=7) fig.subplots_adjust(hspace=.2) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 21 * cm plt.savefig(imgdata, format='PDF', dpi=150, 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 (VISAMP)") plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.") if hasattr(p2vmreduced,'oi_vis_sc_visamp') & (p2vmreduced.polarsplit == False): valmin = np.max([0,np.nanpercentile(p2vmreduced.oi_vis_sc_visamp,1)]) valmax = np.min([1.3,np.nanpercentile(p2vmreduced.oi_vis_sc_visamp,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(p2vmreduced.oi_vis_sc_visamp[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[0]+' - Comb', fontsize=8) axarr[1].imshow(p2vmreduced.oi_vis_sc_visamp[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[1]+' - Comb', fontsize=8) axarr[2].imshow(p2vmreduced.oi_vis_sc_visamp[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[2]+' - Comb', fontsize=8) axarr[3].imshow(p2vmreduced.oi_vis_sc_visamp[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[3]+' - Comb', fontsize=8) axarr[4].imshow(p2vmreduced.oi_vis_sc_visamp[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[4].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[4]+' - Comb', fontsize=8) axarr[5].imshow(p2vmreduced.oi_vis_sc_visamp[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[5].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[5]+' - Comb', fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if hasattr(p2vmreduced,'oi_vis_sc_visamp_s') & (p2vmreduced.polarsplit == True): valmin = np.max([0,np.nanpercentile(p2vmreduced.oi_vis_sc_visamp_s,1)]) valmax = np.min([1.3,np.nanpercentile(p2vmreduced.oi_vis_sc_visamp_s,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(p2vmreduced.oi_vis_sc_visamp_s[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[0]+' - S', fontsize=8) axarr[1].imshow(p2vmreduced.oi_vis_sc_visamp_s[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[1]+' - S', fontsize=8) axarr[2].imshow(p2vmreduced.oi_vis_sc_visamp_s[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[2]+' - S', fontsize=8) axarr[3].imshow(p2vmreduced.oi_vis_sc_visamp_s[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[3]+' - S', fontsize=8) axarr[4].imshow(p2vmreduced.oi_vis_sc_visamp_s[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[4].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[4]+' - S', fontsize=8) axarr[5].imshow(p2vmreduced.oi_vis_sc_visamp_s[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[5].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[5]+' - S', fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) plotTitle(Story,"SC visibility amplitude waterfall (VISAMP)") plotSubtitle(Story,"Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.") 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(p2vmreduced.oi_vis_sc_visamp_p[0::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[0]+' - P', fontsize=8) axarr[1].imshow(p2vmreduced.oi_vis_sc_visamp_p[1::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[1]+' - P', fontsize=8) axarr[2].imshow(p2vmreduced.oi_vis_sc_visamp_p[2::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[2]+' - P', fontsize=8) axarr[3].imshow(p2vmreduced.oi_vis_sc_visamp_p[3::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[3]+' - P', fontsize=8) axarr[4].imshow(p2vmreduced.oi_vis_sc_visamp_p[4::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[4].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[4]+' - P', fontsize=8) axarr[5].imshow(p2vmreduced.oi_vis_sc_visamp_p[5::6,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[5].set_title('SC VISAMP Baseline '+p2vmreduced.basenames[5]+' - P', fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Wavelength average FT VISAMP vs. group delay #============================================================================== plotTitle(Story,"FT visibility amp. vs. group delay (VISAMP, wavelength averaged)") plotSubtitle(Story,"Visibility vs. group delay (microns). Red = S, Orange = P. (VISAMP/GDELAY)") if p2vmreduced.gdelay_in == True: d=() data=[] resamp = nframeplot // 500 # take one value every resamp value for the plot hsize=16*cm vsize=3*cm if (p2vmreduced.polarsplit == False): visamp_wavg = np.zeros((nbase,p2vmreduced.oi_vis_ft_visamp[baseline::6*resamp,:].shape[0]),dtype='d') for baseline in range(0,nbase): visamp_wavg[baseline,:] = np.nanmean(p2vmreduced.oi_vis_ft_visamp[baseline::6*resamp,:], axis=1) mingd = np.nanmin(p2vmreduced.gdelay_ft) maxgd = np.nanmax(p2vmreduced.gdelay_ft) for baseline in range(0,nbase): data=[] yminval = np.nanpercentile(visamp_wavg[baseline,:],2)-0.1 ymaxval = np.nanpercentile(visamp_wavg[baseline,:],98)+0.1 ticky = (ymaxval - yminval)/5. tickx = (maxgd-mingd)/5. data.append(list(zip(p2vmreduced.gdelay_ft[baseline,::resamp], clipdata(visamp_wavg[baseline,:],yminval,ymaxval)))) a = graphscatteraxes(data, mingd, maxgd,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if (p2vmreduced.polarsplit == True): visamp_wavg_s = np.zeros((nbase,p2vmreduced.oi_vis_ft_visamp_s[baseline::6*resamp,:].shape[0]),dtype='d') visamp_wavg_p = np.zeros((nbase,p2vmreduced.oi_vis_ft_visamp_p[baseline::6*resamp,:].shape[0]),dtype='d') for baseline in range(0,nbase): visamp_wavg_s[baseline,:] = np.nanmean(p2vmreduced.oi_vis_ft_visamp_s[baseline::6*resamp,:], axis=1) visamp_wavg_p[baseline,:] = np.nanmean(p2vmreduced.oi_vis_ft_visamp_p[baseline::6*resamp,:], axis=1) mingd = np.nanmin([p2vmreduced.gdelay_ft_s,p2vmreduced.gdelay_ft_p]) maxgd = np.nanmax([p2vmreduced.gdelay_ft_s,p2vmreduced.gdelay_ft_p]) for baseline in range(0,nbase): data=[] yminval = np.nanpercentile([visamp_wavg_s[baseline,:],visamp_wavg_p[baseline,:]],2)-0.1 ymaxval = np.nanpercentile([visamp_wavg_s[baseline,:],visamp_wavg_p[baseline,:]],98)+0.1 ticky = (ymaxval - yminval)/5. tickx = (maxgd-mingd)/5. data.append(list(zip(p2vmreduced.gdelay_ft_s[baseline,::resamp], clipdata(visamp_wavg_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(p2vmreduced.gdelay_ft_p[baseline,::resamp], clipdata(visamp_wavg_p[baseline,:],yminval,ymaxval)))) a = graphscatteraxes(data,mingd,maxgd,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Wavelength average SC VISAMP vs. group delay #============================================================================== plotTitle(Story,"SC visibility amp. vs. group delay (VISAMP wavelength averaged)") plotSubtitle(Story,"Visibility vs. group delay (microns). Red = S, Orange = P. (VISAMP/GDELAY)") if p2vmreduced.gdelay_in == True: d=() data=[] resamp = 1 # take one value every resamp value for the plot hsize=16*cm vsize=3*cm if (p2vmreduced.polarsplit == False): visamp_wavg = np.zeros((nbase,p2vmreduced.oi_vis_sc_visamp[baseline::6*resamp,:].shape[0]),dtype='d') for baseline in range(0,nbase): visamp_wavg[baseline,:] = np.nanmean(p2vmreduced.oi_vis_sc_visamp[baseline::6*resamp,:], axis=1) print("Baseline %s = %s mean VISAMP = %.3f"%(gravitybasenames[baseline],p2vmreduced.basenames[baseline],np.nanmean(visamp_wavg[baseline,:]))) mingd = np.nanmin(p2vmreduced.gdelay_sc) maxgd = np.nanmax(p2vmreduced.gdelay_sc) for baseline in range(0,nbase): data=[] yminval = np.nanpercentile(visamp_wavg[baseline,:],2)-0.1 ymaxval = np.nanpercentile(visamp_wavg[baseline,:],98)+0.1 ticky = (ymaxval - yminval)/5. tickx = (maxgd-mingd)/5. data.append(list(zip(p2vmreduced.gdelay_sc[baseline,::resamp], clipdata(visamp_wavg[baseline,:],yminval,ymaxval)))) a = graphscatteraxes(data, mingd, maxgd,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if (p2vmreduced.polarsplit == True): visamp_wavg_s = np.zeros((nbase,p2vmreduced.oi_vis_sc_visamp_s[baseline::6*resamp,:].shape[0]),dtype='d') visamp_wavg_p = np.zeros((nbase,p2vmreduced.oi_vis_sc_visamp_p[baseline::6*resamp,:].shape[0]),dtype='d') for baseline in range(0,nbase): visamp_wavg_s[baseline,:] = np.nanmean(p2vmreduced.oi_vis_sc_visamp_s[baseline::6*resamp,:], axis=1) visamp_wavg_p[baseline,:] = np.nanmean(p2vmreduced.oi_vis_sc_visamp_p[baseline::6*resamp,:], axis=1) print("Baseline %s = %s S polar mean VISAMP = %.3f P polar mean VISAMP = %.3f"%\ (gravitybasenames[baseline],p2vmreduced.basenames[baseline],\ np.nanmean(visamp_wavg_s[baseline,:]),np.nanmean(visamp_wavg_p[baseline,:]))) mingd = np.nanmin([p2vmreduced.gdelay_sc_s,p2vmreduced.gdelay_sc_p]) maxgd = np.nanmax([p2vmreduced.gdelay_sc_s,p2vmreduced.gdelay_sc_p]) for baseline in range(0,nbase): data=[] yminval = np.nanpercentile([visamp_wavg_s[baseline,:],visamp_wavg_p[baseline,:]],2)-0.1 ymaxval = np.nanpercentile([visamp_wavg_s[baseline,:],visamp_wavg_p[baseline,:]],98)+0.1 ticky = (ymaxval - yminval)/5. tickx = (maxgd-mingd)/5. data.append(list(zip(p2vmreduced.gdelay_sc_s[baseline,::resamp], clipdata(visamp_wavg_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(p2vmreduced.gdelay_sc_p[baseline,::resamp], clipdata(visamp_wavg_p[baseline,:],yminval,ymaxval)))) a = graphscatteraxes(data,mingd,maxgd,None,yminval,ymaxval,None,hsize,vsize,tickx,ticky,2,p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Temporal quality of the PHASE_REF arg(VISDATA_SC * PHASE_REF) #============================================================================== if hasattr(p2vmreduced,'oi_vis_sc_exp_phase') or hasattr(p2vmreduced,'oi_vis_sc_exp_phase_s'): minw = int(0.4*p2vmreduced.nwave_sc) maxw = int(0.6*p2vmreduced.nwave_sc) plotTitle(Story,"arg(VISDATA_SC * PHASE_REF) vs. time") resamp = 1 # take one value every resamp value for the plot if hasattr(p2vmred,'oi_vis_sc_phi_met') & hasattr(p2vmreduced,'oi_vis_sc_visphi') & (p2vmreduced.polarsplit == False): plotSubtitle(Story,"If phase is not flat => DITSHIFT. In rad vs. time (seconds). Mean over wavelengths.") hsize=5*cm vsize=10*cm opd_sc = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') timescale_sc = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') for baseline in range(0,nbase): phi_met = p2vmreduced.oi_vis_sc_phi_met[baseline::nbase] phi_met = np.outer(phi_met - np.mean(phi_met),1.908287/p2vmreduced.wave_sc) # In pipeline: exp_phase = -phi_ft + phi_met or exp_phase = -phi_ft if --use-met=FALSE f1 = (p2vmreduced.oi_vis_sc_visdata[baseline::nbase,:] * np.exp(1.J*(p2vmreduced.oi_vis_sc_exp_phase[baseline::nbase,:]))) opd_sc[baseline,:] = np.angle(np.nanmean(f1[:,minw:maxw],axis=1))# mean over wavelengths timescale_sc = p2vmreduced.time_sc yminval = np.nanmin(opd_sc)-2 ymaxval = np.nanmax(opd_sc)+2 ticky = np.pi/2. d=() data=[] for baseline in range(0,nbase): data.append(list(zip(timescale_sc[0:nframeplot:resamp], clipdata(opd_sc[baseline,0:nframeplot:resamp],yminval,ymaxval)))) if(baseline == 3): a = graphoutaxes( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,p2vmreduced.basenames[baseline]), else: a = graphoutnoaxis( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+p2vmreduced.basenames[baseline]), data = [] d = d + a plotmatrix = [list(d[0:int(nbase/2)]),list(d[int(nbase/2):nbase])] t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize) Story.append(t) Story.append(Spacer(1,2*mm)) if hasattr(p2vmreduced,'oi_vis_sc_visphi_s') & (p2vmreduced.polarsplit == True): plotSubtitle(Story,"If phase is not flat => DITSHIFT. In rad vs. time (s), mean over wavelengths. Red = S, orange = P.") hsize=5*cm vsize=10*cm opd_sc_s = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') opd_sc_p = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') timescale_sc = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') for baseline in range(0,nbase): f1 = (p2vmreduced.oi_vis_sc_visdata_s[baseline::nbase,:] * np.exp(1.J*(p2vmreduced.oi_vis_sc_exp_phase_s[baseline::nbase,:]))) opd_sc_s[baseline,:] = np.angle(np.nanmean(f1[:,minw:maxw],axis=1)) f1 = (p2vmreduced.oi_vis_sc_visdata_p[baseline::nbase,:] * np.exp(1.J*(p2vmreduced.oi_vis_sc_exp_phase_p[baseline::nbase,:]))) opd_sc_p[baseline,:] = np.angle(np.nanmean(f1[:,minw:maxw],axis=1)) timescale_sc = p2vmreduced.time_sc yminval = np.nanmin(opd_sc_s)-2 ymaxval = np.nanmax(opd_sc_s)+2 ticky = np.pi/2. d=() data=[] for baseline in range(0,nbase): data.append(list(zip(timescale_sc[0:nframeplot:resamp], clipdata(opd_sc_s[baseline,0:nframeplot:resamp],yminval,ymaxval)))) data.append(list(zip(timescale_sc[0:nframeplot:resamp], clipdata(opd_sc_p[baseline,0:nframeplot:resamp],yminval,ymaxval)))) #print data if(baseline == 3): a = graphoutaxes( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),None,yminval,ymaxval,None,hsize,vsize,None,ticky,p2vmreduced.basenames[baseline]), else: a = graphoutnoaxis( data, min(timescale_sc[0:nframeplot]), max(timescale_sc[0:nframeplot]),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+p2vmreduced.basenames[baseline]), data = [] d = d + a plotmatrix = [list(d[0:int(nbase/2)]),list(d[int(nbase/2):nbase])] t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize) Story.append(t) Story.append(Spacer(1,2*mm)) Story.append(PageBreak()) #============================================================================== # Fitting quality of the PHASE_REF #============================================================================== plotTitle(Story,"Referencing SC VISDATA with VISDATA_FT (first sample)") plotSubtitle(Story,"phases (rad) versus wavelengths (microns). Red = -PHASE_REF, Orange = SC, Blue=FT") if (p2vmreduced.polarsplit == False): ref_phase = np.angle(np.exp(-1.j*p2vmreduced.oi_vis_sc_exp_phase[0:nbase,:])) vissc_phase = np.angle(p2vmreduced.oi_vis_sc_visdata[0:nbase,:]) visft_phase = np.angle(p2vmreduced.oi_vis_sc_visdataft[0:nbase,:]) else: ref_phase = np.angle(np.exp(-1.j*p2vmreduced.oi_vis_sc_exp_phase_s[0:nbase,:])) vissc_phase = np.angle(p2vmreduced.oi_vis_sc_visdata_s[0:nbase,:]) visft_phase = np.angle(p2vmreduced.oi_vis_sc_visdataft_s[0:nbase,:]) for baseline in range(0,nbase): data = [] data.append(list(zip(p2vmreduced.wave_sc,ref_phase[baseline,:]))) data.append(list(zip(p2vmreduced.wave_sc,np.angle(clean_gdelay_full(np.exp(1.J*vissc_phase[baseline,:]), p2vmreduced.wave_sc))))) # data.append(zip(p2vmreduced.wave_sc,vissc_phase[baseline,:])) data.append(list(zip(p2vmreduced.wave_ft,visft_phase[baseline,:]))) a = graphoutaxes( data, np.nanmin(p2vmreduced.wave_sc), np.nanmax(p2vmreduced.wave_sc),None,-3.15,+3.15,None,16*cm,3*cm,None,1.0,'Baseline '+p2vmreduced.basenames[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 p2vmreduced.polarsplit == True: plotTitle(Story,"FT phasediff, difference between S and P polar. vs. wavelength") plotSubtitle(Story,"In phase (deg) vs. wavelength (microns).") d=() hsize=16*cm vsize=3*cm yminval=-180 ymaxval=+180 ticky = 180 for baseline in range(0,nbase): data=[] tmp = np.angle(np.mean(p2vmreduced.oi_vis_ft_visdata_s[baseline::6,:] * np.conj(p2vmreduced.oi_vis_ft_visdata_p[baseline::6,:]), axis=0))*180./np.pi; data.append(list(zip(p2vmreduced.wave_ft, tmp))) a = graphoutaxes( data, min(p2vmreduced.wave_ft), max(p2vmreduced.wave_ft),None,yminval,ymaxval,None, hsize,vsize,None,ticky,'FT diff. phase S-P Baseline '+p2vmreduced.basenames[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 p2vmreduced.polarsplit == True: plotTitle(Story,"SC phasediff, difference between S and P polar. vs. wavelength") plotSubtitle(Story,"In phase (deg) vs. wavelength (microns).") d=() hsize=16*cm vsize=3*cm yminval=-180 ymaxval=+180 ticky = 180 for baseline in range(0,nbase): data=[] tmp = np.angle(np.mean(p2vmreduced.oi_vis_sc_visdata_s[baseline::6,:] * np.conj(p2vmreduced.oi_vis_sc_visdata_p[baseline::6,:]), axis=0))*180./np.pi; data.append(list(zip(p2vmreduced.wave_sc, tmp))) a = graphoutaxes( data, min(p2vmreduced.wave_sc), max(p2vmreduced.wave_sc),None,yminval,ymaxval,None, hsize,vsize,None,ticky,'SC diff. phase S-P Baseline '+p2vmreduced.basenames[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # ACQ_CAM OI_VIS_ACQ #============================================================================== if p2vmreduced.acq_in is True: plotTitle(Story,"Pupil from ACQ_CAM versus time") plotSubtitle(Story,"XYZ in [pix] from reference, and R in [deg] from vertical") time = p2vmreduced.oi_acq_time[0::ntel] xmin = np.nanmin (time) - 1 xmax = np.nanmax (time) + 1 data = [] for tel in range(0,ntel): ids = (p2vmreduced.oi_acq_pup_n[tel::ntel] != 0) if (np.sum(ids) == 0): continue data.append(list(zip(time[ids], p2vmreduced.oi_acq_pup_n[tel::ntel][ids] + tel*0.05))) val = p2vmreduced.oi_acq_pup_n[p2vmreduced.oi_acq_pup_n!=0] if (len(data)>0): print (xmin) print (xmax) print((np.nanmin (val)-1)) a = graphoutaxes(data, xmin, xmax, None,np.nanmin (val)-1, np.nanmax (val)+2, None, 16*cm,3*cm,None,1.0,'Pupil_Nspot') Story.append(a) Story.append(Spacer(1,14*mm)) data = [] for tel in range(0,ntel): ids = (p2vmreduced.oi_acq_pup_x[tel::ntel] == p2vmreduced.oi_acq_pup_x[tel::ntel]) if (np.sum(ids) == 0): continue data.append(list(zip(time[ids],p2vmreduced.oi_acq_pup_x[tel::ntel][ids]))) if (len(data)>0): a = graphoutaxes(data, xmin, xmax, None,np.nanmin (p2vmreduced.oi_acq_pup_x-0.01), np.nanmax (p2vmreduced.oi_acq_pup_x+0.01), None, 16*cm,3*cm,None,None,'Pupil_X [pix]') Story.append(a) Story.append(Spacer(1,7*mm)) data = [] for tel in range(0,ntel): ids = (p2vmreduced.oi_acq_pup_y[tel::ntel] == p2vmreduced.oi_acq_pup_y[tel::ntel]) if (np.sum(ids) == 0): continue data.append(list(zip(time[ids],p2vmreduced.oi_acq_pup_y[tel::ntel][ids]))) if (len(data)>0): a = graphoutaxes(data, xmin, xmax, None,np.nanmin (p2vmreduced.oi_acq_pup_y-0.01), np.nanmax (p2vmreduced.oi_acq_pup_y+0.01), None, 16*cm,3*cm,None,None,'Pupil_Y [pix]') Story.append(a) Story.append(Spacer(1,7*mm)) data = [] for tel in range(0,ntel): ids = p2vmreduced.oi_acq_pup_z[tel::ntel] == p2vmreduced.oi_acq_pup_z[tel::ntel] if (np.sum(ids) == 0): continue data.append(list(zip(time[ids],p2vmreduced.oi_acq_pup_z[tel::ntel][ids]))) if (len(data)>0): a = graphoutaxes(data, xmin, xmax, None,np.nanmin (p2vmreduced.oi_acq_pup_z-0.1), np.nanmax (p2vmreduced.oi_acq_pup_z+0.1), None, 16*cm,3*cm,None,None,'Pupil_Z [pix]') Story.append(a) Story.append(Spacer(1,14*mm)) data = [] for tel in range(0,ntel): ids = p2vmreduced.oi_acq_pup_r[tel::ntel] == p2vmreduced.oi_acq_pup_r[tel::ntel] if (np.sum(ids) == 0): continue data.append(list(zip(time[ids],p2vmreduced.oi_acq_pup_r[tel::ntel][ids]))) if (len(data)>0): a = graphoutaxes(data, xmin, xmax, None,np.nanmin (p2vmreduced.oi_acq_pup_r-1), np.nanmax (p2vmreduced.oi_acq_pup_r+1), None, 16*cm,3*cm,None,None,'Pupil_R [deg]') Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Metrology TAC algorithm VAMP #============================================================================== if hasattr(p2vmreduced,'met_vamp_tel_sc'): for tel in range(0,4): plotTitle(Story,"OI_VIS_MET - VAMP") plotSubtitle(Story,"TAC algo VAMP (Volts) for Tel %s versus time (s)"%(p2vmreduced.stations[tel],)) hsize=16*cm vsize=1.8*cm resamp = 100 xval = p2vmreduced.oi_vis_met_time[::resamp] xminval = min(p2vmreduced.oi_vis_met_time) xmaxval = max(p2vmreduced.oi_vis_met_time) yminval = 0. ymaxval = 0.15 ticky = (ymaxval - yminval)/2. d = [] for diode in range(0,4): tmp = p2vmreduced.met_vamp_tel_sc[::resamp,tel*4+diode] d.append( graphoutnoxaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))], xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s SC diode %i'%(p2vmreduced.stations[tel],diode+1)) ) tmp = p2vmreduced.met_vamp_tel_ft[::resamp,tel*4+diode] d.append( graphoutnoxaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))], xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s FT diode %i'%(p2vmreduced.stations[tel],diode+1)) ) yminval = 0. ymaxval = 1. ticky = (ymaxval - yminval)/2. # last diode is FC yminval = 0. ymaxval = 3. ticky = (ymaxval - yminval)/2. tmp = p2vmreduced.met_vamp_fc_sc[::resamp,tel] d.append( graphoutnoxaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))],xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s SC -- FC'%(p2vmreduced.stations[tel],)) ) tmp = p2vmreduced.met_vamp_fc_ft[::resamp,tel] d.append( graphoutaxes( [list(zip(xval, clipdata(tmp,yminval,ymaxval)))],xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Tel %s FT -- FC'%(p2vmreduced.stations[tel],)) ) plotmatrix = [ [d[0]],[d[1]],[d[2]],[d[3]],[d[4]],[d[5]],[d[6]],[d[7]],[d[8]],[d[9]] ] t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize+4*mm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # Metrology TAC and DRS algorithm OPD #============================================================================== if hasattr(p2vmreduced,'met_opd_tel'): for tel in range(0,4): plotTitle(Story,"OI_VIS_MET OPD") plotSubtitle(Story,"red: -1 * TAC algo OPD (microns) for Tel %s versus time (s)"%(p2vmreduced.stations[tel],),spacer=None) plotSubtitle(Story,"orange: DRS algo PHASE (microns) for Tel %s versus time (s)"%(p2vmreduced.stations[tel],)) plotSubtitle(Story,"WARNING the TAC opd is reversed (-1*TAC)") hsize=16*cm vsize=3.6*cm resamp = 100 xval = p2vmreduced.oi_vis_met_time[::resamp] xminval = min(p2vmreduced.oi_vis_met_time) xmaxval = max(p2vmreduced.oi_vis_met_time) d = [] for diode in range(0,4): tmp = -1*p2vmreduced.met_opd_tel[::resamp,tel*4+diode] * 1e6 tmp2 = p2vmreduced.oi_vis_met_phase_tel[::resamp,tel*4+diode] * Lambda_met / (2*np.pi) yminval = np.minimum( np.min(tmp), np.min(tmp2)) ymaxval = np.maximum( np.max(tmp), np.max(tmp2)) mean = 0.5*(yminval+ymaxval) if (ymaxval-yminval)<2.0: yminval = mean - 1.0 ymaxval = mean + 1.0 ticky = (ymaxval - yminval)/4. d.append( graphoutnoxaxes( [list(zip(xval,tmp)), list(zip(xval,tmp2))], xminval,xmaxval,yminval,ymaxval,hsize,vsize,ticky,'Tel %s diode %i'%(p2vmreduced.stations[tel],diode+1)) ) tmp = -1*p2vmreduced.met_opd_fc[::resamp,tel] * 1e6 tmp2 = p2vmreduced.oi_vis_met_phase[::resamp,tel] * Lambda_met / (2*np.pi) yminval = np.minimum( np.min(tmp), np.min(tmp2)) ymaxval = np.maximum( np.max(tmp), np.max(tmp2)) ticky = (ymaxval - yminval)/4. d.append( graphoutaxes( [list(zip(xval,tmp)),list(zip(xval,tmp2))],xminval,xmaxval,None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Tel %s -- FC'%(p2vmreduced.stations[tel],)) ) plotmatrix = [ [d[0]],[d[1]],[d[2]],[d[3]],[d[4]] ] t = Table(plotmatrix,colWidths=hsize,rowHeights=vsize+4*mm) t.setStyle(TableStyle([('ALIGN', (1,1), (-1,-1), 'LEFT')])) t.hAlign = 'LEFT' Story.append(t) Story.append(PageBreak()) #============================================================================== # Metrology angles #============================================================================== if hasattr(p2vmreduced,'met_opd_tel'): plotTitle(Story,"MET_TEL differential Tip Tilt and Astig of telescope metrology",spacer=Spacer(1,1*mm)) plotSubtitle(Story,"Average measured TAC algo OPD (microns) versus time (s) "+\ p2vmreduced.stations[0] +" = red, "+\ p2vmreduced.stations[1] +" = orange, "+\ p2vmreduced.stations[2] +" = blue, "+\ p2vmreduced.stations[3] +" = green.") hsize = 16*cm vsize = 4.5*cm resamp = 50 # met_tel = p2vmreduced.oi_vis_met_phase_tel[::resamp,:] * Lambda_met / (2*np.pi) met_tel = p2vmreduced.met_opd_tel[::resamp,:] hsize = 16*cm vsize = 4.5*cm plotSubtitle(Story,"Tip in microns vs. time (seconds).",spacer=Spacer(1,2*mm)) data, yminval, ymaxval = ([], 1e20, -1e20) for tel in range(0,4): phase_met_plot = met_tel[:,tel*4+0] + met_tel[:,tel*4+1] - met_tel[:,tel*4+2] - met_tel[:,tel*4+3] phase_met_plot -= np.nanmedian(phase_met_plot) yminval = np.minimum( yminval, np.nanmin(phase_met_plot)-1e-12 ) ymaxval = np.maximum( ymaxval, np.nanmax(phase_met_plot)+1e-12 ) ticky = (ymaxval-yminval)/5. # in microns data.append( tuple(zip(p2vmreduced.oi_vis_met_time[::resamp], phase_met_plot) )) a = graphoutaxes( data, min(p2vmreduced.oi_vis_met_time),max(p2vmreduced.oi_vis_met_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Tip') Story.append(a) Story.append(Spacer(1,12*mm)) plotSubtitle(Story,"Tilt in microns vs. time (seconds).",spacer=Spacer(1,2*mm)) data, yminval, ymaxval = ([], 1e20, -1e20) for tel in range(0,4): phase_met_plot = met_tel[:,tel*4+0] - met_tel[:,tel*4+1] + met_tel[:,tel*4+2] - met_tel[:,tel*4+3] phase_met_plot -= np.nanmedian(phase_met_plot) yminval = np.minimum( yminval, np.nanmin(phase_met_plot)-1e-12 ) ymaxval = np.maximum( ymaxval, np.nanmax(phase_met_plot)+1e-12 ) ticky = (ymaxval-yminval)/5. # in microns data.append( tuple(zip(p2vmreduced.oi_vis_met_time[::resamp], phase_met_plot) )) a = graphoutaxes( data, min(p2vmreduced.oi_vis_met_time),max(p2vmreduced.oi_vis_met_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Tilt') Story.append(a) Story.append(Spacer(1,12*mm)) plotSubtitle(Story,"Astigmatism in microns vs. time (seconds).",spacer=Spacer(1,2*mm)) data, yminval, ymaxval = ([], 1e20, -1e20) for tel in range(0,4): phase_met_plot = met_tel[:,tel*4+0] - met_tel[:,tel*4+1] - met_tel[:,tel*4+2] + met_tel[:,tel*4+3] phase_met_plot -= np.nanmedian(phase_met_plot) yminval = np.minimum( yminval,np.nanmin(phase_met_plot)-1e-12 ) ymaxval = np.maximum( yminval,np.nanmax(phase_met_plot)+1e-12 ) ticky = (ymaxval-yminval)/5. # in microns data.append( tuple(zip(p2vmreduced.oi_vis_met_time[::resamp], phase_met_plot) )) a = graphoutaxes( data, min(p2vmreduced.oi_vis_met_time),max(p2vmreduced.oi_vis_met_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Astigmatism') Story.append(a) #Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Metrology signals by baseline #============================================================================== plotTitle(Story,"MET_FC Fiber coupler met. signals (OI_VIS_MET)") plotSubtitle(Story,"Microns vs. time (seconds).") #\nBaselines: "+p2vmreduced.basenames[0]+" = red, "+p2vmreduced.basenames[1]+" = orange, "+p2vmreduced.basenames[2]+" = blue, "+p2vmreduced.basenames[3]+" = green, "+p2vmreduced.basenames[4]+" =cyan, "+p2vmreduced.basenames[5]+" =purple.") met = baseline_phases(p2vmreduced.oi_vis_met_phase) nframe_met = met.shape[0] resamp = nframe_met // 500 # take one value every resamp value for the plot ticky = (ymaxval-yminval)/10. # in microns hsize = 16*cm vsize = 3*cm for baseline in range(0,nbase): phase_met_plot = met[::resamp,baseline]*Lambda_met/(2*np.pi) mean_met = np.nanmean(phase_met_plot) phase_met_plot -= mean_met # subtraction of mean value data = [] yminval = np.nanmin(phase_met_plot)-1 ymaxval = np.nanmax(phase_met_plot)+1 ticky = (ymaxval-yminval)/5. # in microns data.append( tuple(zip(p2vmreduced.oi_vis_met_time[::resamp], clipdata(phase_met_plot,yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.oi_vis_met_time),max(p2vmreduced.oi_vis_met_time),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+' Mean = %.3e mic.' % mean_met) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Power spectral density of the METROLOGY unwrapped phase vs. time #============================================================================== plotTitle(Story,"Periodogram of the metrology OPD values (METROLOGY)") plotSubtitle(Story,"Microns RMS vs. Hz.") if welchpresent == True: resamp=1 phasemet = np.zeros((p2vmreduced.oi_vis_met_phase[::resamp].shape[0],4),dtype='d') for tel in range(0,4): phasemet[:,tel] = p2vmreduced.oi_vis_met_phase[::resamp,tel]*Lambda_met/(2*np.pi) PSD_met=[1,2,3,4] metfreq = resamp/(p2vmreduced.oi_vis_met_time[1]-p2vmreduced.oi_vis_met_time[0]) for tel in range(0,4): f_met, PSD_met[tel] = welch(np.unwrap(phasemet[:,tel]), fs=metfreq, detrend='linear', nperseg=1024, scaling='spectrum') yminval = 0 hsize = 16*cm vsize = 4.5*cm for tel in range(0,4): data = [] ymaxval = np.nanmean(np.sqrt(PSD_met[tel]))*5 ticky = (ymaxval-yminval)/5 data.append(tuple(zip(f_met, clipdata(np.sqrt(PSD_met[tel]), yminval,ymaxval)))) a = graphoutaxes( data, min(f_met),max(f_met),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Telescope '+p2vmreduced.stations[tel]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Differential GDELAY and PHASE signals #============================================================================== hsize = 16*cm vsize = 2.7*cm yminval = -2.2 ymaxval = +2.2 ticky = 1. if hasattr(p2vmreduced,'opd_diff') & (p2vmreduced.polarsplit == False): plotTitle(Story,"Differential OPD from group delay and phase") plotSubtitle(Story,"red: GDELAY_SC-GDELAY_FT, mean GD subtracted (microns) vs. time (seconds)",spacer=None) plotSubtitle(Story,"orange: PHASE_SC-PHASE_FT, mean GD subtracted before averaging all wavelengths",spacer=None) plotSubtitle(Story,"blue: PHASE_SC-PHASE_FT at the center of the band (1 single SC channel), mean subtracted") for baseline in range(0,nbase): data = [] data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.gd_diff[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.opd_diff_avg[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.opd_diff[baseline,:,int(p2vmreduced.nwave_sc/2)],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+' - Comb') Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) if hasattr(p2vmreduced,'opd_diff_s') & (p2vmreduced.polarsplit == True): plotTitle(Story,"Differential OPD from group delay and phase (S polar)") plotSubtitle(Story,"red: GDELAY_SC-GDELAY_FT, mean GD subtracted (microns) vs. time (seconds)",spacer=None) plotSubtitle(Story,"orange: PHASE_SC-PHASE_FT, mean GD subtracted before averaging all wavelengths",spacer=None) plotSubtitle(Story,"blue: PHASE_SC-PHASE_FT at the center of the band (1 single SC channel), mean subtracted") for baseline in range(0,nbase): data = [] data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.gd_diff_s[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.opd_diff_avg_s[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.opd_diff_s[baseline,:,int(p2vmreduced.nwave_sc/2)],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+' - S') Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) plotTitle(Story,"Differential OPD from group delay and phase (P polar)") plotSubtitle(Story,"red: GDELAY_SC-GDELAY_FT, mean GD subtracted (microns) vs. time (seconds)",spacer=None) plotSubtitle(Story,"orange: PHASE_SC-PHASE_FT, mean GD subtracted before averaging all wavelengths",spacer=None) plotSubtitle(Story,"blue: PHASE_SC-PHASE_FT at the center of the band (1 single SC channel), mean subtracted") for baseline in range(0,nbase): data = [] data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.gd_diff_p[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.opd_diff_avg_p[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(p2vmreduced.time_sc, clipdata(p2vmreduced.opd_diff_p[baseline,:,int(p2vmreduced.nwave_sc/2)],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.time_sc),max(p2vmreduced.time_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+' - P') Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of p2vmreduced.opd_diff2 (differential OPD) #============================================================================== if hasattr(p2vmreduced,'opd_diff') & (p2vmreduced.polarsplit == False): plotTitle(Story,"Differential OPD SC-FT from PHASE waterfall") plotSubtitle(Story,"In microns (mean subtracted), wavelength in horizontal axis, frame in sequence in vertical axis.") valamp = np.max([np.abs(np.nanpercentile(p2vmreduced.opd_diff,1)),np.abs(np.nanpercentile(p2vmreduced.opd_diff,99))]) valmin = -valamp valmax = +valamp 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)) for baseline in range(0,nbase): im = axarr[baseline].imshow(p2vmreduced.opd_diff[baseline,:,:],vmin=valmin,vmax=valmax,cmap='RdYlBu', interpolation='nearest', aspect='auto') axarr[baseline].set_title('dOPD Baseline '+p2vmreduced.basenames[baseline]+' - 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) widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if hasattr(p2vmreduced,'opd_diff_s') & (p2vmreduced.polarsplit == True): plotTitle(Story,"Differential OPD SC-FT from PHASE waterfall (S polar)") plotSubtitle(Story,"In microns (mean subtracted), wavelength in horizontal axis, frame in sequence in vertical axis.") valamp = np.max([np.abs(np.nanpercentile(p2vmreduced.opd_diff_s,1)),np.abs(np.nanpercentile(p2vmreduced.opd_diff_s,99))]) valmin = -valamp valmax = +valamp 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)) for baseline in range(0,nbase): im = axarr[baseline].imshow(p2vmreduced.opd_diff_s[baseline,:,:],vmin=valmin,vmax=valmax,cmap='RdYlBu', interpolation='nearest', aspect='auto') axarr[baseline].set_title('dOPD Baseline '+p2vmreduced.basenames[baseline]+' - 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) widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) plotTitle(Story,"Differential OPD SC-FT from PHASE waterfall (P polar)") plotSubtitle(Story,"In microns (mean GD_SC-GD_FT subtracted), wavelength in X, frame in sequence in Y.") valamp = np.max([np.abs(np.nanpercentile(p2vmreduced.opd_diff_p,1)),np.abs(np.nanpercentile(p2vmreduced.opd_diff_p,99))]) valmin = -valamp valmax = +valamp 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)) for baseline in range(0,nbase): im = axarr[baseline].imshow(p2vmreduced.opd_diff_p[baseline,:,:],vmin=valmin,vmax=valmax,cmap='RdYlBu', interpolation='nearest', aspect='auto') axarr[baseline].set_title('dOPD Baseline '+p2vmreduced.basenames[baseline]+' - 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) widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=150, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Time sequence of wavelength averaged SC - [FT+MET] differential OPD #============================================================================== if (p2vmreduced.polarsplit == False): hsize=5*cm vsize=10*cm opd_sc = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') opd_ft = np.zeros((nbase,p2vmreduced.time_ft.shape[0]),dtype='d') opd_ft_int = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') opd_met_int = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') timescale_sc = np.zeros((nbase,p2vmreduced.time_sc.shape[0]),dtype='d') unwrapped_met = np.zeros((nbase,p2vmreduced.nframe_met),dtype='d') unwrapped_met_smooth = np.zeros((nbase,p2vmreduced.nframe_met),dtype='d') opd_ft_smooth = np.zeros((nbase,p2vmreduced.nframe_ft),dtype='d') opd_sc_zp = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_abs_zp =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') timescale_sc = p2vmreduced.time_sc timescale_ft = p2vmreduced.time_ft timescale_met = p2vmreduced.oi_vis_met_time dt_sc = timescale_sc[1]-timescale_sc[0] # time between each frame of the sc #dit_sc = p2vmreduced.header['HIERARCH ESO DET2 SEQ1 DIT'] # actual exposure time of SC dt_ft = timescale_ft[1]-timescale_ft[0] # time between each frame of the sc dt_met = timescale_met[1]-timescale_met[0] # time between each frame of the sc opd_sc_abs =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_abs_zp =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_ft_int = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') metrange = np.zeros(nbase) boxsmooth_ft = int(0.5*dt_sc/dt_ft) # Moving average window for boxcar smoothing of the FDDL signals boxsmooth_met = int(0.5*dt_sc/dt_met) # Moving average window for boxcar smoothing of the MET signals for baseline in range(0,nbase): # Group delay instead of unwrapped phase - pkervell 09/02/16 # Science combiner opd_sc[baseline,:] = p2vmreduced.gdelay_sc[baseline,:] # Fringe tracker if hasattr(p2vmreduced,"gdelay_ft_resamp"): opd_ft_int[baseline,:] = p2vmreduced.gdelay_ft_resamp[baseline,:] # from the pipeline #print " FT interpolated OPD from pipeline" else: opd_ft_smooth[baseline,:] = np.convolve(p2vmreduced.gdelay_ft[baseline,:], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') opd_ft_int[baseline,:] = np.interp(timescale_sc+(dt_sc/2.0), timescale_ft+(dt_ft/2.0), opd_ft_smooth[baseline,:]) # Metrology if hasattr(p2vmreduced,"oi_vis_sc_opd_met"): opd_met_int[baseline,:] = p2vmreduced.oi_vis_sc_opd_met[baseline,:] # from the pipeline #print " MET interpolated OPD from pipeline" else: unwrapped_met[baseline,:] = np.unwrap(met[:,baseline])*Lambda_met/(2*np.pi) # computation unwrapped_met_smooth[baseline,:]=np.convolve(unwrapped_met[baseline,:], np.ones((boxsmooth_met,))/boxsmooth_met, mode='same') opd_met_int[baseline,:] = np.interp(timescale_sc+(dt_sc/2.0), timescale_met, unwrapped_met_smooth[baseline,:]) mean_sc = np.nanmean(opd_sc[baseline,:]) mean_ft = np.nanmean(opd_ft_int[baseline,:]) mean_met = np.nanmean(opd_met_int[baseline,:]) # ======================================================== # Computation of the differential quantities # BEWARE OF SIGNS ! # ======================================================== # Without mean value subtractions: opd_sc_abs[baseline,:] = opd_met_int[baseline,:] + opd_ft_int[baseline,:] # With mean value subtractions: opd_sc_zp[baseline,:] = opd_sc[baseline,:] - mean_sc opd_sc_abs_zp[baseline,:] = (opd_met_int[baseline,:]-mean_met) + (opd_ft_int[baseline,:]-mean_ft) metrange[baseline] = np.max(opd_met_int[baseline,:])-np.min(opd_met_int[baseline,:]) metamplitude = np.max(metrange) data=[] for baseline in range(0,nbase): plotTitle(Story, "Astrometric differential OPD vs. time - Baseline "+p2vmreduced.basenames[baseline]) Story.append(Spacer(1,5*mm)) plotSubtitle(Story,"Group delay of SC and FT (mic) vs. time (s). Mean over wavelengths. Red = SC, Orange = FT.") data = [] yminval = np.min([opd_sc[baseline,:],opd_ft_int[baseline,:]]) ymaxval = np.max([opd_sc[baseline,:],opd_ft_int[baseline,:]]) # margin for plot yminval -= 0.1*np.abs(ymaxval-yminval) ymaxval += 0.1*np.abs(ymaxval-yminval) ticky = (ymaxval-yminval)/5 data.append(tuple(zip(timescale_sc, clipdata(opd_sc[baseline,:],yminval,ymaxval)))) data.append(tuple(zip(timescale_sc, clipdata(opd_ft_int[baseline,:],yminval,ymaxval)))) hsize=16*cm vsize=4*cm a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+'; GD_FT and GD_SC; Mean GD_SC-GD_FT: %(mean_gdscft).3f microns'%{"mean_gdscft":p2vmreduced.mean_gd_sc[baseline]-p2vmreduced.mean_gd_ft[baseline]}) Story.append(a) Story.append(Spacer(1,7*mm)) plotSubtitle(Story,"[OPD_SC-OPD_FT] from PHASE vs. time (seconds). Time avg [GD_SC-GD_FT] subtracted.") data = [] resamp = 50 # ymaxval = np.max(np.abs(p2vmreduced.opd_diff_avg[baseline,:])) # yminval = -ymaxval ymaxval = 1. #microns yminval = -ymaxval ticky = (ymaxval-yminval)/6 hsize=16*cm vsize=4*cm data.append( tuple(zip(timescale_sc, clipdata(p2vmreduced.opd_diff_avg[baseline,:],yminval,ymaxval)) )) a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+'; dOPD signal; Mean value: %(meanval).3f microns; RMS: %(rmsval).3f microns '%{"meanval":np.mean(p2vmreduced.opd_diff_avg[baseline,:]),"rmsval":np.std(p2vmreduced.opd_diff_avg[baseline,:])}) Story.append(a) Story.append(Spacer(1,7*mm)) plotSubtitle(Story,"Fiber coupler MET resampled signal vs. time (seconds). Slope subtracted.") data = [] phase_met_plot = np.copy(p2vmreduced.oi_vis_sc_opd_met[baseline,:]) z = np.polyfit(timescale_sc, phase_met_plot, 1) poly = np.poly1d(z) meanmet = np.nanmean(phase_met_plot) #phase_met_plot -= meanmet phase_met_plot -= poly(timescale_sc) ymaxval = np.max(np.abs(phase_met_plot))*1.2 yminval = -ymaxval # ymaxval = np.mean(phase_met_plot) + 1.05*metamplitude ticky = (ymaxval-yminval)/5 hsize=16*cm vsize=4*cm data.append( tuple(zip(timescale_sc, clipdata(phase_met_plot,yminval,ymaxval)) )) a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+'; MET signal; Slope: %(slopemet).3f microns/s; Mean value: %(meanval).3f microns; RMS residual: %(rmsmet).3f microns '%{"slopemet":z[0],"meanval":meanmet,"rmsmet":np.std(phase_met_plot)}) Story.append(a) Story.append(Spacer(1,7*mm)) plotSubtitle(Story,"Total dOPD=[OPD_SC-OPD_FT]+[GD_SC-GD_FT]+MET vs. time (s). Slope subtracted.") data = [] meandopd = np.mean(p2vmreduced.total_dopd_avg[baseline,:]) total_dopd_plot = np.copy(p2vmreduced.total_dopd_avg[baseline,:]) z = np.polyfit(timescale_sc, total_dopd_plot, 1) poly = np.poly1d(z) total_dopd_plot -= poly(timescale_sc) ymaxval = 1. #microns yminval = -ymaxval ticky = (ymaxval-yminval)/6 hsize=16*cm vsize=4*cm data.append( tuple(zip(timescale_sc, clipdata(total_dopd_plot,yminval,ymaxval)) )) a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'Baseline '+p2vmreduced.basenames[baseline]+'; Total dOPD; Slope: %(slopedopd).3f microns/s; Mean value: %(meandopd).3f microns; RMS residual: %(rmsdopd).3f microns '%{"slopedopd":z[0],"meandopd":meandopd,"rmsdopd":np.std(total_dopd_plot)}) Story.append(a) # plotSubtitle(Story,"Fiber coupler MET signal vs. time (seconds). Slope subtracted.") # data = [] # resamp = 50 # phase_met_plot = met[::resamp,baseline]*Lambda_met/(2*np.pi) # z = np.polyfit(p2vmreduced.oi_vis_met_time[::resamp], phase_met_plot, 1) # poly = np.poly1d(z) # phase_met_plot -= poly(p2vmreduced.oi_vis_met_time[::resamp]) # ymaxval = np.max(np.abs(phase_met_plot))*1.05 # yminval = -ymaxval*1.05 # ticky = (ymaxval-yminval)/5 # hsize=16*cm # vsize=4*cm # data.append( tuple(zip(p2vmreduced.oi_vis_met_time[::resamp], clipdata(phase_met_plot,yminval,ymaxval)) )) # a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+p2vmreduced.basenames[baseline]+'; MET-slope signal; RMS: %(rmsval).3f microns '%{"rmsval":np.std(phase_met_plot)}) # Story.append(a) # # Story.append(Spacer(1,7*mm)) # # plotSubtitle(Story,"Astrometric GD_SC-[GD_FT+MET] residual in microns vs. time (s). Linear fit subtracted.") # data = [] # diff_opd = opd_sc[baseline,:]-opd_sc_abs[baseline,:] # z = np.polyfit(timescale_sc, diff_opd, 1) # poly = np.poly1d(z) # opd_slopesub = diff_opd - poly(timescale_sc) # yminval = np.min(opd_slopesub) # ymaxval = np.max(opd_slopesub) # yminval -= 0.1*np.abs(ymaxval-yminval) # ymaxval += 0.1*np.abs(ymaxval-yminval) # ticky = (ymaxval-yminval)/5 # hsize=16*cm # vsize=4*cm # err = np.nanstd(opd_slopesub) # data.append(tuple(zip(timescale_sc, clipdata(opd_slopesub,yminval,ymaxval)))) # a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+p2vmreduced.basenames[baseline]+'; SC-[FT+MET]-linear ; Slope: %(slope).3f microns/s; ZeroPt(t=0s): %(zeropt).3f microns; RMS resid.: %(err).3f microns' % {"slope":z[0],"zeropt":z[1],"err":err} ) # Story.append(a) Story.append(PageBreak()) if hasattr(p2vmreduced,'oi_vis_sc_visphi_s') & (p2vmreduced.polarsplit == True): opd_sc_s = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_p = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') timescale_sc = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') # Metrology and FT OPD interpolated to SC frequency opd_met_int = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_s_zp = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_p_zp = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') unwrapped_met = np.zeros((nbase,p2vmreduced.nframe_met),dtype='d') unwrapped_met_smooth = np.zeros((nbase,p2vmreduced.nframe_met),dtype='d') timescale_sc = p2vmreduced.time_sc timescale_ft = p2vmreduced.time_ft timescale_met = p2vmreduced.oi_vis_met_time dt_sc = timescale_sc[1]-timescale_sc[0] # time between each frame of the sc dt_ft = timescale_ft[1]-timescale_ft[0] # time between each frame of the sc dt_met = timescale_met[1]-timescale_met[0] # time between each frame of the sc opd_sc_abs_s =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_abs_p =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_abs_s_zp =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_sc_abs_p_zp =np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_ft_s_int = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_ft_p_int = np.zeros((nbase,p2vmreduced.nframe_sc),dtype='d') opd_ft_s_smooth = np.zeros((nbase,p2vmreduced.nframe_ft),dtype='d') opd_ft_p_smooth = np.zeros((nbase,p2vmreduced.nframe_ft),dtype='d') boxsmooth_ft = int(0.5*dt_sc/dt_ft) # Moving average window for boxcar smoothing of the FDDL signals boxsmooth_met = int(0.5*dt_sc/dt_met) # Moving average window for boxcar smoothing of the MET signals for baseline in range(0,nbase): # Group delay instead of unwrapped phase - pkervell 09/02/16 # Science combiner opd_sc_s[baseline,:] = p2vmreduced.gdelay_sc_s[baseline,:] opd_sc_p[baseline,:] = p2vmreduced.gdelay_sc_p[baseline,:] # Metrology if hasattr(p2vmreduced,"oi_vis_sc_opd_met"): opd_met_int[baseline,:] = p2vmreduced.oi_vis_sc_opd_met[baseline,:] # from the pipeline else: unwrapped_met[baseline,:] = np.unwrap(met[:,baseline])*Lambda_met/(2*np.pi) unwrapped_met_smooth[baseline,:]=np.convolve(unwrapped_met[baseline,:], np.ones((boxsmooth_met,))/boxsmooth_met, mode='same') opd_met_int[baseline,:] = np.interp(timescale_sc+(dt_sc/2.0), timescale_met, unwrapped_met_smooth[baseline,:]) mean_met = np.nanmean(opd_met_int[baseline,:]) # Fringe tracker if hasattr(p2vmreduced,"gdelay_ft_resamp_s"): opd_ft_s_int[baseline,:] = p2vmreduced.gdelay_ft_resamp_s[baseline,:] opd_ft_p_int[baseline,:] = p2vmreduced.gdelay_ft_resamp_p[baseline,:] else: opd_ft_s_smooth[baseline,:]=np.convolve(p2vmreduced.gdelay_ft_s[baseline,:], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') opd_ft_p_smooth[baseline,:]=np.convolve(p2vmreduced.gdelay_ft_p[baseline,:], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') opd_ft_s_int[baseline,:] = np.interp(timescale_sc+(dt_sc/2.0), timescale_ft+(dt_ft/2.0), opd_ft_s_smooth[baseline,:]) opd_ft_p_int[baseline,:] = np.interp(timescale_sc+(dt_sc/2.0), timescale_ft+(dt_ft/2.0), opd_ft_p_smooth[baseline,:]) # ======================================================== # Computation of the differential quantities # BEWARE OF SIGNS ! # ======================================================== # opd_sc_s_zp[baseline,:] = opd_sc_s[baseline,:] - mean_sc_s # opd_sc_p_zp[baseline,:] = opd_sc_p[baseline,:] - mean_sc_p # opd_sc_abs_s[baseline,:] = opd_met_int[baseline,:] + opd_ft_s_int[baseline,:] # opd_sc_abs_p[baseline,:] = opd_met_int[baseline,:] + opd_ft_p_int[baseline,:] # opd_sc_abs_s_zp[baseline,:] = (opd_met_int[baseline,:]-mean_met) + (opd_ft_s_int[baseline,:]-mean_ft_s) # opd_sc_abs_p_zp[baseline,:] = (opd_met_int[baseline,:]-mean_met) + (opd_ft_p_int[baseline,:]-mean_ft_p) opd_sc_s_zp[baseline,:] = opd_sc_s[baseline,:] opd_sc_p_zp[baseline,:] = opd_sc_p[baseline,:] opd_sc_abs_s[baseline,:] = opd_met_int[baseline,:] + opd_ft_s_int[baseline,:] opd_sc_abs_p[baseline,:] = opd_met_int[baseline,:] + opd_ft_p_int[baseline,:] opd_sc_abs_s_zp[baseline,:] = (opd_met_int[baseline,:]) + (opd_ft_s_int[baseline,:]) opd_sc_abs_p_zp[baseline,:] = (opd_met_int[baseline,:]) + (opd_ft_p_int[baseline,:]) data=[] for baseline in range(0,nbase): plotTitle(Story, "Astrometric differential OPD vs. time - Baseline "+p2vmreduced.basenames[baseline]) Story.append(Spacer(1,5*mm)) plotSubtitle(Story,"Group delay of SC and FT (mic) vs. time (s). Mean over wavelengths. Red = SC S, Orange = SC P, Blue = FT S, Green = FT P.") data = [] yminval = np.min([opd_sc_s[baseline,:],opd_ft_s_int[baseline,:],opd_sc_p[baseline,:],opd_ft_p_int[baseline,:]]) ymaxval = np.max([opd_sc_s[baseline,:],opd_ft_s_int[baseline,:],opd_sc_p[baseline,:],opd_ft_p_int[baseline,:]]) # margin for plot yminval -= 0.1*np.abs(ymaxval-yminval) ymaxval += 0.1*np.abs(ymaxval-yminval) ticky = (ymaxval-yminval)/5 data.append(tuple(zip(timescale_sc, clipdata(opd_sc_s[baseline,:],yminval,ymaxval)))) data.append(tuple(zip(timescale_sc, clipdata(opd_sc_p[baseline,:],yminval,ymaxval)))) data.append(tuple(zip(timescale_sc, clipdata(opd_ft_s_int[baseline,:],yminval,ymaxval)))) data.append(tuple(zip(timescale_sc, clipdata(opd_ft_p_int[baseline,:],yminval,ymaxval)))) hsize=16*cm vsize=4*cm a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ p2vmreduced.basenames[baseline]+'-S; Mean GD_SC-GD_FT: %(mean_gdscfts).3f mic'%\ {"mean_gdscfts":p2vmreduced.mean_gd_sc_s[baseline]-p2vmreduced.mean_gd_ft_s[baseline]}+\ ' '+p2vmreduced.basenames[baseline]+'-P; Mean GD_SC-GD_FT: %(mean_gdscftp).3f mic'%\ {"mean_gdscftp":p2vmreduced.mean_gd_sc_p[baseline]-p2vmreduced.mean_gd_ft_p[baseline]}) Story.append(a) Story.append(Spacer(1,7*mm)) plotSubtitle(Story,"[OPD_SC-OPD_FT] from PHASE vs. time (seconds). Time avg [GD_SC-GD_FT] subtracted.") data = [] resamp = 50 ymaxval = 1. #microns yminval = -ymaxval ticky = (ymaxval-yminval)/6 hsize=16*cm vsize=4*cm data.append( tuple(zip(timescale_sc, clipdata(p2vmreduced.opd_diff_avg_s[baseline,:],yminval,ymaxval)) )) data.append( tuple(zip(timescale_sc, clipdata(p2vmreduced.opd_diff_avg_p[baseline,:],yminval,ymaxval)) )) a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ p2vmreduced.basenames[baseline]+'-S; dOPD; Mean: %(meanval).3f mic; RMS: %(rmsval).3f mic '\ %{"meanval":np.mean(p2vmreduced.opd_diff_avg_s[baseline,:]),"rmsval":np.std(p2vmreduced.opd_diff_avg_s[baseline,:])}+\ p2vmreduced.basenames[baseline]+'-P; dOPD; Mean: %(meanval).3f mic; RMS: %(rmsval).3f mic '\ %{"meanval":np.mean(p2vmreduced.opd_diff_avg_p[baseline,:]),"rmsval":np.std(p2vmreduced.opd_diff_avg_p[baseline,:])} ) Story.append(a) Story.append(Spacer(1,7*mm)) plotSubtitle(Story,"Fiber coupler MET resampled signal vs. time (seconds). Slope subtracted.") data = [] phase_met_plot = np.copy(p2vmreduced.oi_vis_sc_opd_met[baseline,:]) z = np.polyfit(timescale_sc, phase_met_plot, 1) poly = np.poly1d(z) meanmet = np.nanmean(phase_met_plot) #phase_met_plot -= meanmet phase_met_plot -= poly(timescale_sc) ymaxval = np.max(np.abs(phase_met_plot))*1.2 yminval = -ymaxval # ymaxval = np.mean(phase_met_plot) + 1.05*metamplitude ticky = (ymaxval-yminval)/5 hsize=16*cm vsize=4*cm data.append( tuple(zip(timescale_sc, clipdata(phase_met_plot,yminval,ymaxval)) )) a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,p2vmreduced.basenames[baseline]+'; MET; Slope: %(slopemet).3f micr/s; Mean: %(meanval).3f mic; RMS: %(rmsmet).3f microns '%{"slopemet":z[0],"meanval":meanmet,"rmsmet":np.std(phase_met_plot)}) Story.append(a) Story.append(Spacer(1,7*mm)) plotSubtitle(Story,"Total dOPD=[OPD_SC-OPD_FT]+[GD_SC-GD_FT]+MET vs. time (s). Slope subtracted.") data = [] meandopd_s = np.mean(p2vmreduced.total_dopd_avg_s[baseline,:]) meandopd_p = np.mean(p2vmreduced.total_dopd_avg_p[baseline,:]) total_dopd_plot_s = np.copy(p2vmreduced.total_dopd_avg_s[baseline,:]) total_dopd_plot_p = np.copy(p2vmreduced.total_dopd_avg_p[baseline,:]) z_s = np.polyfit(timescale_sc, total_dopd_plot_s, 1) z_p = np.polyfit(timescale_sc, total_dopd_plot_p, 1) poly_s = np.poly1d(z_s) poly_p = np.poly1d(z_p) total_dopd_plot_s -= poly_s(timescale_sc) total_dopd_plot_p -= poly_p(timescale_sc) ymaxval = 1. #microns yminval = -ymaxval ticky = (ymaxval-yminval)/6 hsize=16*cm vsize=4*cm data.append( tuple(zip(timescale_sc, clipdata(total_dopd_plot_s,yminval,ymaxval)) )) data.append( tuple(zip(timescale_sc, clipdata(total_dopd_plot_p,yminval,ymaxval)) )) a = graphoutaxes( data, np.min(timescale_sc), np.max(timescale_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,\ 'B'+p2vmreduced.basenames[baseline]+'-S; Slope: %(slopedopd).3f mic/s; Mean: %(meandopd).3f microns; RMS: %(rmsdopd).3f mic '\ %{"slopedopd":z_s[0],"meandopd":meandopd_s,"rmsdopd":np.std(total_dopd_plot_s)}+\ 'B'+p2vmreduced.basenames[baseline]+'-P; Slope: %(slopedopd).3f mic/s; Mean: %(meandopd).3f microns; RMS: %(rmsdopd).3f mic '\ %{"slopedopd":z_p[0],"meandopd":meandopd_p,"rmsdopd":np.std(total_dopd_plot_p)} ) Story.append(a) Story.append(PageBreak()) #============================================================================== # Table of the dOPD properties #============================================================================== if hasattr(p2vmreduced,'oi_vis_met_phase'): plotTitle(Story,"Differential OPD SC-[FT+MET] overview") Story.append(Spacer(1,10*mm)) nboot = 500 # number of bootstrap samples if p2vmreduced.polarsplit == False: plotSubtitle(Story,"Differential OPD expressed in microns. Bootstrapped error bars with "+str(nboot)+" samples. Error bars computed over the sequence of "+str(p2vmreduced.nframe_sc)+" frames over %.2f seconds." % np.max(p2vmreduced.time_sc)) Story.append(Spacer(1,10*mm)) avgdopd = [] diffopd = [] dopd_boot = np.zeros(nbase) dopd_std = np.zeros(nbase) for baseline in range(0,3): diffopd = p2vmreduced.total_dopd_avg[baseline,:] # opd_sc[baseline,:]-(opd_ft_int[baseline,:]+opd_met_int[baseline,:]) dopd = np.nanmean(diffopd) # Subtraction of slope z = np.polyfit(timescale_sc, diffopd, 1) poly = np.poly1d(z) dopd_slopesub = diffopd - poly(timescale_sc) dopd_std[baseline] = np.nanstd(dopd_slopesub) dopd_boot[baseline] = np.nanmean(boot.ci( dopd_slopesub, np.average, n_samples=nboot, output='errorbar')) avgdopd.append('%(val)5.3f ' % {"val":dopd}+ '\u00B1' +' %(err).3f' % {"err":dopd_boot[baseline]}) data = avgdopd data = np.vstack((p2vmreduced.basenames[0:3],data)) data = np.hstack(([["Baseline"],["Combined (microns)"]],data)) t=Table(data.tolist(),colWidths=4*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,10*mm)) avgdopd = [] diffopd = [] for baseline in range(3,nbase): diffopd = p2vmreduced.total_dopd_avg[baseline,:] # opd_sc[baseline,:]-(opd_ft_int[baseline,:]+opd_met_int[baseline,:]) dopd = np.nanmean(diffopd) # Subtraction of slope z = np.polyfit(timescale_sc, diffopd, 1) poly = np.poly1d(z) dopd_slopesub = diffopd - poly(timescale_sc) dopd_std[baseline] = np.nanstd(dopd_slopesub) dopd_boot[baseline] = np.nanmean(boot.ci( dopd_slopesub, np.average, n_samples=nboot, output='errorbar')) avgdopd.append('%(val).4f ' % {"val":dopd}+ '\u00B1' +' %(err).3f' % {"err":dopd_boot[baseline]}) data = avgdopd data = np.vstack((p2vmreduced.basenames[3:nbase],data)) data = np.hstack(([["Baseline"],["Combined (microns)"]],data)) t=Table(data.tolist(),colWidths=4*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,10*mm)) plotTitle(Story,"Mean dOPD SC-[FT+MET] error bar: %.1f nm" % (np.nanmean(dopd_boot)*1000)) plotTitle(Story,"Mean dOPD SC-[FT+MET] RMS: %.1f nm" % (np.nanmean(dopd_std)*1000)) plotSubtitle(Story,"Linear fit of dOPD is subtracted for error bar and RMS computation.") if hasattr(p2vmreduced,'oi_vis_met_phase') & (p2vmreduced.polarsplit == True): plotSubtitle(Story,"Differential OPD expressed in microns. Bootstrapped error bars with "+str(nboot)+" samples.") plotSubtitle(Story,"Error bars computed over the sequence of "+str(p2vmreduced.nframe_sc)+" frames over %.2f seconds." % np.max(p2vmreduced.time_sc)) Story.append(Spacer(1,10*mm)) avgdopds = [] avgdopdp = [] diffopd_s = [] diffopd_p = [] dopd_s = np.zeros(nbase) dopd_p = np.zeros(nbase) dopd_boot_s = np.zeros(nbase) dopd_boot_p = np.zeros(nbase) dopd_std_s = np.zeros(nbase) dopd_std_p = np.zeros(nbase) for baseline in range(0,nbase): diffopd_s = p2vmreduced.total_dopd_avg_s[baseline,:] dopd_s[baseline] = np.nanmean(diffopd_s) diffopd_p = p2vmreduced.total_dopd_avg_p[baseline,:] dopd_p[baseline] = np.nanmean(diffopd_p) # Subtraction of slope z_s = np.polyfit(timescale_sc, diffopd_s, 1) z_p = np.polyfit(timescale_sc, diffopd_s, 1) poly_s = np.poly1d(z_s) poly_p = np.poly1d(z_p) dopd_slopesub_s = diffopd_s - poly_s(timescale_sc) dopd_slopesub_p = diffopd_p - poly_p(timescale_sc) dopd_std_s[baseline] = np.nanstd(dopd_slopesub_s) dopd_boot_s[baseline] = np.nanmean(boot.ci( dopd_slopesub_s, np.average, n_samples=nboot, output='errorbar')) dopd_std_p[baseline] = np.nanstd(dopd_slopesub_p) dopd_boot_p[baseline] = np.nanmean(boot.ci( dopd_slopesub_p, np.average, n_samples=nboot, output='errorbar')) for baseline in range(0,3): avgdopds.append('%(val)5.3f '%{"val":dopd_s[baseline]}+'\u00B1'+ ' %(err).3f'%{"err":dopd_boot_s[baseline]}) avgdopdp.append('%(val)5.3f '%{"val":dopd_p[baseline]}+'\u00B1'+ ' %(err).3f'%{"err":dopd_boot_p[baseline]}) # average flux over wavelength data = avgdopds data = np.vstack((p2vmreduced.basenames[0:3],data)) data = np.vstack((data,avgdopdp)) data = np.hstack(([["Baseline"],["Polar. S (microns)"],["Polar. P (microns)"]],data)) t=Table(data.tolist(),colWidths=4*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,10*mm)) avgdopds = [] avgdopdp = [] for baseline in range(3,nbase): avgdopds.append('%(val)5.3f '%{"val":dopd_s[baseline]}+'\u00B1'+ ' %(err).3f' % {"err":dopd_boot_s[baseline]}) avgdopdp.append('%(val)5.3f '%{"val":dopd_p[baseline]}+'\u00B1'+ ' %(err).3f' % {"err":dopd_boot_p[baseline]}) # average flux over wavelength data = avgdopds data = np.vstack((p2vmreduced.basenames[3:nbase],data)) data = np.vstack((data,avgdopdp)) data = np.hstack(([["Baseline"],["Polar. S (microns)"],["Polar. P (microns)"]],data)) t=Table(data.tolist(),colWidths=4*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,10*mm)) plotTitle(Story,"Mean dOPD SC-[FT+MET] error bar: %.1f nm" % (np.nanmean([dopd_boot_s,dopd_boot_p])*1000)) plotTitle(Story,"Mean dOPD SC-[FT+MET] RMS: %.1f nm" % (np.nanmean([dopd_std_s,dopd_std_p])*1000)) plotSubtitle(Story,"Linear fit of dOPD is subtracted for error bar and RMS computation.") Story.append(PageBreak()) #============================================================================== # FDDL signals SC #============================================================================== boxsmooth = 50 # Moving average window for boxcar smoothing of the FDDL signals plotTitle(Story,"SC FDDL voltage (FDDL)") plotSubtitle(Story,"In volts vs. time (seconds). Smoothing by moving average of "+str(boxsmooth)+" samples applied for clarity.") if hasattr(p2vmreduced, 'fddl_sc_pos'): nframe_fddl = p2vmreduced.fddl_sc_pos.shape[0] fddl_sc = np.zeros((nframe_fddl,4)) resamp = nframe_fddl // 500 # take one value every resamp value for the plot for tel in range(0,4): fddl_sc[:,tel]=np.convolve(p2vmreduced.fddl_sc_pos[:,tel], np.ones((boxsmooth,))/boxsmooth, mode='same') for tel in range(0,4): amplitude_fddl = np.nanmax(fddl_sc[boxsmooth:nframe_fddl-boxsmooth,tel])-np.nanmin(fddl_sc[boxsmooth:nframe_fddl-boxsmooth,tel]) yminval = np.nanmin(fddl_sc[boxsmooth:nframe_fddl-boxsmooth,tel])-0.1*amplitude_fddl ymaxval = np.nanmax(fddl_sc[boxsmooth:nframe_fddl-boxsmooth,tel])+0.1*amplitude_fddl ticky = (ymaxval-yminval)/10. # in volts hsize = 16*cm vsize = 4.5*cm data = [] data.append( tuple(zip(p2vmreduced.fddl_time[::resamp], clipdata(fddl_sc[::resamp,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.fddl_time),max(p2vmreduced.fddl_time),None,yminval,ymaxval,None,\ hsize,vsize,None,ticky,'SC FDDL Telescope '+p2vmreduced.stations[tel] ) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # FDDL signals FT #============================================================================== boxsmooth = 50 # Moving average window for boxcar smoothing of the FDDL signals plotTitle(Story,"FT FDDL voltage (FDDL)") plotSubtitle(Story,"In volts vs. time (seconds). Smoothing by moving average of "+str(boxsmooth)+" samples applied for clarity.") if hasattr(p2vmreduced, 'fddl_sc_pos'): nframe_fddl = p2vmreduced.fddl_ft_pos.shape[0] fddl_ft = np.zeros((nframe_fddl,4)) for tel in range(0,4): fddl_ft[:,tel]=np.convolve(p2vmreduced.fddl_ft_pos[:,tel], np.ones((boxsmooth,))/boxsmooth, mode='same') for tel in range(0,4): amplitude_fddl = np.nanmax(fddl_ft[boxsmooth:nframe_fddl-boxsmooth,tel])-np.nanmin(fddl_ft[boxsmooth:nframe_fddl-boxsmooth,tel]) yminval = np.nanmin(fddl_ft[boxsmooth:nframe_fddl-boxsmooth,tel])-0.1*amplitude_fddl ymaxval = np.nanmax(fddl_ft[boxsmooth:nframe_fddl-boxsmooth,tel])+0.1*amplitude_fddl ticky = (ymaxval-yminval)/10. # in volts hsize = 16*cm vsize = 4.5*cm data = [] data.append( tuple(zip(p2vmreduced.fddl_time[::resamp], clipdata(fddl_ft[::resamp,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(p2vmreduced.fddl_time),max(p2vmreduced.fddl_time),None,yminval,ymaxval,None,\ hsize,vsize,None,ticky,'FT FDDL Telescope '+p2vmreduced.stations[tel]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Power spectral density of the FDDL volt signals #============================================================================== plotTitle(Story,"Periodogram of the SC FDDL voltage values (FDDL)") if welchpresent == True: fddlfreq = 1./(p2vmreduced.fddl_time[1]-p2vmreduced.fddl_time[0]) plotSubtitle(Story,"In volts RMS vs. Hz. FDDL data frequency = %.2f Hz." % fddlfreq) PSD_met_sc=[1,2,3,4] for tel in range(0,4): f_met_sc, PSD_met_sc[tel] = welch(p2vmreduced.fddl_sc_pos[:,tel], fs=fddlfreq, detrend='linear', nperseg=1024, scaling='spectrum') yminval = 0 hsize = 16*cm vsize = 4.5*cm for tel in range(0,4): data = [] ymaxval = np.nanmean(np.sqrt(PSD_met_sc[tel]))*5 ticky = (ymaxval-yminval)/5 data.append(tuple(zip(f_met_sc, clipdata(np.sqrt(PSD_met_sc[tel]), yminval,ymaxval)))) a = graphoutaxes( data, min(f_met_sc),max(f_met_sc),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'SC FDDL Telescope '+p2vmreduced.stations[tel]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Power spectral density of the FDDL volt signals #============================================================================== plotTitle(Story,"Periodogram of the FT FDDL voltage values (FDDL)") if welchpresent == True: fddlfreq = 1./(p2vmreduced.fddl_time[1]-p2vmreduced.fddl_time[0]) plotSubtitle(Story,"In volts RMS vs. Hz. FDDL data frequency = %.2f Hz." % fddlfreq) PSD_met_ft=[1,2,3,4] for tel in range(0,4): f_met_ft, PSD_met_ft[tel] = welch(p2vmreduced.fddl_ft_pos[:,tel], fs=fddlfreq, detrend='linear', nperseg=1024, scaling='spectrum') yminval = 0 hsize = 16*cm vsize = 4.5*cm for tel in range(0,4): data = [] ymaxval = np.nanmean(np.sqrt(PSD_met_ft[tel]))*5 ticky = (ymaxval-yminval)/5. data.append(tuple(zip(f_met_ft, clipdata(np.sqrt(PSD_met_ft[tel]), yminval,ymaxval)))) a = graphoutaxes( data, min(f_met_ft),max(f_met_ft),None,yminval,ymaxval,None,hsize,vsize,None,ticky,'FT FDDL Telescope '+p2vmreduced.stations[tel]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Create PDF report from Story #============================================================================== print("Create the PDF") gravi_visual_class.TITLE = "GRAVITY P2VMREDUCED Quality Control Report" gravi_visual_class.PAGEINFO = "P2VMREDUCED file: "+filename+".fits" reportname = filename+"-P2VMRED.pdf" doc = SimpleDocTemplate(reportname) doc.build(Story, onFirstPage=myFirstPage, onLaterPages=myLaterPages) print((" "+reportname)) #============================================================================== # MAIN PROGRAM #============================================================================== if __name__ == '__main__': filename='' # filename = '../GRAVITY.2015-11-16T00-06-33_0001' # filename = 'GRAVITY.2016-02-28T06-03-18_p2vmreduced' filename='../../GRAVITY.2016-03-25T00-08-57_p2vmreduced' if len(sys.argv) == 2 : filename=sys.argv[1] if filename == '' : filename=input(" Enter P2VMREDUCED file name (without .fits) : ") p2vmreduced = gravi_visual_class.P2vmreduced(filename+'.fits') produce_p2vmreduced_report(p2vmreduced,filename)