#!/usr/bin/env python3 import warnings warnings.filterwarnings('ignore') import platform, sys, os if platform.uname()[1] == 'wvgoff': # -- VLTI Offline Machine on Paranal sys.path.append('/diska/home/astrov/vltipso/python27/anaconda') sys.path = sys.path[::-1] import numpy as np warnings.simplefilter('ignore', np.RankWarning) import matplotlib matplotlib.use('TkAgg') matplotlib.rcParams['toolbar'] = 'toolbar2' from matplotlib import pyplot as plt plt.ion() from astropy.io import fits import pickle import tkinter, tkinter.filedialog, tkinter.messagebox, tkinter.font import webbrowser import time import scipy.special import scipy.signal # -- Obsolete: __correctP2VMFlux = False #if __correctP2VMFlux: # import p2vmCont Vmodel_info = """Coma-separated destription of the visibility model, all fluxes assume f=1 (primary) - 'ud=1.0' uniform disk diameter, in mas - 'f=1.0' continuum flux of the primary (default is 1.0) - 'f_2.166_0.002=0.1' gaussian emission (centered 2.166um, fwhm=0.002um) - 'fres=0.1' continuum resolved flux - 'fres_2.166_0.002=0.1' gaussian spectral emission (centered 2.166um, fwhm=0.002um) of resolved flux - 'udc=1.0, xc=1.0, yc=1.2, fc=0.1' uniform disk diameter of a companion, at -1.0mas towards east and 1.2 mas towards north. companion flux is 0.1 - 'fc_2.166_0.002=0.1': gaussian line (centered 2.166um, fwhm=0.002um) for the companion flux gaussian feature: 'dg' (size) 'xg', 'yg' (position), 'fg' (flux) and 'fg_...' for emmission / absorption lines. example: ud=0.4, f_2.1665_0.0008=0.3, f=1.0, xc=1.0, yc=1.25, fc_2.1665_0.0008=0.8, fc=0.05, udc=0.8, fres=0, fres_2.090_0.001=0.15 """ def clickVisModel(): toplevel = tkinter.Toplevel() label1 = tkinter.Label(toplevel, text=Vmodel_info, height=0, width=100, anchor='w') label1.pack() graviqlmod = '.graviqlmod' def openUrl(url): webbrowser.open_new(url) def loadGraviMulti(filenames, insname='GRAVITY_SC', wlmin=None, wlmax=None): data = [loadGravi(f, insname=insname) for f in filenames] # check if all observations done for sma object in spectro same mode modeObj = set(['-'.join([d['TARG NAME'],d['SPEC RES'], d['POLA'],d['BASELINE']]) for d in data]) if len(modeObj)!=1: #print modeObj return None # -- averages, stored in first file for o in ['VIS', 'V2', 'T3', 'VISPHI']: for k in list(data[0][o].keys()): # -- for each baseline for i in range(len(data)): # -- for each data file if i==0: flag = np.isnan(data[i][o][k]) data[i][o][k][flag] = 0.0 n = np.ones(len(flag)) n[flag] = 0.0 else: flag = np.isnan(data[i][o][k]) data[0][o][k][~flag] += data[i][o][k][~flag] n[~flag] += 1 data[0][o][k] /= n for o in ['uV2', 'vV2', 'u1T3', 'v1T3', 'u2T3', 'v2T3', 'uVIS', 'vVIS']: for k in list(data[0][o].keys()): # -- for each baseline for i in range(len(data)): # -- for each data file if i==0: data[i][o][k] = data[i][o][k]/len(data) else: data[0][o][k] += data[i][o][k]/len(data) # -- averaging differential phase, not phase if wlmin is None or wlmindata[0]['wl'].max(): wlmax = data[0]['wl'].max() if (wlmax-wlmin)=wlmin)* (np.abs(data[0]['wl'] - 0.5*(wlmin+wlmax)) > 0.33*(wlmax-wlmin))) else: wc = np.where(data[0]['wl']) for o in ['VISPHI']: # for each baseline for k in list(data[0][o].keys()): # for each baseline for i in range(len(data)): #print '->', data[0]['wl'].shape, tmp.shape cc = nanpolyfit(data[0]['wl'][wc]-0.5*(wlmin+wlmax), data[i][o][k][wc], 1) #print filenames[i], cc, tmp.mean() if i==0: data[i][o][k] = (data[i][o][k] - np.polyval(cc, data[0]['wl'] - 0.5*(wlmin+wlmax)))/float(len(data)) else: data[0][o][k] += (data[i][o][k] - np.polyval(cc, data[0]['wl'] - 0.5*(wlmin+wlmax)))/float(len(data)) return data[0] def loadGravi(filename, insname='GRAVITY_SC'): f = fits.open(filename) res = {} insnames = [] #print filename if 'SC' in insname: res['DIAM'] = f[0].header['ESO INS SOBJ DIAMETER'] else: res['DIAM'] = f[0].header['ESO FT ROBJ DIAMETER'] res['TARG NAME'] = f[0].header['ESO INS SOBJ NAME'] res['OB NAME'] = f[0].header['ESO OBS NAME'] res['PI'] = f[0].header['ESO OBS PI-COI NAME'] res['OBS ID'] = f[0].header['ESO OBS ID'] res['PROG ID'] = f[0].header['ESO OBS PROG ID'] res['SPEC RES'] = f[0].header['ESO INS SPEC RES'] res['POLA'] = f[0].header['ESO INS POLA MODE'] res['DIT'] = f[0].header['ESO DET2 SEQ1 DIT'] res['BASELINE']= '-'.join([f[0].header['ESO ISS CONF STATION%d'%i] for i in [1,2,3,4]]) res['PIPELINE'] = f[0].header['ESO PRO REC1 PIPE ID'] if insname=='auto_SC': if ' 0.8' in res['PIPELINE'] or\ ' 0.7' in res['PIPELINE']: insname = 'SPECTRO_SC'+('' if res['POLA']=='COMBINED' else '_P') else: insname = 'GRAVITY_SC'+('' if res['POLA']=='COMBINED' else '_P') elif insname=='auto_FT': if ' 0.8' in res['PIPELINE'] or\ ' 0.7' in res['PIPELINE']: insname = 'SPECTRO_FT'+('' if res['POLA']=='COMBINED' else '_P') else: insname = 'GRAVITY_FT'+('' if res['POLA']=='COMBINED' else '_P') res['INSNAME'] = insname # -- wavelength: ---- for h in f: if 'EXTNAME' in list(h.header.keys()) and \ h.header['EXTNAME'] == 'OI_WAVELENGTH': if insname in h.header['INSNAME']: res['wl'] = h.data['EFF_WAVE']*1e6 insnames.append(h.header['INSNAME']) if not 'wl' in list(res.keys()): print('ERROR: could no find INSNAME="'+insname+'" within:', insnames) return # -- oi array: ---- oiarray = dict(list(zip(f['OI_ARRAY'].data['STA_INDEX'], np.char.strip(f['OI_ARRAY'].data['STA_NAME'])))) # -- V2: ---- res['V2'] = None n = {'all':0.0} for h in f: if 'EXTNAME' in list(h.header.keys()) and \ h.header['EXTNAME'] == 'OI_VIS2': if insname in h.header['INSNAME']: if res['V2'] is None: res['V2'] = {} res['uV2'] = {} res['vV2'] = {} for i in range(6): k = oiarray[h.data['STA_INDEX'][i][0]]+\ oiarray[h.data['STA_INDEX'][i][1]] flag = h.data['FLAG'][i] + np.isnan(h.data['VIS2DATA'][i]) res['V2'][k] = h.data['VIS2DATA'][i].copy() res['V2'][k][flag] = 0.0 res['uV2'][k] = h.data['UCOORD'][i].copy() res['vV2'][k] = h.data['VCOORD'][i].copy() n[k] = np.ones(len(res['V2'][k])) n[k][flag] = 0.0 n['all'] += 1.0 else: for i in range(6): k = oiarray[h.data['STA_INDEX'][i][0]]+\ oiarray[h.data['STA_INDEX'][i][1]] flag = h.data['FLAG'][i] + np.isnan(h.data['VIS2DATA'][i]) res['V2'][k][~flag] += h.data['VIS2DATA'][i][~flag] res['uV2'][k] += h.data['UCOORD'][i].copy() res['vV2'][k] += h.data['VCOORD'][i].copy() n[k][~flag] += 1.0 n['all'] += 1.0 for k in list(res['V2'].keys()): res['V2'][k] /= n[k] res['V2'][k][n[k]==0] = np.nan res['uV2'][k] /= float(n['all']) res['vV2'][k] /= float(n['all']) # -- T3: ---- res['T3'] = None n = {'all':0.0} # average counter for h in f: if 'EXTNAME' in list(h.header.keys()) and \ h.header['EXTNAME'] == 'OI_T3': if insname in h.header['INSNAME']: if res['T3']==None: res['T3'] = {} res['u1T3'],res['v1T3'], res['u2T3'],res['v2T3'] = {}, {}, {}, {} for i in range(4): k = oiarray[h.data['STA_INDEX'][i][0]]+\ oiarray[h.data['STA_INDEX'][i][1]]+\ oiarray[h.data['STA_INDEX'][i][2]] res['T3'][k] = h.data['T3PHI'][i].copy() flag = h.data['FLAG'][i] + np.isnan(res['T3'][k]) res['T3'][k][flag] = 0.0 # no data res['u1T3'][k] = h.data['U1COORD'][i].copy() res['v1T3'][k] = h.data['V1COORD'][i].copy() res['u2T3'][k] = h.data['U2COORD'][i].copy() res['v2T3'][k] = h.data['V2COORD'][i].copy() n[k] = np.ones(len(res['T3'][k])) n[k][flag] = 0.0 # no part of average n['all'] += 1.0 else: for i in range(4): k = oiarray[h.data['STA_INDEX'][i][0]]+\ oiarray[h.data['STA_INDEX'][i][1]]+\ oiarray[h.data['STA_INDEX'][i][2]] flag = h.data['FLAG'][i] + np.isnan(res['T3'][k]) res['T3'][k][~flag] += h.data['T3PHI'][i][~flag] res['u1T3'][k] += h.data['U1COORD'][i].copy() res['v1T3'][k] += h.data['V1COORD'][i].copy() res['u2T3'][k] += h.data['U2COORD'][i].copy() res['v2T3'][k] += h.data['V2COORD'][i].copy() n[k][~flag] += 1.0 n['all'] += 1.0 for k in list(res['T3'].keys()): res['T3'][k] /= n[k] res['T3'][k][n[k]==0] = np.nan res['T3'][k] = (res['T3'][k]+180)%360 - 180 # TODO: unwrap res['u1T3'][k] /= n['all'] res['v1T3'][k] /= n['all'] res['u2T3'][k] /= n['all'] res['v2T3'][k] /= n['all'] # -- diff Phase --- res['VISPHI'] = None res['VIS'] = None n = {'all':0.0} for h in f: if 'EXTNAME' in list(h.header.keys()) and \ h.header['EXTNAME'] == 'OI_VIS': if insname in h.header['INSNAME']: if res['VISPHI'] is None: res['VIS'] = {} res['VISPHI'] = {} res['uVIS'] = {} res['vVIS'] = {} for i in range(6): k = oiarray[h.data['STA_INDEX'][i][0]]+\ oiarray[h.data['STA_INDEX'][i][1]] res['VISPHI'][k] = h.data['VISPHI'][i].copy() res['VIS'][k] = h.data['VISAMP'][i].copy() flag = h.data['FLAG'][i] + np.isnan(h.data['VISPHI'][i]) res['VISPHI'][k][flag] = 0.0 res['VIS'][k][flag] = 0.0 res['uVIS'][k] = h.data['UCOORD'][i].copy() res['vVIS'][k] = h.data['VCOORD'][i].copy() n[k] = np.ones(len(res['VISPHI'][k])) n[k][flag] = 0.0 n['all'] += 1 else: for i in range(6): k = oiarray[h.data['STA_INDEX'][i][0]]+\ oiarray[h.data['STA_INDEX'][i][1]] flag = h.data['FLAG'][i] + np.isnan(h.data['VISPHI'][i]) res['VISPHI'][k][~flag] += h.data['VISPHI'][i][~flag] res['VIS'][k][~flag] += h.data['VISAMP'][i][~flag] res['uVIS'][k] += h.data['UCOORD'][i] res['vVIS'][k] += h.data['VCOORD'][i] n[k][~flag] += 1.0 n['all'] += 1 for k in list(res['VISPHI'].keys()): res['VIS'][k] /= n[k] res['VISPHI'][k] /= n[k] res['VISPHI'][k][n[k]==0] = np.nan res['VIS'][k][n[k]==0] = np.nan res['VISPHI'][k] = (res['VISPHI'][k]+180)%360 - 180 # TODO: unwrap res['uVIS'][k] /= float(n['all']) res['vVIS'][k] /= float(n['all']) # -- Spctr: ---- res['FLUX'] = {} n = 0.0 for h in f: if 'EXTNAME' in list(h.header.keys()) and \ h.header['EXTNAME'] == 'OI_FLUX': if insname in h.header['INSNAME']: if len(res['FLUX'])==0: for i in range(4): k = oiarray[h.data['STA_INDEX'][i]] if __correctP2VMFlux: cont = np.interp(res['wl'], p2vmCont.cc[res['SPEC RES'], res['POLA']]['WL'], p2vmCont.cc[res['SPEC RES'], res['POLA']]['T%d'%(i+1)]) else: cont = 1. res['FLUX'][k] = h.data['FLUX'][i].copy()/cont n += 1 #print res['FLUX'].keys() else: for i in range(4): k = oiarray[h.data['STA_INDEX'][i]] if __correctP2VMFlux: cont = np.interp(res['wl'], p2vmCont.cc[res['SPEC RES'], res['POLA']]['WL'], p2vmCont.cc[res['SPEC RES'], res['POLA']]['T%d'%(i+1)]) else: cont = 1. res['FLUX'][k] += h.data['FLUX'][i]/cont n += 1 for k in list(res['FLUX'].keys()): res['FLUX'][k] /= n f.close() return res def slidingOp(x,y,dx): """ returns 25%, median, 75% """ res = {} for k in ['25%', '75%', 'median', '1sigma', 'mean', '90%']: res[k] = np.zeros(len(x)) for i in range(len(x)): yp = y[np.abs(x-x[i])<=dx/2] #res['25%'][i] = np.percentile(yp, 25) res['median'][i] = np.percentile(yp, 50) #res['mean'][i] = np.mean(yp) res['75%'][i] = np.percentile(yp, 75) res['90%'][i] = np.percentile(yp, 90) res['1sigma'][i] = (np.percentile(yp, 84)-np.percentile(yp, 16))/2 return res __lbda, __tran = None, None def tellTrans(wl, width=2.0): global __lbda, __tran if __lbda is None: wv = '020' # 2mm Water Vapor if os.path.exists('telluric_'+wv+'.pickle'): f = open('telluric_'+wv+'.pickle', 'rb') __lbda, __tran = pickle.load(f) f.close() else: f = fits.open('/Users/amerand/DATA/sky_transmission/'+'transnir'+wv+'.fits') __lbda = f[1].data['WAVELENGTH'][0].copy() __tran = f[1].data['TRANSMISSION'][0].copy() f.close() __lbda = __lbda[::150] __tran = __tran[::150] __tran = __tran[(__lbda<2.5)*(__lbda>1.9)] __lbda = __lbda[(__lbda<2.5)*(__lbda>1.9)] f = fits.open('/Users/amerand/DATA/sky_transmission/'+'emisnir'+wv+'.fits') __l = f[1].data['WAVELENGTH'][0].copy() __e = f[1].data['EMISSION'][0].copy() f.close() __l = __l[::150] __e = __e[::150] __tran += np.interp(__lbda, __l, __e) f = open('telluric_'+wv+'.pickle', 'wb') pickle.dump((__lbda, __tran), f) f.close() res = [] gwl = np.gradient(wl) for i in range(len(wl)): res.append(np.mean(__tran[ np.abs(__lbda-wl[i])0: for f_ in F_: wl = float(f_.split('_')[1]) dwl = float(f_.split('_')[2]) f += params[f_]*np.exp(-4*np.log(2)*(res['wl']-wl)**2/dwl**2) else: # -- case no lines and no flux is given, but a diameter if not 'F' in list(params.keys()) and 'UD' in list(params.keys()): f = 1.0 #ck = filter(lambda x: x=='FC' or (x[:2]=='FC' and x[2:].isdigit()), params.keys()) # -- looks for flux definition # -- companion flux if not 'FC' in list(params.keys()): fc = 1.0 # assume equal companion else: fc = params['FC'] F_ = [k for k in list(params.keys()) if k.startswith('FC_')] if len(F_)>0: for f_ in F_: wl = float(f_.split('_')[1]) dwl = float(f_.split('_')[2]) fc += params[f_]*np.exp(-4*np.log(2)*(res['wl']-wl)**2/dwl**2) # -- gaussian flux if not 'FG' in list(params.keys()): fg = 0.0 # assume null flux else: fg = params['FG'] F_ = [k for k in list(params.keys()) if k.startswith('FG_')] if len(F_)>0: for f_ in F_: wl = float(f_.split('_')[1]) dwl = float(f_.split('_')[2]) fg += params[f_]*np.exp(-4*np.log(2)*(res['wl']-wl)**2/dwl**2) # -- fully resolved component flux if 'FRES' in list(params.keys()): fres = params['FRES'] else: fres = 0.0 F_ = [k for k in list(params.keys()) if k.startswith('FRES_')] if len(F_)>0: for f_ in F_: wl = float(f_.split('_')[1]) dwl = float(f_.split('_')[2]) fres += params[f_]*np.exp(-4*np.log(2)*(res['wl']-wl)**2/dwl**2) if 'V2' in list(r.keys()): for k in list(r['V2'].keys()): B = np.sqrt(r['uV2'][k]**2 + r['vV2'][k]**2) # -- primary V = f + 0j # unresolved object if 'UD' in list(params.keys()): V = f*_Vud(B, params['UD'], res['wl'])*(1.+0j) # -- secondary star if 'XC' in list(params.keys()) and 'YC' in list(params.keys()): if 'UDC' in list(params.keys()): Vc = _Vud(B, params['UDC'], res['wl'])*(1.+0j) else: Vc = 1.0+0.0j phic = -_s*2*np.pi*c*(r['uV2'][k]*params['XC']+r['vV2'][k]*params['YC'])/res['wl'] V += fc*Vc*np.exp(1j*phic) else: fc = 0.0 # -- gaussian if 'DG' in list(params.keys()): Vg = np.exp(-(np.pi*c*params['DG']*B/res['wl'])**2/(4*np.log(2))) else: Vg = 0.0 if 'XG' in list(params.keys()) and 'YG' in list(params.keys()): phig = -_s*2*np.pi*c*(r['uV2'][k]*params['XG']+r['vV2'][k]*params['YG'])/res['wl'] else: phig = 0. V += fg*Vg*np.exp(1j*phig) V /= (f + fres + fc + fg) res['V'][k] = V # -- complex vis: keep it for later res['V2'][k] = np.abs(V)**2 res['VIS'][k] = np.abs(V) res['uV2'][k] = r['uV2'][k] res['vV2'][k] = r['vV2'][k] res['uVIS'][k] = r['uV2'][k] res['vVIS'][k] = r['vV2'][k] # -- normalized Spectrum: if isinstance(f + fres + fc + fg, np.ndarray): res['FLUX'] = f + fres + fc + fg # -- crude normalization res['FLUX'] /= np.median(res['FLUX']) else: res['FLUX'] = np.ones(len(res['wl'])) if 'VISPHI' in list(r.keys()): for k in list(res['V'].keys()): res['VISPHI'][k] = 180*np.angle(res['V'][k])/np.pi for k in list(r['T3'].keys()): T = (k[0:2],k[2:4],k[4:6]) # -- three telescpoes f = [] # -- compute formula for closure phase if T[0]+T[1] in list(res['VISPHI'].keys()): f.append((T[0]+T[1], _s)) else: f.append((T[1]+T[0], - _s)) if T[1]+T[2] in list(res['VISPHI'].keys()): f.append((T[1]+T[2], _s)) else: f.append((T[2]+T[1], - _s)) if T[2]+T[0] in list(res['VISPHI'].keys()): f.append((T[2]+T[0], _s)) else: f.append((T[0]+T[2], - _s)) # -- signed sum res['T3'][k] = 0. for x in f: res['T3'][k] += x[1]*res['VISPHI'][x[0]] res['T3'][k] = (res['T3'][k]+180.)%360. - 180. return res def plotGravi(filename, insname='auto_SC', wlmin=None, wlmax=None, filt=False, onlySpectrum=False, exportFilename='', v2b=False, model='', withdPhi=False, visamp=False): #print '#'*5, exportFilename, '#'*5 top = 0.1 if isinstance(filename, list) or isinstance(filename, tuple): r = loadGraviMulti(filename, insname) if r is None: return False top = 0.06*(min(len(filename)//2, 2)) tmp = os.path.basename(filename[0]) for i,f in enumerate(filename[1:3]): if not withdPhi: tmp += '\n'+os.path.basename(f) else: tmp += ('\n' if i%2==1 else ';')+os.path.basename(f) if len(filename)>3: tmp += '...' filename = tmp else: r = loadGravi(filename, insname) top += 0.05 if r is None: return False if onlySpectrum: plt.close(2) plt.figure(2, figsize=(15/1.5,5/1.5)) elif v2b: plt.close(0) plt.figure(0, figsize=(12/1.5,8/1.5)) else: plt.close(1) if withdPhi: plt.figure(1, figsize=(15/1.5,9/1.5)) else: plt.figure(1, figsize=(11/1.5,9/1.5)) plt.clf() plt.suptitle(filename+'\n'+' | '.join([r['PROG ID'], str(r['OBS ID']), r['OB NAME'], r['INSNAME']]), fontsize=9) # plt.suptitle(' | '.join([r['PROG ID'], str(r['OBS ID']), # r['OB NAME'], r['INSNAME']]), # fontsize=9) if wlmin is None or wlminr['wl'].max(): wlmax = r['wl'].max() w = np.where((r['wl']<=wlmax)*(r['wl']>=wlmin)) r['wl band'] = r['wl'][w] # -- compute Vsibility model rm = Vmodel(r,model, r['TARG NAME']) model = rm['modelstr'] if visamp: vis = 'VIS' else: vis = 'V2' if v2b: for i,k in enumerate(r[vis].keys()): B = np.sqrt(r['u'+vis][k]**2 + r['v'+vis][k]**2) plt.plot(B/r['wl'][w], r[vis][k][w], label=k+' (%5.1fm)'%B, alpha=0.5, marker='.', linestyle='-') plt.plot(B/rm['wl'][w], rm[vis][k][w], color='k', alpha=0.5, marker='.', linestyle='-') plt.grid() plt.legend(loc='upper right', fontsize=7) plt.ylim(0.0, 1.0) plt.xlabel(r'B / M$\lambda$') plt.ylabel('V$^2$') return model # -- decide whether or not we compute differential quantities if len(w[0])=wlmin)* (np.abs(r['wl'] - 0.5*(wlmin+wlmax)) > 0.33*(wlmax-wlmin)) ) if computeDiff: cc = nanpolyfit(r['wl'][wc],tmp[wc],1) r['spectrum band'] = tmp[w] - np.polyval(cc,r['wl'][w]) + 1 plt.ylim(getYlim(tmp[w])) # -- spectral line database # -- https://arxiv.org/pdf/astro-ph/0008213.pdf, table 1 # -- https://iopscience.iop.org/article/10.1086/312501/fulltext/995828.tb2.html lines = {'HeI':[(2.0581, 2.11274, 2.11267, 2.11258, 2.161), 'c', 'dashed'], 'HeII':((2.1885), 'c', 'dotted'), 'MgII':((2.137, 2.143, ), '0.5', 'dashed'), #'MgI':((2.10655, 2.10666), '0.5', 'dotted'), 'HI':[(2.1661, 2.360, 2.367, 2.374, 2.383, 2.392, 2.404, 2.416, 2.431, 2.449), 'b', 'dashed'], 'H2':[(2.122), 'b', 'dotted'], 'NaI':[(2.208, 2.208969), 'y', 'dashed'], 'NIII':[(2.1155, 2.247, 2.251), (0,1,0.5), 'dashed'], 'FeII':((2.089, 2.145, 2.241, 2.368), 'g', 'dashed'), r'$^{12}$C$^{16}$O HB':([2.2935, 2.3227, 2.3535, 2.3829, 2.4142], 'r', 'dashed'), r'$^{13}$C$^{16}$O HB':([2.3448, 2.3739, 2.4037, 2.4341, 2.4971], 'r', 'dotted'), #'AlI':((2.109884,2.116958),(0.2,0.5,0.8), 'dotted'), #'MgI':((2.121393,2.146472,2.386573),'r'), #'ScI':((2.20581,2.20714),'m'), #'SiII':((2.200),'g'), #'CaI':((1.89282, 1.94508, 1.98702, 1.99318, 2.261410, 2.263110, 2.265741),'g') } if onlySpectrum: for k in list(lines.keys()): plt.vlines(lines[k][0], 0, 5*plt.ylim()[1], color=lines[k][1], linestyle=lines[k][2], label=k, alpha=0.3, linewidth=1) # -- https://iopscience.iop.org/article/10.1088/0067-0049/185/2/289/pdf table 6 # -- solar: datal = [(1.97280, 'SiI'), (1.97822, 'CaI'), (1.98204, 'CaI'), (1.98585, 'CaI'), (1.98676, 'CaI'), (1.99226, 'CaI'), (1.99344, 'SiI'), (1.99392, 'CaI'), (1.99673, 'CaI'), (2.03019, 'SiI'), (2.03075, 'SiI'), (2.03494, 'SiI'), (2.03840, 'SiI'), (2.06085, 'SiI'), (2.0635328, 'FeI'), (2.0634902, 'FeI'), (2.07040, 'FeI'), (2.07043, 'SiI'), (2.07226, 'FeI'), (2.08099, 'SiI'), (2.08107, 'FeI'), (2.08465, 'FeI'), (2.09224, 'SiI'), (2.10655, 'MgI'), (2.10666, 'MgI'), (2.10988, 'AlI'), (2.11696, 'AlI'), (2.13601, 'SiI'), (2.17857, 'SiI'), (2.18257, 'SiI'), (2.18853, 'SiI'), (2.20624, 'NaI'), (2.20688, 'SiI'), (2.20897, 'NaI'), (2.22632, 'FeI'), (2.23869, 'FeI'), (2.24794, 'FeI'), (2.25438, 'SiI'), (2.26141, 'CaI'), (2.26260, 'FeI'), (2.26311, 'CaI'), (2.26573, 'CaI'), (2.26720, 'SiI'), (2.28142, 'MgI'), (2.31509, 'FeI'), (2.33548, 'NaI'), (2.33855, 'NaI'),] # -- arcturus: datal = [(1.9727936, 'SiI'), (1.9728928, 'MgI'), (1.9739257, 'MgI'), (1.9782195, 'CaI'), (1.9797279, 'FeI'), (1.9820428, 'CaI'), (1.9852129, 'FeI'), (1.9858512, 'CaI'), (1.9867646, 'CaI'), (1.9922628, 'CaI'), (1.9928778, 'FeI'), (1.9934361, 'SiI'), (1.9939170, 'CaI'), (1.9967282, 'CaI'), (2.0097535, 'FeI'), (2.0286617, 'FeI'), (2.0301933, 'SiI'), (2.0349441, 'SiI'), (2.0355272, 'FeI'), (2.0383973, 'SiI'), (2.0569580, 'FeI'), (2.0608501, 'SiI'), (2.0635321, 'FeI'), (2.0703959, 'FeI'), (2.0704272, 'SiI'), (2.0722606, 'FeI'), (2.0809875, 'SiI'), (2.0810766, 'FeI'), (2.0846491, 'FeI'), (2.0922821, 'SiI'), (2.1065505, 'MgI'), (2.1066642, 'MgI'), (2.1098835, 'AlI'), (2.1169577, 'AlI'), (2.1244261, 'FeI'), (2.1360087, 'SiI'), (2.1785696, 'SiI'), (2.1788888, 'TiI'), (2.1825655, 'SiI'), (2.1885319, 'SiI'), (2.1903353, 'TiI'), (2.2010505, 'TiI'), (2.2057965, 'ScI'), (2.2062447, 'NaI'), (2.2068751, 'SiI'), (2.2089689, 'NaI'), (2.2217295, 'TiI'), (2.2238901, 'TiI'), (2.2263176, 'FeI'), (2.2266250, 'FeI'), (2.2280107, 'TiI'), (2.2316670, 'TiI'), (2.2386901, 'FeI'), (2.2398975, 'FeI'), (2.2450020, 'TiI'), (2.2479405, 'FeI'), (2.2543838, 'SiI'), (2.2614114, 'CaI'), (2.2626002, 'FeI'), (2.2631137, 'CaI'), (2.2632897, 'CaI'), (2.2657359, 'CaI'), (2.2671964, 'SiI'), (2.2814252, 'MgI'), (2.2969596, 'TiI'), (2.3150909, 'FeI'), (2.3354795, 'NaI'), (2.3385518, 'NaI'), (2.3447857, 'TiI'), (2.4572999, 'MgI'), (2.4825687, 'MgI'), (2.4867327, 'MgI'),] # -- organise lines: lines = {} for l in datal: if l[1] in lines: lines[l[1]].append(l[0]) else: lines[l[1]] = [l[0]] if onlySpectrum: for k in lines.keys(): plt.plot(lines[k], 1.1*np.ones(len(lines[k])), '|k') for x in lines[k]: plt.text(x, 1.1, k, ha='center', va='bottom') plt.legend(loc='upper left', fontsize=5) plt.hlines(1, wlmin, wlmax, linestyle='dotted') if not model is '': plt.plot(rm['wl'][w], rm['FLUX'][w], '-r', alpha=0.5, linewidth=1, linestyle='steps') if onlySpectrum: plt.legend(loc='lower center', fontsize=6, ncol=10) plt.xlabel('wavelength (um)') plt.xlim(wlmin, wlmax) plt.ylim(0, 1.05*tmp[w].max()) #plt.show() return model ax.xaxis.grid() if filt: import d4 # -- V2 and visPHI ---- r['VIS band'] = {} r['dVIS band'] = {} r['V2 band'] = {} r['dV2 band'] = {} r['dPHI band'] = {} r['CP band'] = {} r['dCP band'] = {} if not model is '': rm['VIS band'] = {} rm['dVIS band'] = {} rm['V2 band'] = {} rm['dV2 band'] = {} rm['dPHI band'] = {} rm['CP band'] = {} rm['dCP band'] = {} for i,B in enumerate(r[vis].keys()): if withdPhi: axv = plt.subplot(6,3,3+3*i, sharex=ax) else: axv = plt.subplot(6,2,2+2*i, sharex=ax) if i==0: if visamp: plt.title('VISAMP') else: plt.title(r'V$^2$') # -- filter negative measurements: for j in np.where(r[vis][B]<=0)[0]: r[vis][B][j] = np.median(r[vis][B][max(j-7, 0): min(j+7, len(r[vis][B])-1)]) filtV = np.linspace(1,0,8)**0.5 filtO = 16 if r['SPEC RES']=='HIGH' and filt: print('--smoothing--') spr = r[vis][B][w] spr = removenans(rm['wl'][w], spr) # -- sigma clipping #tmp = slidingOp(r['wl'][w], spr, 0.05) #w_ = np.where(np.abs(spr-tmp['median'])>3*tmp['1sigma']) #spr[w_] = tmp['median'][w_] plt.plot(r['wl'][w], d4.filter1D(spr, filtV, order=filtO), '-k', alpha=0.8, linewidth=1) plt.plot(r['wl'][w], r[vis][B][w], '-k', alpha=0.3, linewidth=1,label=B) if not model is '': plt.plot(rm['wl'][w], rm[vis][B][w], '-r', alpha=0.5, linewidth=1, linestyle='steps') else: plt.plot(r['wl'][w], r[vis][B][w], '-k', alpha=0.5, linewidth=1,label=B, linestyle='steps') r[vis+' band'][B] = r[vis][B][w] if computeDiff: cc = nanpolyfit(r['wl'][wc],r[vis][B][wc],1) r['d'+vis+' band'][B] = r[vis][B][w]/np.polyval(cc, r['wl'][w]) if not model is '': plt.plot(rm['wl'][w], rm[vis][B][w], '-r', alpha=0.5, linewidth=1, linestyle='steps') # for k in list(lines.keys()): # plt.vlines(lines[k][0], 0, 2*plt.ylim()[1], color=lines[k][1], # linestyle='dashed') tmp = getYlim(r[vis][B][w]) plt.ylim(max(tmp[0], 0), tmp[1]) plt.legend(loc='upper left', fontsize=7, ncol=1) axv.xaxis.grid() # -- visphi if withdPhi: axp = plt.subplot(6,3,2+3*i, sharex=ax) if i==0: plt.title('Diff Phase (deg)') if computeDiff: # -- fit slope at edges cc = nanpolyfit(r['wl'][wc],r['VISPHI'][B][wc], 1) if not model is '': ccm = nanpolyfit(rm['wl'][wc],rm['VISPHI'][B][wc], 1) else: # -- fit global slope cc = nanpolyfit(r['wl'][w], r['VISPHI'][B][w] , 1) if not model is '': ccm = nanpolyfit(rm['wl'][w], rm['VISPHI'][B][w], 1) # -- store result in another variable r['dPHI band'][B] = r['VISPHI'][B][w] - np.polyval(cc, r['wl'][w]) if not model is '': rm['dPHI band'][B] = rm['VISPHI'][B][w] - np.polyval(ccm, rm['wl'][w]) if r['SPEC RES']=='HIGH' and filt: print('--smoothing--') spr = r['VISPHI'][B][w] - np.polyval(cc, r['wl'][w]) spr = removenans(rm['wl'][w], spr) # -- sigma clipping #tmp = slidingOp(r['wl'][w], spr, 0.05) #w_ = np.where(np.abs(spr-tmp['median'])>3*tmp['1sigma']) #spr[w_] = tmp['median'][w_] spr = d4.filter1D(spr, filtV, order=filtO) plt.plot(r['wl'][w], spr, '-k', linewidth=1, alpha=0.8) plt.plot(r['wl'][w], r['VISPHI'][B][w] - np.polyval(cc, r['wl'][w]), '-k', alpha=0.3, linewidth=1, label=B) plt.ylim(getYlim(r['VISPHI'][B][w] - np.polyval(cc, r['wl'][w]) )) else: plt.plot(r['wl'][w], r['dPHI band'][B], '-k', alpha=0.5, linewidth=1, label=B, linestyle='steps') if not model is '': plt.plot(rm['wl'][w], rm['dPHI band'][B], '-r', alpha=0.5, linewidth=1, linestyle='steps') plt.ylim(getYlim(r['VISPHI'][B][w]- np.polyval(cc, r['wl'][w]))) plt.hlines(0, wlmin, wlmax, linestyle='dotted') # for k in list(lines.keys()): # plt.vlines(lines[k][0], -90, 90, color=lines[k][1], # linestyle='dashed') plt.legend(loc='upper left', fontsize=7, ncol=1) axp.xaxis.grid() axv.set_xlabel('wavelength (um)') if withdPhi: axp.set_xlabel('wavelength (um)') # -- T3 ---- for i,B in enumerate(r['T3'].keys()): if withdPhi: axx = plt.subplot(5,3,4+3*i, sharex=ax) else: axx = plt.subplot(5,2,3+2*i, sharex=ax) if True: tmp = r['T3'][B] else: tmp = (r['T3'][B]+90)%180-90 wr = r['T3'][B][w]<=-90.0 if wr[0].sum(): plt.plot(r['wl'][w][wr], tmp[w][wr], '.r', alpha=0.5, label='<-90') wb = r['T3'][B][w]>=90.0 if wb[0].sum(): plt.plot(r['wl'][w][wb], tmp[w][wb], '.b', alpha=0.5, label='>90') if r['SPEC RES']=='HIGH' and filt: print('--smoothing--') s, c = np.sin(np.pi*r['T3'][B]/180.), np.cos(np.pi*r['T3'][B]/180.) s = removenans(rm['wl'], s) c = removenans(rm['wl'], c) # -- sigma clipping #ts, tc = slidingOp(r['wl'], s, 0.05), slidingOp(r['wl'], c, 0.05) #w_ = np.where(np.abs(s-ts['median'])>3*ts['1sigma']) #s[w_] = ts['median'][w_] #w_ = np.where(np.abs(s-tc['median'])>3*tc['1sigma']) #c[w_] = tc['median'][w_] s = d4.filter1D(s, filtV, order=filtO)[w] c = d4.filter1D(c, filtV, order=filtO)[w] a = np.arctan2(s,c)*180/np.pi a = (a+90)%180-90 plt.plot(r['wl'][w], a, '-k', alpha=0.8, linewidth=1) plt.plot(r['wl'][w], tmp[w], '-k', alpha=0.3, linewidth=1, label=B) else: plt.plot(r['wl'][w], tmp[w], '-k', alpha=0.5, linewidth=1, label=B, linestyle='steps') if not model is '': plt.plot(rm['wl'][w], rm['T3'][B][w], '-r', alpha=0.5, linewidth=1, linestyle='steps') r['CP band'][B] = tmp[w] if computeDiff: cc = nanpolyfit(r['wl'][wc], tmp[wc], 1) r['dCP band'][B] = tmp[w] - np.polyval(cc, r['wl'][w]) if not model is '': ccm = nanpolyfit(rm['wl'][wc], rm['T3'][B][wc], 1) rm['dCP band'][B] = rm['T3'][B][w] - np.polyval(ccm, r['wl'][w]) plt.ylim(getYlim(tmp[w])) plt.hlines(0, wlmin, wlmax, linestyle='dotted') # for k in list(lines.keys()): # plt.vlines(lines[k][0], -150 , 150, color=lines[k][1], # linestyle='dashed') axx.xaxis.grid() plt.legend(loc='upper left', fontsize=7, ncol=1) plt.xlabel('wavelength (um)') plt.xlim(wlmin, wlmax) print('done') if not exportFilename is '': print('exporting to', exportFilename, '(open with cPickle)') f = open(exportFilename, 'wb') pickle.dump(r, f) f.close() #plt.show() return model def removenans(x,y,n=5): x, y = np.array(x), np.array(y) w = np.isfinite(x*y) yp = np.zeros(len(y)) s = scipy.signal.medfilt(y[w], n) wn = np.isnan(x*y) yp[wn] = np.interp(x[wn], x[w], s) yp[w] = y[w] return yp def nanpolyfit(x,y,o): x, y = np.array(x), np.array(y) w = np.isfinite(x*y) try: return np.polyfit(x[w], y[w], o) except: return [0] _gray80 = '#BBBBBB' _gray30 = '#444444' _crimson = '#DC143C' _myblue = '#224488' _myorange = '#886622' _myLightblue = '#44BBFF' _myLightorange = '#FFBB44' class guiPlot(tkinter.Frame): export = False # -- enable the export menu: smooth = False # -- smoothing of high-res data def __init__(self,root, directory=None): self.root = root self.widthProgressBar = 80 self.n = 0 if directory is None or not os.path.exists(directory): self.directory = None self.changeDir() else: self.directory = directory self.mainFrame = None # -- default font self.font = tkinter.font.Font(family='courier', size=12) self.smallfont = tkinter.font.Font(family='courier', size=10) if platform.uname()[0]=='Darwin': # -- Mac OS self.font = tkinter.font.Font(family='Andale Mono', size=14) self.smallfont = tkinter.font.Font(family='Andale Mono', size=12) elif platform.uname()[1]=='wvgoff': # -- Paranal VLTI offline machine self.font = None self.smallfont = None self.makeMainFrame() def quit(self): self.root.destroy() quit() def get_listFiles(self, quiet=False): files = os.listdir(self.directory) mjdobs, tplid = [], [] if not quiet: print(time.asctime()+' Filtering %d FITS files ...'%len(files)) files = [x for x in files if x.endswith('.fits')] # KLUDGE! files = [x for x in files if not x.startswith('r.AMBER')] self.checkDir = len(files) N = self.widthProgressBar for i, f in enumerate(files): n = int(i*N/float(len(files))+1) if not quiet: print('|'+'='*n+' '*(N-n-1)+'|') try: h = fits.open(os.path.join(self.directory, f)) except: continue # -- must be from GRAVITY if 'INSTRUME' in h[0].header and \ h[0].header['INSTRUME']!='GRAVITY': h.close() continue if 'MJD-OBS' in h[0].header: mjdobs.append(h[0].header['MJD-OBS']) else: mjdobs.append(0) if 'ESO TPL ID' in h[0].header and 'ESO PRO CATG' in h[0].header: tplid.append('/'.join([h[0].header['ESO TPL ID'], h[0].header['ESO PRO CATG']])) else: tplid.append('') h.close() if not quiet: print('\033[F', end=' ') mjdobs = np.array(mjdobs) #print tplid w = np.where([i.startswith('GRAVITY') and # -- GRAVITY template '_obs_' in i and # -- an observation (not a calibration) not '_SKY' in i and # -- not a sky 'VIS' in i # -- visibilities (not dark, etc.) for i in tplid]) # -- debug: #print '-->', tplid files = list(np.array(files)[w][np.argsort(mjdobs[w])]) return files def makeMainFrame(self): self.mainFrame = tkinter.Frame.__init__(self, self.root, bg=_gray30) self.waveFrame = tkinter.Frame(self.mainFrame, bg=_gray30) self.actFrame = tkinter.Frame(self.mainFrame, bg=_gray30) self.modelFrame = tkinter.Frame(self.mainFrame, bg=_gray30) if self.export: self.exportFrame = tkinter.Frame(self.mainFrame, bg=_gray30) self.listFrame = None self.root.title('GRAVIQL '+os.path.abspath(self.directory)) bo = {'fill':'both', 'padx':1, 'pady':1, 'side':'left'} b = tkinter.Button(self.actFrame, text='Sp, CP, dPhi, V', font=self.font, command = self.quickViewAll) b.pack(**bo); b.config(bg=_gray80, fg=_myblue) b = tkinter.Button(self.actFrame, text='CP, V', font=self.font, command = self.quickViewCPV2) b.pack(**bo); b.config(bg=_gray80, fg=_myblue) b = tkinter.Button(self.actFrame, text='V(B)', font=self.font, command = self.quickViewV2) b.pack(**bo); b.config(bg=_gray80, fg=_myblue) b = tkinter.Button(self.actFrame, text='Spectrum', font=self.font, command= self.quickViewSpctr) b.pack(**bo); b.config(bg=_gray80, fg=_myblue) self.vis = tkinter.StringVar() self.vis.set('VIS') # default value b = tkinter.OptionMenu(self.actFrame, self.vis, 'VIS', 'V2') b.pack(**bo); b.config(bg=_gray80, fg=_myorange) self.spectro = tkinter.StringVar() self.spectro.set('SC') # default value b = tkinter.OptionMenu(self.actFrame, self.spectro, 'SC', 'FT') b.pack(**bo); b.config(bg=_gray80, fg=_myorange) url = 'https://github.com/amerand/GRAVIQL' #b = Tkinter.Label(self.actFrame, text=url, font=self.font, # justify='center', anchor='center') b = tkinter.Button(self.actFrame, text=url, command=lambda u=url:openUrl(u)) b.pack(**bo); b.config(bg=_gray80, fg='#224488') b = tkinter.Button(self.actFrame, text='Reload Files', font=self.font, command= self.makeFileFrame) b.pack(**bo); b.config(bg=_gray80, fg=_myorange) b = tkinter.Button(self.actFrame, text='Change Directory', font=self.font, command= self.changeDir) b.pack(**bo); b.config(bg=_gray80, fg=_myorange) b = tkinter.Button(self.actFrame, text='QUIT', font=self.font, command= self.quit) b.pack(**bo); b.config(bg=_crimson, fg='#000000') modes = [('Full', 'None None'), #('High SNR', '2.05 2.42'), #('Brd 1.945', '%.4f %.4f'%(1.9451-0.02, 1.9451+0.02)), ('HeI 2.058', '2.038 2.078'), ('FeII 2.089', '2.074 2.104'), # Warm and dense CSE https://arxiv.org/pdf/1408.6658.pdf sec 3.2.1 ('HeI 2.112', '2.0926 2.1326'), ('NIII 2.115', '%.4f %.4f'%(2.1155-0.02, 2.1155+0.02)), # http://adsabs.harvard.edu/abs/1996ApJS..107..281H ('MgII 2.140','2.130 2.150'), # doublet ('Brg 2.166', '2.136 2.196'), ('HeII 2.188', '%.4f %.4f'%(2.1885-0.02, 2.1885+0.02)), # http://adsabs.harvard.edu/abs/1996ApJS..107..281H ('NaI 2.206', '2.198 2.218'), ('NIII 2.249', '2.237 2.261'), ('CO bands', '2.28 None')] bo = {'fill':'both', 'side':'left', 'padx':1, 'pady':1} self.wlrange = tkinter.StringVar() self.wlrange.set('None None') # default sepactral range self.plot_opt={} self.setPlotRange() # -- create a button for each spectral range: for text, mode in modes: b = tkinter.Radiobutton(self.waveFrame, text=text, variable=self.wlrange, value=mode, font=self.smallfont, indicatoron=0) b.pack(**bo) b.config(selectcolor=_myorange, fg=_gray80, bg=_gray30) # -- back from the right now: bo = {'fill':'both', 'side':'right', 'padx':0, 'pady':1} # -- drop down menu for instrument (SC of FT) # self.spectro = tkinter.StringVar() # self.spectro.set('SC') # default value # b = tkinter.OptionMenu(self.waveFrame, self.spectro, 'SC', 'FT') # b.pack(**bo); b.config(bg=_gray80, fg=_myorange) b = tkinter.Label(self.waveFrame, text='um', bg=_gray30, fg=_gray80, font=self.font) b.pack(**bo) b = tkinter.Entry(self.waveFrame, textvariable=self.wlrange, width=12, font=self.font) b.pack(**bo) b = tkinter.Label(self.waveFrame, text='WL range', bg=_gray30, fg=_gray80, font=self.font) b.pack(**bo) # -- input for visibility model self.modelStr = tkinter.StringVar() self.modelStr.set('') # default value b = tkinter.Entry(self.modelFrame, textvariable=self.modelStr, width=115, font=self.font) b.pack(**bo) # b = Tkinter.Label(self.modelFrame, text='Vis. Model:', # bg=_gray30, fg=_gray80, font=self.font) b = tkinter.Button(self.modelFrame, text='Vis. Model:', bg=_gray30, fg=_gray80, font=self.font, command=clickVisModel) b.pack(**bo) if self.export: self.exportStr = tkinter.StringVar() self.exportStr.set('filename.dpy') # default value b = tkinter.Entry(self.exportFrame, textvariable=self.exportStr, width=100, font=self.font) b.pack(**bo) b = tkinter.Label(self.exportFrame, text='export to binary file:', bg=_gray30, fg=_gray80, font=self.font) b.pack(**bo) self.makeFileFrame() return def tick(self): #print ['tic', 'tac'][self.n%2], time.asctime() self.n += 1 tmp = len([f for f in os.listdir(self.directory) if f.endswith('fits')]) if tmp!=self.checkDir: if tkinter.messagebox.askyesno('GRAVIQL', 'the directory has changed, do you want to update the file list?'): self.makeFileFrame() self.listFrame.after(60000, self.tick) return def changeDir(self): self.root.update() self.directory = tkinter.filedialog.askdirectory(initialdir=self.directory) self.root.title('GRAVIQL '+self.directory) self.makeFileFrame() self.root.update() return def makeFileFrame(self): if not self.listFrame is None: # -- reload for widget in self.listFrame.winfo_children(): widget.destroy() self.listFrame.destroy() self.listFrame.pack_forget() # -- temporary message: #c = Tkinter.Label(self.listwaveFrame, text='Reloading ...') #c.pack(fill='both') #self.actFrame.pack(anchor='nw', fill='x') #self.waveFrame.pack(anchor='nw', fill='x') # -- list all files self.filename = None self.checkList = {} self.listFiles = self.get_listFiles() _format = '%3s %-13s %14s %7d %2s/%3s %ss %11s %3s"/%2sms %3.0f%%(%3.0fms) %3.0f%%(%1.0f) %7s %16s %5s %5s' _legend = ' Object Prog ID Contain. Mode DIT Baseline See/T0 @V FT(T0@K) SC(nB) corr. Date-Obs LST Pipe' #print '' #print legend container = None if len(self.listFiles)==0: self.listFrame = tkinter.Frame(self.mainFrame, bg=_gray30) c = tkinter.Label(self.listFrame, text='--- no relevant files to display ---') c.pack(fill='both') self.actFrame.pack(anchor='nw', fill='x') self.waveFrame.pack(anchor='nw', fill='x') self.modelFrame.pack(anchor='nw', fill='x') if self.export: self.exportFrame.pack(anchor='nw', fill='x') self.listFrame.pack() return else: self.listFrame = tkinter.Frame(self.mainFrame, bg=_gray30) self.scrollbar = tkinter.Scrollbar(self.listFrame, orient='vertical') c = tkinter.Label(self.listFrame, text=_legend, bg=_gray30, fg=_gray80, font=self.font) c.pack(fill='both') self.listBox = tkinter.Listbox(self.listFrame, selectmode='multiple', yscrollcommand=self.scrollbar.set, height=min(len(self.listFiles), 35), width=len(_legend)+5) self.scrollbar.config(command=self.listBox.yview) self.scrollbar.pack(side='right', fill='y') BG = [ (_myblue, _gray80), (_myorange, _gray80) ] ABG = [(_myLightblue, _gray30), (_myLightorange, _gray30) ] FG = [(_gray80, _myblue), (_gray80, _myorange) ] ic = -1 container = -1 self.listAllFiles = [] print('\033[F', end=' ') print(time.asctime()+' Loading %d relevant files...'%len(self.listFiles)) N = self.widthProgressBar for _i, fi in enumerate(self.listFiles): n = int(_i*N/float(len(self.listFiles))+1) print('|'+'='*n+' '*(N-n-1)+'|') key = os.path.join(self.directory, fi) f = fits.open(key) if 'INSTRUME' in f[0].header and \ f[0].header['INSTRUME']!='GRAVITY': f.close() continue self.listAllFiles.append(key) self.checkList[key] = tkinter.IntVar() if not 'ESO OBS CONTAINER ID' in list(f[0].header.keys()): cont = -1 else: cont = f[0].header['ESO OBS CONTAINER ID'] if container != cont: ic+=1 container = cont seeing = 0.5*(f[0].header['ESO ISS AMBI FWHM START']+ f[0].header['ESO ISS AMBI FWHM END']) tau0 = 500*(f[0].header['ESO ISS AMBI TAU0 START']+ f[0].header['ESO ISS AMBI TAU0 END']) # -- QC fringe tracker: FT = [] SC = [] T0 = [] for b in ['12','13','14','23','24','34']: _key = 'ESO QC TAU0 OPDC'+b if _key in list(f[0].header.keys()): T0.append(1000*f[0].header[_key]) _key = 'ESO QC ACCEPTED_RATIO_FT%s_P1'%b if _key in list(f[0].header.keys()): FT.append(f[0].header[_key]) _key = 'ESO QC ACCEPTED_RATIO_FT%s_P2'%b if _key in list(f[0].header.keys()): FT.append(f[0].header[_key]) _key = 'ESO QC ACCEPTED_RATIO_SC%s_P1'%b if _key in list(f[0].header.keys()): SC.append(f[0].header[_key]) _key = 'ESO QC ACCEPTED_RATIO_SC%s_P2'%b if _key in list(f[0].header.keys()): SC.append(f[0].header[_key]) dit = f[0].header['ESO DET2 SEQ1 DIT'] baseline = '-'.join([f[0].header['ESO ISS CONF STATION%d'%i] for i in [1,2,3,4]]) try: lst = f[0].header['LST']/3600. except: lst = 0.0 lst ='%02d:%02.0f'%(int(lst), 60*(lst%1)) dit = f[0].header['ESO DET2 SEQ1 DIT'] dit = '%2.0f'%dit if dit>1 else ('%2.1f'%dit)[1:] tau0 = '%2.0f'%tau0 if tau0>1 else ('%2.1f'%tau0)[1:] seeing = '%3.1f'%seeing if seeing>1 else ('%4.2f'%seeing)[1:] # -- reduction parameters: p = {f[0].header[k]:int(k.split('PARAM')[1].split()[0]) for k in list(f[0].header.keys()) if 'ESO PRO REC1 PARAM' in k and 'NAME' in k} p = {k:f[0].header['ESO PRO REC1 PARAM%d VALUE'%p[k]] for k in list(p.keys())} text = _format%('CAL' if 'Calibrator' in f[0].header['ESO TPL NAME'] else 'SCI', f[0].header['ESO INS SOBJ NAME'][:12]+('*' if 'calibrated' in fi else ' '), f[0].header['ESO OBS PROG ID'], container, f[0].header['ESO INS SPEC RES'][:2], f[0].header['ESO INS POLA MODE'][:3], dit, baseline, seeing, tau0, np.median(FT), max(-1, np.median(T0)), np.median(SC), np.sum(np.array(SC)>0)/float(len(SC))*6., '', #p['vis-correction'], f[0].header['DATE-OBS'][:-3], lst, f[0].header['ESO PRO REC1 PIPE ID'].split()[-1]) f.close() self.listBox.insert('end', text) # c = Tkinter.Checkbutton(self.listFrame, text=text, font=self.font, # variable=self.checkList[key], # onvalue=1, offvalue=0) # c.pack(fill='both', side='top') # c.config(justify='left', borderwidth=0, padx=0, pady=0, anchor='e', # foreground =FG[ic%2][0] if 'sciraw' in key else FG[ic%2][1], # selectcolor=BG[ic%2][0] if 'sciraw' in key else BG[ic%2][1], # fg=FG[ic%2][0] if 'sciraw' in key else FG[ic%2][1], # bg=BG[ic%2][0] if 'sciraw' in key else BG[ic%2][1], # activebackground=ABG[ic%2][0] if 'sciraw' in key else BG[ic%2][1], # highlightbackground=ABG[ic%2][0] if 'sciraw' in key else BG[ic%2][1], # #relief='flat', offrelief='flat', overrelief='flat', # highlightthickness=0,) #colorsT = [('46', '36'), ('43', '33')] # if 'sciraw' in key: # print '\033[%sm'%colorsT[ic%2][0]+text+'\033[0m' # else: # print '\033[%sm'%colorsT[ic%2][1]+text+'\033[0m' print('\033[F', end=' ') print('\033[F', end=' ') print('') print('') if len(self.listFiles)>0: self.listBox.pack(fill='both', expand=1) self.listBox.config(font=self.font, highlightbackground=FG[0][0]) self.actFrame.pack(anchor='nw', fill='x') self.waveFrame.pack(anchor='nw', fill='x') self.modelFrame.pack(anchor='nw', fill='x') if self.export: self.exportFrame.pack(anchor='nw', fill='x') self.listFrame.pack() self.listFrame.after(60000, self.tick) return def setFileList(self): self.filename = [] #for k in self.checkList.keys(): # if self.checkList[k].get(): # self.filename.append(k) items = self.listBox.curselection() if len(items)>0 and not isinstance(items[0], int): items = list(map(int, items)) self.filename = [self.listAllFiles[i] for i in items] return def setPlotRange(self): try: if self.wlrange.get().split()[0]=='None': self.plot_opt['wlmin'] = None else: self.plot_opt['wlmin'] = float(self.wlrange.get().split()[0]) if self.wlrange.get().split()[1]=='None': self.plot_opt['wlmax'] = None else: self.plot_opt['wlmax'] = float(self.wlrange.get().split()[1]) except: tkinter.messagebox.showerror('ERROR', 'range format is "wlmin wlmax" in um') return False return True def quickViewV2(self): self.setFileList() if not self.setPlotRange(): return if self.filename is None: tkinter.messagebox.showerror('ERROR', 'no file selected') return test = plotGravi(self.filename, v2b=True, model=self.modelStr.get(), insname='auto_'+self.spectro.get(), filt=self.smooth, exportFilename='' if not self.export else self.exportStr.get(), visamp=self.vis.get()=='VIS', **self.plot_opt) #print 'TEST:', test if isinstance(test, str): self.modelStr.set(test) #tkMessageBox.showerror('ERROR', 'Incorrect file selection') return def quickViewAll(self): self.setFileList() if not self.setPlotRange(): return if self.filename is None: tkinter.messagebox.showerror('ERROR', 'no file selected') return test = plotGravi(self.filename, insname='auto_'+self.spectro.get(), model=self.modelStr.get(), filt=self.smooth, exportFilename='' if not self.export else self.exportStr.get(), withdPhi = True, visamp=self.vis.get()=='VIS', **self.plot_opt) #print 'TEST:', test if isinstance(test, str): self.modelStr.set(test) #tkMessageBox.showerror('ERROR', 'Incorrect file selection') return def quickViewCPV2(self): self.setFileList() if not self.setPlotRange(): return if self.filename is None: tkinter.messagebox.showerror('ERROR', 'no file selected') return test = plotGravi(self.filename, insname='auto_'+self.spectro.get(), model=self.modelStr.get(), filt=self.smooth, exportFilename='' if not self.export else self.exportStr.get(), withdPhi = False, visamp=self.vis.get()=='VIS', **self.plot_opt) #print 'TEST:', test if isinstance(test, str): self.modelStr.set(test) #tkMessageBox.showerror('ERROR', 'Incorrect file selection') return def quickViewSpctr(self): self.setFileList() if not self.setPlotRange(): return if self.filename is None: tkinter.messagebox.showerror('ERROR', 'no file selected') return test = plotGravi(self.filename, onlySpectrum=True, model=self.modelStr.get(), insname='auto_'+self.spectro.get(), filt=self.smooth, exportFilename='' if not self.export else self.exportStr.get(), visamp=self.vis.get()=='VIS', **self.plot_opt) #print 'TEST:', test if isinstance(test, str): self.modelStr.set(test) #tkMessageBox.showerror('ERROR', 'Incorrect file selection') return if __name__=='__main__': root = tkinter.Tk() root.config(bg=_gray30) if len(sys.argv)==1: directory = tkinter.filedialog.askdirectory() else: directory = sys.argv[1] guiPlot(root, directory).pack() root.mainloop()