#!/usr/bin/env python3 """ Author: Ole Koenig (ole.koenig@fau.de) This software is part of the High-energy Lightcurve Generator (HILIGT) software. http://xmmuls.esac.esa.int/hiligt/ This code plots a lightcurve and hardness ratios by polling the HILIGT servers through the available terminal client. Requirements: * Python3 with pandas, matplotlib, astropy * The scripts source.py and client.py must be soft-linked into this folder Usage: see ./lightcurve.py --help Example: ./lightcurve.py 187.2779 2.0523 or ./lightcurve.py 187.2779 2.0523 --mission="XMMpnt,XMMslew" --label="3C273" --legend="0,upper right" --scale=linear --band="soft" --nh=3.0E20 --specparam=2.0 Notes: The INTEGRAL bands are currently not implemented Changelog: 2020-07-08: Rewrite to use terminal client and not front-end CSV file, add hardness ratio plots """ import pandas as pd import os import matplotlib.pyplot as plt from matplotlib.lines import Line2D # from matplotlib.ticker import FormatStrFormatter # to increase ticks import numpy as np import matplotlib.cm as cm from argparse import ArgumentParser from argparse import RawTextHelpFormatter from astropy import units as u from astropy.time import Time from astropy.coordinates import Angle from client import ULS_mission from source import Source # ====== Parameters ================================================= plotsize_x = 10 # [cm] plotsize_y = 10 # [cm] upper_limit_arrow = 5 # [%] # =================================================================== allMissions = ['Vela5B', 'Ginga', 'Asca', 'Ariel5', 'XMMpnt', 'XMMslew', 'XMMstacked', 'RosatSurvey', 'RosatPointedPSPC', 'RosatPointedHRI', 'Integral', 'ExosatLE', 'ExosatME', 'Einstein', 'SwiftXRT', 'Heao1', 'Uhuru'] def main(): args = arguments() mission = args.mission if type(args.mission) == list else list(args.mission.split(",")) source = Source(args.ra, args.dec, args.label, args.model, args.specparam, args.nh, args.ulsig, args.usecat) lc = Lightcurve(source, mission) if args.noquery is False: lc.queryHILIGT() data = lc.loadData() lc.plotLightcurve(data, args.band, args.scale, args.legend, args.outfile) lc.plotHardnessratio(data) def arguments(): parser = ArgumentParser(description="This is a helper script in order to plot the lightcurves from the\n" + "High-Energy Lightcurve Generator (HILIGT).", formatter_class=RawTextHelpFormatter) # Parsed arguments for terminal client parser.add_argument("ra", nargs='?', help="right ascencion", default=0.0) parser.add_argument("dec", nargs='?', help="declination") parser.add_argument("--label", help="label for source (cannot contain whitespace)", default="") parser.add_argument("--mission", help="selected mission, choose from {}".format(",".join(allMissions)), default=allMissions) parser.add_argument("--model", help="spectral model", choices=["bbody", "plaw"], default="plaw") parser.add_argument("--specparam", help="index/temperature of spectral model", default="2.0") parser.add_argument("--nh", help="absorption column (cm^-2)", choices=["1.0E20", "3.0E20", "1.0E21"], default="3.0E20") parser.add_argument("--ulsig", help="upper limit significance", choices=["1", "2", "3"], default="2") parser.add_argument("--usecat", help="enable/disable use of catalogue values", choices=["yes", "no"], default="yes") parser.add_argument("--noquery", help="Set to disable HILIGT web call, can be used if JSON already downloaded", action='store_true') # Arguments for plotting parser.add_argument('--band', type=str, default="all", help=("Which band to plot: 'soft'=0.2-2keV, 'hard'=2-12keV, 'total'=0.2-12keV\n" "'all' will produce a multipanel plot")) parser.add_argument('--scale', type=str, choices=["linear", "log"], default="log", help="Linear or logarithmic scale") parser.add_argument('--legend', type=str, default="0,best", help="Set legend position in format panel,'position' with panel={0,1,2}") parser.add_argument('-o', '--outfile', type=str, help="Output filename; ending can be pdf,png,jpg") args = parser.parse_args() return args class Lightcurve: def __init__(self, source, mission=allMissions): """ Initializes all parameters used for the HILIGT query :param source: Source object containing ra, dec, label, model, specparam, nh, ulsig, usecat """ self.ra = source.ra() self.dec = source.dec() self.label = source.label() self.model = source.model() self.specparam = source.specparam() self.nh = source.nh() self.ulsig = source.ulsig() self.usecat = source.usecat() self.mission = mission ra = R"{0:.0f}h {1:.0f}m {2:.2f}s".format(*Angle(self.ra, unit=u.deg).hms) dec = R"{0:.0f}$^\circ$ {1:.0f}$'$ {2:.2f}$''$".format(*Angle(self.dec, unit=u.deg).dms) if self.label != "": self.title = R"{}: $\alpha$ = {}, $\delta$ = {}".format(self.label, ra, dec) else: self.title = R"$\alpha$ = {}, $\delta$ = {}".format(ra, dec) self.colormap = cm.tab20(range(len(self.mission))) # e.g., Set1, plasma, Spectral, tab20 (default) def queryHILIGT(self): """ Performs the HILIGT web query and writes the output to a file "hiligt.json". The output can be loaded with the loadData. :return: A list of the resulting JSON strings """ results = [] with open("hiligt.json", "w") as fp: for mission in self.mission: print("Querying "+mission+"...") result = ULS_mission(mission, self.ra, self.dec, self.label, self.model, self.specparam, self.nh, self.ulsig, self.usecat, "JSON", verbosity=0) results.append(result) fp.write(result.data.decode('utf-8') + "\n") def loadData(self): """ Load the HILIGT JSON output files. :return: Pandas dataframe """ # Modify the JSON file because only the individual mission queries can be loaded easily # - I'm not particularly happy with this solution -> should be better done in queryHILIGT() os.system("tr -d '\n\r' < hiligt.json > hiligt_mod.json") # delete the newlines os.system("sed -i 's/\]\[/,/g' hiligt_mod.json") # replace ][ with , to read all JSON entries at once data = pd.read_json("hiligt_mod.json", convert_dates=["_start_date", "_end_date"]) # (in the long run MJD should be added in the HILIGT backend) data["_mjd"] = Time(data['_start_date']).mjd # For convenience append the name of the flux band (how can this be done # more effectively - there must be a way to not create so many lists! bandlabels = [] for index, row in data.iterrows(): # for convenience, label the bands with "soft, hard, total" (surely this can be done better) if (row['_elow'] == 0.2) & (row['_ehigh'] == 2.0): bandlabels.append("soft") elif (row['_elow'] == 2.0) & (row['_ehigh'] == 12.0): bandlabels.append("hard") elif (row['_elow'] == 0.2) & (row['_ehigh'] == 12.0): bandlabels.append("total") else: bandlabels.append("unknown") data["_band"] = bandlabels os.remove("hiligt_mod.json") return data def _getUpperLimitArrowSize(self, data, size_percent, fluxbands): """ If several flux bands are plotted the upper limit error bars should all have the same size. However, the y-axis in the three flux band windows may have different sizes. This function goes through the pandas data frame to determine the minimum and maximum flux values. :param size_percent: Size of the upper limit error [default: 5%] :action: max(flux values)-min(flux values):math:`\cdot` , size_percent :returns: upplimsize """ upplimsize = [] for idx, band in enumerate(fluxbands): fluxmin = min(data[(data['_band'] == band) & (data['_crate_flux'] != " ")]['_crate_flux']) fluxmax = max(data[(data['_band'] == band) & (data['_crate_flux'] != " ")]['_crate_flux']) upplimsize.append((fluxmax-fluxmin)*size_percent/100) return upplimsize def plotLightcurve(self, data, bands="all", scale="log", legend="0,best", outfile=None): """ Plots the data with user defined parameters (which bands, which scale, where is the legend) using matplotlib. :param data: Pandas dataframe :param bands: Can be soft (0.2-2keV), hard (2-12keV), total (0.2-12keV), all" (default: all) :param scale: Can be log or linear (default: log) :param legend: Matplotlib position of legend with "panel,position" :param outfile: Output filename """ if bands is None or bands == "all": fluxbands = ["soft", "hard", "total"] else: fluxbands = list(bands.split(",")) # Set number of panels if len(fluxbands) > 1: fig, ax = plt.subplots(nrows=len(fluxbands), sharex=True, sharey=False, figsize=(plotsize_x, plotsize_y)) fig.subplots_adjust(hspace=0) else: ax = [None] # Hack such that axis indexing works fig, ax[0] = plt.subplots(figsize=(plotsize_x, plotsize_y)) # Second axis to display time in years ax2 = ax[0].twiny() if scale is None or scale == "log": # Set log scale ('clip' to not allow negative values, for minor ticks set subsy=[2,3,4,5,6,7,8,9] [ax[idx].set_yscale("log", nonpositive='clip') for idx in range(len(fluxbands))] if scale == "linear": upplimsize = self._getUpperLimitArrowSize(data, upper_limit_arrow, fluxbands) used_missions = [] used_colors = [] for mission, c in zip(self.mission, self.colormap): for axIdx, band in enumerate(fluxbands): # select data of current band only banddata = data[(data['_band'] == band) & (data['_mission'] == missionDict[mission]) & (data['_status'] != "NoDataFound")].copy() if len(banddata) == 0: # Continue if a mission doesn't have data at this band continue datapoints = banddata[banddata['_crate'].notnull()] ax[axIdx].errorbar(datapoints['_mjd'], datapoints['_crate_flux'], yerr=datapoints['_crate_flux_err'], fmt='.', color=c, markeredgewidth=2, capthick=2, elinewidth=2, linestyle='None') upperlimits = banddata[(banddata['_ul'].notnull()) & (banddata['_ul'] > 0)] if len(upperlimits) > 0: if scale == "linear": ax[axIdx].errorbar(upperlimits['_mjd'], upperlimits['_ul_flux'], yerr=upplimsize[axIdx], color=c, uplims=True, linestyle="None") else: # A constant has a length of (exp(flux)-1)/scaling on log scale ax[axIdx].errorbar(upperlimits['_mjd'], upperlimits['_ul_flux'], yerr=(np.exp(list(upperlimits['_ul_flux']))-1)/upper_limit_arrow, color=c, uplims=True, linestyle="None") # For legend plotting later if mission not in used_missions: used_missions.append(mission) used_colors.append(c) # HACK to display years if axIdx == 0: ax2.plot(datapoints['_start_date'], datapoints['_crate_flux'], alpha=0) ax[axIdx].set_ylabel(R"%s keV" % fluxDict[band]) # ax[axIdx].grid(True,alpha=0.1) # ATTENTION: autofmt_xdate could kills some labels # fig.autofmt_xdate() # should fix the tic label positions ax[0].set_title(self.title) ax[len(fluxbands)-1].set_xlabel(R"MJD (d)") ax2.set_xlabel(r"Date") legend_panel_pos = int(legend.split(",")[0]) legend_pos = legend.split(",")[1] custom_lines = [Line2D([0], [0], color=c, lw=2) for c in used_colors] ax[legend_panel_pos].legend(custom_lines, used_missions, loc=legend_pos) if outfile is None: outfile = f"hiligt_lc_{self.label}_{self.ra}_{self.dec}.pdf" fig.savefig(outfile, bbox_inches='tight') print("Saving to", outfile) def plotHardnessratio(self, data): """ Plot of the hardness ratio (2-12 flux - 0.2-2 flux ) / (2-12 flux + 0.2-2 flux). Only works for sources where sufficient data in all bands is available. """ fig, ax = plt.subplots(figsize=(10, 7)) ax.set_title(R"Hardness Ratio of " + self.title) ax.set_ylabel(R"(H-S)/(H+S) : hard={} keV, soft={} keV".format(fluxDict["hard"], fluxDict["soft"])) ax.set_xlabel(R"MJD (d)") ax2 = ax.twiny() ax2.set_xlabel(R"Date") used_missions = [] used_colors = [] for mission, c in zip(self.mission, self.colormap): softdata = data[(data['_band'] == "soft") & (data['_mission'] == missionDict[mission]) & (data['_status'] != "NoDataFound") & (data['_crate'].notnull())].copy() harddata = data[(data['_band'] == "hard") & (data['_mission'] == missionDict[mission]) & (data['_status'] != "NoDataFound") & (data['_crate'].notnull())].copy() # Continue, if a mission doesn't have data in both bands if len(softdata) == 0 or len(harddata) == 0 or len(softdata) != len(harddata): continue mjds = np.array(softdata['_mjd'], dtype=np.float64) hard = np.array(harddata['_crate_flux'], dtype=np.float64) dhard = np.array(harddata['_crate_flux_err'], dtype=np.float64) soft = np.array(softdata['_crate_flux'], dtype=np.float64) dsoft = np.array(softdata['_crate_flux_err'], dtype=np.float64) hr = (hard-soft)/(hard+soft) dhr = np.sqrt((2*soft/(hard+soft)**2)**2 * dhard**2 + (2*hard/(hard+soft)**2)**2 * dsoft**2) ax.errorbar(mjds, hr, yerr=dhr, linestyle='None', fmt=".") ax2.plot(softdata['_start_date'], hr, alpha=0) if mission not in used_missions: used_missions.append(mission) used_colors.append(c) custom_lines = [Line2D([0], [0], color=c, lw=2) for c in used_colors] ax.legend(custom_lines, used_missions) outfile = f"hiligt_hr_{self.label}_{self.ra}_{self.dec}.pdf" fig.savefig(outfile, bbox_inches='tight') print("Saving to", outfile) missionDict = dict(Vela5B="Vela5B", Ginga="GINGA-LAC", Asca="ASCA", Ariel5="ARIEL5", XMMpnt="XMM-Newton pointed", XMMslew="XMM-Newton slew", XMMstacked="XMM-Newton stacked", RosatSurvey="ROSAT-Survey", RosatPointedPSPC="ROSAT-pointed-PSPC", RosatPointedHRI="ROSAT-pointed-HRI", Integral="Integral", ExosatLE="EXOSAT-LE", ExosatME="EXOSAT-ME", Einstein="EINSTEIN", SwiftXRT="Swift-XRT", Heao1="HEAO-1", Uhuru="Uhuru") fluxDict = dict(soft="Flux 0.2 - 2", hard="Flux 2 - 12", total="Flux 0.2 - 12") # for y-labels of plot if __name__ == "__main__": main()