#!/usr/bin/env python # -*- coding: utf-8 -*- """Measure North position angle and plate-scale on acquisition camera. Measure the plate scale and error in orientation of the acquisition camera by fitting Gaussians on the FT and SC star and comparing to the information in the FITS headers. The default output consists in one line listing the plate scales for the for GRAVITY ports (left to right on the acquisition camera), followed by the error made in derotation for the four ports. To apply this correction, subtract it from either kmirrorOffset, rotationOffset or whatever gets into INS.DROTOFF[1-4]. By default, a plot of the time evolution of these parameters will be dsiplayed (unless -p0 is passed). The input consists in one or several raw object GRAVITY files obtained in dual-field mode (DPR.TYPE=OBJECT,DUAL). These files may have been acquired on distinct astrophysical objects. By default, the script will "./GRAVI*.fits". The input files will be filtered according to --keywords_equal, by default only files with "ESO DPR TYPE"=="OBJECT,DUAL" and "ESO DPR CATG"=="SCIENCE" will be kept. If "ESO INS SOBJ X" and "ESO INS SOBJ Y" do not reflect the correct separation (partcicularly the case when offsets are used in exposure template), the correct separation for each pair may be given with --rho. By default, the script will trust the first separation for each pair (see --first_offsets), so that a mapping sequence will work assuming the offsets where 0 for the first exposure. The astrophysical objects must be such that the targets that were positioned on the FT and SC fibres during acquisition are the brightest in there neighbourhood (the size of this neighbourhood can be specified with --window). It is recommended to use one or more dark frames. By default, the script will look for master dark frames with the appropriate DIT as gvacq_MasterSkyField.fits and gvacq_MasterSkyField_2.8.fits (simple full-frame images used by the real time software) in: - the current direcory; - ../common_calibration; - ${INS_ROOT}; - in the gravi_acqscale package. If such dark frames are not available or not suitable, normal dark frames (e.g. recorded as part of a P2VM) can be supplied using the --dark option (this option can be repeated to provide more than one sky frame). Such dark frames and sky frames must be cut the same way as the object frames (i.e. the value of the FITS cards DET1.FRAMES.NX, DET1.FRAMES.NY, DET1.FRAM[1-4].STRX and DET1.FRAM[1-4].STRY must match). The master dark frames mentioned above do not suffer from this limitation as they are not cut. The script first guesses the location of the two stars in each telescope using INS.SOBJ.X and INS.SOBJ.Y for the separation and IN.DROTOFF[1-4] for the position angle. It then fits Gaussians on the two stars. It finally estimates the plate scale and a correction to the field rotation that should be applied. The default behaviour (measuring and plotting plate-scales and rotations): run_gravi_acqscale.py To save an accumulated dark frame: run_gravi_acqscale.py -d dark1.fits ... -d darkN.fits -D master.fits To output the pixel positions instead of the plate-scales and rotations, use the --print_what option. Various options can be used to set the fit window, verbosity level, correct the ephemerides after the fact etc. """ try: import pyfits except: from astropy.io import fits as pyfits import numpy as np from matplotlib import pyplot as plt, dates as mdates from matplotlib.lines import Line2D from dateutil.tz import tzutc from datetime import timedelta from scipy import ndimage from astropy import time as astrotime import argparse import glob from . import parser_actions from . import io def mean_rms(data, axis=None, mean_fcn=np.mean): m = mean_fcn(data, axis=axis) r = np.sqrt(mean_fcn(np.square(data-m), axis=axis)) return m, r def format_results(label, data, uncertainties, fmt_str): out = label for i in np.arange(len(data)): out += ' ' + fmt_str.format(data[i], uncertainties[i]) return out def main(): parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter, prefix_chars='-+') parser.add_argument('fname', nargs='*') parser.add_argument("-w", "--window", help="full width of window for peak fitting [30]", type=int, default=30) parser.add_argument("-g", "--group", help="number of frames to average together [DET1.NDIT/4]", type=int) parser.add_argument("-r", "--rho", help="separation of binary (mas) [from FITS header]. May be " "prefixed with FT and SC taget names and used several times: --rho=IRS16C:S2:1200. " "The prefixed values are used for the specific pair (or the reverse) while a value " "without prefix is a global default.", type=str, action=parser_actions.Rho) # parser.add_argument("-t", "--theta", help="position angle of binary on sky (East of North), [from " # "FITS header]. May be prefixed with FT and SC target names and used several times: " # "--theta=IRS16C:S2:150. " # "The prefixed values are used for the specific pair (or the reverse) while a value " # "without prefix is a global default.", type=str, action=RhoAction) parser.add_argument("-f", "+f", "--first_offsets", "--no-first_offsets", action=parser_actions.Toggle, nargs=0, default=True, help= "trust the first offsets found for each FT/SC pair. " "This is useful if the first file for a pair " "corresponds to the real position while other files " "were taken at non-zero offset in the exposure template. " "Use +f or --no-first-offsets to disable. [Enabled]") parser.add_argument("-d", "--dark_file", help="dark file (can be used " "several times to accumulate several dark frames). " "Unnecessary if a master dark is available.", action='append') parser.add_argument("-D", "--dark_outfile", help="file to save preprocessed dark") parser.add_argument("-p", "--plot_what", help= "what to plot. 0: nothing, [1]: final fit results, 2: " "final results+individual fit residuals", type=int, default=1) parser.add_argument("-P", "--print_what", help= "what to print. 0: nothing, [1]: plate scales and PA " "errors (average), 2: pixel positions (per file), " "3: pixel offsets (per file), 4: fiber offsets (per file)", type=int, default=1) parser.add_argument("-v", "--verbose", help="verbosity level.", type=int, default=0) parser.add_argument("-k", "--keywords_equal", help= "keywords should exist and match the value. May be " "repeated. " "[--keywords_equal='ESO DPR TYPE:OBJECT,DUAL' " "--keywords_equal='ESO DPR CATG:SCIENCE']", action=parser_actions.Keyword, type=str) parser.add_argument("-x", "--exclude_object", action=parser_actions.Object, type=str, help="exclude FT:SC doublet. " "[-xCalibration_unit_fiber_1:Calibration_unit_fiber_2]") parser.add_argument("-i", "--include_object", action=parser_actions.Object, type=str, help="include FT:SC doublet. May be specified more " "than once. If specified, all other doublets will be " "excluded. --exclude_object has precedence over " "--include_object. [None]") parser.add_argument("-c", "--commoncalib-dir", dest="commoncalibdir", default=None, help="Add this directory in the search path for calibration files," "in addition of current [../common_calibration]") args = parser.parse_args() rho=args.rho if rho is None: rho=dict() if args.first_offsets: first_offsets = dict() else: first_offsets = None verbose=args.verbose if args.dark_file is None: dark=0. else: dark_file=args.dark_file infile=fits.open(dark_file.pop()) prihdr=infile[0].header dark=infile['IMAGING_DATA_ACQ'] darkhdr=dark.header dark=dark.data for f in dark_file: dark=np.append(dark, fits.open(f)['IMAGING_DATA_ACQ'].data, axis=0) NDIT=dark.shape[0] dark=np.median(dark, axis=0) if args.dark_outfile is not None: pri=fits.PrimaryHDU(header=prihdr) img=fits.ImageHDU(dark[None,:,:], header=darkhdr) pri.header.set('HIERARCH ESO DET1 NDIT', NDIT, 'number of DITs accumulated together') newfits=fits.HDUList([pri, img]) newfits.writeto(args.dark_outfile) ny=dark.shape[0]//4 dark=dark[:ny,:] if args.commoncalibdir is not None: load_master_darks(args.commoncalibdir) if args.keywords_equal is None: args.keywords_equal = {"ESO DPR TYPE": "OBJECT,DUAL"} if args.exclude_object is None: args.exclude_object =set( ["Calibration_unit_fiber_1:Calibration_unit_fiber_2"]) # If file list is emtpy, use "./GRAVI*.fits" if len(args.fname) is 0: args.fname=sorted(glob.glob("./GRAVI*.fits")) filter_files=io.Filter(include_object=args.include_object, exclude_object=args.exclude_object, keywords_equal=args.keywords_equal) args.fname=list(filter(filter_files, args.fname)) if len(args.fname) == 0: raise RuntimeError("No valid file found!") x1s_mean=np.zeros((len(args.fname),4)) x1s_rms=np.zeros((len(args.fname),4)) x2s_mean=np.zeros((len(args.fname),4)) x2s_rms=np.zeros((len(args.fname),4)) y1s_mean=np.zeros((len(args.fname),4)) y1s_rms=np.zeros((len(args.fname),4)) y2s_mean=np.zeros((len(args.fname),4)) y2s_rms=np.zeros((len(args.fname),4)) NPA_mean=np.zeros((len(args.fname),4)) NPA_rms=np.zeros((len(args.fname),4)) scales_mean=np.zeros((len(args.fname),4)) scales_rms=np.zeros((len(args.fname),4)) RelPA_mean=np.zeros((len(args.fname),4)) RelPA_rms=np.zeros((len(args.fname),4)) RelPS_mean=np.zeros((len(args.fname),4)) RelPS_rms=np.zeros((len(args.fname),4)) MJDs=np.zeros(len(args.fname)) source_names=np.ndarray(len(args.fname), dtype=np.object_) all_markers=Line2D.filled_markers cur_marker=0 n_markers=len(all_markers) markers=dict() for f in np.arange(len(args.fname)): fitsfile=io.File(args.fname[f], rho=rho, first_offsets=first_offsets, group=args.group, dark=dark, window=args.window, plot_fit=args.plot_what) MJDs[f]=fitsfile.hdulist[0].header["MJD-OBS"] source_names[f]=fitsfile.key if fitsfile.key not in list(markers.keys()): markers[fitsfile.key]=all_markers[cur_marker] cur_marker = (cur_marker+1)%n_markers if verbose: print("*******************************************************************************") print(("File name : {}".format(fitsfile.filename))) print(("FT target : {}".format(fitsfile.FTname))) print(("SC target : {}".format(fitsfile.SCname))) print(("Position angle: {}°".format(fitsfile.theta_in))) print(("Separation : {} mas".format(fitsfile.rho_in))) print(("Configuration :{}".format(fitsfile.config))) x1s=fitsfile.field_ft_x y1s=fitsfile.field_ft_y x2s=fitsfile.field_sc_x y2s=fitsfile.field_sc_y apparent_sep=np.sqrt((x2s-x1s)**2+(y2s-y1s)**2) apparent_PA=np.arctan2(x2s-x1s, y2s-y1s)*180./np.pi x1s_mean[f,:], x1s_rms[f,:]=mean_rms(x1s, axis=0, mean_fcn=np.median) y1s_mean[f,:], y1s_rms[f,:]=mean_rms(y1s, axis=0, mean_fcn=np.median) x2s_mean[f,:], x2s_rms[f,:]=mean_rms(x2s, axis=0, mean_fcn=np.median) y2s_mean[f,:], y2s_rms[f,:]=mean_rms(y2s, axis=0, mean_fcn=np.median) APA_mean, APA_rms = mean_rms(apparent_PA, axis=0, mean_fcn=np.median) psep_mean, psep_rms=mean_rms(apparent_sep, axis=0, mean_fcn=np.median) North_PA=(apparent_PA-fitsfile.approx_PA-fitsfile.offangles+90.)%360.-90. NPA_mean[f,:], NPA_rms[f,:] = mean_rms(North_PA, axis=0, mean_fcn=np.median) scales=fitsfile.rho_in/apparent_sep scales_mean[f,:], scales_rms[f,:] = mean_rms(scales, axis=0, mean_fcn=np.median) RelPA=North_PA-North_PA[:,2][:,np.newaxis] RelPA_mean[f,:], RelPA_rms[f,:] = mean_rms(RelPA, axis=0, mean_fcn=np.median) RelPS=scales/scales[:,2][:,np.newaxis] RelPS_mean[f,:], RelPS_rms[f,:] = mean_rms(RelPS, axis=0, mean_fcn=np.median) if verbose: print((format_results('FT obj X (pix):', x1s_mean[f,:], x1s_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('FT obj Y (pix):', y1s_mean[f,:], y1s_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('SC obj X (pix):', x2s_mean[f,:], x2s_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('SC obj Y (pix):', y2s_mean[f,:], y2s_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('Apparent PA(°):', APA_mean%360, APA_rms, '{:8.3f}±{:.4f}'))) print((format_results('App. pix. sep.:', psep_mean, psep_rms, '{:8.3f}±{:.4f}'))) print((format_results('North PA :', NPA_mean[f,:], NPA_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('P scales :', scales_mean[f,:], scales_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('PA(x)-PA(GV3) :', RelPA_mean[f,:], RelPA_rms[f,:], '{:8.3f}±{:.4f}'))) print((format_results('PS(x)/PS(GV3) :', RelPS_mean[f,:], RelPS_rms[f,:], '{:8.3f}±{:.4f}'))) if args.print_what is 2: txt = '{} {:8s} {:8s} '.format(fitsfile.filename, fitsfile.FTname, fitsfile.SCname) for i in np.arange(4): txt += '{:8.3f} {:8.3f} {:8.3f} {:8.3f} '.format(x1s_mean[f,i], y1s_mean[f,i], x2s_mean[f,i], y2s_mean[f,i]) print (txt) elif args.print_what is 3: txt = '{} {:8s} {:8s} '.format(fitsfile.filename, fitsfile.FTname, fitsfile.SCname) for i in np.arange(4): txt += '{:8.3f} {:8.3f} '.format(x2s_mean[f,i]-x1s_mean[f,i], y2s_mean[f,i]-y1s_mean[f,i]) print (txt) elif args.print_what is 4: if args.first_offsets is not None: ddx=fitsfile.dx_in - first_offsets[fitsfile.key][0] ddy=fitsfile.dy_in - first_offsets[fitsfile.key][1] drho=np.sqrt(ddx*ddx+ddy*ddy)/scales_mean[f,:] dth=np.arctan2(ddx, ddy)*180./np.pi PA=fitsfile.approx_PA+dth-first_offsets[fitsfile.key][2] ddx_pix = drho*np.sin(PA/180.*np.pi) ddy_pix = drho*np.cos(PA/180.*np.pi) else: ddx_pix = np.zeros(4) ddy_pix = np.zeros(4) txt = '{} {:8s} {:8s} '.format(fitsfile.filename, fitsfile.FTname, fitsfile.SCname) for i in np.arange(4): txt += '{:8.3f} {:8.3f} '.format((x2s_mean[f,i]-x1s_mean[f,i])-(fitsfile.fiber_sc_x[i]-fitsfile.fiber_ft_x[i])+ddx_pix[i], (y2s_mean[f,i]-y1s_mean[f,i])-(fitsfile.fiber_sc_y[i]-fitsfile.fiber_ft_y[i])+ddy_pix[i]) print (txt) # if args.plot_what >= 1: # colors=['b', 'g', 'r', 'c'] # markers=['+', 'o', 's', 'x'] # plt.figure(figsize=(16, 5)) # plt.subplot(1, 2, 1) # ngroups = North_PA.shape[0] # for p in np.arange(4): # plt.plot(North_PA[:,p], colors[p]+markers[p], label='GV'+str(p+1)) # plt.plot([0, ngroups], [NPA_mean[f,p], NPA_mean[f,p]], colors[p]+'-') # plt.plot([0, ngroups], [NPA_mean[f,p], NPA_mean[f,p]]-NPA_rms[f,p], colors[p]+'--') # plt.plot([0, ngroups], [NPA_mean[f,p], NPA_mean[f,p]]+NPA_rms[f,p], colors[p]+'--') # plt.legend(numpoints=1, fontsize='x-small') # plt.title("North PA") # plt.subplot(1, 2, 2) # for p in np.arange(4): # plt.plot(scales[:,p], colors[p]+markers[p], label='GV'+str(p+1)) # plt.plot([0, ngroups], [scales_mean[f,p], scales_mean[f,p]], colors[p]+'-') # plt.plot([0, ngroups], [scales_mean[f,p], scales_mean[f,p]]-scales_rms[f,p], colors[p]+'--') # plt.plot([0, ngroups], [scales_mean[f,p], scales_mean[f,p]]+scales_rms[f,p], colors[p]+'--') # plt.title("Pixel scales") # plt.legend(numpoints=1, fontsize='x-small') # plt.show() x1s_mean_mean=np.mean(x1s_mean, axis=0) x1s_rms_mean=np.mean(x1s_rms, axis=0) x1s_mean_rms=np.sqrt(np.mean(np.square(x1s_mean-x1s_mean_mean), axis=0)) y1s_mean_mean=np.mean(y1s_mean, axis=0) y1s_rms_mean=np.mean(y1s_rms, axis=0) y1s_mean_rms=np.sqrt(np.mean(np.square(y1s_mean-y1s_mean_mean), axis=0)) x2s_mean_mean=np.mean(x2s_mean, axis=0) x2s_rms_mean=np.mean(x2s_rms, axis=0) x2s_mean_rms=np.sqrt(np.mean(np.square(x2s_mean-x2s_mean_mean), axis=0)) y2s_mean_mean=np.mean(y2s_mean, axis=0) y2s_rms_mean=np.mean(y2s_rms, axis=0) y2s_mean_rms=np.sqrt(np.mean(np.square(y2s_mean-y2s_mean_mean), axis=0)) NPA_mean_mean=np.mean(NPA_mean, axis=0) NPA_rms_mean=np.mean(NPA_rms, axis=0) NPA_mean_rms=np.sqrt(np.mean(np.square(NPA_mean-NPA_mean_mean), axis=0)) PS_mean_mean=np.mean(scales_mean, axis=0) PS_rms_mean=np.mean(scales_rms, axis=0) PS_mean_rms=np.sqrt(np.mean(np.square(scales_mean-PS_mean_mean), axis=0)) RPA_mean_mean=np.mean(RelPA_mean, axis=0) RPA_rms_mean=np.mean(RelPA_rms, axis=0) RPA_mean_rms=np.sqrt(np.mean(np.square(RelPA_mean-RPA_mean_mean), axis=0)) RPS_mean_mean=np.mean(RelPS_mean, axis=0) RPS_rms_mean=np.mean(RelPS_rms, axis=0) RPS_mean_rms=np.sqrt(np.mean(np.square(RelPS_mean-RPS_mean_mean), axis=0)) if verbose: print("*******************************************************************************") print("*******************************************************************************") print("Final results : mean(mean(quantity))±rms(mean(quantity); mean(rms(quantity))") print(" PA: Position angle of North on acquisition camera,") print(" clockwise from top, for encoder position 0") print(" PS: plate scale, mas/pix") print(" RelPA: PA(GVx)-PA(GV3)") print(" RelPS: PS(GVx)/PS(GV3)") format_string='{:8.4f}±{:.4f} {:8.4f}±{:.4f} {:8.4f}±{:.4f} {:8.4f}±{:.4f}' format_string_bis='{:15.4f} {:15.4f} {:15.4f} {:15.4f}' print(('FTx (mean±rms): '+format_string.format(x1s_mean_mean[0], x1s_mean_rms[0], x1s_mean_mean[1], x1s_mean_rms[1], x1s_mean_mean[2], x1s_mean_rms[2], x1s_mean_mean[3], x1s_mean_rms[3]))) print(('mean(rms(FTx)): '+format_string_bis.format(x1s_rms_mean[0], x1s_rms_mean[1], x1s_rms_mean[2], x1s_rms_mean[3]))) print(('FTy (mean±rms): '+format_string.format(y1s_mean_mean[0], y1s_mean_rms[0], y1s_mean_mean[1], y1s_mean_rms[1], y1s_mean_mean[2], y1s_mean_rms[2], y1s_mean_mean[3], y1s_mean_rms[3]))) print(('mean(rms(FTy)): '+format_string_bis.format(y1s_rms_mean[0], y1s_rms_mean[1], y1s_rms_mean[2], y1s_rms_mean[3]))) print(('SCx (mean±rms): '+format_string.format(x2s_mean_mean[0], x2s_mean_rms[0], x2s_mean_mean[1], x2s_mean_rms[1], x2s_mean_mean[2], x2s_mean_rms[2], x2s_mean_mean[3], x2s_mean_rms[3]))) print(('mean(rms(SCx)): '+format_string_bis.format(x2s_rms_mean[0], x2s_rms_mean[1], x2s_rms_mean[2], x2s_rms_mean[3]))) print(('SCy (mean±rms): '+format_string.format(y2s_mean_mean[0], y2s_mean_rms[0], y2s_mean_mean[1], y2s_mean_rms[1], y2s_mean_mean[2], y2s_mean_rms[2], y2s_mean_mean[3], y2s_mean_rms[3]))) print(('mean(rms(SCy)): '+format_string_bis.format(y2s_rms_mean[0], y2s_rms_mean[1], y2s_rms_mean[2], y2s_rms_mean[3]))) print(('PA (mean±rms) : '+format_string.format(NPA_mean_mean[0], NPA_mean_rms[0], NPA_mean_mean[1], NPA_mean_rms[1], NPA_mean_mean[2], NPA_mean_rms[2], NPA_mean_mean[3], NPA_mean_rms[3]))) print(('mean(rms(PA)) : '+format_string_bis.format(NPA_rms_mean[0], NPA_rms_mean[1], NPA_rms_mean[2], NPA_rms_mean[3]))) print(('PS (mean±rms) : '+format_string.format(PS_mean_mean[0], PS_mean_rms[0], PS_mean_mean[1], PS_mean_rms[1], PS_mean_mean[2], PS_mean_rms[2], PS_mean_mean[3], PS_mean_rms[3]))) print(('mean(rms(PS)) : '+format_string_bis.format(PS_rms_mean[0], PS_rms_mean[1], PS_rms_mean[2], PS_rms_mean[3]))) print(('RelPA mean±rms: '+format_string.format(RPA_mean_mean[0], RPA_mean_rms[0], RPA_mean_mean[1], RPA_mean_rms[1], RPA_mean_mean[2], RPA_mean_rms[2], RPA_mean_mean[3], RPA_mean_rms[3]))) print(('mean(rms(RPA)): '+format_string_bis.format(RPA_rms_mean[0], RPA_rms_mean[1], RPA_rms_mean[2], RPA_rms_mean[3]))) print(('RelPS mean±rms: '+format_string.format(RPS_mean_mean[0], RPS_mean_rms[0], RPS_mean_mean[1], RPS_mean_rms[1], RPS_mean_mean[2], RPS_mean_rms[2], RPS_mean_mean[3], RPS_mean_rms[3]))) print(('mean(rms(RPS)): '+format_string_bis.format(RPS_rms_mean[0], RPS_rms_mean[1], RPS_rms_mean[2], RPS_rms_mean[3]))) if args.print_what is 1: print(('{:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} '.format(PS_mean_mean[0], PS_mean_mean[1], PS_mean_mean[2], PS_mean_mean[3], NPA_mean_mean[0], NPA_mean_mean[1], NPA_mean_mean[2], NPA_mean_mean[3]))) nfiles=len(args.fname) if args.plot_what >= 1: times=astrotime.Time(MJDs, format='mjd').datetime mtime=np.min(times) Mtime=np.max(times) if len(times) is 1: dtime=timedelta(0, 0, 0, 1) else: dtime=(Mtime-mtime)/len(times) mtime -= dtime Mtime += dtime colors=['b', 'g', 'r', 'c'] # markers=['+', 'o', 's', 'x'] fig=plt.figure(figsize=(16, 5)) ax1=plt.subplot(1, 2, 1) ngroups = North_PA.shape[0] for p in np.arange(4): ax1.errorbar(times, NPA_mean[:,p], yerr=NPA_rms[:,p], ecolor=colors[p]) ax1.plot(times, NPA_mean[:,p], colors[p]+'-', label='GV'+str(p+1)) for names in list(markers.keys()): if p == 3: label=names else: label=None indices=np.where(source_names==names) ax1.plot(times[indices], NPA_mean[indices, p].flatten(), colors[p]+markers[names], label=label) cur_marker = (cur_marker+1)%n_markers ax1.plot([mtime, Mtime], [NPA_mean_mean[p], NPA_mean_mean[p]], colors[p]+'-') ax1.plot([mtime, Mtime], [NPA_mean_mean[p], NPA_mean_mean[p]]-NPA_mean_rms[p], colors[p]+'--') ax1.plot([mtime, Mtime], [NPA_mean_mean[p], NPA_mean_mean[p]]+NPA_mean_rms[p], colors[p]+'--') plt.legend(numpoints=1, fontsize='x-small', fancybox=True, framealpha=0.5) plt.title("North PA") ax1.fmt_xdata = mdates.DateFormatter('%Y:%m:%dT%H:%M:%S.%f', tz=tzutc) ax2=plt.subplot(1, 2, 2, sharex=ax1) for p in np.arange(4): ax2.errorbar(times, scales_mean[:,p], yerr=scales_rms[:,p], ecolor=colors[p]) ax2.plot(times, scales_mean[:,p], colors[p]+'-', label='GV'+str(p+1)) for names in list(markers.keys()): if p == 3: label=names else: label=None indices=np.where(source_names==names) ax2.plot(times[indices], scales_mean[indices, p].flatten(), colors[p]+markers[names], label=label) cur_marker = (cur_marker+1)%n_markers ax2.plot([mtime, Mtime], [PS_mean_mean[p], PS_mean_mean[p]], colors[p]+'-') ax2.plot([mtime, Mtime], [PS_mean_mean[p], PS_mean_mean[p]]-PS_mean_rms[p], colors[p]+'--') ax2.plot([mtime, Mtime], [PS_mean_mean[p], PS_mean_mean[p]]+PS_mean_rms[p], colors[p]+'--') fig.autofmt_xdate() plt.title("Pixel scales") plt.legend(numpoints=1, fontsize='x-small', fancybox=True, framealpha=0.5) plt.show() if (__name__ == "__main__"): main()