# -*- coding: utf-8 -*- """ Created on Mon Sep 14 11:41:26 2015 @author: kervella """ try: import pyfits except: from astropy.io import fits as pyfits import numpy as np import datetime import math import matplotlib.pyplot as plt #import os #import gc # =============================================================================== # Python plotting utilities for reports # =============================================================================== from reportlab.platypus import SimpleDocTemplate, Paragraph, Spacer, Table from reportlab.platypus import TableStyle, Image, PageBreak from reportlab.graphics.shapes import Drawing, String, Rect, Path from reportlab.graphics.charts.legends import LineLegend from reportlab.lib import colors from reportlab.rl_config import defaultPageSize from reportlab.lib.units import inch, cm, mm from reportlab.lib.styles import getSampleStyleSheet from reportlab.rl_config import defaultPageSize from reportlab.graphics.charts.lineplots import LinePlot from reportlab.graphics.charts.barcharts import VerticalBarChart from reportlab.graphics.charts.lineplots import LinePlot, ScatterPlot from reportlab.graphics.shapes import Drawing, String, Rect, Path # 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 import io from io import StringIO # Global variable base_list=["12", "13", "14", "23", "24", "34"] tel_list = np.array([[0,1],[0,2],[0,3],[1,2],[1,3],[2,3]]) nbase = 6 ntel = 4 styles = getSampleStyleSheet() # These can be modified by gravi_visual_class.TITLE # and gravi_visual_class.PAGEINFO TITLE="default title" PAGEINFO="default page info" PAGE_HEIGHT=defaultPageSize[1] PAGE_WIDTH=defaultPageSize[0] def round_to_n(x, n=1): if not x: return 0 power = -int(math.floor(math.log10(abs(x)))) + (n - 1) factor = (10 ** power) return round(x * factor) / factor def computeVisAmp(vis,flux): visamp = np.abs(vis.field('VISDATA'))**2 - vis.field('VISERR').real - vis.field('VISERR').imag for base in range(nbase): f1f2 = flux.field('FLUX')[tel_list[base][0]::ntel,:] * flux.field('FLUX')[tel_list[base][1]::ntel,:] visamp[base::nbase] /= f1f2 return np.sqrt(np.abs(visamp)) 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.nan_to_num(np.asarray(datalist)) minval = np.nan_to_num(minval)+0 maxval = np.nan_to_num(maxval)+0 try: low_values_indices = array_np < minval # Where values are low array_np[low_values_indices] = minval # All low values set to ymin except: print ('cannot clipdata') array_np[:] = minval try: high_values_indices = array_np > maxval # Where values are high array_np[high_values_indices] = maxval # All high values set to ymax except: print ('cannot clipdata') array_np[:] = maxval return array_np def transbarchartout(data,xlabels,ymin,ymax,stepy,sizex,sizey,colorchoice0,colorchoice1,title=''): drawing = Drawing(sizex,sizey) bc = VerticalBarChart() bc.x = sizey*.1 bc.y = sizey*.1 bc.height = sizey*.8 bc.width = sizex*.8 #print data bc.data = data bc.strokeColor = colors.black bc.bars[0].fillColor = colorchoice0 bc.bars[1].fillColor = colorchoice1 bc.valueAxis.valueMin = ymin bc.valueAxis.valueMax = ymax bc.valueAxis.valueStep = stepy bc.valueAxis.visibleGrid = 1 bc.valueAxis.labels.fontSize = 8 bc.categoryAxis.labels.boxAnchor = 'ne' bc.categoryAxis.labels.dx = -2*mm #bc.categoryAxis.labels.dy = -2*mm bc.categoryAxis.labels.angle = 90 bc.categoryAxis.categoryNames = xlabels bc.categoryAxis.labels.fontSize = 8 drawing.add(bc) drawing.add(String(0.1*sizex, 0.8*sizey, title, fontSize=10)) return drawing def graphoutnoxaxes(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 = 0 lp.xValueAxis.visibleLabels = 0 #lp.xValueAxis.labelTextFormat = '%.2e' #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 = '%.2e' 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 graphoutaxes(data,xmin,xmax,formx,ymin,ymax,formy,sizex,sizey,tickx,ticky,title): formx = '%.2f' if formx is None else formx formy = '%.2f' if formy is None else formy tickx = (xmax-xmin)/5 if tickx is None else tickx ticky = (ymax-ymin)/5 if ticky is None else ticky if (tickx <= 0 or ticky<=0): raise ValueError("Negative or null tick") 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.valueStep = tickx lp.xValueAxis.visibleTicks = 1 lp.xValueAxis.visibleLabels = 1 lp.xValueAxis.labelTextFormat = formx #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 = formy lp.yValueAxis.labels.fontSize = 8 lp.xValueAxis.visibleGrid = 1 lp.yValueAxis.visibleGrid = 1 # Title of sub plot drawing.add(lp) if title is not None: 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): tickx = (xmax-xmin)/5 ticky = (ymax-ymin)/5 if ticky is None else ticky 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.valueStep = tickx lp.xValueAxis.visibleTicks = 0 lp.xValueAxis.visibleLabels = 0 lp.yValueAxis.valueMin = ymin lp.yValueAxis.valueMax = ymax lp.yValueAxis.valueStep = ticky lp.yValueAxis.visibleTicks = 0 lp.yValueAxis.visibleLabels = 0 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 graphscatteraxes(data,xmin,xmax,formx,ymin,ymax,formy,sizex,sizey,tickx,ticky,dotsize,title): dotsize = 2 if dotsize is None else dotsize formx = '%.2f' if formx is None else formx formy = '%.2f' if formy is None else formy tickx = (xmax-xmin)/5 if tickx is None else tickx ticky = (ymax-ymin)/5 if ticky is None else ticky drawing = Drawing(sizex,sizey) sp = ScatterPlot() sp.x = 0 sp.y = 0 sp.height = sizey sp.width = sizex sp.data = data sp.joinedLines = False sp.lines[0].symbol.size = dotsize sp.lines[0].strokeColor = colors.red sp.lines[0].symbol.kind = 'Circle' sp.lines[1].symbol.size = dotsize sp.lines[1].strokeColor = colors.orange sp.lines[1].symbol.kind = 'Circle' sp.lines[2].symbol.size = dotsize sp.lines[2].strokeColor = colors.blue sp.lines[2].symbol.kind = 'Circle' sp.lines[3].symbol.size = dotsize sp.lines[3].strokeColor = colors.green sp.lines[3].symbol.kind = 'Circle' sp.lines[4].symbol.size = dotsize sp.lines[4].strokeColor = colors.cyan sp.lines[4].symbol.kind = 'Circle' sp.lines[5].symbol.size = dotsize sp.lines[5].strokeColor = colors.purple sp.lines[5].symbol.kind = 'Circle' sp.lineLabelFormat = None sp.xLabel='' sp.yLabel='' sp.strokeColor = colors.black sp.xValueAxis.valueMin = xmin sp.xValueAxis.valueMax = xmax sp.xValueAxis.valueStep = tickx sp.xValueAxis.visibleTicks = 1 sp.xValueAxis.visibleLabels = 1 sp.xValueAxis.labels.fontSize = 8 sp.yValueAxis.valueMin = ymin sp.yValueAxis.valueMax = ymax sp.yValueAxis.valueStep = ticky sp.yValueAxis.visibleTicks = 1 sp.yValueAxis.visibleLabels = 1 sp.yValueAxis.labels.fontSize = 8 sp.xValueAxis.visibleGrid = 1 sp.yValueAxis.visibleGrid = 1 sp.xValueAxis.labelTextFormat = formx sp.yValueAxis.labelTextFormat = formy # Title of sub plot drawing.add(sp) if title is not None: drawing.add(String(0.1*sizex, 0.8*sizey, title, fontSize=8)) return drawing def graphaxesplt(data,xmin,xmax,ymin,ymax,sizex,sizey,ticky,title): imgdata = io.StringIO() plt.close('all') figwidth = sizex figheight = sizey plt.figure(figsize=(figwidth,figheight)) colorlist = ['red', 'orange', 'blue', 'green', 'cyan', 'purple'] for base in range(0,nbase): if data.polarsplit == False: plt.plot(data[:,0], data[:,1], color=colorlist[base], markersize=3, markeredgecolor='gray') plt.set_xlim([xmin,xmax]) plt.set_ylim([ymin,ymax]) return imgdata def graphoutdarkaxes(xvalues,yvalues,xmin,xmax,ymin,ymax,sizex,sizey,title): drawing = Drawing(sizex,sizey) lp = LinePlot() lp.x = 0 lp.y = 0 lp.height = sizey lp.width = sizex # Patch to remove 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.strokeColor = colors.black lp.xValueAxis.valueMin = xmin lp.xValueAxis.valueMax = xmax lp.xValueAxis.visibleTicks = 1 lp.xValueAxis.visibleLabels = 1 lp.xValueAxis.labels.angle = 0 # lp.xValueAxis.labelTextFormat = '%2.1f' # lp.xValueAxis.labelAxisMode = 'wavelength (mic)' lp.xValueAxis.labels.fontSize = 10 lp.yValueAxis.valueMin = ymin lp.yValueAxis.valueMax = ymax lp.yValueAxis.visibleTicks = 1 lp.yValueAxis.visibleLabels = 1 lp.yValueAxis.labelTextFormat = '%2.1f' lp.yValueAxis.labels.fontSize = 10 lp.xValueAxis.visibleGrid = 1 lp.yValueAxis.visibleGrid = 1 # Title of sub plot drawing.add(lp) drawing.add(String(0.1*sizex, 0.85*sizey, title, fontSize=10)) return drawing def stepbarchartout(data,xlabels,ymin,ymax,stepy,sizex,sizey,colorchoice): drawing = Drawing(sizex,sizey) bc = VerticalBarChart() bc.x = sizey*.1 bc.y = sizey*.1 bc.height = sizey*.8 bc.width = sizex*.8 bc.data = data bc.strokeColor = colors.black bc.bars[0].fillColor = colorchoice bc.valueAxis.valueMin = ymin bc.valueAxis.valueMax = ymax bc.valueAxis.valueStep = stepy bc.valueAxis.visibleGrid = 1 bc.valueAxis.labels.fontSize = 8 bc.categoryAxis.labels.boxAnchor = 'ne' bc.categoryAxis.labels.dx = -2*mm bc.categoryAxis.labels.angle = 90 bc.categoryAxis.categoryNames = xlabels bc.categoryAxis.labels.fontSize = 8 drawing.add(bc) return drawing def plotTitle(Story,title,spacer=Spacer(1,-2*mm)): print (title) Story.append(Paragraph(title, style = styles["Heading2"])) if spacer is not None: Story.append(spacer) def plotSubtitle(Story,title,spacer=Spacer(1,5*mm)): Story.append(Paragraph(title, style = styles["Normal"])) if spacer is not None: Story.append(spacer) 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() #============================================================================== # Page about general information #============================================================================== def plotSummary (Story, filename, product, onTarget=False): if 'HIERARCH ESO DPR CATG' in product.header: catg = product.header['HIERARCH ESO DPR CATG'] elif 'HIERARCH ESO PRO CATG' in product.header: catg = product.header['HIERARCH ESO PRO CATG'] else: catg='' objname = get_key_withdefault(product.header,'HIERARCH ESO FT ROBJ NAME','') objalpha = get_key_withdefault(product.header,'HIERARCH ESO FT ROBJ ALPHA','') objdelta = get_key_withdefault(product.header,'HIERARCH ESO FT ROBJ DELTA','') objmag = str(get_key_withdefault(product.header,'HIERARCH ESO FT ROBJ MAG','')) sobjname = get_key_withdefault(product.header,'HIERARCH ESO INS SOBJ NAME','') sobjdalpha = get_key_withdefault(product.header,'HIERARCH ESO INS SOBJ X','0.0') sobjddelta = get_key_withdefault(product.header,'HIERARCH ESO INS SOBJ Y','0.0') sobjmag = str(get_key_withdefault(product.header,'HIERARCH ESO INS SOBJ MAG','0.0')) if hasattr(product,'time_ft'): ftfreq = 1./np.mean(np.diff(product.time_ft)) elif (product.header['HIERARCH ESO DET3 SEQ1 DIT'] == 0.00085): ftfreq = 909 elif (product.header['HIERARCH ESO DET3 SEQ1 DIT'] == 0.00084): ftfreq = 909 elif (product.header['HIERARCH ESO DET3 SEQ1 DIT'] == 0.003): ftfreq = 303 else: ftfreq = 0 if hasattr(product,'stations'): config = 'GV1='+product.telnames[0]+'='+product.stations[0]+ \ ', GV2='+product.telnames[1]+'='+product.stations[1]+ \ ', GV3='+product.telnames[2]+'='+product.stations[2]+ \ ', GV4='+product.telnames[3]+'='+product.stations[3] else: config = '' if 'HIERARCH ESO ISS AMBI FWHM START' in product.header and 'HIERARCH ESO QC TAU0 OPDC12' in product.header: # mean tau0 from GRAVITY FT gravitau0 = np.nanmean([product.header['HIERARCH ESO QC TAU0 OPDC12'], product.header['HIERARCH ESO QC TAU0 OPDC13'], product.header['HIERARCH ESO QC TAU0 OPDC14'], product.header['HIERARCH ESO QC TAU0 OPDC23'], product.header['HIERARCH ESO QC TAU0 OPDC24'], product.header['HIERARCH ESO QC TAU0 OPDC34']]) ambi = '{0:.2f} arcsec, {1:.1f} m/s, {2:.1f} ms, {3:.1f} ms'.format(product.header['HIERARCH ESO ISS AMBI FWHM START'],\ product.header['HIERARCH ESO ISS AMBI WINDSP'], product.header['HIERARCH ESO ISS AMBI TAU0 START']*1e3,\ gravitau0*1e3) trackratio = '{0:.1f}%, {1:.1f}%, {2:.1f}%, {3:.1f}%, {4:.1f}%, {5:.1f}%'.format(product.header['HIERARCH ESO QC TRACKING_RATIO_FT12'], product.header['HIERARCH ESO QC TRACKING_RATIO_FT13'], product.header['HIERARCH ESO QC TRACKING_RATIO_FT14'], product.header['HIERARCH ESO QC TRACKING_RATIO_FT23'], product.header['HIERARCH ESO QC TRACKING_RATIO_FT24'], product.header['HIERARCH ESO QC TRACKING_RATIO_FT34']) lambda0 = 2190 # lambda0 for the FT to compute rms in nm trackresiduals = 'B12: {0:.0f}, B13: {1:.0f}, B14: {2:.0f}, B23: {3:.0f}, B24: {4:.0f}, B34: {5:.0f}'.format(product.header['HIERARCH ESO QC PHASE_FT12 RMS']*lambda0/(2*np.pi), product.header['HIERARCH ESO QC PHASE_FT13 RMS']*lambda0/(2*np.pi), product.header['HIERARCH ESO QC PHASE_FT14 RMS']*lambda0/(2*np.pi), product.header['HIERARCH ESO QC PHASE_FT23 RMS']*lambda0/(2*np.pi), product.header['HIERARCH ESO QC PHASE_FT24 RMS']*lambda0/(2*np.pi), product.header['HIERARCH ESO QC PHASE_FT34 RMS']*lambda0/(2*np.pi)) else: ambi = '' trackratio = '' trackresiduals = '' keywords = ["file name:", #"OB start date:", "Observing date:", "Processing/report date:", "Product category:", "SC & Polar setup:", "SC NDIT x DIT:", "FT DIT freq.:"] keyvals = [filename+'.fits', #product.header['HIERARCH ESO OBS START'], product.header['DATE-OBS'], product.header['DATE']+' '+datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S"), catg, product.header['HIERARCH ESO INS SPEC RES']+', '+product.header['HIERARCH ESO INS POLA MODE']+', '+ product.header['HIERARCH ESO FT POLA MODE'], str(product.header['HIERARCH ESO DET2 NDIT'])+' x '+str(product.header['HIERARCH ESO DET2 SEQ1 DIT'])+' s', str(product.header['HIERARCH ESO DET3 SEQ1 DIT'])+' s, '+'%.0f'%ftfreq +' Hz'] if onTarget is True: keywords += [ "FT object name", "FT object RA Dec Kmag", "SC object name", "SC object dRA dDec Kmag", "Telescope stations", "Seeing Wind T0(V) T0(K)", "FT tracking ratio", "FT OPD residuals (nm)"] keyvals += [ objname, objalpha +' '+ objdelta +' K='+ objmag , sobjname, '%(val).3f'%{"val":float(sobjdalpha)/1000.} +'" '+ '%(val).3f'%{"val":float(sobjddelta)/1000.} +'" K='+ sobjmag, config, ambi, trackratio, trackresiduals] textmatrix = list(zip(keywords,keyvals)) t = Table(textmatrix,colWidths=[5*cm, 11*cm],rowHeights=5*mm) t.setStyle(TableStyle([('INNERGRID', (0,0), (-1,-1), 0.25, colors.black), ('BOX', (0,0), (-1,-1), 2, colors.black),])) Story.append(t) #============================================================================== # Page about reduction parameter #============================================================================== def plotQCSummary (Story, p2vmreduced, parameters): plotTitle(Story,"Summary of QC parameters") plotSubtitle(Story,"Parameters and values") # Get the name of the parameters print (parameters) names = p2vmreduced.header['*ESO QC DISP*'] keys = [] vals = [] for key, value in names.items(): keys.append(key) vals.append(value) hsize = [9*cm, 5.5*cm] vsize = 5*mm textmatrix = list(zip(keys,vals)) 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) #============================================================================== # Page about reduction parameter #============================================================================== def plotReductionSummary (Story, p2vmreduced): plotTitle(Story,"Summary of reduction parameters") plotSubtitle(Story,"Parameters and value for the recipe called") recnum = len(p2vmreduced.header['*REC* PIPE ID']) print(('Recipe ID = %i'%recnum)) recnum = 'REC%i'%recnum # Get the name of the recipe and all the options keywords = ['Pipeline version','Recipe'] keyvals = [p2vmreduced.header['HIERARCH ESO PRO '+recnum+' PIPE ID'], p2vmreduced.header['HIERARCH ESO PRO '+recnum+' ID']] names = p2vmreduced.header['*'+recnum+' PARAM* NAME'] values = p2vmreduced.header['*'+recnum+' PARAM* VALUE'] for i in range(len(names)): keywords.append(names[i]) keyvals.append(values[i]) hsize = [4.5*cm, 10*cm] vsize = 5*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) # Get all check warning Story.append(Spacer(1,5*mm)) plotSubtitle(Story,"Warning in QC CHECKS") keywords = ['Check'] keyvals = ['Message'] checks = p2vmreduced.header['*QC CHECK *'] for c in list(checks.items()): keywords.append(c[0]) keyvals.append(c[1]) hsize = [6*cm, 8.5*cm,] vsize = 5*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) # Get all files Story.append(Spacer(1,5*mm)) plotSubtitle(Story,"Files for the recipe") keywords = ['CATG'] keyvals = ['NAME'] names = p2vmreduced.header['*'+recnum+' RAW* CATG'] values = p2vmreduced.header['*'+recnum+' RAW* NAME'] for i in range(len(names)): keywords.append(names[i]) keyvals.append(values[i]) names = p2vmreduced.header['*'+recnum+' CAL* CATG'] values = p2vmreduced.header['*'+recnum+' CAL* NAME'] for i in range(len(names)): keywords.append(names[i]) keyvals.append(values[i]) hsize = [6*cm, 8.5*cm,] vsize = 5*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) #============================================================================== # Functions #============================================================================== def get_base_number(regnames,reg): base_order=["12","13","14","23","24","34"] base_inv=["21","31","41","32","42","43"] for b in np.arange(6): if base_order[b] == regnames[reg][0:2] or base_inv[b] == regnames[reg][0:2]: return b return -1 def get_reg_number(regnames,base,pola): # pola = 'S','P' or 'C' base_order=["12","13","14","23","24","34"] base_inv=["21","31","41","32","42","43"] noutputs = regnames.shape[0] regbase = [] regabcd = [] regpola = [] for output in range(0,noutputs): regbase.append(''.join(list(str(regnames[output]))[0:2])) # baseline number regabcd.append(list(str(regnames[output]))[3]) # ABCD output regpola.append(list(str(regnames[output]))[5]) # polarization reg = [] for output in range(0,noutputs): if ((regbase[output] == base_order[base]) or (regbase[output] == base_inv[base])): if (regabcd[output] == "A"): if (regpola[output] == pola): reg.append(output) if (regabcd[output] == "B"): if (regpola[output] == pola): reg.append(output) if (regabcd[output] == "C"): if (regpola[output] == pola): reg.append(output) if (regabcd[output] == "D"): if (regpola[output] == pola): reg.append(output) return reg def get_key_withdefault(dict, key, default): if key in list(dict.keys()): return dict[key] if key in ["HIERARCH "+c for c in list(dict.keys())]: return dict[key] else: return default def create_array_from_list(data,axis=None): out = np.array(data[0]) for d in data[1:]: out = np.append(out,d,axis=axis) return out def baseline_phases(telphases): M_matrix = np.array([-1.0,+1.0, 0.0, 0.0, # baseline 12 -1.0, 0.0,+1.0, 0.0, # baseline 13 -1.0, 0.0, 0.0,+1.0, # baseline 14 0.0,-1.0,+1.0, 0.0, # baseline 23 0.0,-1.0, 0.0,+1.0, # baseline 24 0.0, 0.0,-1.0,+1.0]) # baseline 34 M_matrix = M_matrix.reshape((6,4)) met = np.dot(M_matrix,telphases.T).T return met # mean angle of a numpy array of angles in radians or in degrees (if deg=True) def mean_angle(angles,axis=None,deg=False): if deg == True: angles = angles % 360 # Bring the angles between 0 and 360 complex_angles = np.cos(angles*np.pi/180.) + np.sin(angles*np.pi/180.)*1j meanangle = np.angle(np.nanmean(complex_angles,axis=axis))*180./np.pi else: angles = angles % (2*np.pi) # Bring the angles between 0 and 2pi complex_angles = np.cos(angles) + np.sin(angles)*1j meanangle = np.angle(np.nanmean(complex_angles,axis=axis)) return meanangle # standard deviation of a numpy array of angles in radians or in degrees (if deg=True) def std_angle(angles,axis=None,deg=False): if deg == True: angles = angles % 360 # Bring the angles between 0 and 360 #uw_angles = np.unwrap(angles*np.pi/180.,axis=axis) #stdangle = np.nanstd(uw_angles,axis=axis)*180./np.pi complex_angles = np.cos(angles*np.pi/180.) + np.sin(angles*np.pi/180.)*1j stdangle = np.nanstd(np.angle(complex_angles),axis=axis)*180./np.pi # return degrees else: angles = angles % (2*np.pi) # Bring the angles between 0 and 2pi #uw_angles = np.unwrap(angles,axis=axis) #stdangle = np.nanstd(uw_angles,axis=axis) complex_angles = np.cos(angles) + np.sin(angles)*1j stdangle = np.nanstd(np.angle(complex_angles),axis=axis) # return radians return stdangle def clean_gdelay_is(visdata, wave): nw = wave.size delta = (2*np.pi) * np.mean(1./wave[1:nw-1] - 1./wave[0:nw-2]) tmp = visdata gd = np.angle(np.mean(tmp[1:nw-1] * np.conj(tmp[0:nw-2]))) / delta tmp = tmp * np.exp(-2.j*np.pi * gd / wave) return tmp def clean_gdelay_fft(visdata, wave): nw = wave.size zp = 100 delta = 1./ (1./wave[nw-1] - 1./wave[0]) / zp tmp = visdata psd = np.abs(np.fft.fftn(tmp,s=[nw*zp])) mnx = np.argmax(psd) gd = mnx*delta if mnxnp.abs(threshold): for i in range(0,clean_step+1): # linear interpolation over the phase jump spring_phase[time+i] = spring_phase[time-1] + np.sign(bump[time])*2*np.pi*((1.0*i)/(1.0*clean_step)) time = time + clean_step + 1 # skip the bump region else: spring_phase[time]=spring_phase[time-1] time = time + 1 phase = phase - spring_phase return phase # Returns the weighted average of a list ignoring the nan and null values def nanaverage(dat,wei): #import pandas mdat = np.ma.masked_array(dat,np.isnan(dat)) #mdat2 = np.ma.masked_array(dat,pandas.isnull(mdat)) mwei = np.ma.masked_array(wei,np.isnan(dat)) #mwei2 = np.ma.masked_array(mwei,pandas.isnull(mdat)) #mm = np.average(mdat2,weights=mwei2) mm = np.average(mdat,weights=mwei) return mm # =============================================================================== # Python classes to read GRAVITY files at different stages of the data processing # =============================================================================== # # Oifits, Oilist, P2vm, Wave, Preproc, P2vmreduced, Rawdata, Badpix, Profile, Dark, Argon # class Oifits: def __init__(self, filename): # open fits file self.filename = filename gravi_oifits=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_oifits[0].header.copy() self.exptime_sc = self.header['HIERARCH ESO DET2 SEQ1 DIT'] # Mode of the instrument if 'SINGLE' in str(self.header['HIERARCH ESO PRO CATG']): self.single = True else: self.single = False if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True # get the FDDL, OI_ARRAY, OI_WAVELENGTH and OI_VIS_MET extentions for hdu in gravi_oifits : # OI_TARGET if hdu.name == 'OI_TARGET': # test if OI_TARGET is present in header self.oi_target = gravi_oifits['OI_TARGET'].data if hdu.name == 'OI_ARRAY' : oi_array=hdu.data.copy() if hdu.name == 'ARRAY_GEOMETRY' : # TBC, not present (2015-09) oi_array=hdu.data.copy() # get the OI_VIS, OI_FLUX and T3 extensions in combined polarization mode if self.polarsplit == False: for hdu in gravi_oifits : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_SC' in str(hdu.header['INSNAME'])): oi_vis_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('_FT' in str(hdu.header['INSNAME'])): oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') : oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_SC' in str(hdu.header['INSNAME'])): oi_vis2_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS2') & ('_FT' in str(hdu.header['INSNAME'])): oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') : oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('_SC' in str(hdu.header['INSNAME'])): oi_flux_sc=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc.field('TIME')[::4]/3600. # in decimal hours #headerflux = hdu.header.copy() columnflux = hdu.columns.names # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('_FT' in str(hdu.header['INSNAME'])): oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') : oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_T3') & ('_SC' in str(hdu.header['INSNAME'])): oi_t3_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('_FT' in str(hdu.header['INSNAME'])): oi_t3_ft=hdu.data.copy() if (hdu.name == 'OI_T3_FT') : oi_t3_ft=hdu.data.copy() # get the OI_VIS, OI_FLUX and T3 extensions in split polarization mode (S considered =P1, P considered =P2) if self.polarsplit == True: for hdu in gravi_oifits : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_flux_sc_p=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc_p.field('TIME')[::4]/3600. # in decimal hours #headerflux = hdu.header.copy() columnflux = hdu.columns.names if (hdu.name == 'OI_FLUX') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_flux_sc_s=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_vis_sc_p=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_vis_sc_s=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_vis2_sc_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_vis2_sc_s=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_t3_sc_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_t3_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() # Old format if (hdu.name == 'OI_VIS_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() gravi_oifits.close() del gravi_oifits # gc.collect() # Wavelength scales self.wave_sc=oi_wave_sc.field('EFF_WAVE')*1.e6 # wavelength scale in microns self.minwave_sc=np.min(self.wave_sc) self.maxwave_sc=np.max(self.wave_sc) self.wave_ft=oi_wave_ft.field('EFF_WAVE')*1.e6 self.minwave_ft=np.min(self.wave_sc) self.maxwave_ft=np.max(self.wave_sc) self.nwave_sc = self.wave_sc.shape[0] self.nwave_ft = self.wave_ft.shape[0] # OI_ARRAY parameters self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') self.stations = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION4','S4'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION3','S3'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION2','S2'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION1','S1')] self.telnames = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T4NAME','T4'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T3NAME','T3'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T2NAME','T2'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T1NAME','T1')] # baseline names self.basenames = [self.stations[0]+self.stations[1], self.stations[0]+self.stations[2], self.stations[0]+self.stations[3], self.stations[1]+self.stations[2], self.stations[1]+self.stations[3], self.stations[2]+self.stations[3]] # Time scale of the different records in the OIFITS file self.time = self.oi_flux_sc_time # Number of records self.nrecords = self.time.shape[0] # OI_FLUX table parameters if self.polarsplit == False: if 'VISFLUX' in columnflux: self.oi_flux_sc = oi_flux_sc.field('VISFLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('VISFLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc = oi_flux_sc.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft = oi_flux_ft.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_ft) else: self.oi_flux_sc = oi_flux_sc.field('FLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('FLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc = oi_flux_sc.field('FLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft = oi_flux_ft.field('FLUXERR').reshape(self.nrecords,4,self.nwave_ft) if self.polarsplit == True: if 'VISFLUX' in columnflux: self.oi_flux_sc_s = oi_flux_sc_s.field('VISFLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('VISFLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('VISFLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('VISFLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_s = oi_flux_sc_s.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_s = oi_flux_ft_s.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_p = oi_flux_sc_p.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_p = oi_flux_ft_p.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_ft) else: self.oi_flux_sc_s = oi_flux_sc_s.field('FLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('FLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('FLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('FLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_s = oi_flux_sc_s.field('FLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_s = oi_flux_ft_s.field('FLUXERR').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_p = oi_flux_sc_p.field('FLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_p = oi_flux_ft_p.field('FLUXERR').reshape(self.nrecords,4,self.nwave_ft) # OI_VIS and T3 table parameters if self.polarsplit == False: self.oi_vis_sc_time = oi_vis_sc.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_sc_mjd = oi_vis_sc.field('MJD')[::6] self.oi_vis_sc_target_id = oi_vis_sc.field('TARGET_ID')[::6] self.oi_vis_sc_int_time = oi_vis_sc.field('INT_TIME') # in seconds # Complex visibilities self.oi_vis_sc_visdata = oi_vis_sc.field('VISDATA') self.oi_vis_sc_viserr = oi_vis_sc.field('VISERR') self.oi_vis_sc_visamp = oi_vis_sc.field('VISAMP') self.oi_vis_sc_visamperr = oi_vis_sc.field('VISAMPERR') self.oi_vis_sc_visphi = oi_vis_sc.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr = oi_vis_sc.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_flag = oi_vis_sc.field('FLAG') self.oi_vis_ft_time = oi_vis_ft.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_ft_mjd = oi_vis_ft.field('MJD')[::6] self.oi_vis_ft_target_id = oi_vis_ft.field('TARGET_ID')[::6] self.oi_vis_ft_int_time = oi_vis_ft.field('INT_TIME') # in seconds self.oi_vis_ft_visdata = oi_vis_ft.field('VISDATA') self.oi_vis_ft_viserr = oi_vis_ft.field('VISERR') self.oi_vis_ft_visamp = oi_vis_ft.field('VISAMP') self.oi_vis_ft_visamperr = oi_vis_ft.field('VISAMPERR') self.oi_vis_ft_visphi = oi_vis_ft.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr = oi_vis_ft.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) # u,v coordinates different for SC and FT nbaseline= 6 self.ucoord_ft = oi_vis_ft.field('UCOORD') self.vcoord_ft = oi_vis_ft.field('VCOORD') self.projbase_ft = np.zeros(nbaseline) self.azimbase_ft = np.zeros(nbaseline) self.spfreq_ft = np.zeros((nbaseline,self.nwave_ft)) for baseline in range(0,nbaseline): self.projbase_ft[baseline] = np.sqrt(np.square(self.ucoord_ft[baseline])+np.square(self.vcoord_ft[baseline])) self.azimbase_ft[baseline] = np.arctan2(self.vcoord_ft[baseline], self.ucoord_ft[baseline]) * 180./np.pi self.spfreq_ft[baseline,:] = self.projbase_ft[baseline]/(self.wave_ft/1.E6)/(180.*3600./np.pi) # B/lambda cycles/arcsec self.ucoord_sc = oi_vis_sc.field('UCOORD') self.vcoord_sc = oi_vis_sc.field('VCOORD') self.projbase_sc = np.zeros(nbaseline) self.azimbase_sc = np.zeros(nbaseline) self.spfreq_sc = np.zeros((nbaseline,self.nwave_sc)) for baseline in range(0,nbaseline): self.projbase_sc[baseline] = np.sqrt(np.square(self.ucoord_sc[baseline])+np.square(self.vcoord_sc[baseline])) self.azimbase_sc[baseline] = np.arctan2(self.vcoord_sc[baseline], self.ucoord_sc[baseline]) * 180./np.pi self.spfreq_sc[baseline,:] = self.projbase_sc[baseline]/(self.wave_sc/1.E6)/(180.*3600./np.pi) # B/lambda cycles/arcsec # Station index common self.sta_index = oi_vis_sc.field('STA_INDEX') # Squared visibilities self.oi_vis2_sc_vis2data = oi_vis2_sc.field('VIS2DATA') self.oi_vis2_sc_vis2err = oi_vis2_sc.field('VIS2ERR') self.oi_vis2_sc_ucoord = oi_vis2_sc.field('UCOORD') self.oi_vis2_sc_vcoord = oi_vis2_sc.field('VCOORD') self.oi_vis2_sc_sta_index = oi_vis2_sc.field('STA_INDEX') self.oi_vis2_sc_flag = oi_vis2_sc.field('FLAG') self.oi_vis2_ft_vis2data = oi_vis2_ft.field('VIS2DATA') self.oi_vis2_ft_vis2err = oi_vis2_ft.field('VIS2ERR') self.oi_vis2_ft_ucoord = oi_vis2_ft.field('UCOORD') self.oi_vis2_ft_vcoord = oi_vis2_ft.field('VCOORD') self.oi_vis2_ft_sta_index = oi_vis2_ft.field('STA_INDEX') self.oi_vis2_ft_flag = oi_vis2_ft.field('FLAG') # Closure quantities self.t3_u1coord = oi_t3_sc.field('U1COORD') self.t3_v1coord = oi_t3_sc.field('V1COORD') self.t3_baseline1 = np.sqrt(np.square(self.t3_u1coord)+np.square(self.t3_v1coord)) self.t3_u2coord = oi_t3_sc.field('U2COORD') self.t3_v2coord = oi_t3_sc.field('V2COORD') self.t3_baseline2 = np.sqrt(np.square(self.t3_u2coord)+np.square(self.t3_v2coord)) self.t3_u3coord = self.t3_u1coord+self.t3_u2coord self.t3_v3coord = self.t3_v1coord+self.t3_v2coord self.t3_baseline3 = np.sqrt(np.square(self.t3_u3coord)+np.square(self.t3_v3coord)) self.t3_baseavg = np.nanmean([self.t3_baseline1,self.t3_baseline2,self.t3_baseline3],axis=0) self.t3_basemax = np.nanmax([self.t3_baseline1,self.t3_baseline2,self.t3_baseline3],axis=0) self.t3_spfreq_ft = np.zeros((4,self.nwave_ft)) self.t3_spfreq_sc = np.zeros((4,self.nwave_sc)) for triplet in range(0,4): self.t3_spfreq_ft[triplet,:] = self.t3_basemax[triplet]/(self.wave_ft/1.E6)/(180.*3600./np.pi) self.t3_spfreq_sc[triplet,:] = self.t3_basemax[triplet]/(self.wave_sc/1.E6)/(180.*3600./np.pi) self.t3phi_sc = oi_t3_sc.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_sc = oi_t3_sc.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_ft = oi_t3_ft.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_ft = oi_t3_ft.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3amp_sc = oi_t3_sc.field('T3AMP') self.t3amperr_sc = oi_t3_sc.field('T3AMPERR') self.t3amp_ft = oi_t3_ft.field('T3AMP') self.t3amperr_ft = oi_t3_ft.field('T3AMPERR') self.t3_flag = oi_t3_ft.field('FLAG') # Station indices for the closure quantities self.t3_sc_staindex = oi_t3_sc.field('STA_INDEX') self.t3_ft_staindex = oi_t3_ft.field('STA_INDEX') if self.polarsplit == True: self.oi_vis_sc_time = oi_vis_sc_s.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_sc_mjd = oi_vis_sc_s.field('MJD')[::6] self.oi_vis_sc_target_id = oi_vis_sc_s.field('TARGET_ID') # same for s and p polar self.oi_vis_sc_int_time = oi_vis_sc_s.field('INT_TIME') # in seconds # Complex visibilities self.oi_vis_sc_visdata_s = oi_vis_sc_s.field('VISDATA') self.oi_vis_sc_viserr_s = oi_vis_sc_s.field('VISERR') self.oi_vis_sc_visamp_s = oi_vis_sc_s.field('VISAMP') self.oi_vis_sc_visamperr_s = oi_vis_sc_s.field('VISAMPERR') self.oi_vis_sc_visphi_s = oi_vis_sc_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr_s = oi_vis_sc_s.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_flag_s = oi_vis_sc_s.field('FLAG') self.oi_vis_sc_visdata_p = oi_vis_sc_p.field('VISDATA') self.oi_vis_sc_viserr_p = oi_vis_sc_p.field('VISERR') self.oi_vis_sc_visamp_p = oi_vis_sc_p.field('VISAMP') self.oi_vis_sc_visamperr_p = oi_vis_sc_p.field('VISAMPERR') self.oi_vis_sc_visphi_p = oi_vis_sc_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr_p = oi_vis_sc_p.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_time = oi_vis_ft_s.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_ft_mjd = oi_vis_ft_s.field('MJD')[::6] self.oi_vis_ft_target_id = oi_vis_ft_s.field('TARGET_ID') # same for s and p polar self.oi_vis_ft_int_time = oi_vis_ft_s.field('INT_TIME') # in seconds self.oi_vis_ft_visdata_s = oi_vis_ft_s.field('VISDATA') self.oi_vis_ft_viserr_s = oi_vis_ft_s.field('VISERR') self.oi_vis_ft_visamp_s = oi_vis_ft_s.field('VISAMP') self.oi_vis_ft_visamperr_s = oi_vis_ft_s.field('VISAMPERR') self.oi_vis_ft_visphi_s = oi_vis_ft_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr_s = oi_vis_ft_s.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_sta_index_s = oi_vis_ft_s.field('STA_INDEX') self.oi_vis_ft_visdata_p = oi_vis_ft_p.field('VISDATA') self.oi_vis_ft_viserr_p = oi_vis_ft_p.field('VISERR') self.oi_vis_ft_visamp_p = oi_vis_ft_p.field('VISAMP') self.oi_vis_ft_visamperr_p = oi_vis_ft_p.field('VISAMPERR') self.oi_vis_ft_visphi_p = oi_vis_ft_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr_p = oi_vis_ft_p.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) # u,v coordinates common to both polarizations, different for SC and FT nbaseline= 6 self.ucoord_ft = oi_vis_ft_p.field('UCOORD') self.vcoord_ft = oi_vis_ft_p.field('VCOORD') self.projbase_ft = np.zeros(nbaseline) self.azimbase_ft = np.zeros(nbaseline) self.spfreq_ft = np.zeros((nbaseline,self.nwave_ft)) for baseline in range(0,nbaseline): self.projbase_ft[baseline] = np.sqrt(np.square(self.ucoord_ft[baseline])+np.square(self.vcoord_ft[baseline])) self.azimbase_ft[baseline] = np.arctan2(self.vcoord_ft[baseline], self.ucoord_ft[baseline]) * 180./np.pi self.spfreq_ft[baseline,:] = self.projbase_ft[baseline]/(self.wave_ft/1.E6)/(180.*3600./np.pi) # B/lambda cycles/arcsec self.ucoord_sc = oi_vis_sc_p.field('UCOORD') self.vcoord_sc = oi_vis_sc_p.field('VCOORD') self.projbase_sc = np.zeros(nbaseline) self.azimbase_sc = np.zeros(nbaseline) self.spfreq_sc = np.zeros((nbaseline,self.nwave_sc)) for baseline in range(0,nbaseline): self.projbase_sc[baseline] = np.sqrt(np.square(self.ucoord_sc[baseline])+np.square(self.vcoord_sc[baseline])) self.azimbase_sc[baseline] = np.arctan2(self.vcoord_sc[baseline], self.ucoord_sc[baseline]) * 180./np.pi self.spfreq_sc[baseline,:] = self.projbase_sc[baseline]/(self.wave_sc/1.E6)/(180.*3600./np.pi) # B/lambda cycles/arcsec # Station index common to both polarizations, SC and FT self.sta_index = oi_vis_sc_s.field('STA_INDEX') # Squared visibilities self.oi_vis2_sc_vis2data_s = oi_vis2_sc_s.field('VIS2DATA') self.oi_vis2_sc_vis2err_s = oi_vis2_sc_s.field('VIS2ERR') self.oi_vis2_sc_sta_index_s = oi_vis2_sc_s.field('STA_INDEX') self.oi_vis2_sc_flag_s = oi_vis2_sc_s.field('FLAG') self.oi_vis2_sc_vis2data_p = oi_vis2_sc_p.field('VIS2DATA') self.oi_vis2_sc_vis2err_p = oi_vis2_sc_p.field('VIS2ERR') self.oi_vis2_sc_sta_index_p = oi_vis2_sc_p.field('STA_INDEX') self.oi_vis2_sc_flag_p = oi_vis2_sc_p.field('FLAG') self.oi_vis2_ft_vis2data_s = oi_vis2_ft_s.field('VIS2DATA') self.oi_vis2_ft_vis2err_s = oi_vis2_ft_s.field('VIS2ERR') self.oi_vis2_ft_flag_s = oi_vis2_ft_s.field('FLAG') self.oi_vis2_ft_vis2data_p = oi_vis2_ft_p.field('VIS2DATA') self.oi_vis2_ft_vis2err_p = oi_vis2_ft_p.field('VIS2ERR') self.oi_vis2_ft_flag_p = oi_vis2_ft_p.field('FLAG') # Closure quantities self.t3_u1coord = oi_t3_sc_s.field('U1COORD') self.t3_v1coord = oi_t3_sc_s.field('V1COORD') self.t3_baseline1 = np.sqrt(np.square(self.t3_u1coord)+np.square(self.t3_v1coord)) self.t3_u2coord = oi_t3_sc_s.field('U2COORD') self.t3_v2coord = oi_t3_sc_s.field('V2COORD') self.t3_baseline2 = np.sqrt(np.square(self.t3_u2coord)+np.square(self.t3_v2coord)) self.t3_u3coord = self.t3_u1coord+self.t3_u2coord self.t3_v3coord = self.t3_v1coord+self.t3_v2coord self.t3_baseline3 = np.sqrt(np.square(self.t3_u3coord)+np.square(self.t3_v3coord)) self.t3_baseavg = np.nanmean([self.t3_baseline1,self.t3_baseline2,self.t3_baseline3],axis=0) self.t3_basemax = np.nanmax([self.t3_baseline1,self.t3_baseline2,self.t3_baseline3],axis=0) self.t3_spfreq_ft = np.zeros((4,self.nwave_ft)) self.t3_spfreq_sc = np.zeros((4,self.nwave_sc)) for triplet in range(0,4): self.t3_spfreq_ft[triplet,:] = self.t3_basemax[triplet]/(self.wave_ft/1.E6)/(180.*3600./np.pi) self.t3_spfreq_sc[triplet,:] = self.t3_basemax[triplet]/(self.wave_sc/1.E6)/(180.*3600./np.pi) self.t3phi_sc_s = oi_t3_sc_s.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_sc_s = oi_t3_sc_s.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_ft_s = oi_t3_ft_s.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_ft_s = oi_t3_ft_s.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_sc_p = oi_t3_sc_p.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_sc_p = oi_t3_sc_p.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_ft_p = oi_t3_ft_p.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_ft_p = oi_t3_ft_p.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3amp_sc_s = oi_t3_sc_s.field('T3AMP') self.t3amperr_sc_s = oi_t3_sc_s.field('T3AMPERR') self.t3amp_ft_s = oi_t3_ft_s.field('T3AMP') self.t3amperr_ft_s = oi_t3_ft_s.field('T3AMPERR') self.t3_flag_s = oi_t3_ft_s.field('FLAG') self.t3amp_sc_p = oi_t3_sc_p.field('T3AMP') self.t3amperr_sc_p = oi_t3_sc_p.field('T3AMPERR') self.t3amp_ft_p = oi_t3_ft_p.field('T3AMP') self.t3amperr_ft_p = oi_t3_ft_p.field('T3AMPERR') self.t3_flag_p = oi_t3_ft_p.field('FLAG') # Station indices for closure quantities self.t3_sc_staindex = oi_t3_sc_s.field('STA_INDEX') self.t3_ft_staindex = oi_t3_ft_s.field('STA_INDEX') class Oilist: def __init__(self, filelist, fits_keywords): self.oi = [] filelist_out = [] timescale_out = [] oi_out = [] nfiles_in = len(filelist) nfiles = 0 for filename in filelist: gravi_file = pyfits.open(filename+'.fits',memmap=True) header = gravi_file[0].header.copy() gravi_file.close() del gravi_file # Final OIFITS data files only key_names = list(fits_keywords.keys()) type_ok = True for strname in key_names: type_ok *= (strname in header) if type_ok: keys_ok = True # for name in key_names: # keys_ok *= (header[name] in fits_keywords[name]) # print keys_ok if keys_ok == True: nfiles = nfiles + 1 filelist_out.append(filename) print(' Loading OIFITS file %i/%i: %s.fits'%(nfiles,nfiles_in,filename)) a = Oifits(filename+'.fits') oi_out.append(a) timescale_out.append(float(a.header['MJD-OBS'])) # timescale of the measurements # TO BE CHANGED: (SEVERAL CAN BE PRESENT IN EACH FILE) del a del header # gc.collect() # Sort the output structure as a function of MJD of observation sorted_timescale = np.array(np.argsort(timescale_out)) # indices of sorted timescale self.filelist = [] self.timescale = [] for i in range(0,sorted_timescale.shape[0]): self.oi.append(oi_out[sorted_timescale[i]]) self.filelist.append(filelist_out[sorted_timescale[i]]) self.timescale.append(timescale_out[sorted_timescale[i]]) self.nfiles = nfiles class P2vm: def __init__(self, filename): # open fits file self.filename = filename gravi_p2vm=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_p2vm[0].header #print self.header['HIERARCH ESO PRO REC1 PIPE ID'] if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True if "COMBINED" in str(self.header['HIERARCH ESO INS POLA MODE']): self.polarsplitSC = False else: self.polarsplitSC = True #print self.polarsplit self.ports_sc = np.array(gravi_p2vm['IMAGING_DETECTOR_SC'].data['PORTS'])-1 # telescopes from 0 to 3 self.ports_ft = np.array(gravi_p2vm['IMAGING_DETECTOR_FT'].data['PORTS'])-1 self.regname_sc=gravi_p2vm['P2VM_SC'].data['REGNAME'] self.regname_ft=gravi_p2vm['P2VM_FT'].data['REGNAME'] self.transmission_sc=gravi_p2vm['P2VM_SC'].data['TRANSMISSION'] self.transmission_ft=gravi_p2vm['P2VM_FT'].data['TRANSMISSION'] self.transmission_met=gravi_p2vm['P2VM_MET'].data['TRANSMISSION'] self.nwave_sc=self.transmission_sc.shape[2] self.nwave_ft=self.transmission_ft.shape[2] self.nregion_sc=self.transmission_sc.shape[0] self.nregion_ft=self.transmission_ft.shape[0] self.coherence_sc=gravi_p2vm['P2VM_SC'].data['COHERENCE'] self.coherence_ft=gravi_p2vm['P2VM_FT'].data['COHERENCE'] self.coherence_met=gravi_p2vm['P2VM_MET'].data['COHERENCE'] self.phase_sc=gravi_p2vm['P2VM_SC'].data['PHASE'] self.phase_ft=gravi_p2vm['P2VM_FT'].data['PHASE'] self.phase_met=gravi_p2vm['P2VM_MET'].data['PHASE'] self.cmatrix_sc=np.array(gravi_p2vm['P2VM_SC'].data['C MATRIX']).reshape(self.nregion_sc,6,self.nwave_sc) self.cmatrix_ft=np.array(gravi_p2vm['P2VM_FT'].data['C MATRIX']).reshape(self.nregion_ft,6,self.nwave_ft) # get the OI_WAVELENGTH extention # these tables are present twice in the P2VM files (one for each polarization) but contain the same data if (self.polarsplit == False): for hdu in gravi_p2vm : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC' in str(hdu.header['INSNAME'])): self.wave_sc=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_SC') : self.wave_sc=hdu.data['EFF_WAVE']*1.e6 # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT' in str(hdu.header['INSNAME'])): self.wave_ft=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_FT') : self.wave_ft=hdu.data['EFF_WAVE']*1.e6 if (self.polarsplit == True): for hdu in gravi_p2vm : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC_P1' in str(hdu.header['INSNAME'])): self.wave_sc=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_SC') : self.wave_sc=hdu.data['EFF_WAVE']*1.e6 # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT_P1' in str(hdu.header['INSNAME'])): self.wave_ft=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_FT') : self.wave_ft=hdu.data['EFF_WAVE']*1.e6 self.minwave_sc = np.min(self.wave_sc) self.maxwave_sc = np.max(self.wave_sc) self.minwave_ft = np.min(self.wave_ft) self.maxwave_ft = np.max(self.wave_ft) gravi_p2vm.close() class P2vmlist: def __init__(self, filelist, fits_keywords): self.p2vm = [] filelist_out = [] timescale_out = [] p2vm_out = [] nfiles_in = len(filelist) nfiles = 0 for filename in filelist: gravi_file = pyfits.open(filename+'.fits',memmap=True) header = gravi_file[0].header.copy() gravi_file.close() del gravi_file # P2VM data files only with proper keywords key_names = list(fits_keywords.keys()) type_ok = True for strname in key_names: type_ok *= (strname in header) if type_ok: keys_ok = True for name in key_names: keys_ok *= (header[name] in fits_keywords[name]) if keys_ok == True: nfiles = nfiles + 1 filelist_out.append(filename) print(' Loading P2VM file %i/%i: %s.fits'%(nfiles,nfiles_in,filename)) a = P2vm(filename+'.fits') p2vm_out.append(a) timescale_out.append(float(a.header['MJD-OBS'])) # timescale of the measurements # TO BE CHANGED: (SEVERAL CAN BE PRESENT IN EACH FILE) del a del header # gc.collect() # Sort the output structure as a function of MJD of observation sorted_timescale = np.array(np.argsort(timescale_out)) # indices of sorted timescale self.filelist = [] self.timescale = [] for i in range(0,sorted_timescale.shape[0]): self.p2vm.append(p2vm_out[sorted_timescale[i]]) self.filelist.append(filelist_out[sorted_timescale[i]]) self.timescale.append(timescale_out[sorted_timescale[i]]) self.nfiles = nfiles class Wave: def __init__(self, filename): # open fits file self.filename = filename gravi_wave=pyfits.open(filename, mode='readonly') # get the WAVE_DATA_SC extention self.valid = 0 for hdu in gravi_wave : #print hdu.name if hdu.name=='WAVE_DATA_SC' : self.valid=1 if self.valid == 0 : print("No WAVE_DATA_SC table in this file") return None self.valid = 0 for hdu in gravi_wave : #print hdu.name if hdu.name=='WAVE_DATA_FT' : self.valid=1 if self.valid == 0 : print("No WAVE_DATA_FT table in this file") return None # Main header self.header = gravi_wave[0].header if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True self.wavefiber_in = False for hdu in gravi_wave : if hdu.name=='WAVE_DATA_SC' : h = hdu.header self.startx = h['STARTX'] if 'STARTX' in h else h['HIERARCH ESO PRO PROFILE STARTX'] self.nwave_sc= hdu.data.field('DATA1').size self.nregion_sc= h['TFIELDS'] if hdu.name=='WAVE_DATA_FT' : h = hdu.header self.nwave_ft= hdu.data.field('DATA1').size self.nregion_ft= h['TFIELDS'] if hdu.name=='WAVE_FIBRE_SC' : h = hdu.header self.wavefiber_in = True if self.polarsplit == True: self.nwave_sc= hdu.data.field('BASE_12_S').size wfs = np.zeros((6,self.nwave_sc)) wfp = np.zeros((6,self.nwave_sc)) wfs[0,:]=hdu.data['BASE_12_S']*1E6 wfs[1,:]=hdu.data['BASE_13_S']*1E6 wfs[2,:]=hdu.data['BASE_14_S']*1E6 wfs[3,:]=hdu.data['BASE_23_S']*1E6 wfs[4,:]=hdu.data['BASE_24_S']*1E6 wfs[5,:]=hdu.data['BASE_34_S']*1E6 wfp[0,:]=hdu.data['BASE_12_P']*1E6 wfp[1,:]=hdu.data['BASE_13_P']*1E6 wfp[2,:]=hdu.data['BASE_14_P']*1E6 wfp[3,:]=hdu.data['BASE_23_P']*1E6 wfp[4,:]=hdu.data['BASE_24_P']*1E6 wfp[5,:]=hdu.data['BASE_34_P']*1E6 self.wavefiber_s = wfs self.wavefiber_p = wfp if self.polarsplit == False: self.nwave_sc= hdu.data.field('BASE_12_C').size wfc = np.zeros((6,self.nwave_sc)) wfc[0,:]=hdu.data['BASE_12_C']*1E6 wfc[1,:]=hdu.data['BASE_13_C']*1E6 wfc[2,:]=hdu.data['BASE_14_C']*1E6 wfc[3,:]=hdu.data['BASE_23_C']*1E6 wfc[4,:]=hdu.data['BASE_24_C']*1E6 wfc[5,:]=hdu.data['BASE_34_C']*1E6 self.wavefiber_c = wfc #self.wave_sc = wave_sc.reshape(self.nregion_sc,self.nwave_sc) #self.wave_ft = wave_ft.reshape(self.nregion_ft,self.nwave_ft) # Reshaping of wavelength vectors (output in microns) wavetable_sc = np.zeros([self.nregion_sc,self.nwave_sc],dtype='d') wave_sc = gravi_wave['WAVE_DATA_SC'].data for output in range(1,self.nregion_sc+1): wavetable_sc[output-1,:] = wave_sc.field('DATA'+str(output))[0].reshape(1,self.nwave_sc) self.wave_sc = wavetable_sc*1E6 wavetable_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') wave_ft = gravi_wave['WAVE_DATA_FT'].data for output in range(1,self.nregion_ft+1): wavetable_ft[output-1,:] = wave_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft) self.wave_ft = wavetable_ft*1E6 self.minwave_ft = np.min(self.wave_ft) self.maxwave_ft = np.max(self.wave_ft) imaging_det_sc = gravi_wave['IMAGING_DETECTOR_SC'].data self.regname_sc = imaging_det_sc.field('REGNAME') self.center_sc = (imaging_det_sc.field('CENTER')) imaging_det_ft = gravi_wave['IMAGING_DETECTOR_FT'].data self.regname_ft = imaging_det_ft.field('REGNAME') self.center_ft = (imaging_det_ft.field('CENTER')) if self.polarsplit == True: self.base_list_sc = [self.regname_sc[i*8][0:2] for i in range(0,6)] self.base_list_ft = [self.regname_ft[i*8][0:2] for i in range(0,6)] if self.polarsplit == False: self.base_list_sc = [self.regname_sc[i*4][0:2] for i in range(0,6)] self.base_list_ft = [self.regname_ft[i*4][0:2] for i in range(0,6)] # get the TEST_WAVE engineering image cube extention: map of the SC spectral extraction + map of wavelengths # 0: interpolated wavelength map # 1: profile map # 2: non-interpolated wavelength map (not available for older processed data) self.test_wave = gravi_wave['TEST_WAVE'].data self.test_wave[0,:,:] *= 1E6 # to get microns in the wavelength tables self.test_wave[2,:,:] *= 1E6 # to get microns in the wavelength tables ind=np.where(self.test_wave[0,:,:] > 1) self.minwave_sc = np.min(self.test_wave[0,ind[0],ind[1]]) self.maxwave_sc = np.max(self.test_wave[0,ind[0],ind[1]]) #self.minwave_sc = np.min(np.select(self.test_wave[0,:,:] > 1.0,self.test_wave[0,:,:])) #waself.maxwave_sc = np.max(np.select(self.test_wave[0,:,:] > 1.0,self.test_wave[0,:,:])) gravi_wave.close() class Preproc: def __init__(self, filename): # open fits file self.filename = filename gravi_preproc=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_preproc[0].header # Single mode of the instrument (no metrology data) #if 'SINGLE' in str(self.header['HIERARCH ESO DO CATG']): if 'SINGLE' in str(self.header['HIERARCH ESO PRO REC1 RAW1 CATG']): self.single = True else: self.single = False if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False self.nregion = 24 else: self.polarsplit = True self.nregion = 48 fddl = None # in case the FDDL data are not in the file detector_sc = None detector_ft = None oi_array = None spectrum_sc = None spectrum_ft = None metrology = None self.fddlpresent = False self.gdelay_in = False # get the FDDL, OI_ARRAY, OI_WAVELENGTH, OI_VIS_MET and OPDC extentions for hdu in gravi_preproc : if hdu.name == 'ARRAY_GEOMETRY' : # TBC oi_array=hdu.data if hdu.name == 'FDDL': fddl = hdu.data self.fddlpresent = True if hdu.name == 'IMAGING_DETECTOR_SC' : detector_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_FT' : detector_ft=hdu.data if hdu.name == 'METROLOGY' : metrology=hdu.data if hdu.name == 'SPECTRUM_DATA_FT' : spectrum_ft=hdu.data if hdu.name == 'SPECTRUM_DATA_SC' : spectrum_sc=hdu.data if hdu.name == 'OPDC' : opdc=hdu.data self.gdelay_in = True if self.polarsplit == False: for hdu in gravi_preproc : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if self.polarsplit == True: for hdu in gravi_preproc : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() # Wavelength scales self.wave_sc=oi_wave_sc.field('EFF_WAVE')*1.e6 # wavelength scale in microns self.minwave_sc=np.min(self.wave_sc) self.maxwave_sc=np.max(self.wave_sc) self.wave_ft=oi_wave_ft.field('EFF_WAVE')*1.e6 self.minwave_ft=np.min(self.wave_sc) self.maxwave_ft=np.max(self.wave_sc) self.nwave_sc = self.wave_sc.shape[0] self.nwave_ft = self.wave_ft.shape[0] # OPDC parameters if self.gdelay_in == True: self.opdc_time = opdc.field('TIME')/1e6 # in seconds self.opdc_piezo_offset = opdc.field('PIEZO_DL_OFFSET') # in VOLTS self.opdc_vltidl_offset = opdc.field('VLTI_DL_OFFSET') self.opdc_state = opdc.field('STATE') # IMAGING_DETECTOR parameters if detector_sc is not None: self.regname_sc = detector_sc.field('REGNAME') self.ports_sc = detector_sc.field('PORTS') else: self.regname_sc = '' if detector_ft is not None: self.regname_ft = detector_ft.field('REGNAME') self.ports_ft = detector_ft.field('PORTS') else: self.regname_ft = '' # OI_ARRAY parameters if oi_array is not None: self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # Output spectra if spectrum_sc is not None: self.spectrum_sc = [] # initialization self.time_sc = spectrum_sc.field('TIME')/1E6 # in seconds self.spectrum_sc = np.zeros((self.nregion, self.time_sc.shape[0], self.nwave_sc)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_sc[i-1,:,:] = spectrum_sc.field('DATA'+str(i)) if spectrum_ft is not None: self.spectrum_ft = [] # initialization self.time_ft = spectrum_ft.field('TIME')/1E6 # in seconds self.spectrum_ft = np.zeros((self.nregion , self.time_ft.shape[0], self.nwave_ft)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_ft[i-1,:,:] = spectrum_ft.field('DATA'+str(i)) # METROLOGY parameters self.time_met = None self.metrology_volt = None if metrology is not None: self.time_met = metrology.field('TIME')/1e6 # in seconds self.metrology_volt = metrology.field('VOLT') self.nframe_met = self.time_met.shape[0] # Number of frames if 'HIERARCH ESO DET2 NDIT' in self.header: self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] else: self.nframe_sc = 0 self.nframe_ft = self.spectrum_ft.shape[1] gravi_preproc.close() class P2vmreduced: #------------------------------- #Columns in the OI_VIS of the SC # #TARGET_ID : id listed in OI_TARGET #TIME [usec] : timestamp #MJD [day] : #INT_TIME [s] : integration time of this frame #VISDATA [e,e] : complex coherent flux spectra of SC in this frame #VISERR [e,e] : complex error on the coherent flux spectra #UCOORD [m] : uv-plane of this SC frame #VCOORD [m] : uv-plane of this SC frame #STA_INDEX : station index in the OI_ARRAY # #PHASE_MET_FC [rad] : unwrap FT-SC phase as computed by the DRS algorithm #PHASE_MET_TEL [rad] : unwrap FT-SC phase as computed by the DRS algorithm, mean of 4 diodes #OPD_MET_FC [m] : unwrap SC-FT delay as computed by the TAC algorihm #OPD_MET_TEL [m] : unwrap SC-FT delay as computed by the TAC algorihm, 4 diodes # #VISDATA_FT [e,e] : VISDATA spectra of FT (integrated in this SC frame) #PHASE_REF [rad] : reference phase, actually -1*arg{VISDATA_FT}, re-interpolated in the SC wavelength. # #NFRAME_FT : number of averaged FT frame in this SC frame #NFRAME_MET : number of averaged MET frame in this SC frame # #GDELAY [m] : real-time GD #SNR : real-time SNR #GDELAY_BOOT [m] : best GD estimate, accounting closing triangles #SNR_BOOT : best SNR estimate, accounting closing triangles # #V_FACTOR_FT : measured visibility loss on the FT #V_FACTOR : predicted visibility loss of this SC frame (re-interpolation of V_FACTOR_FT on the SC wavelengths) # #FRINGEDET_RATIO : fraction of FT frame accepted in this SC frame #REJECTION_FLAG : this frame is accepted/rejected # #-------------------------------- #Columns in the OI_FLUX of the SC # #TARGET_ID : id listed in OI_TARGET #TIME [usec] : timestamp #MJD [day] : #INT_TIME [s] : integration time of this frame #FLUX [e] : flux spectra #FLUXERR [e] : error on flux spectra #STA_INDEX : station index in the OI_ARRAY #TOTALFLUX_SC [e] : total flux of SC in this SC frame (integrated over spectrum) #TOTALFLUX_FT [e] : total flux of FT in this SC frame (integrated over spectrum) # #------------------------- #Columns in the OI_VIS_MET # #PHASE_TEL [rad] : 4 diodes x 4 beams phases at telescope, unwrap by pipeline algorithm (FT-SC) #PHASE_FC [rad] : 4 beams phases at combiner, unwrap by pipeline algorithm (FT-SC) #OPD_TEL [m] : 4 diodes x 4 beams OPD at telescope, unwrap by TAC algorithm (SC-FT) #OPD_FC [m] : 4 beams OPDs at telescope, unwrap by TAC algorithm (SC-FT) #FLAG_FC, FLAG_TEL: flags computed by TAC algorithm #VAMP_FC_FT, VAMP_FC_SC, VAMP_TEL_FT, VAMP_TEL_SC: Volt amplitude # #------------------------- #Columns in the OPDC table # #TIME [usec] : timestamp #PIEZO_DL_OFFSET #VLTI_DL_OFFSET #STATE : global fringe tracking state #STEPS : target phase modulation per baseline (scrambled), in units of pi/8 # #---------------------------------------- #HEADER information about target position # #FT.ROBJ.ALPHA #FT.ROBJ.DELTA #FT.ROBJ.EPOCH #FT.ROBJ.EQUINOX #FT.ROBJ.PARALLAX #FT.ROBJ.PMA #FT.ROBJ.PMD #FT.ROBJ.MAG #FT.ROBJ.DIAMETER #FT.ROBJ.NAME # #INS.SOBJ.X #INS.SOBJ.Y #INS.SOBJ.DIAMETER #INS.SOBJ.MAG #INS.SOBJ.NAME # #And during a SWAP: # #FT.ROBJ.ALPHA/DELTA += INS.SOBJ.X/Y #INS.SOBJ.XY *= -1 #INS.SOBJ.DIAMETER <=> INS.ROBJ.DIAMETER #INS.SOBJ.MAG <=> INS.ROBJ.MAG #INS.SOBJ.NAME <=> FT.ROBJ.NAME # #------------------- #Changes 2016-02-01: # #OI_VIS_MET #rename: PHI -> PHASE_FC [rad], 4 beams #rename: PHI_TEL -> PHASE_TEL, 4 diodes x 4 beams #change: OPD_TEL and OPD_FC to add the zero # #OI_VIS of the SC #rename: PHI_MET -> PHASE_MET_FC [rad], scalar #rename: PHI_MET_TEL -> PHASE_MET_TEL [rad], scalar, mean of four diodes #add: OPD_MET_FC [m] scalar #add: OPD_MET_TEL [m] four diodes #rename: EXP_PHASE -> PHASE_REF #change: GDELAY and GDELAY_BOOT -> now in [m] # def __init__(self, filename): # open fits file self.filename = filename gravi_p2vmreduced=pyfits.open(filename, mode='readonly') nbase = 6 # Main FITS header self.header = gravi_p2vmreduced[0].header # Single mode of the instrument (no OI_VIS_MET data) #if 'SINGLE' in str(self.header['HIERARCH ESO DO CATG']): if 'SINGLE' in str(self.header['HIERARCH ESO PRO REC1 RAW1 CATG']): self.single = True else: self.single = False if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True fddl = None # in case the FDDL data are not in the file oi_vis_met = None oi_flux_sc = None oi_flux_sc_s = None oi_vis_sc = None oi_vis_sc_s = None self.gdelay_in=False self.vfactor_in=False # get the FDDL, OI_ARRAY extentions for hdu in gravi_p2vmreduced : if hdu.name == 'OI_ARRAY' : oi_array=hdu.data if hdu.name == 'ARRAY_GEOMETRY' : # TBC oi_array=hdu.data if hdu.name == 'FDDL': fddl = hdu.data if hdu.name == 'OI_VIS_MET' : oi_vis_met=hdu.data if hdu.name == 'OPDC' : opdc=hdu.data self.gdelay_in = True # ACQ CAM if 'OI_VIS_ACQ' in gravi_p2vmreduced: print ('OI_VIS_ACQ is in') self.acq_in = True data = gravi_p2vmreduced['OI_VIS_ACQ'].data self.oi_acq_time = data.field('TIME').copy() / 1e6 self.oi_acq_pup_x = data.field('PUPIL_X').copy() self.oi_acq_pup_y = data.field('PUPIL_Y').copy() self.oi_acq_pup_z = data.field('PUPIL_Z').copy() self.oi_acq_pup_r = data.field('PUPIL_R').copy() self.oi_acq_pup_n = data.field('PUPIL_NSPOT').copy() else: self.acq_in = False # IP numbers self.ipnames = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF INPUT1',1), get_key_withdefault(self.header,'HIERARCH ESO ISS CONF INPUT2',3), get_key_withdefault(self.header,'HIERARCH ESO ISS CONF INPUT3',5), get_key_withdefault(self.header,'HIERARCH ESO ISS CONF INPUT4',7)] # Match of GRAVITY inputs to IPs: # IP7-GRAV1, IP5-GRAV2, IP3-GRAV3, IP1-GRAV1 self.stations = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION4','S4'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION3','S3'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION2','S2'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION1','S1')] self.telnames = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T4NAME','T4'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T3NAME','T3'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T2NAME','T2'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T1NAME','T1')] # baseline names self.basenames = [self.stations[0]+self.stations[1], self.stations[0]+self.stations[2], self.stations[0]+self.stations[3], self.stations[1]+self.stations[2], self.stations[1]+self.stations[3], self.stations[2]+self.stations[3]] # get the OI_VIS, OI_FLUX and T3 extensions in combined polarization mode if self.polarsplit == False: for hdu in gravi_p2vmreduced : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_SC' in str(hdu.header['INSNAME'])): oi_vis_sc=hdu.data.copy() #header_oi_vis_sc = hdu.header.copy() column_oi_vis_sc = hdu.columns.names # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('_FT' in str(hdu.header['INSNAME'])): oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') : oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_SC' in str(hdu.header['INSNAME'])): oi_vis2_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS2') & ('_FT' in str(hdu.header['INSNAME'])): oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') : oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('_SC' in str(hdu.header['INSNAME'])): oi_flux_sc=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc.field('TIME')[::4]/3600. # in decimal hours #headerflux = hdu.header.copy() columnflux = hdu.columns.names # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('_FT' in str(hdu.header['INSNAME'])): oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') : oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_T3') & ('_SC' in str(hdu.header['INSNAME'])): oi_t3_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('_FT' in str(hdu.header['INSNAME'])): oi_t3_ft=hdu.data.copy() if (hdu.name == 'OI_T3_FT') : oi_t3_ft=hdu.data.copy() # get the OI_VIS, OI_FLUX and T3 extensions in split polarization mode (S considered =P1, P considered =P2) if self.polarsplit == True: for hdu in gravi_p2vmreduced : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_vis_sc_p=hdu.data.copy() header_oi_vis_sc_p = hdu.header.copy() column_oi_vis_sc_p = hdu.columns.names if (hdu.name == 'OI_VIS') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_vis_sc_s=hdu.data.copy() header_oi_vis_sc_s = hdu.header.copy() column_oi_vis_sc_s = hdu.columns.names if (hdu.name == 'OI_VIS2') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_vis2_sc_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_vis2_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_flux_sc_p=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc_p.field('TIME')[::4]/3600. # in decimal hours headerflux = hdu.header.copy() columnflux = hdu.columns.names if (hdu.name == 'OI_FLUX') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_flux_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_t3_sc_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('_SC_P2' in str(hdu.header['INSNAME'])): oi_t3_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() # Wavelength scales self.wave_sc = oi_wave_sc.field('EFF_WAVE')*1.e6 # wavelength scale in microns self.minwave_sc=np.min(self.wave_sc) self.maxwave_sc=np.max(self.wave_sc) self.nwave_sc = self.wave_sc.shape[0] self.wave_ft = oi_wave_ft.field('EFF_WAVE')*1.e6 self.minwave_ft=np.min(self.wave_sc) self.maxwave_ft=np.max(self.wave_sc) self.nwave_ft = self.wave_ft.shape[0] # OI_ARRAY parameters self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # OPDC parameters if self.gdelay_in == True: # print("Read OPDC data") self.opdc_time = opdc.field('TIME')/1e6 # in seconds self.opdc_piezo_offset = opdc.field('PIEZO_DL_OFFSET') # in VOLTS self.opdc_vltidl_offset = opdc.field('VLTI_DL_OFFSET') self.opdc_state = opdc.field('STATE') #nrow = len(self.opdc_time) try: self.opdc_kalman_piezo = opdc.field('KALMAN_PIEZO') # Time, telescope self.opdc_opd = opdc.field('OPD') # Time, baseline self.opdc_kalman_opd = opdc.field('KALMAN_OPD') # Time, baseline except KeyError: print("No Kalman OPD table") # OI_VIS_MET parameters if oi_vis_met is not None: if len(oi_vis_met.field('PHASE_FC_TAC').shape) == 1: print("Load the MET and TAC data in new format") self.oi_vis_met_time = oi_vis_met.field('TIME')[::4]/1e6 # in seconds nrow = len(self.oi_vis_met_time) try: self.oi_vis_met_phase = oi_vis_met.field('PHASE_FC_TAC').reshape((nrow,4)) self.oi_vis_met_phase_tel = oi_vis_met.field('PHASE_TEL_TAC').reshape((nrow,16)) except: self.oi_vis_met_phase = oi_vis_met.field('PHASE_FC').reshape((nrow,4)) self.oi_vis_met_phase_tel = oi_vis_met.field('PHASE_TEL').reshape((nrow,16)) self.met_opd_tel = oi_vis_met.field('OPD_TEL').reshape((nrow,16)) self.met_opd_fc = oi_vis_met.field('OPD_FC').reshape((nrow,4)) self.met_flag_tel = oi_vis_met.field('FLAG_TEL').reshape((nrow,16)) self.met_flag_fc = oi_vis_met.field('FLAG_FC').reshape((nrow,4)) self.met_vamp_tel_sc = oi_vis_met.field('VAMP_TEL_SC').reshape((nrow,16)) self.met_vamp_tel_ft = oi_vis_met.field('VAMP_TEL_FT').reshape((nrow,16)) self.met_vamp_fc_sc = oi_vis_met.field('VAMP_FC_SC').reshape((nrow,4)) self.met_vamp_fc_ft = oi_vis_met.field('VAMP_FC_FT').reshape((nrow,4)) else: print("Load the MET and TAC data in old format") self.oi_vis_met_time = oi_vis_met.field('TIME')/1e6 # in seconds try: if 'PHASE_FC_TAC' in oi_vis_met.columns.names: self.oi_vis_met_phase = oi_vis_met.field('PHASE_FC_TAC') if 'PHASE_TEL_TAC' in oi_vis_met.columns.names: self.oi_vis_met_phase_tel = oi_vis_met.field('PHASE_TEL_TAC') except: if 'PHASE_FC' in oi_vis_met.columns.names: self.oi_vis_met_phase = oi_vis_met.field('PHASE_FC') if 'PHASE_TEL' in oi_vis_met.columns.names: self.oi_vis_met_phase_tel = oi_vis_met.field('PHASE_TEL') if 'OPD_TEL' in oi_vis_met.columns.names: self.met_opd_tel = oi_vis_met.field('OPD_TEL') if 'OPD_FC' in oi_vis_met.columns.names: self.met_opd_fc = oi_vis_met.field('OPD_FC') if 'FLAG_TEL' in oi_vis_met.columns.names: self.met_flag_tel = oi_vis_met.field('FLAG_TEL') if 'FLAG_FC' in oi_vis_met.columns.names: self.met_flag_fc = oi_vis_met.field('FLAG_FC') if 'VAMP_TEL_SC' in oi_vis_met.columns.names: self.met_vamp_tel_sc = oi_vis_met.field('VAMP_TEL_SC') if 'VAMP_TEL_FT' in oi_vis_met.columns.names: self.met_vamp_tel_ft = oi_vis_met.field('VAMP_TEL_FT') if 'VAMP_FC_SC' in oi_vis_met.columns.names: self.met_vamp_fc_sc = oi_vis_met.field('VAMP_FC_SC') if 'VAMP_FC_FT' in oi_vis_met.columns.names: self.met_vamp_fc_ft = oi_vis_met.field('VAMP_FC_FT') # DELAY to account for the missing first frame of the SC # SCdelay = self.header['HIERARCH ESO INS TIM1 PERIOD'] > now included in DRS # OI_FLUX table parameters if (oi_flux_sc is not None) & (self.polarsplit == False): print("Load FLUX data") self.oi_flux_sc_time = oi_flux_sc.field('TIME')[::4]/1e6 self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/1e6 if 'VISFLUX' in columnflux: self.oi_flux_sc = oi_flux_sc.field('VISFLUX').reshape(self.oi_flux_sc_time.shape[0],4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('VISFLUX').reshape(self.oi_flux_ft_time.shape[0],4,self.nwave_ft) else: self.oi_flux_sc = oi_flux_sc.field('FLUX').reshape(self.oi_flux_sc_time.shape[0],4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('FLUX').reshape(self.oi_flux_ft_time.shape[0],4,self.nwave_ft) if 'TOTALFLUX_SC' in columnflux: self.totalflux_sc = oi_flux_sc.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time.shape[0],4) self.totalflux_ft = oi_flux_sc.field('TOTALFLUX_FT').reshape(self.oi_flux_sc_time.shape[0],4) if (oi_flux_sc_s is not None) & (self.polarsplit == True): self.oi_flux_sc_time_s = oi_flux_sc_s.field('TIME')[::4]/1e6 self.oi_flux_ft_time_s = oi_flux_ft_s.field('TIME')[::4]/1e6 self.oi_flux_sc_time_p = oi_flux_sc_p.field('TIME')[::4]/1e6 self.oi_flux_ft_time_p = oi_flux_ft_p.field('TIME')[::4]/1e6 if 'VISFLUX' in columnflux: self.oi_flux_sc_s = oi_flux_sc_s.field('VISFLUX').reshape(self.oi_flux_sc_time_s.shape[0],4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('VISFLUX').reshape(self.oi_flux_ft_time_s.shape[0],4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('VISFLUX').reshape(self.oi_flux_sc_time_p.shape[0],4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('VISFLUX').reshape(self.oi_flux_ft_time_p.shape[0],4,self.nwave_ft) else: self.oi_flux_sc_s = oi_flux_sc_s.field('FLUX').reshape(self.oi_flux_sc_time_s.shape[0],4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('FLUX').reshape(self.oi_flux_ft_time_s.shape[0],4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('FLUX').reshape(self.oi_flux_sc_time_p.shape[0],4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('FLUX').reshape(self.oi_flux_ft_time_p.shape[0],4,self.nwave_ft) if 'TOTALFLUX_SC' in columnflux: self.totalflux_sc_s = oi_flux_sc_s.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time_s.shape[0],4) self.totalflux_ft_s = oi_flux_sc_s.field('TOTALFLUX_FT').reshape(self.oi_flux_sc_time_s.shape[0],4) self.totalflux_sc_p = oi_flux_sc_p.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time_p.shape[0],4) self.totalflux_ft_p = oi_flux_sc_p.field('TOTALFLUX_FT').reshape(self.oi_flux_sc_time_p.shape[0],4) # Number of frames self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] if self.polarsplit == False: self.nframe_ft = self.oi_flux_ft_time.shape[0] if self.polarsplit == True: self.nframe_ft = self.oi_flux_ft_time_s.shape[0] self.nframe_met = self.oi_vis_met_time.shape[0] # OI_VIS table parameters if (oi_vis_sc is not None) & (self.polarsplit == False): print("Load VIS data") self.time_sc = oi_vis_sc.field('TIME')[::nbase]/1e6 if 'EXP_PHASE' in oi_vis_sc.columns.names: self.oi_vis_sc_exp_phase = oi_vis_sc.field('EXP_PHASE') if 'PHASE_REF' in oi_vis_sc.columns.names: self.oi_vis_sc_exp_phase = oi_vis_sc.field('PHASE_REF') if 'PHI_MET' in oi_vis_sc.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc.field('PHI_MET') if 'PHASE_MET_FC' in oi_vis_sc.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc.field('PHASE_MET_FC') if 'PHASE_MET_TEL' in oi_vis_sc.columns.names: self.oi_vis_sc_phi_met_tel = oi_vis_sc.field('PHASE_MET_TEL') if 'OPD_MET_FC' in oi_vis_sc.columns.names: self.oi_vis_sc_opd_met = oi_vis_sc.field('OPD_MET_FC').reshape(self.nframe_sc,nbase).T*1e6 # in microns if 'OPD_MET_TEL' in oi_vis_sc.columns.names: self.oi_vis_sc_opd_met_tel = oi_vis_sc.field('OPD_MET_TEL')*1e6 # in microns self.oi_vis_sc_visdata = oi_vis_sc.field('VISDATA') self.oi_vis_sc_viserr = oi_vis_sc.field('VISERR') if 'VISAMP' in oi_vis_sc.columns.names: self.oi_vis_sc_visamp = oi_vis_sc.field('VISAMP') self.oi_vis_sc_visphi = oi_vis_sc.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) else: self.oi_vis_sc_visamp = computeVisAmp(oi_vis_sc, oi_flux_sc) self.oi_vis_sc_visphi = np.angle(oi_vis_sc.field('VISDATA'))*np.pi/180. self.oi_vis_sc_ucoord = oi_vis_sc.field('UCOORD') self.oi_vis_sc_vcoord = oi_vis_sc.field('VCOORD') self.oi_vis_sc_sta_index = oi_vis_sc.field('STA_INDEX') if 'VISDATA_FT' in oi_vis_sc.columns.names: self.oi_vis_sc_visdataft = oi_vis_sc.field('VISDATA_FT') if 'VISVAR_FT' in oi_vis_sc.columns.names: self.oi_vis_sc_visvarft = oi_vis_sc.field('VISVAR_FT') self.oi_vis_sc_vispowerft = oi_vis_sc.field('VISPOWER_FT') self.time_ft = oi_vis_ft.field('TIME')[::nbase]/1e6 self.oi_vis_ft_visdata = oi_vis_ft.field('VISDATA') self.oi_vis_ft_viserr = oi_vis_ft.field('VISERR') if 'VISAMP' in oi_vis_ft.columns.names: self.oi_vis_ft_visamp = oi_vis_ft.field('VISAMP') self.oi_vis_ft_visphi = oi_vis_ft.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) else: self.oi_vis_ft_visamp = computeVisAmp(oi_vis_ft, oi_flux_ft) self.oi_vis_ft_visphi = np.angle(oi_vis_ft.field('VISDATA'))*np.pi/180. self.oi_vis_ft_ucoord = oi_vis_ft.field('UCOORD') self.oi_vis_ft_vcoord = oi_vis_ft.field('VCOORD') self.oi_vis_ft_sta_index = oi_vis_ft.field('STA_INDEX') self.oi_vis_ft_target_phase = oi_vis_ft.field('TARGET_PHASE').reshape(self.nframe_ft,nbase) self.oi_vis_ft_state = oi_vis_ft.field('STATE').reshape(self.nframe_ft,nbase) self.oi_vis_ft_reject_flag = oi_vis_ft.field('REJECTION_FLAG').reshape(self.nframe_ft,nbase) # Quadratic VIS2 quantities are COMPUTED and not present in the p2vmreduced file self.oi_vis2_sc = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) for baseline in range(0,nbase): self.oi_vis2_sc[baseline,:,:] = self.oi_vis_sc_visamp[baseline::nbase,:]**2 * np.sign(self.oi_vis_sc_visamp[baseline::nbase,:]) self.oi_vis2_ft = np.zeros((nbase,self.nframe_ft,self.nwave_ft)) for baseline in range(0,nbase): self.oi_vis2_ft[baseline,:,:] = self.oi_vis_ft_visamp[baseline::nbase,:]**2 * np.sign(self.oi_vis_ft_visamp[baseline::nbase,:]) if self.gdelay_in == True: # Group delay in microns self.gdelay_sc = oi_vis_sc.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 if 'GDELAY_FT' in oi_vis_sc.columns.names: self.gdelay_ft_resamp = oi_vis_sc.field('GDELAY_FT').reshape(self.nframe_sc,nbase).T*1e6 # GD of FT at SC frequency self.gdelay_ft = oi_vis_ft.field('GDELAY').reshape(self.nframe_ft,nbase).T*1e6 self.gdelay_boot_ft = oi_vis_ft.field('GDELAY_BOOT').reshape(self.nframe_ft,nbase).T*1e6 # Signal to noise of fringes self.snr_sc = oi_vis_sc.field('SNR').reshape(self.nframe_sc,nbase).T self.snr_ft = oi_vis_ft.field('SNR').reshape(self.nframe_ft,nbase).T self.snr_boot_ft = oi_vis_ft.field('SNR_BOOT').reshape(self.nframe_ft,nbase).T if 'V_FACTOR' in column_oi_vis_sc: print("Load VFACTOR data") self.vfactor_in = True self.oi_vis_vfactor_ft = oi_vis_sc.field('V_FACTOR_FT').reshape(self.nframe_sc,nbase,self.nwave_ft) self.oi_vis_vfactor = oi_vis_sc.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) self.oi_vis_reject_flag = oi_vis_sc.field('REJECTION_FLAG').reshape(self.nframe_sc,nbase) self.oi_vis_fringe_det_ratio = oi_vis_sc.field('FRINGEDET_RATIO').reshape(self.nframe_sc,nbase) if 'P_FACTOR' in column_oi_vis_sc: self.oi_vis_pfactor = oi_vis_sc.field('P_FACTOR') if (oi_vis_sc_s is not None) & (self.polarsplit == True): self.time_sc = oi_vis_sc_s.field('TIME')[::nbase]/1e6 if 'EXP_PHASE' in oi_vis_sc_s.columns.names: self.oi_vis_sc_exp_phase_s = oi_vis_sc_s.field('EXP_PHASE') self.oi_vis_sc_exp_phase_p = oi_vis_sc_p.field('EXP_PHASE') if 'PHASE_REF' in oi_vis_sc_s.columns.names: self.oi_vis_sc_exp_phase_s = oi_vis_sc_s.field('PHASE_REF') self.oi_vis_sc_exp_phase_p = oi_vis_sc_p.field('PHASE_REF') if 'PHI_MET' in oi_vis_sc_s.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc_s.field('PHI_MET') # self.oi_vis_sc_phi_met_p = oi_vis_sc_p.field('PHI_MET') if 'PHASE_MET_FC' in oi_vis_sc_s.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc_s.field('PHASE_MET_FC') # self.oi_vis_sc_phi_met_p = oi_vis_sc_p.field('PHASE_MET_FC') if 'PHASE_MET_TEL' in oi_vis_sc_s.columns.names: self.oi_vis_sc_phi_met_tel = oi_vis_sc_s.field('PHASE_MET_TEL') #self.oi_vis_sc_phi_met_tel_p = oi_vis_sc_p.field('PHASE_MET_TEL') if 'OPD_MET_FC' in oi_vis_sc_s.columns.names: self.oi_vis_sc_opd_met = oi_vis_sc_s.field('OPD_MET_FC').reshape(self.nframe_sc,nbase).T*1e6 # in microns #self.oi_vis_sc_opd_met_p = oi_vis_sc_p.field('OPD_MET_FC').reshape(self.nframe_sc,nbase).T*1e6 # in microns if 'OPD_MET_TEL' in oi_vis_sc_s.columns.names: self.oi_vis_sc_opd_met_tel = oi_vis_sc_s.field('OPD_MET_TEL')*1e6 # in microns #self.oi_vis_sc_opd_met_tel_p = oi_vis_sc_p.field('OPD_MET_TEL') self.oi_vis_sc_visdata_s = oi_vis_sc_s.field('VISDATA') self.oi_vis_sc_viserr_s = oi_vis_sc_s.field('VISERR') if 'VISAMP' in oi_vis_sc_s.columns.names: self.oi_vis_sc_visamp_s = oi_vis_sc_s.field('VISAMP') self.oi_vis_sc_visphi_s = oi_vis_sc_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) else: self.oi_vis_sc_visamp_s = computeVisAmp(oi_vis_sc_s, oi_flux_sc_s) self.oi_vis_sc_visphi_s = np.angle(oi_vis_sc_s.field('VISDATA'))*np.pi/180. self.oi_vis_sc_ucoord_s = oi_vis_sc_s.field('UCOORD') self.oi_vis_sc_vcoord_s = oi_vis_sc_s.field('VCOORD') self.oi_vis_sc_sta_index_s = oi_vis_sc_s.field('STA_INDEX') if 'VISDATA_FT' in oi_vis_sc_s.columns.names: self.oi_vis_sc_visdataft_s = oi_vis_sc_s.field('VISDATA_FT') self.oi_vis_sc_visdataft_p = oi_vis_sc_p.field('VISDATA_FT') if 'VISVAR_FT' in oi_vis_sc_s.columns.names: self.oi_vis_sc_visvarft_s = oi_vis_sc_s.field('VISVAR_FT') self.oi_vis_sc_visvarft_p = oi_vis_sc_p.field('VISVAR_FT') self.oi_vis_sc_vispowerft_s = oi_vis_sc_s.field('VISPOWER_FT') self.oi_vis_sc_vispowerft_p = oi_vis_sc_p.field('VISPOWER_FT') self.oi_vis_sc_visdata_p = oi_vis_sc_p.field('VISDATA') self.oi_vis_sc_viserr_p = oi_vis_sc_p.field('VISERR') if 'VISAMP' in oi_vis_sc_p.columns.names: self.oi_vis_sc_visamp_p = oi_vis_sc_p.field('VISAMP') self.oi_vis_sc_visphi_p = oi_vis_sc_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) else: self.oi_vis_sc_visamp_p = computeVisAmp(oi_vis_sc_p, oi_flux_sc_p) self.oi_vis_sc_visphi_p = np.angle(oi_vis_sc_p.field('VISDATA'))*np.pi/180. self.oi_vis_sc_ucoord_p = oi_vis_sc_p.field('UCOORD') self.oi_vis_sc_vcoord_p = oi_vis_sc_p.field('VCOORD') self.oi_vis_sc_sta_index_p = oi_vis_sc_p.field('STA_INDEX') self.time_ft = oi_vis_ft_s.field('TIME')[::nbase]/1e6 self.oi_vis_ft_visdata_s = oi_vis_ft_s.field('VISDATA') self.oi_vis_ft_viserr_s = oi_vis_ft_s.field('VISERR') if 'VISAMP' in oi_vis_ft_s.columns.names: self.oi_vis_ft_visamp_s = oi_vis_ft_s.field('VISAMP') self.oi_vis_ft_visphi_s = oi_vis_ft_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) else: self.oi_vis_ft_visamp_s = computeVisAmp(oi_vis_ft_s, oi_flux_ft_s) self.oi_vis_ft_visphi_s = np.angle(oi_vis_ft_s.field('VISDATA'))*np.pi/180. self.oi_vis_ft_ucoord_s = oi_vis_ft_s.field('UCOORD') self.oi_vis_ft_vcoord_s = oi_vis_ft_s.field('VCOORD') self.oi_vis_ft_sta_index_s = oi_vis_ft_s.field('STA_INDEX') self.oi_vis_ft_visdata_p = oi_vis_ft_p.field('VISDATA') self.oi_vis_ft_viserr_p = oi_vis_ft_p.field('VISERR') if 'VISAMP' in oi_vis_ft_p.columns.names: self.oi_vis_ft_visamp_p = oi_vis_ft_p.field('VISAMP') self.oi_vis_ft_visphi_p = oi_vis_ft_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) else: self.oi_vis_ft_visamp_p = computeVisAmp(oi_vis_ft_p, oi_flux_ft_p) self.oi_vis_ft_visphi_p = np.angle(oi_vis_ft_p.field('VISDATA'))*np.pi/180. self.oi_vis_ft_ucoord_p = oi_vis_ft_p.field('UCOORD') self.oi_vis_ft_vcoord_p = oi_vis_ft_p.field('VCOORD') self.oi_vis_ft_sta_index_p = oi_vis_ft_p.field('STA_INDEX') self.oi_vis_ft_target_phase_s = oi_vis_ft_s.field('TARGET_PHASE').reshape(self.nframe_ft,nbase) self.oi_vis_ft_state_s = oi_vis_ft_s.field('STATE').reshape(self.nframe_ft,nbase) self.oi_vis_ft_reject_flag_s = oi_vis_ft_s.field('REJECTION_FLAG').reshape(self.nframe_ft,nbase) self.oi_vis_ft_target_phase_p = oi_vis_ft_p.field('TARGET_PHASE').reshape(self.nframe_ft,nbase) self.oi_vis_ft_state_p = oi_vis_ft_p.field('STATE').reshape(self.nframe_ft,nbase) self.oi_vis_ft_reject_flag_p = oi_vis_ft_p.field('REJECTION_FLAG').reshape(self.nframe_ft,nbase) if self.gdelay_in == True: # Group delay in microns self.gdelay_sc_s = oi_vis_sc_s.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 self.gdelay_sc_p = oi_vis_sc_p.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 if 'GDELAY_FT' in oi_vis_sc_s.columns.names: self.gdelay_ft_resamp_s = oi_vis_sc_s.field('GDELAY_FT').reshape(self.nframe_sc,nbase).T*1e6 # GD of FT at SC frequency self.gdelay_ft_resamp_p = oi_vis_sc_p.field('GDELAY_FT').reshape(self.nframe_sc,nbase).T*1e6 # GD of FT at SC frequency self.gdelay_ft_s = oi_vis_ft_s.field('GDELAY').reshape(self.nframe_ft,nbase).T*1e6 # at FT frequency self.gdelay_ft_p = oi_vis_ft_p.field('GDELAY').reshape(self.nframe_ft,nbase).T*1e6 # at FT frequency # self.gdelay_boot_sc_s = oi_vis_sc_s.field('GDELAY_BOOT').reshape(self.nframe_sc,nbase).T*coh_length_sc # self.gdelay_boot_sc_p = oi_vis_sc_p.field('GDELAY_BOOT').reshape(self.nframe_sc,nbase).T*coh_length_sc self.gdelay_boot_ft_s = oi_vis_ft_s.field('GDELAY_BOOT').reshape(self.nframe_ft,nbase).T*1e6 self.gdelay_boot_ft_p = oi_vis_ft_p.field('GDELAY_BOOT').reshape(self.nframe_ft,nbase).T*1e6 # Signal to noise of fringes self.snr_sc_s = oi_vis_sc_s.field('SNR').reshape(self.nframe_sc,nbase).T self.snr_sc_p = oi_vis_sc_p.field('SNR').reshape(self.nframe_sc,nbase).T self.snr_ft_s = oi_vis_ft_s.field('SNR').reshape(self.nframe_ft,nbase).T self.snr_ft_p = oi_vis_ft_p.field('SNR').reshape(self.nframe_ft,nbase).T # self.snr_boot_sc_s = oi_vis_sc_s.field('SNR_BOOT').reshape(self.nframe_sc,nbase).T # self.snr_boot_sc_p = oi_vis_sc_p.field('SNR_BOOT').reshape(self.nframe_sc,nbase).T self.snr_boot_ft_s = oi_vis_ft_s.field('SNR_BOOT').reshape(self.nframe_ft,nbase).T self.snr_boot_ft_p = oi_vis_ft_p.field('SNR_BOOT').reshape(self.nframe_ft,nbase).T if 'V_FACTOR' in column_oi_vis_sc_p: self.vfactor_in = True self.oi_vis_vfactor_ft_s = oi_vis_sc_s.field('V_FACTOR_FT').reshape(self.nframe_sc,nbase,self.nwave_ft) # time_sc, baseline, wavelength_ft self.oi_vis_vfactor_s = oi_vis_sc_s.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) self.oi_vis_vfactor_ft_p = oi_vis_sc_p.field('V_FACTOR_FT').reshape(self.nframe_sc,nbase,self.nwave_ft) self.oi_vis_vfactor_p = oi_vis_sc_p.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) self.oi_vis_reject_flag_s = oi_vis_sc_s.field('REJECTION_FLAG').reshape(self.nframe_sc,nbase) self.oi_vis_reject_flag_p = oi_vis_sc_p.field('REJECTION_FLAG').reshape(self.nframe_sc,nbase) self.oi_vis_fringe_det_ratio_s = oi_vis_sc_s.field('FRINGEDET_RATIO').reshape(self.nframe_sc,nbase) self.oi_vis_fringe_det_ratio_p = oi_vis_sc_p.field('FRINGEDET_RATIO').reshape(self.nframe_sc,nbase) if 'P_FACTOR' in column_oi_vis_sc_p: self.oi_vis_pfactor_s = oi_vis_sc_s.field('P_FACTOR') self.oi_vis_pfactor_p = oi_vis_sc_p.field('P_FACTOR') # Quadratic VIS2 quantities are COMPUTED and not present in the p2vmreduced files self.oi_vis2_sc_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) self.oi_vis2_sc_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) for baseline in range(0,nbase): self.oi_vis2_sc_s[baseline,:,:] = self.oi_vis_sc_visamp_s[baseline::nbase,:]**2 * np.sign(self.oi_vis_sc_visamp_s[baseline::nbase,:]) self.oi_vis2_sc_p[baseline,:,:] = self.oi_vis_sc_visamp_p[baseline::nbase,:]**2 * np.sign(self.oi_vis_sc_visamp_p[baseline::nbase,:]) self.oi_vis2_ft_s = np.zeros((nbase,self.nframe_ft,self.nwave_ft)) self.oi_vis2_ft_p = np.zeros((nbase,self.nframe_ft,self.nwave_ft)) for baseline in range(0,nbase): self.oi_vis2_ft_s[baseline,:,:] = self.oi_vis_ft_visamp_s[baseline::nbase,:]**2 * np.sign(self.oi_vis_ft_visamp_s[baseline::nbase,:]) self.oi_vis2_ft_p[baseline,:,:] = self.oi_vis_ft_visamp_p[baseline::nbase,:]**2 * np.sign(self.oi_vis_ft_visamp_p[baseline::nbase,:]) #============================================================================== # Computation of differential OPD quantities #============================================================================== if (self.gdelay_in == True) & (self.polarsplit == False): self.mean_gd_sc = np.zeros(nbase) self.mean_gd_ft = np.zeros(nbase) self.gd_diff = np.zeros((nbase,self.nframe_sc)) self.opd_diff_avg = np.zeros((nbase,self.nframe_sc)) self.opd_diff = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.total_dopd_avg = np.zeros((nbase,self.nframe_sc)) self.total_dopd = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC for baseline in range(0,nbase): # time average of GD over the sequence self.mean_gd_sc[baseline] = np.nanmean(self.gdelay_sc[baseline,:]) self.mean_gd_ft[baseline] = np.nanmean(self.gdelay_ft_resamp[baseline,:]) # Differential group delay in microns self.gd_diff[baseline,:] = self.gdelay_sc[baseline,:] - self.gdelay_ft_resamp[baseline,:] - (self.mean_gd_sc[baseline] - self.mean_gd_ft[baseline]) # Wavelength average of differential phase, in microns, once the mean GD has been subtracted visdata_sc = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_sc[baseline]/self.wave_sc)) for data in self.oi_vis_sc_visdata[baseline::nbase,:] ]) visdata_ft = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_ft[baseline]/self.wave_ft)) for data in self.oi_vis_sc_visdataft[baseline::nbase,:] ]) ph_diff_avg = visdata_sc * np.conj(visdata_ft) self.opd_diff_avg[baseline,:] = np.angle ( ph_diff_avg * np.conj(np.mean(ph_diff_avg)) ) / (2.*np.pi) * np.mean(self.wave_sc) # Differential phase in microns, function of wavelength visdata_sc = self.oi_vis_sc_visdata[baseline::nbase,:] * np.exp(1.j*self.oi_vis_sc_exp_phase[baseline::nbase,:]) for wave in range(0,self.nwave_sc): # self.opd_diff[baseline,:,wave] = np.angle( visdata_sc[:,wave] * np.conj(visdata_sc[:,wave]) ) / (2.*np.pi) * self.wave_sc[wave] self.opd_diff[baseline,:,wave] = np.angle( visdata_sc[:,wave] * np.conj(np.mean(visdata_sc[:,wave])) ) / (2.*np.pi) * self.wave_sc[wave] # Construction of the total dOPD SC-FT including metrology and GD self.total_dopd_avg[baseline,:] = self.opd_diff_avg[baseline,:] + (self.mean_gd_sc[baseline] - self.mean_gd_ft[baseline]) + self.oi_vis_sc_opd_met[baseline,:] # Construction of the total dOPD SC-FT including metrology and GD as a function of wavelength for wave in range(0,self.nwave_sc): self.total_dopd[baseline,:, wave] = self.opd_diff[baseline,:,wave] + (self.mean_gd_sc[baseline] - self.mean_gd_ft[baseline]) + self.oi_vis_sc_opd_met[baseline,:] if (self.gdelay_in == True) & (self.polarsplit == True): self.mean_gd_sc_s = np.zeros(nbase) self.mean_gd_sc_p = np.zeros(nbase) self.mean_gd_ft_s = np.zeros(nbase) self.mean_gd_ft_p = np.zeros(nbase) self.gd_diff_s = np.zeros((nbase,self.nframe_sc)) self.gd_diff_p = np.zeros((nbase,self.nframe_sc)) self.opd_diff_avg_s = np.zeros((nbase,self.nframe_sc)) self.opd_diff_avg_p = np.zeros((nbase,self.nframe_sc)) self.opd_diff_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.opd_diff_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.total_dopd_avg_s = np.zeros((nbase,self.nframe_sc)) self.total_dopd_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.total_dopd_avg_p = np.zeros((nbase,self.nframe_sc)) self.total_dopd_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC for baseline in range(0,nbase): self.mean_gd_sc_s[baseline] = np.nanmean(self.gdelay_sc_s[baseline,:]) self.mean_gd_sc_p[baseline] = np.nanmean(self.gdelay_sc_p[baseline,:]) self.mean_gd_ft_s[baseline] = np.nanmean(self.gdelay_ft_resamp_s[baseline,:]) self.mean_gd_ft_p[baseline] = np.nanmean(self.gdelay_ft_resamp_p[baseline,:]) # Differential group delay in microns self.gd_diff_s[baseline,:] = self.gdelay_sc_s[baseline,:] - self.gdelay_ft_resamp_s[baseline,:] - (self.mean_gd_sc_s[baseline] - self.mean_gd_ft_s[baseline]) self.gd_diff_p[baseline,:] = self.gdelay_sc_p[baseline,:] - self.gdelay_ft_resamp_p[baseline,:] - (self.mean_gd_sc_p[baseline] - self.mean_gd_ft_p[baseline]) # Differential phase average in the band, in microns, once the mean_gd has been removed visdata_sc_s = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_sc_s[baseline]/self.wave_sc)) for data in self.oi_vis_sc_visdata_s[baseline::nbase,:] ]) visdata_sc_p = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_sc_p[baseline]/self.wave_sc)) for data in self.oi_vis_sc_visdata_p[baseline::nbase,:] ]) visdata_ft_s = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_ft_s[baseline]/self.wave_ft)) for data in self.oi_vis_sc_visdataft_s[baseline::nbase,:] ]) visdata_ft_p = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_ft_p[baseline]/self.wave_ft)) for data in self.oi_vis_sc_visdataft_p[baseline::nbase,:] ]) ph_diff_avg_s = visdata_sc_s * np.conj(visdata_ft_s); ph_diff_avg_p = visdata_sc_p * np.conj(visdata_ft_p); self.opd_diff_avg_s[baseline,:] = np.angle ( ph_diff_avg_s * np.conj(np.mean(ph_diff_avg_s)) ) / (2.*np.pi) * np.mean(self.wave_sc) self.opd_diff_avg_p[baseline,:] = np.angle ( ph_diff_avg_p * np.conj(np.mean(ph_diff_avg_p)) ) / (2.*np.pi) * np.mean(self.wave_sc) # Differential phase in microns, function of wavelength visdata_sc_s = self.oi_vis_sc_visdata_s[baseline::nbase,:] * np.exp (1.j*self.oi_vis_sc_exp_phase_s[baseline::nbase,:]) visdata_sc_p = self.oi_vis_sc_visdata_p[baseline::nbase,:] * np.exp (1.j*self.oi_vis_sc_exp_phase_p[baseline::nbase,:]) for wave in range(0,self.nwave_sc): self.opd_diff_s[baseline,:,wave] = np.angle ( visdata_sc_s[:,wave] * np.conj(np.mean(visdata_sc_s[:,wave])) ) / (2.*np.pi) * self.wave_sc[wave] self.opd_diff_p[baseline,:,wave] = np.angle ( visdata_sc_p[:,wave] * np.conj(np.mean(visdata_sc_p[:,wave])) ) / (2.*np.pi) * self.wave_sc[wave] # Construction of the total dOPD SC-FT including metrology and GD self.total_dopd_avg_s[baseline,:] = self.opd_diff_avg_s[baseline,:] + (self.mean_gd_sc_s[baseline] - self.mean_gd_ft_s[baseline]) + self.oi_vis_sc_opd_met[baseline,:] self.total_dopd_avg_p[baseline,:] = self.opd_diff_avg_p[baseline,:] + (self.mean_gd_sc_p[baseline] - self.mean_gd_ft_p[baseline]) + self.oi_vis_sc_opd_met[baseline,:] # Construction of the total dOPD SC-FT including metrology and GD as a function of wavelength for wave in range(0,self.nwave_sc): self.total_dopd_s[baseline,:, wave] = self.opd_diff_s[baseline,:,wave] + (self.mean_gd_sc_s[baseline] - self.mean_gd_ft_s[baseline]) + self.oi_vis_sc_opd_met[baseline,:] self.total_dopd_p[baseline,:, wave] = self.opd_diff_p[baseline,:,wave] + (self.mean_gd_sc_p[baseline] - self.mean_gd_ft_p[baseline]) + self.oi_vis_sc_opd_met[baseline,:] # #============================================================================== # # Computation of the group delay signals # #============================================================================== # # # Smoothing of the FT group delay signals # smooth_gd_ft = 30 # number of samples for smoothing # # if (self.polarsplit == False): # nwave_sc = self.nwave_sc # nwave_ft = self.nwave_ft # # # SC group delay signals # visdata_sc = np.zeros((nbase,self.nframe_sc,self.nwave_sc),dtype=complex) # for base in range(0,nbase): # for wave in range(0,self.nwave_sc): # visdata_sc[base,:,wave] = self.oi_vis_sc_visdata[base::nbase,wave] # # lambda_mean_sc = np.mean(self.wave_sc) # delta_lambda_sc = self.wave_sc[1]-self.wave_sc[0] # self.group_delay_sc = np.zeros((nbase,self.nframe_sc)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_sc): # incstep_sc = visdata_sc[base,frame,0:nwave_sc-1]*np.conjugate(visdata_sc[base,frame,1:nwave_sc]) # self.group_delay_sc[base,frame] = np.angle(np.sum(incstep_sc))/(2*np.pi)*(lambda_mean_sc**2)/delta_lambda_sc # # # FT group delay signals # visdata_ft = np.zeros((nbase,self.nframe_ft,self.nwave_ft),dtype=complex) # for base in range(0,nbase): # for wave in range(0,self.nwave_ft): # visdata_ft[base,:,wave] = self.oi_vis_ft_visdata[base::nbase,wave] # # lambda_mean_ft = np.mean(self.wave_ft) # delta_lambda_ft = self.wave_ft[1]-self.wave_ft[0] # incstep_ft = np.zeros((nbase,self.nframe_ft),dtype=complex) # self.group_delay_ft = np.zeros((nbase,self.nframe_ft)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # incstep_ft[base,frame] = np.sum(visdata_ft[base,frame,0:nwave_ft-1]*np.conjugate(visdata_ft[base,frame,1:nwave_ft])) # # for base in range(0,nbase): # smoothing of the inter-spectra incstep # incstep_ft[base,:] = np.convolve(incstep_ft[base,:], np.ones((smooth_gd_ft,))/smooth_gd_ft, mode='same') # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # self.group_delay_ft[base,frame] = np.angle(incstep_ft[base,frame])/(2*np.pi)*(lambda_mean_ft**2)/delta_lambda_ft # # if (self.polarsplit == True): # nwave_sc = self.nwave_sc # nwave_ft = self.nwave_ft # # # SC group delay signals # visdata_sc_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc),dtype=complex) # visdata_sc_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc),dtype=complex) # for base in range(0,nbase): # for wave in range(0,self.nwave_sc): # visdata_sc_s[base,:,wave] = self.oi_vis_sc_visdata_s[base::nbase,wave] # visdata_sc_p[base,:,wave] = self.oi_vis_sc_visdata_p[base::nbase,wave] # # lambda_mean_sc = np.mean(self.wave_sc) # delta_lambda_sc = self.wave_sc[1]-self.wave_sc[0] # coh_length_sc = (lambda_mean_sc**2)/delta_lambda_sc # # self.group_delay_sc_s = np.zeros((nbase,self.nframe_sc)) # self.group_delay_sc_p = np.zeros((nbase,self.nframe_sc)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_sc): # incstep_sc_s = visdata_sc_s[base,frame,0:nwave_sc-1]*np.conjugate(visdata_sc_s[base,frame,1:nwave_sc]) # incstep_sc_p = visdata_sc_p[base,frame,0:nwave_sc-1]*np.conjugate(visdata_sc_p[base,frame,1:nwave_sc]) # self.group_delay_sc_s[base,frame] = np.angle(np.sum(incstep_sc_s))/(2*np.pi)*coh_length_sc # self.group_delay_sc_p[base,frame] = np.angle(np.sum(incstep_sc_p))/(2*np.pi)*coh_length_sc # # # FT group delay signals # visdata_ft_s = np.zeros((nbase,self.nframe_ft,self.nwave_ft),dtype=complex) # visdata_ft_p = np.zeros((nbase,self.nframe_ft,self.nwave_ft),dtype=complex) # # for base in range(0,nbase): # for wave in range(0,self.nwave_ft): # visdata_ft_s[base,:,wave] = self.oi_vis_ft_visdata_s[base::nbase,wave] # visdata_ft_p[base,:,wave] = self.oi_vis_ft_visdata_p[base::nbase,wave] # # lambda_mean_ft = np.mean(self.wave_ft) # delta_lambda_ft = self.wave_ft[1]-self.wave_ft[0] # coh_length_ft = (lambda_mean_ft**2)/delta_lambda_ft # # incstep_ft_s = np.zeros((nbase,self.nframe_ft),dtype=complex) # incstep_ft_p = np.zeros((nbase,self.nframe_ft),dtype=complex) # self.group_delay_ft_s = np.zeros((nbase,self.nframe_ft)) # self.group_delay_ft_p = np.zeros((nbase,self.nframe_ft)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # incstep_ft_s[base,frame] = np.sum(visdata_ft_s[base,frame,0:nwave_ft-1]*np.conjugate(visdata_ft_s[base,frame,1:nwave_ft])) # incstep_ft_p[base,frame] = np.sum(visdata_ft_p[base,frame,0:nwave_ft-1]*np.conjugate(visdata_ft_p[base,frame,1:nwave_ft])) # # for base in range(0,nbase): # smoothing of the inter-spectra incstep # incstep_ft_s[base,:] = np.convolve(incstep_ft_s[base,:], np.ones((smooth_gd_ft,))/smooth_gd_ft, mode='same') # incstep_ft_p[base,:] = np.convolve(incstep_ft_p[base,:], np.ones((smooth_gd_ft,))/smooth_gd_ft, mode='same') # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # self.group_delay_ft_s[base,frame] = np.angle(incstep_ft_s[base,frame])/(2*np.pi)*coh_length_ft # self.group_delay_ft_p[base,frame] = np.angle(incstep_ft_p[base,frame])/(2*np.pi)*coh_length_ft # #============================================================================== # # Computation of the K factors # #============================================================================== # # k2factor = np.zeros((nbase,self.nframe_sc)) # dt_sc = self.sc_time[nbase]-self.sc_time[0] # dt_ft = self.ft_time[nbase]-self.ft_time[0] # # if self.polarsplit == True: # # FT frames before the first SC integration (5s) # skip_ft = (self.sc_time[0]-dt_sc/2) // dt_ft # visdata_ft_s = np.mean(self.oi_vis_ft_visdata_s.reshape(self.ft_time.shape[0],6,self.nwave_ft),axis=2) # time, baseline, wavelength # visdata_ft_p = np.mean(self.oi_vis_ft_visdata_p.reshape(self.ft_time.shape[0],6,self.nwave_ft),axis=2) # time, baseline, wavelength # visdata_ft = np.mean([visdata_ft_s, visdata_ft_p],axis=0) # else: # dt_sc = self.oi_vis_sc_time[nbase]-self.oi_vis_sc_time_s[0] # dt_ft = self.oi_vis_ft_time[nbase]-self.oi_vis_ft_time_s[0] # # FT frames before the first SC integration (5s) # skip_ft = (self.oi_vis_sc_time[0]-dt_sc/2) // dt_ft # visdata_ft = np.mean(self.oi_vis_ft_visdata.reshape(self.ft_time.shape[0],6,self.nwave_ft),axis=2) # time, baseline, wavelength # # # Ratio of the inter-frame time between SC and FT # ratio_sc_ft = dt_sc/dt_ft # # # for baseline in range(0,nbase): # for frame_sc in range(0,self.nframe_sc): # minindex_ft = int(ratio_sc_ft * frame_sc)+skip_ft # maxindex_ft = minindex_ft + int(ratio_sc_ft) # # bias_coh = np.var(np.abs(visdata_ft[minindex_ft:maxindex_ft,baseline])) # v2_coh = np.mean(np.abs(np.square(visdata_ft[minindex_ft:maxindex_ft,baseline]))) # v2_incoh = np.square(np.abs(np.mean(visdata_ft[minindex_ft:maxindex_ft,baseline]))) # k2factor[baseline, frame_sc] = v2_coh/v2_incoh # # self.kfactor = np.sqrt(k2factor) # # print self.kfactor gravi_p2vmreduced.close() class Rawdata: def __init__(self, filename): # open fits file self.filename = filename gravi_raw = pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_raw[0].header # Single mode of the instrument (no metrology data) # if 'SINGLE' in str(self.header['HIERARCH ESO DO CATG']): # self.single = True # else: # self.single = False # Raw data types if 'HIERARCH ESO DPR TYPE' in self.header: dtype = str(self.header['HIERARCH ESO DPR TYPE']) self.catg = dtype # default case # if 'DARK' in dtype: self.catg = 'DARK' # if 'FLAT' in dtype: self.catg = 'FLAT' # if 'WAVE' in dtype: self.catg = 'WAVE' # if 'P2VM' in dtype: self.catg = 'P2VM' else: self.catg = 'UNKNOWN' # Interferometer configuration self.stations = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION4','S4'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION3','S3'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION2','S2'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF STATION1','S1')] self.telnames = [get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T4NAME','T4'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T3NAME','T3'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T2NAME','T2'),\ get_key_withdefault(self.header,'HIERARCH ESO ISS CONF T1NAME','T1')] # baseline names self.basenames = [self.stations[0]+self.stations[1], self.stations[0]+self.stations[2], self.stations[0]+self.stations[3], self.stations[1]+self.stations[2], self.stations[1]+self.stations[3], self.stations[2]+self.stations[3]] # Polarization setup if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True fddl = None # in case the FDDL data are not in the file metrology = None detector_sc = None detector_ft = None self.gdelay_in = False # Number of SC frames self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] self.exptime_sc = self.header['HIERARCH ESO DET2 SEQ1 DIT'] from astropy.time import Time def mjd_date(datestring): datestring2 = datestring[:10]+' '+datestring[12:] t = Time([datestring2],format='iso',scale='utc') return t.mjd[0] # get the ACQ, ARRAY_GEOMETRY, FDDL, OPDC, imaging_detector, imaging_data and METROLOGY extentions for hdu in gravi_raw : if hdu.name == 'ARRAY_GEOMETRY' : # array geometry oi_array=hdu.data if hdu.name == 'IMAGING_DATA_ACQ' : # acquisition camera frames self.acq = hdu.data # direct assignation #if hdu.name == 'OPDC': # opdc = hdu.data if hdu.name == 'FDDL': fddl = hdu.data if hdu.name == 'IMAGING_DETECTOR_SC' : detector_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_FT' : detector_ft=hdu.data # IMAGING_DATA SC frames direct assignation if hdu.name == 'IMAGING_DATA_SC' : self.data_sc = hdu.data self.time_sc = np.array(list(range(0,self.nframe_sc))) * self.header['HIERARCH ESO INS TIM1 PERIOD'] + \ mjd_date(self.header['HIERARCH ESO INS TIM1 START'])*86400. # reference in the middle of the inter-frame if hdu.name == 'IMAGING_DATA_FT' : imaging_ft=hdu.data if hdu.name == 'METROLOGY' : metrology=hdu.data if hdu.name == 'OPDC' : opdc=hdu.data self.gdelay_in = True # OPDC parameters if self.gdelay_in == True: self.opdc_time = opdc.field('TIME')/1.0e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.opdc_piezo_offset = opdc.field('PIEZO_DL_OFFSET') # in VOLTS self.opdc_vltidl_offset = opdc.field('VLTI_DL_OFFSET') self.opdc_state = opdc.field('STATE') # OI_ARRAY parameters self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1.0e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # IMAGING_DETECTOR parameters if detector_sc is not None: self.regname_sc = detector_sc.field('REGNAME') else: self.regname_sc = '' if detector_ft is not None: self.regname_ft = detector_ft.field('REGNAME') else: self.regname_ft = '' # IMAGING_DATA_FT self.time_ft = imaging_ft.field('TIME')/1.0e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.exptime_ft = self.header['HIERARCH ESO DET3 SEQ1 DIT'] data = imaging_ft.field('PIX') if self.polarsplit == True: self.data_ft = np.zeros((self.time_ft.shape[0],48,5)) for i in range(0,24): for j in range(0,5): self.data_ft[:,2*i,j] = data[:,j,i] for j in range(5,10): self.data_ft[:,2*i+1,j-5] = data[:,j,i] # raw.data_ft[time,output,wave_channel] else: self.data_ft = np.zeros((self.time_ft.shape[0],24,5)) for i in range(0,24): for j in range(0,5): self.data_ft[:,i,j] = data[:,j,i] self.nframe_ft = self.time_ft.shape[0] # number of FT frames # print self.data_ft # METROLOGY parameters if metrology is not None: self.metrology_time = metrology.field('TIME')/1e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.metrology_volt = metrology.field('VOLT') self.metrology_powerlaser = metrology.field('POWER_LASER') self.metrology_lambdalaser = metrology.field('LAMBDA_LASER') # Data structure consistency checks # ================================= self.datacheck = [] # verification of the number of extensions self.nextensions = len(gravi_raw[:]) if(self.nextensions != 13): self.datacheck.append("# Bad number of FITS extensions") print(" *** ERROR: Bad number of FITS extensions (13 expexted, %i found)." % (self.nextensions)) # time coverage of the RMN data vs SC data mintimes = [self.opdc_time[0], self.fddl_time[0], self.time_ft[0], self.metrology_time[0]] maxtimes = [self.opdc_time[-1], self.fddl_time[-1], self.time_ft[-1], self.metrology_time[-1]] if (self.time_sc[0] < np.max(mintimes)) or (self.time_sc[-1] > np.min(maxtimes)): self.datacheck.append("# SC time range not covered by RMN data") print(" *** ERROR: SC time range not covered by RMN data (s) [OPDC, FDDL, FT, MET]") print(" Min times : ", mintimes) print(" Max times : ", maxtimes) print(" SC times : ", self.time_sc[0], self.time_sc[-1]) # verification of the length of reflective memory data rangemax = 0.020 # maximum acceptable range in seconds self.opdc_time_range = self.opdc_time[-1] - self.opdc_time[0] self.fddl_time_range = self.fddl_time[-1] - self.fddl_time[0] self.ft_time_range = self.time_ft[-1] - self.time_ft[0] met_time = imaging_ft.field('TIME')/1e6 self.metrology_time_range = met_time[-1] - met_time[0] ranges = [self.opdc_time_range, self.fddl_time_range, self.ft_time_range, self.metrology_time_range] if (np.abs(np.max(ranges) - np.min(ranges)) > rangemax): self.datacheck.append("# Inconsistent time range of RMN data") print(" *** ERROR: Inconsistent time range of RMN data (s) [OPDC, FDDL, FT, MET]:") print(" ", ranges) # verification of the continuity of time scales of reflective memory data stepoffsetmax = 0.20 # maximum acceptable change in time step over scales (relative) self.opdc_time_steps = self.opdc_time[1:len(self.opdc_time)-1] - self.opdc_time[0:len(self.opdc_time)-2] print("OPDC min step, max step, mean step:", min(self.opdc_time_steps), max(self.opdc_time_steps), np.mean(self.opdc_time_steps)) print("OPDC max step position:", np.where(self.opdc_time_steps==max(self.opdc_time_steps))) self.fddl_time_steps = self.fddl_time[1:len(self.fddl_time)-1] - self.fddl_time[0:len(self.fddl_time)-2] self.time_ft_steps = self.time_ft[1:len(self.time_ft)-1] - self.time_ft[0:len(self.time_ft)-2] self.metrology_time_steps = met_time[1:len(met_time)-1] - met_time[0:len(met_time)-2] stepranges = [np.ptp(self.opdc_time_steps)/np.mean(self.opdc_time_steps), np.ptp(self.fddl_time_steps)/np.mean(self.fddl_time_steps), np.ptp(self.time_ft_steps)/np.mean(self.time_ft_steps), np.ptp(self.metrology_time_steps)/np.mean(self.metrology_time_steps)] if np.max(stepranges) > stepoffsetmax: self.datacheck.append("# Variable RMN time steps") print(" *** ERROR: Variable RMN time steps (rel. to mean) [OPDC, FDDL, FT, MET]: ") print(" ", stepranges) if self.datacheck == []: self.datacheck = "# Good" gravi_raw.close() class Badpix: def __init__(self, filename): # open fits file self.filename = filename gravi_badpix=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_badpix[0].header for hdu in gravi_badpix : if hdu.name=='IMAGING_DATA_SC' : h = hdu.header self.ncol_sc= h['NAXIS1'] self.nrow_sc= h['NAXIS2'] if hdu.name=='IMAGING_DATA_FT' : h = hdu.header self.nwave_ft= 5 # Polarization setup if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False self.nregion_ft= 24 else: self.polarsplit = True self.nregion_ft = 48 # Image of the badpix map of the SC detector self.badpiximage_sc = gravi_badpix['IMAGING_DATA_SC'].data.squeeze() # squeeze to remove the spurious dimension self.badpix_exptime_sc = 0.0 badpixtable_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') badpix_ft = gravi_badpix['IMAGING_DATA_FT'].data if 'PIX' in badpix_ft.columns.names: print ('New format for BADPIX_FT (PIX)') for output in range(0,self.nregion_ft): ystart = output/24 badpixtable_ft[output,:] = badpix_ft.field('PIX')[0,ystart:ystart+self.nwave_ft,output%24] else: for output in range(1,self.nregion_ft+1): badpixtable_ft[output-1,:] = badpix_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft) self.badpixtable_ft = badpixtable_ft gravi_badpix.close() class Profile: def __init__(self, filename): # open fits file self.filename = filename gravi_profile=pyfits.open(filename, mode='readonly') # Main header self.header = gravi_profile[0].header if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False self.nchannels = 24 else: self.polarsplit = True self.nchannels = 48 imgaging_det_sc = gravi_profile['IMAGING_DETECTOR_SC'].data self.regname_sc = imgaging_det_sc.field('REGNAME') self.center_sc = (imgaging_det_sc.field('CENTER')) imgaging_det_ft = gravi_profile['IMAGING_DETECTOR_FT'].data self.regname_ft = imgaging_det_ft.field('REGNAME') self.center_ft = (imgaging_det_ft.field('CENTER')) # get the PROFILE_DATA useful keywords h = gravi_profile['PROFILE_DATA'].header self.startx = h['STARTX'] if 'STARTX' in h else h['HIERARCH ESO PRO PROFILE STARTX'] self.nx = len(gravi_profile['PROFILE_DATA'].data.field('DATA1')) del h self.imaging_data_sc = np.squeeze(gravi_profile['IMAGING_DATA_SC'].data) self.imaging_data_ft = np.squeeze(gravi_profile['IMAGING_DATA_FT'].data.field('PIX')) # DOES NOT WORK ON MY MAC FOR UNKNOWN REASON # prof = gravi_profile['PROFILE_DATA'].data # print ' *** Reading profiles...' # self.profile_data = np.squeeze(prof['DATA1']) # #print profile_data.shape # for i in range(1,self.nchannels): # self.profile_data += np.squeeze(prof['DATA'+str(i+1)]) # self.profile_data = [] #np.zeros((self.nchannels)) # prof = gravi_profile['PROFILE_DATA'].data # print ' *** Reading profiles...' # for i in range(0,self.nchannels): # self.profile_data.append(np.squeeze(prof['DATA'+str(i+1)])) # print ' *** Reading profiles: OK !' gravi_profile.close() del gravi_profile class Disp: def __init__(self, filename): # open fits file self.filename = filename gravi_disp=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_disp[0].header self.lbdmet = 1.908287 # DISP_MODEL table self.lin_fddl_sc = gravi_disp['DISP_MODEL'].data.field('LIN_FDDL_SC') self.lin_fddl_ft = gravi_disp['DISP_MODEL'].data.field('LIN_FDDL_FT') self.nmean = gravi_disp['DISP_MODEL'].data.field('N_MEAN') self.ndiff = gravi_disp['DISP_MODEL'].data.field('N_DIFF') self.wave0 = gravi_disp['DISP_MODEL'].data.field('WAVE0')[0] * 1e6 self.beta = gravi_disp['DISP_MODEL'].data.field('BETA') self.gamma = gravi_disp['DISP_MODEL'].data.field('GAMMA') # DISP versus WAVE_TH in [um] self.nmean_waveth = gravi_disp['DISP_WAVETH'].data.field('N_MEAN').T self.ndiff_waveth = gravi_disp['DISP_WAVETH'].data.field('N_DIFF').T self.beta_waveth = gravi_disp['DISP_WAVETH'].data.field('BETA').T self.gamma_waveth = gravi_disp['DISP_WAVETH'].data.field('GAMMA').T self.waveth = gravi_disp['DISP_WAVETH'].data.field('WAVE_TH') * 1e6 # DISP versus WAVE in [um] self.beta_wave = gravi_disp['DISP_WAVE'].data.field('BETA').T self.gamma_wave = gravi_disp['DISP_WAVE'].data.field('GAMMA').T self.wave = gravi_disp['DISP_WAVE'].data.field('EFF_WAVE') * 1e6 gravi_disp.close() class Dispvis: def __init__(self, filename): # open fits file self.filename = filename hdulist=pyfits.open(filename, mode='readonly') # Main FITS header self.header = hdulist[0].header # Get if SPLIT (11) or COMBINED (10) pola = self.header['HIERARCH ESO INS POLA MODE'] print(('INS.POLA.MODE = '+pola)) extver = 11 if pola == 'SPLIT' else 10 # Get sizes self.nobs = hdulist['OI_VIS',extver].data.field('VISDATA').shape[0]/nbase self.nlbd = hdulist['OI_VIS',extver].data.field('VISDATA').shape[1] # Get data self.visphase = hdulist['OI_VIS',extver].data.field('PHASE').reshape(self.nobs,nbase) self.visphi = hdulist['OI_VIS',extver].data.field('VISPHI').reshape(self.nobs,nbase,self.nlbd) self.visamp = hdulist['OI_VIS',extver].data.field('VISAMP').reshape(self.nobs,nbase,self.nlbd) self.visdata = hdulist['OI_VIS',extver].data.field('VISDATA').reshape(self.nobs,nbase,self.nlbd) self.gdelay = hdulist['OI_VIS',extver].data.field('GDELAY').reshape(self.nobs,nbase) # OPD_MET_FC per base met_tel = hdulist['OI_FLUX',extver].data.field('OPD_MET_FC').reshape(self.nobs,ntel) self.opd_met = met_tel[:,tel_list[:,0]] - met_tel[:,tel_list[:,1]] # Lock date per beam self.lkdt_met = hdulist['OI_FLUX',extver].data.field('LKDT_MET_FC').reshape(self.nobs,ntel) # Dispersion of GDELAY self.gd_coeff = np.median(np.diff (self.gdelay, axis=0) / np.diff (self.opd_met, axis=0)) print(("gd_coeff = %.5f"%self.gd_coeff)) # Difference from expected GD delta = (self.gdelay - self.opd_met * self.gd_coeff) self.gdelay_delta = delta - np.mean (delta, axis=0)[None,:] class Dark: def __init__(self, filename): # open fits file self.filename = filename gravi_dark=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_dark[0].header # Polarization setup if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False self.nregion_ft= 24 else: self.polarsplit = True self.nregion_ft = 48 for hdu in gravi_dark : if hdu.name=='IMAGING_DATA_SC' : h = hdu.header self.ncol_sc= h['NAXIS1'] self.ncol_sc= h['NAXIS2'] if hdu.name=='IMAGING_DATA_FT' : h = hdu.header self.nwave_ft= 5 # Image of the dark of the SC detector self.darkimage_sc = gravi_dark['IMAGING_DATA_SC'].data.squeeze() # squeeze to remove the spurious dimension self.darkerror_sc = (gravi_dark['IMAGING_ERR_SC'].data.squeeze()) self.dark_exptime_sc = 0.0 # Image of the dark of the FT detector dark_ft = gravi_dark['IMAGING_DATA_FT'].data darktable_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') if 'PIX' in dark_ft.columns.names: print ('New format for BADPIX_FT (PIX)') for output in range(0,self.nregion_ft): ystart = output/24 darktable_ft[output,:] = dark_ft.field('PIX')[0,ystart:ystart+5,output%24] else: for output in range(1,self.nregion_ft+1): darktable_ft[output-1,:] = dark_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft) self.darktable_ft = np.copy(darktable_ft) # Image of the dark stdev of the FT detector dark_ft = gravi_dark['IMAGING_ERR_FT'].data darkerror_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') if 'PIX' in dark_ft.columns.names: print ('New format for BADPIX_FT (PIX)') for output in range(0,self.nregion_ft): ystart = output/24 darkerror_ft[output,:] = dark_ft.field('PIX')[0,ystart:ystart+5,output%24] else: for output in range(1,self.nregion_ft+1): darkerror_ft[output-1,:] = (np.abs(dark_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft))) self.darkerror_ft = np.copy(darkerror_ft) gravi_dark.close() class PhaseFTSC: def __init__(self, filename): # open fits file self.filename = filename gravi_phaseftsc=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_phaseftsc[0].header spectrum_sc = None spectrum_ft = None detector_sc = None detector_ft = None # get the FDDL, OI_ARRAY, OI_WAVELENGTH and OI_VIS_MET extentions for hdu in gravi_phaseftsc : if hdu.name == 'SPECTRUM_DATA_FT' : spectrum_ft=hdu.data if hdu.name == 'SPECTRUM_DATA_SC' : spectrum_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_SC' : detector_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_FT' : detector_ft=hdu.data if hdu.name == 'FDDL': fddl = hdu.data self.fddlpresent = True # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # # Pipeline FT and SC average phase per baseline # ext7 = gravi_phaseftsc[7].data # ext8 = gravi_phaseftsc[8].data # self.met_phase_sc = ext7.field('PHASE') # baseline, time_sc (expressed in radians) # self.met_phase_ft = ext8.field('PHASE') # baseline, time_ft (expressed in radians) # Metrology OPL SC-FT tables per telescope from pipeline ext11 = gravi_phaseftsc[11].data self.time_met = ext11.field('TIME')/1E6 # in seconds self.nframe_met = self.time_met.shape[0] self.met_dopl = np.zeros((4,self.nframe_met)) for tel in range(0,4): self.met_dopl[tel,:] = ext11.field('OPL'+str(tel+1))*1e6 # tel, time_met IN MICRONS # IMAGING_DETECTOR parameters if detector_sc is not None: self.regname_sc = detector_sc.field('REGNAME') self.ports_sc = detector_sc.field('PORTS') else: self.regname_sc = '' if detector_ft is not None: self.regname_ft = detector_ft.field('REGNAME') self.ports_ft = detector_ft.field('PORTS') else: self.regname_ft = '' self.nregion = self.regname_sc.shape[0] self.nwave_sc = spectrum_sc.field('DATA1').shape[1] self.nwave_ft = spectrum_ft.field('DATA1').shape[1] # Extracted non resampled spectra if spectrum_sc is not None: self.spectrum_sc = [] # initialization self.time_sc = spectrum_sc.field('TIME')/1E6 # in seconds self.spectrum_sc = np.zeros((self.nregion, self.time_sc.shape[0], self.nwave_sc)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_sc[i-1,:,:] = np.mean(spectrum_sc.field('DATA'+str(i)), axis=2) if spectrum_ft is not None: self.spectrum_ft = [] # initialization self.time_ft = spectrum_ft.field('TIME')/1E6 # in seconds self.spectrum_ft = np.zeros((self.nregion, self.time_ft.shape[0], self.nwave_ft)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_ft[i-1,:,:] = np.mean(spectrum_ft.field('DATA'+str(i)), axis=2) # Number of frames if 'HIERARCH ESO DET2 NDIT' in self.header: self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] else: self.nframe_sc = 0 self.nframe_ft = self.spectrum_ft.shape[1] gravi_phaseftsc.close() class Argon: def __init__(self, filename): # open fits file self.filename = filename gravi_argon = pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_argon[0].header self.frames = np.array(gravi_argon['IMAGING_DATA_SC'].data) class P2vmreducedAstroList: def __init__(self, filelist, fits_keywords): self.headers = [] self.filename = [] nfiles_in = len(filelist) nfiles = 0 for filename in filelist: gravi_file = pyfits.open(filename+".fits",memmap=True) header = gravi_file[0].header.copy() gravi_file.close() del gravi_file # Final OIFITS data files only key_names = list(fits_keywords.keys()) type_ok = True for strname in key_names: type_ok *= (strname in header) if type_ok == False: continue keys_ok = True for name in key_names: keys_ok *= (header[name] in fits_keywords[name]) if keys_ok == False: continue print(' Loading P2VMREDUCED file %i/%i: %s.fits'%(nfiles,nfiles_in,filename)) hdulist = pyfits.open(filename+".fits") self.headers.append(hdulist[0].header) self.filename.append(filename+".fits") # 10 = GRAVITY_SC, 11 = GRAVITY_SC_P1, 12 = GRAVITY_SC_P2 # 20 = GRAVITY_FT, 21 = GRAVITY_FT_P1, 22 = GRAVITY_FT_P2 extv = 11 # Append the data of the new file to the already loaded tables # FILEID is a new column to store the mapping with headers if nfiles == 0: nrows = hdulist['OI_VIS',extv].data.shape[0] cols = pyfits.ColDefs([pyfits.Column(name='FILEID', format='I',array=np.zeros(nrows))]) oi_vis = pyfits.BinTableHDU.from_columns(hdulist['OI_VIS',extv].columns + cols) else: nrows1 = oi_vis.data.shape[0] nrows2 = hdulist['OI_VIS',extv].data.shape[0] oi_vis = pyfits.BinTableHDU.from_columns(oi_vis.columns, nrows=nrows1+nrows2) for colname in hdulist['OI_VIS',extv].columns.names: oi_vis.data[colname][nrows1:] = hdulist['OI_VIS',extv].data[colname] oi_vis.data['FILEID'][nrows1:] = nfiles nfiles = nfiles + 1 # hdulist.close() ideally we should find a way to close the previous files # End loop on files # Set this big OI_VIS to the class self.oi_vis = oi_vis # # Common plotting system #