#! /usr/bin/env python3 # -*- coding: iso-8859-15 -*- """ Created on Tue Sep 15 16:59:06 2015 @author: kervella """ import matplotlib as mpl mpl.use('Agg') # necessary to be able to generate the reports even without X server import matplotlib.pyplot as plt import os from . import gravi_visual_class import datetime import numpy as np def generate_oi_filelist(directory): filelist = [] for file in os.listdir(directory): if file.endswith(".fits"): filelist.append(directory+'/'+os.path.splitext(file)[0]) return filelist #============================================================================== # Preparation of the report using reportlab #============================================================================== def produce_series_report(oilist,reportname): from reportlab.platypus import SimpleDocTemplate, Paragraph, Spacer, Table, TableStyle, Image, PageBreak from reportlab.lib import colors from reportlab.graphics.shapes import Drawing, String, Rect, Path from reportlab.graphics.charts.legends import LineLegend #from reportlab.graphics.charts.legends import Legend #from reportlab.platypus import * from reportlab.lib.styles import getSampleStyleSheet from reportlab.rl_config import defaultPageSize #from reportlab.lib.pagesizes import A4 from reportlab.lib.units import inch, cm, mm # from reportlab.pdfgen import canvas from reportlab.graphics.charts.lineplots import LinePlot # from reportlab.graphics.charts.barcharts import VerticalBarChart # import sys #import PIL import io # from cStringIO import StringIO # For insertion of matplotlib images (or others)in platypus, through cStringIO from pdfrw import PdfReader from pdfrw.buildxobj import pagexobj from pdfrw.toreportlab import makerl from reportlab.platypus import Flowable from reportlab.lib.enums import TA_LEFT, TA_CENTER, TA_RIGHT PAGE_HEIGHT=defaultPageSize[1] PAGE_WIDTH=defaultPageSize[0] styles = getSampleStyleSheet() class PdfImage(Flowable): """PdfImage wraps the first page from a PDF file as a Flowable which can be included into a ReportLab Platypus document. Based on the vectorpdf extension in rst2pdf (http://code.google.com/p/rst2pdf/)""" def __init__(self, filename_or_object, width=None, height=None, kind='direct'): from reportlab.lib.units import inch # If using StringIO buffer, set pointer to begining if hasattr(filename_or_object, 'read'): filename_or_object.seek(0) page = PdfReader(filename_or_object, decompress=False).pages[0] self.xobj = pagexobj(page) self.imageWidth = width self.imageHeight = height x1, y1, x2, y2 = self.xobj.BBox self._w, self._h = x2 - x1, y2 - y1 if not self.imageWidth: self.imageWidth = self._w if not self.imageHeight: self.imageHeight = self._h self.__ratio = float(self.imageWidth)/self.imageHeight if kind in ['direct','absolute'] or width==None or height==None: self.drawWidth = width or self.imageWidth self.drawHeight = height or self.imageHeight elif kind in ['bound','proportional']: factor = min(float(width)/self._w,float(height)/self._h) self.drawWidth = self._w*factor self.drawHeight = self._h*factor def wrap(self, aW, aH): return self.drawWidth, self.drawHeight def drawOn(self, canv, x, y, _sW=0): if _sW > 0 and hasattr(self, 'hAlign'): a = self.hAlign if a in ('CENTER', 'CENTRE', TA_CENTER): x += 0.5*_sW elif a in ('RIGHT', TA_RIGHT): x += _sW elif a not in ('LEFT', TA_LEFT): raise ValueError("Bad hAlign value " + str(a)) xobj = self.xobj xobj_name = makerl(canv._doc, xobj) xscale = self.drawWidth/self._w yscale = self.drawHeight/self._h x -= xobj.BBox[0] * xscale y -= xobj.BBox[1] * yscale canv.saveState() canv.translate(x, y) canv.scale(xscale, yscale) canv.doForm(xobj_name) canv.restoreState() Title = "GRAVITY OIFITS SERIES Quality Control Report" pageinfo = "SERIES first file: "+oilist.filelist[0]+".fits" def myFirstPage(canvas, doc): canvas.saveState() canvas.setFont('Helvetica',16) canvas.drawCentredString(PAGE_WIDTH/2.0, PAGE_HEIGHT-60, Title) canvas.setFont('Helvetica',9) canvas.drawString(inch, 0.75 * inch,"Page 1 / %s" % pageinfo) canvas.restoreState() def myLaterPages(canvas, doc): canvas.saveState() canvas.setFont('Helvetica', 9) canvas.drawString(inch, 0.75 * inch,"Page %d / %s" % (doc.page, pageinfo)) canvas.restoreState() def clipdata(datalist,minval,maxval): array_np = np.asarray(datalist) low_values_indices = array_np < minval # Where values are low array_np[low_values_indices] = minval # All low values set to ymin high_values_indices = array_np > maxval # Where values are high array_np[high_values_indices] = maxval # All high values set to ymax return array_np def graphoutaxes(data,xmin,xmax,ymin,ymax,sizex,sizey,ticky,title): drawing = Drawing(sizex,sizey) lp = LinePlot() lp.x = 0 lp.y = 0 lp.height = sizey lp.width = sizex lp.data = data lp.joinedLines = 1 lp.lines[0].strokeColor = colors.red lp.lines[1].strokeColor = colors.orange lp.lines[2].strokeColor = colors.blue lp.lines[3].strokeColor = colors.green lp.lines[4].strokeColor = colors.cyan lp.lines[5].strokeColor = colors.purple lp.strokeColor = colors.black lp.xValueAxis.valueMin = xmin lp.xValueAxis.valueMax = xmax lp.xValueAxis.visibleTicks = 1 lp.xValueAxis.visibleLabels = 1 lp.xValueAxis.labelTextFormat = '%.2f' #lp.xValueAxis.labelAxisMode = 'wavelength (mic)' lp.xValueAxis.labels.fontSize = 8 #lp.xValueAxis._text = "Time (s)" lp.yValueAxis.valueMin = ymin lp.yValueAxis.valueMax = ymax lp.yValueAxis.valueStep = ticky lp.yValueAxis.visibleTicks = 1 lp.yValueAxis.visibleLabels = 1 lp.yValueAxis.labelTextFormat = '%.3e' lp.yValueAxis.labels.fontSize = 8 lp.xValueAxis.visibleGrid = 1 lp.yValueAxis.visibleGrid = 1 # Title of sub plot drawing.add(lp) drawing.add(String(0.1*sizex, 0.8*sizey, title, fontSize=8)) return drawing def graphoutnoaxis(data,xmin,xmax,ymin,ymax,sizex,sizey,ticky,title): drawing = Drawing(sizex,sizey) lp = LinePlot() lp.x = 0 lp.y = 0 lp.height = sizey lp.width = sizex # # Patch for values outside limits # array_np = np.asarray(yvalues) # low_values_indices = array_np < ymin # Where values are low # array_np[low_values_indices] = ymin # All low values set to ymin # high_values_indices = array_np > ymax # Where values are high # array_np[high_values_indices] = ymax # All high values set to ymax # data = [tuple(zip(xvalues, array_np) )] lp.data = data lp.joinedLines = 1 lp.lines[0].strokeColor = colors.red lp.lines[1].strokeColor = colors.orange lp.lines[2].strokeColor = colors.blue lp.lines[3].strokeColor = colors.green lp.lines[4].strokeColor = colors.cyan lp.lines[5].strokeColor = colors.purple lp.strokeColor = colors.black lp.xValueAxis.valueMin = xmin lp.xValueAxis.valueMax = xmax lp.xValueAxis.visibleTicks = 0 lp.xValueAxis.visibleLabels = 0 # lp.xValueAxis.labelTextFormat = '%2.1f' # lp.xValueAxis.labels.fontSize = 8 lp.yValueAxis.valueMin = ymin lp.yValueAxis.valueMax = ymax lp.yValueAxis.valueStep = ticky lp.yValueAxis.visibleTicks = 0 lp.yValueAxis.visibleLabels = 0 # lp.xValueAxis.labelTextFormat = '%2.1f' # lp.yValueAxis.labels.fontSize = 8 lp.xValueAxis.visibleGrid = 1 lp.yValueAxis.visibleGrid = 1 # Title of sub plot drawing.add(lp) drawing.add(String(0.1*sizex, 0.8*sizey, title, fontSize=8)) return drawing doc = SimpleDocTemplate(reportname) Story = [Spacer(1,3*cm)] now = datetime.datetime.now() nbase = 6 ntel = 4 nfiles = oilist.nfiles nwave_sc = oilist.oi[0].nwave_sc nwave_ft = oilist.oi[0].nwave_ft base_list=["12", "13", "14", "23", "24", "34"] # style = styles["Heading2"] # p = Paragraph("General information", style) # Story.append(p) style = styles["Normal"] keywords = ["First file in series: ", "Processing date of file 0:", "Last file in series:", "Processing date of file %i:" % oilist.nfiles, "Number of files:", "Start MJD date:", "Total duration:", "Report date", "Data type:", "SC & Polar setup:", "SC NDIT x DIT:", "FT DIT freq.:", "FT object name", "FT object RA Dec Kmag", "Telescope stations" ] header0 = oilist.oi[0].header # header of the first file # The time scale of the series is set to start at zero timescale = (np.array(oilist.timescale)-oilist.timescale[0])*24. # in decimal hours # for i in range(0,nfiles): # timescale.append(float(oilist.oi[i].oi_flux_sc_time)) # timescale = np.array(timescale) if 'HIERARCH ESO DPR CATG' in oilist.oi[0].header: catg = header0['HIERARCH ESO DPR CATG'] else: catg='' if 'HIERARCH ESO PRO CATG' in oilist.oi[0].header: catg = header0['HIERARCH ESO PRO CATG'] else: catg='' if oilist.oi[0].single==True: singledual = 'SINGLE' else: singledual = 'DUAL' if 'HIERARCH ESO TPL NDIT SKY' in oilist.oi[0].header: skyframe = '*** SKY ***' else: skyframe = '*** OBJECT ***' if 'HIERARCH ESO FT ROBJ NAME' in oilist.oi[0].header: objname = header0['HIERARCH ESO FT ROBJ NAME'] objalpha = header0['HIERARCH ESO FT ROBJ ALPHA'] objdelta = header0['HIERARCH ESO FT ROBJ DELTA'] objmag = str(header0['HIERARCH ESO FT ROBJ MAG']) else: objname = '' objalpha = '' objdelta = '' objmag = '' ftdit = header0['HIERARCH ESO DET3 SEQ1 DIT'] ftfreq = 1./ftdit keyvals = [oilist.filelist[0]+'.fits', oilist.oi[0].header['DATE'], oilist.filelist[nfiles-1]+'.fits', oilist.oi[nfiles-1].header['DATE'], oilist.nfiles, oilist.timescale[0], "%.4f hours" % (timescale[nfiles-1]-timescale[0]), now.strftime("%Y-%m-%dT%H:%M:%S"), catg+', '+ skyframe, #tech, header0['HIERARCH ESO INS SPEC RES']+', '+ header0['HIERARCH ESO FT POLA MODE'], str(header0['HIERARCH ESO DET2 NDIT'])+' x '+str(header0['HIERARCH ESO DET2 SEQ1 DIT'])+' s', str(ftdit)+' s, '+'%(val).0f'%{"val":ftfreq} +' Hz', objname, objalpha +' '+ objdelta +' K='+ objmag , 'T1='+header0['HIERARCH ESO ISS CONF STATION1']+', T2='+ header0['HIERARCH ESO ISS CONF STATION2']+', T3='+ header0['HIERARCH ESO ISS CONF STATION3']+', T4='+ header0['HIERARCH ESO ISS CONF STATION4'] ] hsize = [5*cm, 11*cm] vsize = 7*mm textmatrix = list(zip(keywords,keyvals)) t = Table(textmatrix,colWidths=hsize,rowHeights=vsize) t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black), ('BOX', (0,0), (-1,-1), 2, colors.black),])) Story.append(t) Story.append(Spacer(1,2*cm)) style = styles["Heading2"] p = Paragraph("Organization of the report:", style) Story.append(p) style = styles["Heading3"] p = Paragraph("- Photometric quantities of SC, then of FT", style) Story.append(p) p = Paragraph("- Visibility amplitudes of SC, then of FT", style) Story.append(p) p = Paragraph("- Referenced phases of SC, then phases of FT", style) Story.append(p) p = Paragraph("- Statistics of interferometric quantities", style) Story.append(p) Story.append(PageBreak()) #============================================================================== # Flux signals #============================================================================== style = styles["Heading2"] p = Paragraph("Wavelength averaged flux values for SC vs. time (OI_FLUX)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Photo-e-/s vs. seconds. Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.", style) Story.append(p) Story.append(Spacer(1,5*mm)) if hasattr(oilist.oi[0],'oi_flux_sc') & (oilist.oi[0].polarsplit == False): flux_sc = np.zeros((nfiles,ntel,oilist.oi[0].nwave_sc)) # time, tel, wavelength for tel in range(0,ntel): for frame in range(0,nfiles): for wave in range(0,nwave_sc): flux_sc[frame,tel,wave] = oilist.oi[frame].oi_flux_sc[:,tel,wave] avgflux_sc = np.nanmean(flux_sc,axis=2) yminval = 0 ymaxval = np.nanmax(avgflux_sc)+0.1*np.abs(np.nanmax(avgflux_sc)) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 9*cm data = [] for tel in range(0,ntel): data.append( tuple(zip(timescale, clipdata(avgflux_sc[:,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(timescale),max(timescale),yminval,ymaxval,hsize,vsize,ticky,'Combined-polarization') Story.append(a) if hasattr(oilist.oi[0],'oi_flux_sc_s') & (oilist.oi[0].polarsplit == True): flux_sc_s = np.zeros((nfiles,ntel,nwave_sc)) # time, tel, wavelength flux_sc_p = np.zeros((nfiles,ntel,nwave_sc)) for tel in range(0,ntel): for frame in range(0,nfiles): for wave in range(0,nwave_sc): flux_sc_s[frame,tel,wave] = oilist.oi[frame].oi_flux_sc_s[:,tel,wave] flux_sc_p[frame,tel,wave] = oilist.oi[frame].oi_flux_sc_p[:,tel,wave] avgflux_sc_s = np.nanmean(flux_sc_s,axis=2) avgflux_sc_p = np.nanmean(flux_sc_p,axis=2) yminval = 0 ymaxval = np.nanmax(avgflux_sc_s)+0.1*np.abs(np.nanmax(avgflux_sc_s)) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 4.5*cm data = [] for tel in range(0,ntel): data.append( tuple(zip(timescale, clipdata(avgflux_sc_s[:,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(timescale),max(timescale),yminval,ymaxval,hsize,vsize,ticky,'S-polarization') Story.append(a) Story.append(Spacer(1,7*mm)) yminval = 0 ymaxval = np.nanmax(avgflux_sc_p)+0.1*np.abs(np.nanmax(avgflux_sc_p)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,ntel): data.append( tuple(zip(timescale, clipdata(avgflux_sc_p[:,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(timescale),max(timescale),yminval,ymaxval,hsize,vsize,ticky,'P-polarization') Story.append(a) Story.append(Spacer(1,5*mm)) style = styles["Heading2"] p = Paragraph("Wavelength averaged flux values for FT vs. time (OI_FLUX)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Photo-e-/s vs. seconds. Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.", style) Story.append(p) Story.append(Spacer(1,5*mm)) if hasattr(oilist.oi[0],'oi_flux_ft') & (oilist.oi[0].polarsplit == False): flux_ft = np.zeros((nfiles,ntel,nwave_ft)) # time, tel, wavelength for tel in range(0,ntel): for frame in range(0,nfiles): for wave in range(0,nwave_ft): flux_ft[frame,tel,wave] = oilist.oi[frame].oi_flux_ft[:,tel,wave] avgflux_ft = np.nanmean(flux_ft,axis=2) yminval = np.nanmin(avgflux_ft)-0.1*np.abs(np.nanmax(avgflux_ft)) ymaxval = np.nanmax(avgflux_ft)+0.1*np.abs(np.nanmax(avgflux_ft)) ticky = (ymaxval-yminval)/10 # in photo-e- hsize = 16*cm vsize = 9*cm data = [] for tel in range(0,ntel): data.append( tuple(zip(timescale, clipdata(avgflux_ft[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(timescale),max(timescale),yminval,ymaxval,hsize,vsize,ticky,'Combined-polarization') Story.append(a) if hasattr(oilist.oi[0],'oi_flux_ft_s') & (oilist.oi[0].polarsplit == True): flux_ft_s = np.zeros((nfiles,ntel,nwave_ft)) # time, tel, wavelength flux_ft_p = np.zeros((nfiles,ntel,nwave_ft)) for tel in range(0,ntel): for frame in range(0,nfiles): for wave in range(0,nwave_ft): flux_ft_s[frame,tel,wave] = oilist.oi[frame].oi_flux_ft_s[:,tel,wave] flux_ft_p[frame,tel,wave] = oilist.oi[frame].oi_flux_ft_p[:,tel,wave] avgflux_ft_s = np.nanmean(flux_ft_s,axis=2) avgflux_ft_p = np.nanmean(flux_ft_p,axis=2) yminval = np.nanmin(avgflux_ft_s)-0.1*np.abs(np.nanmax(avgflux_ft_s)) ymaxval = np.nanmax(avgflux_ft_s)+0.1*np.abs(np.nanmax(avgflux_ft_s)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,ntel): data.append( tuple(zip(timescale, clipdata(avgflux_ft_s[:,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(timescale),max(timescale),yminval,ymaxval,hsize,vsize,ticky,'S-polarization') Story.append(a) Story.append(Spacer(1,7*mm)) yminval = np.nanmin(avgflux_ft_p)-0.1*np.abs(np.nanmax(avgflux_ft_p)) ymaxval = np.nanmax(avgflux_ft_p)+0.1*np.abs(np.nanmax(avgflux_ft_p)) ticky = (ymaxval-yminval)/10 # in photo-e- data = [] for tel in range(0,ntel): data.append( tuple(zip(timescale, clipdata(avgflux_ft_p[:,tel],yminval,ymaxval)) )) a = graphoutaxes( data, min(timescale),max(timescale),yminval,ymaxval,hsize,vsize,ticky,'P-polarization') Story.append(a) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of SC OIFLUX for the four telescopes #============================================================================== style = styles["Heading2"] p = Paragraph("SC flux relative variation per telescope (OIFLUX SC)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Wavelength channel number in horizontal axis, number of file (time) in vertical axis. The color scale is in percentage of the median value (given in photo-electrons/s).", style) Story.append(p) if oilist.oi[0].polarsplit == False: medianflux = np.zeros((ntel)) mediansub = np.zeros((nfiles,ntel,nwave_sc)) for tel in range(0,ntel): medianflux[tel] = np.nanmedian(flux_sc[:,tel,:]) mediansub[:,tel,:] = flux_sc[:,tel,:] - medianflux[tel] valmin = 0#np.min([np.nanpercentile(100.*flux_sc[:,tel,:]/medianflux[tel],2) for tel in range(0,ntel)]) valmax = np.max([np.nanpercentile(100.*flux_sc[:,tel,:]/medianflux[tel],98) for tel in range(0,ntel)]) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) for tel in range(0,ntel): im = axarr[tel].imshow(100.*flux_sc[:,tel,:]/medianflux[tel],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel].set_title('SC OIFLUX Telescope %i - Comb - Median = %.1e e-/s' % (tel+1, medianflux[tel]), 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) if oilist.oi[0].polarsplit == True: medianfluxs = np.zeros((ntel)) medianfluxp = np.zeros((ntel)) for tel in range(0,ntel): medianfluxs[tel] = np.nanmedian(flux_sc_s[:,tel,:]) medianfluxp[tel] = np.nanmedian(flux_sc_p[:,tel,:]) valmin = np.min([np.nanpercentile(100.*flux_sc_s[:,tel,:]/medianfluxs[tel],2) for tel in range(0,ntel)]) valmax = np.max([np.nanpercentile(100.*flux_sc_s[:,tel,:]/medianfluxs[tel],98) for tel in range(0,ntel)]) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) for tel in range(0,ntel): im = axarr[tel].imshow(100.*flux_sc_s[:,tel,:]/medianfluxs[tel],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel].set_title('SC OIFLUX Telescope %i - S - Median = %.1e e-/s' % (tel+1, medianfluxs[tel]), 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) valmin = np.min([np.nanpercentile(100.*flux_sc_p[:,tel,:]/medianfluxp[tel],2) for tel in range(0,ntel)]) valmax = np.max([np.nanpercentile(100.*flux_sc_p[:,tel,:]/medianfluxp[tel],98) for tel in range(0,ntel)]) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) for tel in range(0,ntel): im = axarr[tel].imshow(100.*flux_sc_p[:,tel,:]/medianfluxp[tel],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel].set_title('SC OIFLUX Telescope %i - P - Median = %.1e e-/s' % (tel+1, medianfluxs[tel]), 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) style = styles["Heading2"] p = Paragraph("SC flux difference S-P waterfall (OIFLUX SC)", style) Story.append(p) style = styles["Normal"] p = Paragraph("In photo-electrons/s, wavelength channel number in horizontal axis, number of frame (time) in vertical axis. The color scale is in percentage of the median.", style) Story.append(p) meansub = np.zeros((nfiles,ntel,nwave_sc)) flux_diff_sp = np.zeros((nfiles,ntel,nwave_sc)) medianflux = np.zeros((ntel)) for tel in range(0,ntel): flux_diff_sp[:,tel,:] = flux_sc_s[:,tel,:] - flux_sc_p[:,tel,:] medianflux[tel] = np.nanmedian(flux_diff_sp[:,tel,:]) valmin = np.min([np.nanpercentile(100.*flux_diff_sp[:,tel,:]/medianflux[tel],2) for tel in range(0,ntel)]) valmax = np.max([np.nanpercentile(100.*flux_diff_sp[:,tel,:]/medianflux[tel],98) for tel in range(0,ntel)]) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) for tel in range(0,ntel): im = axarr[tel].imshow(100.*flux_diff_sp[:,tel,:]/medianflux[tel],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel].set_title('SC OIFLUX Tel. %i - [S-P] - Median = %.1e e-/s' % (tel+1, medianflux[tel]), 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of FT OIFLUX for the four telescopes #============================================================================== style = styles["Heading2"] p = Paragraph("FT flux waterfall view per telescope (OIFLUX FT)", style) Story.append(p) style = styles["Normal"] p = Paragraph("In photo-electrons/s, wavelength channel number in horizontal axis, number of frame (time) in vertical axis. Median subtracted.", style) Story.append(p) if oilist.oi[0].polarsplit == False: mediansub = np.zeros((nfiles,ntel,nwave_ft)) medianflux = np.zeros((ntel)) for tel in range(0,ntel): medianflux[tel] = np.nanmedian(flux_ft[:,tel,:]) mediansub[:,tel,:] = flux_ft[:,tel,:] - medianflux[tel] valmin = np.nanpercentile(mediansub,2) valmax = np.nanpercentile(mediansub,98) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(mediansub[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('FT Tel. 1 - Comb - Median = %.1e' % medianflux[0], fontsize=8) axarr[1].imshow(mediansub[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('FT Tel. 2 - Comb - Median = %.1e' % medianflux[1], fontsize=8) axarr[2].imshow(mediansub[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('FT Tel. 3 - Comb - Median = %.1e' % medianflux[2], fontsize=8) axarr[3].imshow(mediansub[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('FT Tel. 4 - Comb - Median = %.1e' % medianflux[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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if oilist.oi[0].polarsplit == True: mediansubs = np.zeros((nfiles,ntel,nwave_ft)) mediansubsp = np.zeros((nfiles,ntel,nwave_ft)) for tel in range(0,4): mediansubs[:,tel,:] = flux_ft_s[:,tel,:] - np.nanmedian(flux_ft_s[:,tel,:]) mediansubsp[:,tel,:] = flux_ft_s[:,tel,:] -flux_ft_p[:,tel,:] mediansubsp[:,tel,:] -= np.nanmean(mediansubsp[:,tel,:]) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=4, ncols=2, sharex=True, sharey=True, figsize=(5,8)) for tel in range(0,4): valmin = np.nanpercentile(mediansubs,2) valmax = np.nanpercentile(mediansubs,98) im = axarr[tel,0].imshow(mediansubs[:,tel,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel,0].set_title('FT Tel.'+str(tel+1)+'-S - Median = %.1e' % np.nanmedian(flux_ft_s[:,tel,:]), fontsize=7) valmin = np.nanpercentile(mediansubsp,1) valmax = np.nanpercentile(mediansubsp,99) im = axarr[tel,1].imshow(mediansubsp[:,tel,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[tel,1].set_title('FT Tel. '+str(tel+1)+'-[S-P] - Median = %.1e' % np.nanmedian(flux_ft_p[:,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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Average photometric spectra #============================================================================== style = styles["Heading2"] p = Paragraph("Average spectrum of SC (OI_FLUX)", style) Story.append(p) style = styles["Normal"] p = Paragraph("In photo-e-/s vs. wavelength (microns). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.", style) Story.append(p) Story.append(Spacer(1,5*mm)) if hasattr(oilist.oi[0],'oi_flux_sc') & (oilist.oi[0].polarsplit == False): avgflux = np.zeros((ntel,nwave_sc),dtype='d') for tel in range(0,4): avgflux[tel,:]=np.nanmean(flux_sc[:,tel,:],axis=0) # average flux over time 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(oilist.oi[0].wave_sc, clipdata(avgflux[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,'Combined-polarization') Story.append(a) if hasattr(oilist.oi[0],'oi_flux_sc_s') & (oilist.oi[0].polarsplit == True): avgfluxs = np.zeros((ntel,nwave_sc),dtype='d') avgfluxp = np.zeros((ntel,nwave_sc),dtype='d') for tel in range(0,4): avgfluxs[tel,:]=np.nanmean(flux_sc_s[:,tel,:],axis=0) # average flux over time avgfluxp[tel,:]=np.nanmean(flux_sc_p[:,tel,:],axis=0) # average flux over time 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- hsize = 16*cm vsize = 4.5*cm data = [] for tel in range(0,4): data.append( tuple(zip(oilist.oi[0].wave_sc, clipdata(avgfluxs[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,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(oilist.oi[0].wave_sc, clipdata(avgfluxp[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,'P-polarization') Story.append(a) Story.append(Spacer(1,5*mm)) style = styles["Heading2"] p = Paragraph("Average spectrum of FT (OI_FLUX)", style) Story.append(p) style = styles["Normal"] p = Paragraph("In photo-e-/s vs. wavelength (microns). Tel1 = red, Tel2 = orange, Tel3 = blue, Tel4 = green.", style) Story.append(p) Story.append(Spacer(1,5*mm)) if hasattr(oilist.oi[0],'oi_flux_ft_time') & (oilist.oi[0].polarsplit == False): avgflux = np.zeros((4,nwave_ft),dtype='d') for tel in range(0,4): avgflux[tel,:]=np.nanmean(flux_ft[:,tel,:],axis=0) # average flux over time hsize = 16*cm vsize = 9*cm 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- data = [] for tel in range(0,4): data.append( tuple(zip(oilist.oi[0].wave_ft, clipdata(avgflux[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,'SCombined-polarization') Story.append(a) if hasattr(oilist.oi[0],'oi_flux_ft_time_s') & (oilist.oi[0].polarsplit == True): avgfluxs = np.zeros((ntel,oilist.oi[0].nwave_ft),dtype='d') avgfluxp = np.zeros((ntel,oilist.oi[0].nwave_ft),dtype='d') for tel in range(0,4): avgfluxs[tel,:]=np.nanmean(oilist.oi[0].oi_flux_ft_s[:,tel,:],axis=0) # average flux over time avgfluxp[tel,:]=np.nanmean(oilist.oi[0].oi_flux_ft_p[:,tel,:],axis=0) # average flux over time hsize = 16*cm vsize = 4.5*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(oilist.oi[0].wave_ft, clipdata(avgfluxs[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,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(oilist.oi[0].wave_ft, clipdata(avgfluxp[tel,:],yminval,ymaxval)) )) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,'P-polarization') Story.append(a) Story.append(PageBreak()) # ------------------------------------ Science combiner and fringe tracker interferometric quantities ------------------------------------ #============================================================================== # Time average spectrum of SC VISAMP #============================================================================== style = styles["Heading2"] p = Paragraph("Time averaged SC visibility amp. vs. wavelength (VISAMP SC)", style) Story.append(p) yminval = 0. ymaxval = 2. ticky = 0.2 # visibility tick data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visamp') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In visibility vs. wavelength (microns).", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp = np.zeros((nbase,nwave_sc),dtype='d') visamp_sc = np.zeros((nfiles,nbase,nwave_sc)) # time, tel, wavelength for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_sc): visamp_sc[frame,baseline,wave] = oilist.oi[frame].oi_vis_sc_visamp[baseline,wave] for baseline in range(0,nbase): data = [] visamp[baseline,:] = np.nanmean(visamp_sc[:,baseline,:], axis=0) yminval = np.nanmin(visamp[baseline,:])-0.05 ymaxval = np.nanmax(visamp[baseline,:])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(oilist.oi[0].wave_sc, clipdata(visamp[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_sc), max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visamp_s') & (oilist.oi[0].polarsplit == True): visamp_sc_s = np.zeros((nfiles,nbase,nwave_sc)) # time, tel, wavelength visamp_sc_p = np.zeros((nfiles,nbase,nwave_sc)) for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_sc): visamp_sc_s[frame,baseline,wave] = oilist.oi[frame].oi_vis_sc_visamp_s[baseline,wave] visamp_sc_p[frame,baseline,wave] = oilist.oi[frame].oi_vis_sc_visamp_p[baseline,wave] style = styles["Normal"] p = Paragraph("In visibility vs. wavelength (microns). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp_s = np.zeros((nbase,nwave_sc),dtype='d') visamp_p = np.zeros((nbase,nwave_sc),dtype='d') for baseline in range(0,nbase): visamp_s[baseline,:] = np.nanmean(visamp_sc_s[:,baseline,:], axis=0) visamp_p[baseline,:] = np.nanmean(visamp_sc_p[:,baseline,:], axis=0) for baseline in range(0,nbase): yminval = np.nanmin([visamp_s[baseline,:],visamp_p[baseline,:]])-0.05 ymaxval = np.nanmax([visamp_s[baseline,:],visamp_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(oilist.oi[0].wave_sc, clipdata(visamp_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(oilist.oi[0].wave_sc, clipdata(visamp_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_sc), max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Wavelength average SC VISAMP vs. time #============================================================================== style = styles["Heading2"] p = Paragraph("Wavelength averaged SC visibility amp. vs. time (VISAMP SC)", style) Story.append(p) data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visamp') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In visibility vs. time (decimal hours).", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp_sc_wavg = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visamp_sc_wavg[baseline,:] = np.nanmean(visamp_sc[:,baseline,:], axis=1) for baseline in range(0,nbase): data=[] yminval = np.nanmin(visamp_sc_wavg[baseline,:])-0.05 ymaxval = np.nanmax(visamp_sc_wavg[baseline,:])-0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(timescale, clipdata(visamp_sc_wavg[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,\ 'Baseline '+base_list[baseline]+' Mean = %.4e RMS = %.2e ' % \ (np.nanmean(visamp_sc_wavg[baseline,:]),np.std(visamp_sc_wavg[baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visamp_s') & (oilist.oi[0].polarsplit == True): style = styles["Normal"] p = Paragraph("In visibility vs. time (decimal hours). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp_sc_wavg_s = np.zeros((nbase,nfiles),dtype='d') visamp_sc_wavg_p = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visamp_sc_wavg_s[baseline,:] = np.nanmean(visamp_sc_s[:,baseline,:], axis=1) visamp_sc_wavg_p[baseline,:] = np.nanmean(visamp_sc_p[:,baseline,:], axis=1) for baseline in range(0,nbase): yminval = np.nanmin([visamp_sc_wavg_s[baseline,:],visamp_sc_wavg_p[baseline,:]])-0.05 ymaxval = np.nanmax([visamp_sc_wavg_s[baseline,:],visamp_sc_wavg_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(timescale, clipdata(visamp_sc_wavg_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(timescale, clipdata(visamp_sc_wavg_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,\ 'Baseline '+base_list[baseline]+' Mean S = %.4e RMS S = %.2e Mean P = %.4e RMS P = %.2e ' % \ (np.nanmean(visamp_sc_wavg_s[baseline,:]),np.std(visamp_sc_wavg_s[baseline,:]),\ np.nanmean(visamp_sc_wavg_p[baseline,:]),np.std(visamp_sc_wavg_p[baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of SC VISAMP #============================================================================== style = styles["Heading2"] p = Paragraph("SC visibility amplitude waterfall view (VISAMP SC)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Wavelength channel number in horizontal axis, number of file in sequence in vertical axis.", style) Story.append(p) if hasattr(oilist.oi[0],'oi_vis_sc_visamp') & (oilist.oi[0].polarsplit == False): valmin = np.nanpercentile(visamp_sc,1) valmax = np.nanpercentile(visamp_sc,99) valmin -= 0.2*np.abs(valmax-valmin) valmax += 0.2*np.abs(valmax-valmin) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(visamp_sc[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC VISAMP Baseline 12 - Comb', fontsize=8) axarr[1].imshow(visamp_sc[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC VISAMP Baseline 13 - Comb', fontsize=8) axarr[2].imshow(visamp_sc[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC VISAMP Baseline 14 - Comb', fontsize=8) axarr[3].imshow(visamp_sc[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC VISAMP Baseline 23 - Comb', fontsize=8) axarr[4].imshow(visamp_sc[:,4,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[4].set_title('SC VISAMP Baseline 24 - Comb', fontsize=8) axarr[5].imshow(visamp_sc[:,5,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[5].set_title('SC VISAMP Baseline 34 - Comb', fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if hasattr(oilist.oi[0],'oi_vis_sc_visamp_s') & (oilist.oi[0].polarsplit == True): valmin = np.nanpercentile(visamp_sc_s,1) valmax = np.nanpercentile(visamp_sc_s,99) valmin -= 0.2*np.abs(valmax-valmin) valmax += 0.2*np.abs(valmax-valmin) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(visamp_sc_s[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC VISAMP Baseline 12 - S', fontsize=8) axarr[1].imshow(visamp_sc_s[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC VISAMP Baseline 13 - S', fontsize=8) axarr[2].imshow(visamp_sc_s[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC VISAMP Baseline 14 - S', fontsize=8) axarr[3].imshow(visamp_sc_s[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC VISAMP Baseline 23 - S', fontsize=8) axarr[4].imshow(visamp_sc_s[:,4,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[4].set_title('SC VISAMP Baseline 24 - S', fontsize=8) axarr[5].imshow(visamp_sc_s[:,5,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[5].set_title('SC VISAMP Baseline 34 - S', fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) style = styles["Heading2"] p = Paragraph("SC visibility amplitude waterfall view (VISAMP SC)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.", style) Story.append(p) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(visamp_sc_p[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[0].set_title('SC VISAMP Baseline 12 - P', fontsize=8) axarr[1].imshow(visamp_sc_p[:,1,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[1].set_title('SC VISAMP Baseline 13 - P', fontsize=8) axarr[2].imshow(visamp_sc_p[:,2,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[2].set_title('SC VISAMP Baseline 14 - P', fontsize=8) axarr[3].imshow(visamp_sc_p[:,3,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[3].set_title('SC VISAMP Baseline 23 - P', fontsize=8) axarr[4].imshow(visamp_sc_p[:,4,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[4].set_title('SC VISAMP Baseline 24 - P', fontsize=8) axarr[5].imshow(visamp_sc_p[:,5,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[5].set_title('SC VISAMP Baseline 34 - P', fontsize=8) fig.subplots_adjust(hspace=.3) # Color bar plot fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Time average spectrum of FT VISAMP #============================================================================== style = styles["Heading2"] p = Paragraph("Time averaged FT visibility amp. vs. wavelength (VISAMP FT)", style) Story.append(p) yminval = 0. ymaxval = 2. ticky = 0.2 # visibility tick data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visamp') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In visibility vs. wavelength (microns).", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp = np.zeros((nbase,nwave_ft),dtype='d') visamp_ft = np.zeros((nfiles,nbase,nwave_ft)) # time, tel, wavelength for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_ft): visamp_ft[frame,baseline,wave] = oilist.oi[frame].oi_vis_ft_visamp[baseline,wave] for baseline in range(0,nbase): data = [] visamp[baseline,:] = np.nanmean(visamp_ft[:,baseline,:], axis=0) yminval = np.nanmin(visamp[baseline,:])-0.05 ymaxval = np.nanmax(visamp[baseline,:])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(oilist.oi[0].wave_ft, clipdata(visamp[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_ft), max(oilist.oi[0].wave_ft),yminval,ymaxval,hsize,vsize,ticky, \ 'Baseline '+base_list[baseline]+' Mean = %.4e RMS = %.2e ' % \ (np.nanmean(visamp[baseline,:]),np.std(visamp[baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visamp_s') & (oilist.oi[0].polarsplit == True): visamp_ft_s = np.zeros((nfiles,nbase,nwave_ft)) # time, tel, wavelength visamp_ft_p = np.zeros((nfiles,nbase,nwave_ft)) for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_ft): visamp_ft_s[frame,baseline,wave] = oilist.oi[frame].oi_vis_ft_visamp_s[baseline,wave] visamp_ft_p[frame,baseline,wave] = oilist.oi[frame].oi_vis_ft_visamp_p[baseline,wave] style = styles["Normal"] p = Paragraph("In visibility vs. wavelength (microns). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp_s = np.zeros((nbase,nwave_ft),dtype='d') visamp_p = np.zeros((nbase,nwave_ft),dtype='d') for baseline in range(0,nbase): visamp_s[baseline,:] = np.nanmean(visamp_ft_s[:,baseline,:], axis=0) visamp_p[baseline,:] = np.nanmean(visamp_ft_p[:,baseline,:], axis=0) for baseline in range(0,nbase): yminval = np.nanmin([visamp_s[baseline,:],visamp_p[baseline,:]])-0.05 ymaxval = np.nanmax([visamp_s[baseline,:],visamp_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(oilist.oi[0].wave_ft, clipdata(visamp_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(oilist.oi[0].wave_ft, clipdata(visamp_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_ft), max(oilist.oi[0].wave_ft),yminval,ymaxval,hsize,vsize,ticky,\ 'Baseline '+base_list[baseline]+' Mean S = %.4e RMS S = %.2e Mean P = %.4e RMS P = %.2e ' % \ (np.nanmean(visamp_s[baseline,:]),np.std(visamp_s[baseline,:]),np.nanmean(visamp_p[baseline,:]),np.std(visamp_p[baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Wavelength average FT VISAMP vs. time #============================================================================== style = styles["Heading2"] p = Paragraph("Wavelength averaged FT visibility amp. vs. time (VISAMP SC)", style) Story.append(p) data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visamp') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In visibility vs. time (decimal hours).", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp_ft_wavg = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visamp_ft_wavg[baseline,:] = np.nanmean(visamp_ft[:,baseline,:], axis=1) for baseline in range(0,nbase): data=[] yminval = np.nanmin(visamp_ft_wavg[baseline,:])-0.05 ymaxval = np.nanmax(visamp_ft_wavg[baseline,:])-0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(timescale, clipdata(visamp_ft_wavg[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visamp_s') & (oilist.oi[0].polarsplit == True): style = styles["Normal"] p = Paragraph("In visibility vs. time (decimal hours). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visamp_ft_wavg_s = np.zeros((nbase,nfiles),dtype='d') visamp_ft_wavg_p = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visamp_ft_wavg_s[baseline,:] = np.nanmean(visamp_ft_s[:,baseline,:], axis=1) visamp_ft_wavg_p[baseline,:] = np.nanmean(visamp_ft_p[:,baseline,:], axis=1) for baseline in range(0,nbase): yminval = np.nanmin([visamp_ft_wavg_s[baseline,:],visamp_ft_wavg_p[baseline,:]])-0.05 ymaxval = np.nanmax([visamp_ft_wavg_s[baseline,:],visamp_ft_wavg_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(timescale, clipdata(visamp_ft_wavg_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(timescale, clipdata(visamp_ft_wavg_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of FT VISAMP #============================================================================== style = styles["Heading2"] p = Paragraph("FT visibility amplitude waterfall view (VISAMP FT)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Wavelength channel number in horizontal axis, number of file in sequence in vertical axis.", style) Story.append(p) if hasattr(oilist.oi[0],'oi_vis_ft_visamp') & (oilist.oi[0].polarsplit == False): valmin = np.nanpercentile(visamp_ft,2) valmax = np.nanpercentile(visamp_ft,98) valmin -= 0.2*np.abs(valmax-valmin) valmax += 0.2*np.abs(valmax-valmin) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) im = axarr[0].imshow(visamp_ft[:,0,:],vmin=valmin,vmax=valmax,cmap='cubehelix', interpolation='nearest', aspect='auto') for baseline in range(0,nbase): im = axarr[baseline].imshow(visamp_ft[:,baseline,:],vmin=valmin,vmax=valmax,\ cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline].set_title('FT VISAMP Baseline '+base_list[baseline]+' - Comb', fontsize=6) 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if hasattr(oilist.oi[0],'oi_vis_ft_visamp_s') & (oilist.oi[0].polarsplit == True): valmin = np.nanpercentile([visamp_ft_s,visamp_ft_p],2) valmax = np.nanpercentile([visamp_ft_s,visamp_ft_p],98) valmin -= 0.2*np.abs(valmax-valmin) valmax += 0.2*np.abs(valmax-valmin) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=2, sharex=True, figsize=(5,7)) for baseline in range(0,nbase): im = axarr[baseline,0].imshow(visamp_ft_s[:,baseline,:],vmin=valmin,vmax=valmax,\ cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,0].set_title('FT VISAMP Baseline '+base_list[baseline]+' - S', fontsize=6) axarr[baseline,1].imshow(visamp_ft_s[:,baseline,:],vmin=valmin,vmax=valmax,\ cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline,1].set_title('FT VISAMP Baseline '+base_list[baseline]+' - P', fontsize=6) 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Time average spectrum of SC VISPHI #============================================================================== style = styles["Heading2"] p = Paragraph("Time averaged SC phase vs. wavelength (VISPHI SC)", style) Story.append(p) yminval = 0. ymaxval = 2. ticky = 0.2 # visibility tick data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visphi') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In OPD (microns) vs. wavelength (microns). Mean value subtracted.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi = np.zeros((nbase,nwave_sc),dtype='d') visphi_sc = np.zeros((nfiles,nbase,nwave_sc)) # time, tel, wavelength for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_sc): visphi_sc[frame,baseline,wave] = oilist.oi[frame].oi_vis_sc_visphi[baseline,wave] visphi_sc = np.unwrap(visphi_sc,axis=0) # unwrap in time for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_sc): visphi_sc[frame,baseline,wave] *= oilist.oi[frame].wave_sc[wave]/(2*np.pi) # get OPD instead of phase for baseline in range(0,nbase): data = [] visphi[baseline,:] = np.nanmean(visphi_sc[:,baseline,:], axis=0)-np.nanmean(visphi_sc[:,baseline,:]) yminval = np.nanmin(visphi[baseline,:])-0.05 ymaxval = np.nanmax(visphi[baseline,:])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(oilist.oi[0].wave_sc, clipdata(visphi[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_sc),\ max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]+' Mean = %.6e mic' %\ (np.nanmean(visphi_sc[:,baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visphi_s') & (oilist.oi[0].polarsplit == True): visphi_sc_s = np.zeros((nfiles,nbase,nwave_sc)) # time, tel, wavelength visphi_sc_p = np.zeros((nfiles,nbase,nwave_sc)) for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_sc): visphi_sc_s[frame,baseline,wave] = oilist.oi[frame].oi_vis_sc_visphi_s[baseline,wave] visphi_sc_p[frame,baseline,wave] = oilist.oi[frame].oi_vis_sc_visphi_p[baseline,wave] visphi_sc_s = np.unwrap(visphi_sc_s,axis=0) # unwrap in time visphi_sc_p = np.unwrap(visphi_sc_p,axis=0) # unwrap in time for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_sc): visphi_sc_s[frame,baseline,wave] *= oilist.oi[frame].wave_sc[wave]/(2*np.pi) # get OPD instead of phase visphi_sc_p[frame,baseline,wave] *= oilist.oi[frame].wave_sc[wave]/(2*np.pi) style = styles["Normal"] p = Paragraph("In OPD (microns) vs. wavelength (microns). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi_s = np.zeros((nbase,nwave_sc),dtype='d') visphi_p = np.zeros((nbase,nwave_sc),dtype='d') for baseline in range(0,nbase): visphi_s[baseline,:] = np.nanmean(visphi_sc_s[:,baseline,:], axis=0)-np.nanmean(visphi_sc_s[:,baseline,:]) visphi_p[baseline,:] = np.nanmean(visphi_sc_p[:,baseline,:], axis=0)-np.nanmean(visphi_sc_p[:,baseline,:]) for baseline in range(0,nbase): yminval = np.nanmin([visphi_s[baseline,:],visphi_p[baseline,:]])-0.05 ymaxval = np.nanmax([visphi_s[baseline,:],visphi_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(oilist.oi[0].wave_sc, clipdata(visphi_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(oilist.oi[0].wave_sc, clipdata(visphi_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_sc), max(oilist.oi[0].wave_sc),yminval,ymaxval,hsize,vsize,ticky, \ 'Baseline '+base_list[baseline]+' Mean S = %.6e mic Mean P = %.6e mic' % \ (np.nanmean(visphi_sc_s[:,baseline,:]),np.nanmean(visphi_sc_p[:,baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Wavelength average SC VISPHI vs. time #============================================================================== style = styles["Heading2"] p = Paragraph("SC REFERENCED phase vs. time (VISPHI SC)", style) Story.append(p) data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visphi') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("OPD (microns) vs. time (decimal hours). The mean values for S and P over the sequence are subtracted.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi_sc_wavg = np.zeros((nfiles,nbase),dtype='d') print(visphi_sc.shape) for baseline in range(0,nbase): visphi_sc_wavg[:,baseline] = np.nanmean(visphi_sc[:,baseline,:], axis=1) for baseline in range(0,nbase): data=[] yminval = np.nanpercentile(visphi_sc_wavg[baseline,:]-np.nanmean(visphi_sc_wavg[baseline,:]),2)-0.05 ymaxval = np.nanpercentile(visphi_sc_wavg[baseline,:]-np.nanmean(visphi_sc_wavg[baseline,:]),98)+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(timescale, clipdata(visphi_sc_wavg[baseline,:]-np.nanmean(visphi_sc_wavg[baseline,:]),yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,\ base_list[baseline]+" Mean = %.3e microns"%np.nanmean(visphi_sc_wavg[baseline,:])) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_sc_visphi_s') & (oilist.oi[0].polarsplit == True): style = styles["Normal"] p = Paragraph("OPD (microns) vs. time (decimal hours). Red = S, Orange = P. Mean over sequence subracted.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi_sc_wavg_s = np.zeros((nbase,nfiles),dtype='d') visphi_sc_wavg_p = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visphi_sc_wavg_s[baseline,:] = np.nanmean(visphi_sc_s[:,baseline,:], axis=1) visphi_sc_wavg_p[baseline,:] = np.nanmean(visphi_sc_p[:,baseline,:], axis=1) yminval = np.nanpercentile([visphi_sc_wavg_s[baseline,:]-np.nanmean(visphi_sc_wavg_s[baseline,:]),\ visphi_sc_wavg_p[baseline,:]-np.nanmean(visphi_sc_wavg_p[baseline,:])],2)-0.05 ymaxval = np.nanpercentile([visphi_sc_wavg_s[baseline,:]-np.nanmean(visphi_sc_wavg_s[baseline,:]),\ visphi_sc_wavg_p[baseline,:]-np.nanmean(visphi_sc_wavg_p[baseline,:])],98)+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(timescale, clipdata(visphi_sc_wavg_s[baseline,:]-np.nanmean(visphi_sc_wavg_s[baseline,:]),yminval,ymaxval)))) data.append(list(zip(timescale, clipdata(visphi_sc_wavg_p[baseline,:]-np.nanmean(visphi_sc_wavg_p[baseline,:]),yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,\ "Base "+base_list[baseline]+" Mean S = %.6e mic RMS S = %.3e mic Mean P = %.6e mic RMS P = %.3e mic"%\ (np.nanmean(visphi_sc_s[:,baseline,:]),np.std(visphi_sc_wavg_s[baseline,:]),\ np.nanmean(visphi_sc_p[:,baseline,:]),np.std(visphi_sc_wavg_p[baseline,:]))) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Image (waterfall) display of SC VISPHI #============================================================================== style = styles["Heading2"] p = Paragraph("SC REFERENCED phase waterfall view (VISPHI SC)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Wavelength channel number in horizontal axis, number of file in sequence in vertical axis. Color scale in microns (median subtracted).", style) Story.append(p) valmin = np.zeros(nbase) valmax = np.zeros(nbase) if hasattr(oilist.oi[0],'oi_vis_sc_visphi') & (oilist.oi[0].polarsplit == False): imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) valmin = np.nanpercentile(visphi_sc-np.nanmedian(visphi_sc),2) valmax = np.nanpercentile(visphi_sc-np.nanmedian(visphi_sc),98) for baseline in range(0,nbase): im = axarr[baseline].imshow(visphi_sc[:,baseline,:],vmin=valmin,vmax=valmax, cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline].set_title('SC visphi Baseline '+base_list[baseline]+'-Comb Median = %.3e mic' %\ np.nanmean(visphi_sc[:,baseline,:]), 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.colorbar(im, fraction=0.15,aspect=15,orientation='horizontal',shrink=0.7,pad=0.05) #fig.tight_layout() #plt.colorbar() widthfig = 20 * cm heightfig= 22 * cm plt.savefig(imgdata, format='PDF', dpi=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) if hasattr(oilist.oi[0],'oi_vis_sc_visphi_s') & (oilist.oi[0].polarsplit == True): imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) visphi_sc_s_medsub = np.copy(visphi_sc_s) for baseline in range(0,nbase): visphi_sc_s_medsub[:,baseline,:] = visphi_sc_s[:,baseline,:]-np.nanmedian(visphi_sc_s[:,baseline,:]) valmin = np.nanpercentile(visphi_sc_s_medsub,2) valmax = np.nanpercentile(visphi_sc_s_medsub,98) for baseline in range(0,nbase): im = axarr[baseline].imshow(visphi_sc_s_medsub[:,baseline,:],vmin=valmin,vmax=valmax, \ cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline].set_title('SC visphi Baseline '+base_list[baseline]+'-S Median = %.3e mic' % \ np.nanmedian(visphi_sc_s[:,baseline,:]), 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) style = styles["Heading2"] p = Paragraph("SC REFERENCED phase waterfall view (VISPHI SC)", style) Story.append(p) style = styles["Normal"] p = Paragraph("Wavelength channel number in horizontal axis, number of frame in sequence in vertical axis.", style) Story.append(p) imgdata = io.StringIO() plt.close('all') plt.rc('xtick', labelsize=6) plt.rc('ytick', labelsize=6) fig, axarr = plt.subplots(nrows=6, ncols=1, sharex=True, figsize=(5,7)) visphi_sc_p_medsub = np.copy(visphi_sc_p) for baseline in range(0,nbase): visphi_sc_p_medsub[:,baseline,:] = visphi_sc_p[:,baseline,:]-np.nanmedian(visphi_sc_p[:,baseline,:]) valmin = np.percentile(visphi_sc_p_medsub,2) valmax = np.nanpercentile(visphi_sc_p_medsub,98) for baseline in range(0,nbase): im = axarr[baseline].imshow(visphi_sc_p_medsub[:,baseline,:],vmin=valmin,vmax=valmax, \ cmap='cubehelix', interpolation='nearest', aspect='auto') axarr[baseline].set_title('SC visphi Baseline '+base_list[baseline]+'-P Median = %.3e mic' % \ np.nanmedian(visphi_sc_p[:,baseline,:]), 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=250, bbox_inches='tight') pi1 = PdfImage(imgdata, width=widthfig, height=heightfig, kind='bound') plt.close('all') Story.append(pi1) Story.append(PageBreak()) #============================================================================== # Time average spectrum of FT VISPHI #============================================================================== style = styles["Heading2"] p = Paragraph("Time averaged FT phase vs. wavelength (VISPHI FT)", style) Story.append(p) yminval = 0. ymaxval = 2. ticky = 0.2 # visibility tick data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visphi') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In OPD (microns) vs. wavelength (microns).", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi = np.zeros((nbase,nwave_ft),dtype='d') visphi_ft = np.zeros((nfiles,nbase,nwave_ft)) # time, tel, wavelength for baseline in range(0,nbase): for wave in range(0,nwave_ft): for frame in range(0,nfiles): visphi_ft[frame,baseline,wave] = oilist.oi[frame].oi_vis_ft_visphi[baseline,wave] visphi_ft = np.unwrap(visphi_ft,axis=0) # unwrap in time for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_ft): visphi_ft[frame,baseline,wave] *= oilist.oi[frame].wave_ft[wave]/(2*np.pi) # get OPD for baseline in range(0,nbase): data = [] visphi[baseline,:] = np.nanmean(visphi_ft[:,baseline,:], axis=0) yminval = np.nanmin(visphi[baseline,:])-0.05 ymaxval = np.nanmax(visphi[baseline,:])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(oilist.oi[0].wave_ft, clipdata(visphi[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_ft), max(oilist.oi[0].wave_ft),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visphi_s') & (oilist.oi[0].polarsplit == True): visphi_ft_s = np.zeros((nfiles,nbase,nwave_ft)) # time, tel, wavelength visphi_ft_p = np.zeros((nfiles,nbase,nwave_ft)) for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_ft): visphi_ft_s[frame,baseline,wave] = oilist.oi[frame].oi_vis_ft_visphi_s[baseline,wave] visphi_ft_p[frame,baseline,wave] = oilist.oi[frame].oi_vis_ft_visphi_p[baseline,wave] visphi_ft_s = np.unwrap(visphi_ft_s,axis=0) # unwrap in time visphi_ft_p = np.unwrap(visphi_ft_p,axis=0) # unwrap in time for baseline in range(0,nbase): for frame in range(0,nfiles): for wave in range(0,nwave_ft): visphi_ft_s[frame,baseline,wave] *= oilist.oi[frame].wave_ft[wave]/(2*np.pi) # get OPD instead of phase visphi_ft_p[frame,baseline,wave] *= oilist.oi[frame].wave_ft[wave]/(2*np.pi) style = styles["Normal"] p = Paragraph("In OPD (microns) vs. wavelength (microns). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi_s = np.zeros((nbase,nwave_ft),dtype='d') visphi_p = np.zeros((nbase,nwave_ft),dtype='d') for baseline in range(0,nbase): visphi_s[baseline,:] = np.nanmean(visphi_ft_s[:,baseline,:], axis=0) visphi_p[baseline,:] = np.nanmean(visphi_ft_p[:,baseline,:], axis=0) for baseline in range(0,nbase): yminval = np.nanmin([visphi_s[baseline,:],visphi_p[baseline,:]])-0.05 ymaxval = np.nanmax([visphi_s[baseline,:],visphi_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(oilist.oi[0].wave_ft, clipdata(visphi_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(oilist.oi[0].wave_ft, clipdata(visphi_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(oilist.oi[0].wave_ft), max(oilist.oi[0].wave_ft),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) #============================================================================== # Wavelength average FT VISPHI vs. time #============================================================================== style = styles["Heading2"] p = Paragraph("Wavelength averaged FT phase vs. time (VISPHI FT)", style) Story.append(p) data=[] Story.append(Spacer(1,2*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visphi') & (oilist.oi[0].polarsplit == False): style = styles["Normal"] p = Paragraph("In OPD (microns) vs. time (decimal hours).", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi_ft_wavg = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visphi_ft_wavg[baseline,:] = np.nanmean(visphi_ft[:,baseline,:], axis=1) for baseline in range(0,nbase): data=[] yminval = np.nanmin(visphi_ft_wavg[baseline,:])-0.05 ymaxval = np.nanmax(visphi_ft_wavg[baseline,:])-0.05 ticky = (ymaxval - yminval)/5. # visibility tick data.append(list(zip(timescale, clipdata(visphi_ft_wavg[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) if hasattr(oilist.oi[0],'oi_vis_ft_visphi_s') & (oilist.oi[0].polarsplit == True): style = styles["Normal"] p = Paragraph("In OPD (microns) vs. time (decimal hours). Red = polarization S, Orange = polarization P.", style) Story.append(p) Story.append(Spacer(1,5*mm)) hsize=16*cm vsize=3*cm visphi_ft_wavg_s = np.zeros((nbase,nfiles),dtype='d') visphi_ft_wavg_p = np.zeros((nbase,nfiles),dtype='d') for baseline in range(0,nbase): visphi_ft_wavg_s[baseline,:] = np.nanmean(visphi_ft_s[:,baseline,:], axis=1) visphi_ft_wavg_p[baseline,:] = np.nanmean(visphi_ft_p[:,baseline,:], axis=1) for baseline in range(0,nbase): yminval = np.nanmin([visphi_ft_wavg_s[baseline,:],visphi_ft_wavg_p[baseline,:]])-0.05 ymaxval = np.nanmax([visphi_ft_wavg_s[baseline,:],visphi_ft_wavg_p[baseline,:]])+0.05 ticky = (ymaxval - yminval)/5. # visibility tick data = [] data.append(list(zip(timescale, clipdata(visphi_ft_wavg_s[baseline,:],yminval,ymaxval)))) data.append(list(zip(timescale, clipdata(visphi_ft_wavg_p[baseline,:],yminval,ymaxval)))) a = graphoutaxes( data, min(timescale), max(timescale),yminval,ymaxval,hsize,vsize,ticky,'Baseline '+base_list[baseline]) Story.append(a) Story.append(Spacer(1,7*mm)) Story.append(PageBreak()) if False: #============================================================================== # Table of statistical values #============================================================================== style = styles["Heading2"] p = Paragraph("Statistics of interferometric quantities", style) Story.append(p) style = styles["Normal"] p = Paragraph("Averaged quantities over wavelengths.", style) Story.append(p) Story.append(Spacer(1,2*cm)) style = styles["Normal"] keywords = ["Total duration:", "Report date", "Baselines", "SC S visibilities mean", "SC S visibilities RMS", "SC P visibilities mean", "SC P visibilities RMS", "SC REF S phase mean (OPD)", "SC REF S phase RMS (OPD)", "SC REF P phase mean (OPD)", "SC REF P phase RMS (OPD)", "FT S visibilities mean", "FT S visibilities RMS", "FT P visibilities mean", "FT P visibilities RMS", "FT S phase mean (OPD)", "FT S phase RMS (OPD)", "FT P phase mean (OPD)", "FT P phase RMS (OPD)" ] m_v_ft_s = [np.nanmean(visamp_ft_wavg_s[baseline,:]) for baseline in range(0,nbase)] s_v_ft_s = [np.std(visamp_ft_wavg_s[baseline,:]) for baseline in range(0,nbase)] m_v_ft_p = [np.nanmean(visamp_ft_wavg_p[baseline,:]) for baseline in range(0,nbase)] s_v_ft_p = [np.std(visamp_ft_wavg_p[baseline,:]) for baseline in range(0,nbase)] m_p_ft_s = [np.nanmean(visphi_ft_wavg_s[baseline,:]) for baseline in range(0,nbase)] s_p_ft_s = [np.std(visphi_ft_wavg_s[baseline,:]) for baseline in range(0,nbase)] m_p_ft_p = [np.nanmean(visphi_ft_wavg_p[baseline,:]) for baseline in range(0,nbase)] s_p_ft_p = [np.std(visphi_ft_wavg_p[baseline,:]) for baseline in range(0,nbase)] m_v_sc_s = [np.nanmean(visamp_sc_wavg_s[baseline,:]) for baseline in range(0,nbase)] s_v_sc_s = [np.std(visamp_sc_wavg_s[baseline,:]) for baseline in range(0,nbase)] m_v_sc_p = [np.nanmean(visamp_sc_wavg_p[baseline,:]) for baseline in range(0,nbase)] s_v_sc_p = [np.std(visamp_sc_wavg_p[baseline,:]) for baseline in range(0,nbase)] m_p_sc_s = [np.nanmean(visphi_sc_wavg_s[baseline,:]) for baseline in range(0,nbase)] s_p_sc_s = [np.std(visphi_sc_wavg_s[baseline,:]) for baseline in range(0,nbase)] m_p_sc_p = [np.nanmean(visphi_sc_wavg_p[baseline,:]) for baseline in range(0,nbase)] s_p_sc_p = [np.std(visphi_sc_wavg_p[baseline,:]) for baseline in range(0,nbase)] keyvals = ["%.4f hours" % (timescale[nfiles-1]-timescale[0]), now.strftime("%Y-%m-%dT%H:%M:%S"), ' '.join(["%s" % base_list[i] for i in range(0,nbase)]), ' '.join(["%.4f" % m_v_sc_s[i] for i in range(0,nbase)]), ' '.join(["%.4f" % s_v_sc_s[i] for i in range(0,nbase)]), ' '.join(["%.4f" % m_v_sc_p[i] for i in range(0,nbase)]), ' '.join(["%.4f" % s_v_sc_p[i] for i in range(0,nbase)]), (' '.join(["%.4f" % (m_p_sc_s[i]/1000) for i in range(0,nbase)]))+' (mm)', (' '.join(["%.2f" % (s_p_sc_s[i]*1000) for i in range(0,nbase)]))+' (nm)', (' '.join(["%.4f" % (m_p_sc_p[i]/1000) for i in range(0,nbase)]))+' (mm)', (' '.join(["%.2f" % (s_p_sc_p[i]*1000) for i in range(0,nbase)]))+' (nm)', ' '.join(["%.4f" % m_v_ft_s[i] for i in range(0,nbase)]), ' '.join(["%.4f" % s_v_ft_s[i] for i in range(0,nbase)]), ' '.join(["%.4f" % m_v_ft_p[i] for i in range(0,nbase)]), ' '.join(["%.4f" % s_v_ft_p[i] for i in range(0,nbase)]), (' '.join(["%.2f" % (m_p_ft_s[i]*1000) for i in range(0,nbase)]))+' (nm)', (' '.join(["%.2f" % (s_p_ft_s[i]*1000) for i in range(0,nbase)]))+' (nm)', (' '.join(["%.2f" % (m_p_ft_p[i]*1000) for i in range(0,nbase)]))+' (nm)', (' '.join(["%.2f" % (s_p_ft_p[i]*1000) for i in range(0,nbase)]))+' (nm)' ] hsize = [5*cm, 11*cm] vsize = 7*mm textmatrix = list(zip(keywords,keyvals)) t = Table(textmatrix,colWidths=hsize,rowHeights=vsize) t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black), ('BOX', (0,0), (-1,-1), 2, colors.black),])) Story.append(t) Story.append(PageBreak()) doc.build(Story, onFirstPage=myFirstPage, onLaterPages=myLaterPages) print((" "+reportname)) #============================================================================== # MAIN PROGRAM #============================================================================== if __name__ == '__main__': # %reset -f is used to clean all the variables from the memory #directory = '../2015-09-13-stability' #os.getcwd() # > phases in radians ! #reportname = '2015-09-13-stability-series-full.pdf' # > phases in radians ! #directory = '../2015-09-14-stability' #os.getcwd() # > phases in radians ! #reportname = '2015-09-14-stability-series-full.pdf' # > phases in radians ! #directory = '../2015-09-16-stab-v2' #os.getcwd() #reportname = '2015-09-16-stability-v2-000-967.pdf' #reportname = '2015-09-16-stability-v2-000-500.pdf' #reportname = '2015-09-16-stability-v2-600-967.pdf' directory = '.' reportname = 'OIFITS_series.pdf' filelist = generate_oi_filelist(directory) #filelist.sort() # process the files in alphabetical order # Key words in the header to be checked for file selection # fits_keys = {'HIERARCH ESO PRO CATG': ['VIS_DUAL_SCI_RAW', 'VIS_SINGLE_SCI_RAW','VIS_SINGLE_CAL_RAW', 'VIS_DUAL_CAL_RAW'],\ # 'HIERARCH ESO FT POLA MODE': ['COMBINED','SPLIT']} fits_keys = {'HIERARCH ESO PRO CATG': ['VIS_DUAL_SCI_RAW', 'VIS_SINGLE_SCI_RAW','VIS_SINGLE_CAL_RAW', 'VIS_DUAL_CAL_RAW']} oilist = gravi_visual_class.Oilist(filelist[:],fits_keys) nfiles = oilist.nfiles #print oilist.oi[1].nwave_sc produce_series_report(oilist,reportname)