#! /usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Sep 14 11:41:26 2015 @author: kervella """ from astropy.io import fits 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.graphics.charts.axes import LogYValueAxis 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, ParagraphStyle 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 try: fluxdata=flux.field('FLUX') except KeyError: fluxdata=flux.field('FLUXDATA') for base in range(nbase): f1f2 = fluxdata[tel_list[base][0]::ntel,:] * fluxdata[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,log_y=False): 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)" if log_y: lp.yValueAxis = LogYValueAxis() lp.yValueAxis.labelTextFormat = '%.1e' 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,log_y=False): formx = '%.2g' if formx is None else formx formy = '%.2g' 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)" if log_y: lp.yValueAxis = LogYValueAxis() lp.yValueAxis.labelTextFormat = '%.1e' lp.xValueAxis.visibleGrid = 0 lp.yValueAxis.visibleGrid = 0 else: lp.xValueAxis.visibleGrid = 1 lp.yValueAxis.visibleGrid = 1 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 # 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,log_y=False): 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 if log_y: lp.yValueAxis = LogYValueAxis() lp.yValueAxis.labelTextFormat = '%.1e' lp.xValueAxis.visibleGrid = 0 lp.yValueAxis.visibleGrid = 0 else: lp.xValueAxis.visibleGrid = 1 lp.yValueAxis.visibleGrid = 1 lp.yValueAxis.valueMin = ymin lp.yValueAxis.valueMax = ymax lp.yValueAxis.valueStep = ticky lp.yValueAxis.visibleTicks = 0 lp.yValueAxis.visibleLabels = 0 # 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) def plotSubtitleSpace(Story,title): subStyle=ParagraphStyle( name = 'subtitle', parent=styles["Normal"], spaceBefore=5*mm, spaceAfter=3*mm, ) Story.append(Paragraph(title, style = subStyle)) # 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')) sobjoffx = str(get_key_withdefault(product.header,'HIERARCH ESO INS SOBJ OFFX','0.0')) sobjoffy = str(get_key_withdefault(product.header,'HIERARCH ESO INS SOBJ OFFY','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 offx offy", "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 +' '+ sobjoffx +' '+ sobjoffy, '%(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") plotSubtitleSpace(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 plotSubtitleSpace(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 plotSubtitleSpace(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 range(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, Astroreduced, Rawdata, Badpix, Profile, Dark, Argon # class Oifits: def __init__(self, filename): # open fits file self.filename = filename gravi_oifits=fits.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 'VISFLUX' in columnflux: FLUXDATA = 'VISFLUX' FLUXERR = 'VISFLUXERR' elif 'FLUX' in columnflux: FLUXDATA = 'FLUX' FLUXERR = 'FLUXERR' elif 'FLUXDATA' in columnflux: FLUXDATA = 'FLUXDATA' FLUXERR = 'FLUXERR' else: raise Exception('no FLUXDATA column in OI_FLUX') if self.polarsplit == False: self.oi_flux_sc = oi_flux_sc.field(FLUXDATA).reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field(FLUXDATA).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: self.oi_flux_sc_s = oi_flux_sc_s.field(FLUXDATA).reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field(FLUXDATA).reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field(FLUXDATA).reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field(FLUXDATA).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 = fits.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 # ======================== # TODO: Astroreduced class # ======================== class Astroreduced: def __init__(self, filename): # open fits file self.filename = filename gravi_astroreduced = fits.open(filename, mode='readonly') nbase = 6 # Main FITS header self.header = gravi_astroreduced[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=True self.vfactor_in=True # get the extentions for hdu in gravi_astroreduced : if hdu.name == 'OI_ARRAY' : oi_array=hdu.data if hdu.name == 'OI_VIS_MET' : oi_vis_met = hdu.data metcolumns = hdu.columns.names if hdu.name == 'IMAGING_DATA_ACQ' : oi_imaging_data_acq=hdu.data # ACQ CAM if 'OI_VIS_ACQ' in gravi_astroreduced: print (' OI_VIS_ACQ is in') self.acq_in = True data = gravi_astroreduced['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() self.oi_acq_pup_u = data.field('PUPIL_U').copy() self.oi_acq_pup_v = data.field('PUPIL_V').copy() self.oi_acq_pup_w = data.field('PUPIL_W').copy() self.oi_acq_opd_pup = data.field('OPD_PUPIL').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 extensions in combined polarization mode if self.polarsplit == False: for hdu in gravi_astroreduced : header = hdu.header.copy() if 'INSNAME' in header: if (hdu.name == 'OI_WAVELENGTH') & ('_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('_SC' in str(hdu.header['INSNAME'])): oi_vis_sc=hdu.data.copy() column_oi_vis_sc = hdu.columns.names 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 # get the OI_VIS, OI_FLUX extensions in split polarization mode (P considered =P1, S considered =P2) if self.polarsplit == True: for hdu in gravi_astroreduced : header = hdu.header.copy() if 'INSNAME' in header: # The two polarizations have the same wavelength scale if (hdu.name == 'OI_WAVELENGTH') & ('_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc = hdu.data.copy() # if (hdu.name == 'OI_WAVELENGTH') & ('_SC_P2' in str(hdu.header['INSNAME'])): # oi_wave_sc_s=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_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() # Wavelength scale 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] # 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') # OI_VIS_MET parameters if oi_vis_met is not None: 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) if 'OPD_TELFC_CORR_XY' in metcolumns: self.met_opd_telfc_corr_xy = oi_vis_met.field('OPD_TELFC_CORR_XY').reshape((nrow,16)) if 'OPD_FC_CORR' in metcolumns: self.met_opd_fc_corr = oi_vis_met.field('OPD_FC_CORR').reshape((nrow,4)) # OI_FLUX table parameters if 'VISFLUX' in columnflux: FLUXDATA = 'VISFLUX' FLUXERR = 'VISFLUXERR' elif 'FLUX' in columnflux: FLUXDATA = 'FLUX' FLUXERR = 'FLUXERR' elif 'FLUXDATA' in columnflux: FLUXDATA = 'FLUXDATA' FLUXERR = 'FLUXERR' else: raise Exception('no FLUXDATA column in OI_FLUX') 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_sc = oi_flux_sc.field(FLUXDATA).reshape(self.oi_flux_sc_time.shape[0],4,self.nwave_sc) if 'TOTALFLUX_SC' in columnflux: self.totalflux_sc = oi_flux_sc.field('TOTALFLUX_SC').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_sc_time_p = oi_flux_sc_p.field('TIME')[::4]/1e6 self.oi_flux_sc_s = oi_flux_sc_s.field(FLUXDATA).reshape(self.oi_flux_sc_time_s.shape[0],4,self.nwave_sc) self.oi_flux_sc_p = oi_flux_sc_p.field(FLUXDATA).reshape(self.oi_flux_sc_time_p.shape[0],4,self.nwave_sc) 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_sc_p = oi_flux_sc_p.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time_p.shape[0],4) # Number of frames self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] if oi_vis_met is not None: 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') # 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,:]) if self.gdelay_in == True: # Group delay in microns self.gdelay_sc = oi_vis_sc.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 # Signal to noise of fringes self.snr_sc = oi_vis_sc.field('SNR').reshape(self.nframe_sc,nbase).T if ('V_FACTOR' in column_oi_vis_sc): print("Load VFACTOR data") self.vfactor_in = True 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') 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_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 # 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_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 if 'V_FACTOR' in column_oi_vis_sc_p: self.vfactor_in = True self.oi_vis_vfactor_s = oi_vis_sc_s.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) 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,:]) gravi_astroreduced.close() # ============================ # TODO: end Astroreduced class # ============================ class P2vm: def __init__(self, filename): # open fits file self.filename = filename gravi_p2vm=fits.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'] try : 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) except : 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 = fits.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=fits.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=fits.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=fits.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: try: phase_fc=oi_vis_met.field('PHASE_FC_TAC') except: phase_fc=oi_vis_met.field('PHASE_FC') try: phase_tel=oi_vis_met.field('PHASE_TEL_DRS') except: phase_tel=oi_vis_met.field('PHASE_TEL') if len(phase_fc.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) self.oi_vis_met_phase = phase_fc.reshape((nrow,4)) self.oi_vis_met_phase_tel = 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 'VISFLUX' in columnflux: FLUXDATA = 'VISFLUX' FLUXERR = 'VISFLUXERR' elif 'FLUX' in columnflux: FLUXDATA = 'FLUX' FLUXERR = 'FLUXERR' elif 'FLUXDATA' in columnflux: FLUXDATA = 'FLUXDATA' FLUXERR = 'FLUXERR' else: raise Exception('no FLUXDATA column in OI_FLUX') 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 self.oi_flux_sc = oi_flux_sc.field(FLUXDATA).reshape(self.oi_flux_sc_time.shape[0],4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field(FLUXDATA).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 self.oi_flux_sc_s = oi_flux_sc_s.field(FLUXDATA).reshape(self.oi_flux_sc_time_s.shape[0],4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field(FLUXDATA).reshape(self.oi_flux_ft_time_s.shape[0],4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field(FLUXDATA).reshape(self.oi_flux_sc_time_p.shape[0],4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field(FLUXDATA).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 = fits.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 != 21): 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=fits.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 = int(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=fits.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=fits.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=fits.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=fits.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 = int(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 = int(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=fits.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 = fits.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 = fits.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 = fits.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 = fits.ColDefs([fits.Column(name='FILEID', format='I',array=np.zeros(nrows))]) oi_vis = fits.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 = fits.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 #