diff --git a/examples/make_pointings.py b/examples/make_pointings.py index 08cc4b3..3004983 100755 --- a/examples/make_pointings.py +++ b/examples/make_pointings.py @@ -72,8 +72,7 @@ def get_positions(xmin, xmax, ymin, ymax, step): # is invariant of latitude. xstep = step / math.cos(gammalib.deg2rad * y) nx = int((xmax-xmin)/xstep+1.5) - rescale = (xmax-xmin)/(xstep*(nx-1)) - xstep *= rescale + xstep = (xmax-xmin)/float(nx-1) # Set x offset. For every second row the x position is displace by # half a step size. @@ -427,9 +426,9 @@ def set_extgal(separation=3.0, caldb='prod2', lst=True): # Set patch obsdef.extend(set_patch(tmin, lmin=-90.0, lmax=90.0, bmin=+5.0, bmax=+88.0, - separation=separation, hours=1000, + separation=separation, hours=500, run_duration=25, site='Automatic', caldb=caldb, lst=lst, - autodec=10.0)) + autodec=-10.0)) # Return observation definition return obsdef @@ -501,6 +500,9 @@ def set_lmc(hours=250.0, caldb='prod2', lst=True): # Initialise observation definition obsdef = [] + # Initialise start time (1-1-2021) + tmin = 7671.0 * 86400.0 + # Set LMC centre centre = gammalib.GSkyDir() centre.radec_deg(80.0, -69.0) @@ -523,7 +525,7 @@ def set_lmc(hours=250.0, caldb='prod2', lst=True): lat = pnt.b_deg() # Set positions and duration - obs = {'lon': lon, 'lat': lat, 'duration': duration} + obs = {'lon': lon, 'lat': lat, 'tmin': tmin, 'duration': duration} # Set IRF irf = set_irf('South', obs, caldb, lst=lst) @@ -535,6 +537,9 @@ def set_lmc(hours=250.0, caldb='prod2', lst=True): # Append observation obsdef.append(obs) + # Update start time for next pointing + tmin = set_tmin_for_next_pointing(tmin, duration) + # Dump statistics print('Number of pointings: %d (%.2f s)' % (n_pnt,duration)) diff --git a/examples/show_pointings.py b/examples/show_pointings.py new file mode 100755 index 0000000..f662bf1 --- /dev/null +++ b/examples/show_pointings.py @@ -0,0 +1,186 @@ +#! /usr/bin/env python +# ========================================================================== +# Display pointings in observation definition XML file +# +# Copyright (C) 2017 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import sys +import gammalib +import cscripts +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except (ImportError, RuntimeError): + print('This script needs the "matplotlib" module') + sys.exit() + + +# =============================== # +# Extract pointings from XML file # +# =============================== # +def get_pointings(filename): + """ + Extract pointings from XML file + + Parameters + ---------- + filename : str + File name of observation definition XML file + + Returns + ------- + pnt : list of dict + Pointings + """ + # Initialise pointings + pnt = [] + + # Open XML file + xml = gammalib.GXml(filename) + + # Get observation list + obs = xml.element('observation_list') + + # Get number of observations + nobs = obs.elements('observation') + + # Loop over observations + for i in range(nobs): + + # Get observation + run = obs.element('observation', i) + + # Get pointing parameter + npars = run.elements('parameter') + ra = None + dec = None + south = False + for k in range(npars): + par = run.element('parameter', k) + if par.attribute('name') == 'Pointing': + ra = float(par.attribute('ra')) + dec = float(par.attribute('dec')) + if par.attribute('name') == 'Calibration': + if 'South' in par.attribute('response'): + south = True + else: + south = False + if ra != None: + p = gammalib.GSkyDir() + p.radec_deg(ra, dec) + entry = {'l': p.l_deg(), 'b': p.b_deg(), 'ra': ra, 'dec': dec, 'south': south} + pnt.append(entry) + + # Return pointings + return pnt + + +# ================ # +# Plot information # +# ================ # +def plot_pointings(pnt, plotfile): + """ + Plot information + + Parameters + ---------- + pnt : list of dict + Pointings + plotfile : str + Plot filename + """ + # Create figure + plt.figure() + + # Setup figure + ax = plt.gca() + ax.cla() + ax.set_xlim((180, -180)) + ax.set_ylim((-90, 90)) + + # Loop over pointings + for p in pnt: + + # Get longitude + l = p['l'] + if l > 180.0: + l = l - 360.0 + + # Set color + if p['south']: + color='r' + else: + color='b' + + # Set circle + circle = plt.Circle((l, p['b']), 4.0, color=color, fill=False) + + # Add circle + ax.add_artist(circle) + + # Plot title and labels + plt.xlabel('Galactic longitude (deg)') + plt.ylabel('Galactic latitude (deg)') + plt.title('Pointings') + + # Show plots or save it into file + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return + + +# ============== # +# Show pointings # +# ============== # +def show_pointings(): + """ + Show pointings + """ + # Set usage string + usage = 'show_pointings.py [-p plotfile] file' + + # Set default options + options = [{'option': '-p', 'value': ''}] + + # Get arguments and options from command line arguments + args, options = cscripts.ioutils.get_args_options(options, usage) + + # Extract script parameters from options + plotfile = options[0]['value'] + + # Get pointings + pnt = get_pointings(args[0]) + + # Plot pointings + plot_pointings(pnt, plotfile) + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Show pointings + show_pointings() diff --git a/test/test_examples.py.in b/test/test_examples.py.in index 6c6cab8..8dde648 100755 --- a/test/test_examples.py.in +++ b/test/test_examples.py.in @@ -134,6 +134,8 @@ class Test(gammalib.GPythonTestSuite): 'Test ts_distribution.py') self.append(self.test_show_obs, 'Test show_obs.py') + self.append(self.test_show_pointings, + 'Test show_pointings.py') # Return return @@ -460,6 +462,20 @@ class Test(gammalib.GPythonTestSuite): # Return return + # Test show_pointings + def test_show_pointings(self): + """ + Test show_pointings + """ + # Set arguments + args = '-p example_show_pointings.eps %s/obs_unbinned.xml' % self._datadir + + # Execute script + self._execute_python('show_pointings', args=args) + + # Return + return + # ======================== # # Main routine entry point #