Add unit test of make_pointings.py script

Jürgen Knödlseder
1 parent 5947039a
Showing 2 changed files with 195 additions and 64 deletions
examples/make_pointings.py
test/test_examples.py
 #! /usr/bin/env python
 # ==========================================================================
-# This script generates pointing patterns for CTA observations.
+# Generates pointing patterns for CTA observations
 #
-# Usage:
-#   ./make_pointings.py
-#
-# ==========================================================================
-#
-# Copyright (C) 2015 Juergen Knoedlseder
+# Copyright (C) 2015-2016 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
@@ -33,10 +28,30 @@ import gammalib
 # ============================================================ #
 def get_positions(xmin, xmax, ymin, ymax, step):
     """
+    Get pointing positions for an interleaved pattern
+
     Get positions for a patch of interleaved pointings on the sky.
     The function does respect the latitude dependence of the
     pointing distance. The step parameter is for a latitude of
     zero.
+
+    Parameters
+    ----------
+    xmin : float
+        Minimum longitude (deg)
+    xmax : float
+        Maximum longitude (deg)
+    ymin : float
+        Minimum latitude (deg)
+    ymax : float
+        Maximum latitude (deg)
+    step : float
+        Pointing step size (deg)
+
+    Returns
+    -------
+    positions : list of dict
+        List of dictionaries with x and y positions
     """
     # Initialise positions
     positions = []
@@ -46,7 +61,7 @@ def get_positions(xmin, xmax, ymin, ymax, step):
     ny      = int((ymax-ymin)/ystep+1.5)
     rescale = (ymax-ymin)/(ystep*(ny-1))
     ystep  *= rescale
-    print("Number of rows ....: %d" % (ny))
+    print('Number of rows ....: %d' % (ny))
     # Loop over rows
     y = ymin
@@ -63,7 +78,6 @@ def get_positions(xmin, xmax, ymin, ymax, step):
             x = xmin
         else:
             x = xmin + 0.5 * xstep
-        #print(row, y, nx)
         # Append pointings
         for pnt in range(nx):
@@ -86,7 +100,33 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
               separation=3.0, hours=100.0,
               site=None, lst=True, autodec=0.0):
     """
-    Setup pointing patch on the sky.
+    Setup pointing patch on the sky
+
+    Parameters
+    ----------
+    lmin : float, optional
+        Minimum Galactic Longitude (deg)
+    lmax : float, optional
+        Maximum Galactic Longitude (deg)
+    bmin : float, optional
+        Minimum Galactic Latitude (deg)
+    bmax : float, optional
+        Maximum Galactic Latitude (deg)
+    separation : float, optional
+        Pointing separation (deg)
+    hours : float, optional
+        Total observation duration (h)
+    site : str, optional
+        Array site (one of 'Automatic', 'South', 'North')
+    lst : bool, optional
+        Use LSTs
+    autodec : float, optional
+        Declination for switching automatically between South and North site
+
+    Returns
+    -------
+    obsdef : list of dict
+        List of pointing definitions
     """
     # Initialise observation definition
     obsdef = []
@@ -112,30 +152,30 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
         if site != None:
             # Handle automatic site switching
-            if site == "Automatic":
+            if site == 'Automatic':
                 pos = gammalib.GSkyDir()
                 pos.lb_deg(pnt['x'], pnt['y'])
                 dec = pos.dec_deg()
                 if (dec >= autodec):
-                    obs_site = "North"
+                    obs_site = 'North'
                 else:
-                    obs_site = "South"
+                    obs_site = 'South'
             else:
                 obs_site = site
             # Set site IRF
-            caldb = "prod2"
-            if obs_site == "North":
+            caldb = 'prod2'
+            if obs_site == 'North':
                 if lst:
-                    irf = "North_50h"
+                    irf = 'North_50h'
                 else:
-                    irf = "North_50h"
+                    irf = 'North_50h'
                 exposure_north += duration
             else:
                 if lst:
-                    irf = "South_50h"
+                    irf = 'South_50h'
                 else:
-                    irf = "South_50h"
+                    irf = 'South_50h'
                 exposure_south += duration
             # Set IRF information
@@ -146,11 +186,11 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
         obsdef.append(obs)
     # Dump statistics
-    print("Number of pointings: %d (%.2f sec)" % (n_pnt,duration))
-    print("South array .......: %.2f hours (%.2f%%)" %
+    print('Number of pointings: %d (%.2f s)' % (n_pnt,duration))
+    print('South array .......: %.2f hours (%.2f%%)' %
           (exposure_south/3600.0,
            100.0*exposure_south/(exposure_south+exposure_north)))
-    print("North array .......: %.2f hours (%.2f%%)" %
+    print('North array .......: %.2f hours (%.2f%%)' %
           (exposure_north/3600.0,
            100.0*exposure_north/(exposure_south+exposure_north)))
@@ -163,7 +203,23 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
 # =========================== #
 def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
     """
-    Setup Galactic plane survey.
+    Setup Galactic plane survey
+
+    Parameters
+    ----------
+    separation : float, optional
+        Pointing separation (deg)
+    bmin : float, optional
+        Minimum Galactic Latitude (deg)
+    bmax : float, optional
+        Maximum Galactic Latitude (deg)
+    lst : bool, optional
+        Use LSTs
+
+    Returns
+    -------
+    obsdef : list of dict
+        List of pointing definitions
     """
     # Initialise observation definition
     obsdef = []
@@ -171,27 +227,27 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
     # Add inner region South
     obsdef.extend(set_patch(lmin=-60.0, lmax=60.0, bmin=bmin, bmax=bmax,
                             separation=separation, hours=780,
-                            site="South", lst=lst))
+                            site='South', lst=lst))
     # Add Vela & Carina region
     obsdef.extend(set_patch(lmin=240.0, lmax=300.0, bmin=bmin, bmax=bmax,
                             separation=separation, hours=180,
-                            site="South", lst=lst))
+                            site='South', lst=lst))
     # Add 210-240 region
     obsdef.extend(set_patch(lmin=210.0, lmax=240.0, bmin=bmin, bmax=bmax,
                             separation=separation, hours=60,
-                            site="South", lst=lst))
+                            site='South', lst=lst))
     # Add Cygnus, Perseus
     obsdef.extend(set_patch(lmin=60.0, lmax=150.0, bmin=bmin, bmax=bmax,
                             separation=separation, hours=450,
-                            site="North", lst=lst))
+                            site='North', lst=lst))
     # Add Anticentre
     obsdef.extend(set_patch(lmin=150.0, lmax=210.0, bmin=bmin, bmax=bmax,
                             separation=separation, hours=150,
-                            site="North", lst=lst))
+                            site='North', lst=lst))
     # Return observation definition
     return obsdef
@@ -203,6 +259,18 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
 def set_extgal(separation=3.0, lst=True):
     """
     Setup Extragalactic survey.
+
+    Parameters
+    ----------
+    separation : float, optional
+        Pointing separation (deg)
+    lst : bool, optional
+        Use LSTs
+
+    Returns
+    -------
+    obsdef : list of dict
+        List of pointing definitions
     """
     # Initialise observation definition
     obsdef = []
@@ -210,7 +278,7 @@ def set_extgal(separation=3.0, lst=True):
     # Set patch
     obsdef.extend(set_patch(lmin=-90.0, lmax=90.0, bmin=+5.0, bmax=+88.0,
                             separation=separation, hours=1000,
-                            site="Automatic", lst=lst, autodec=10.0))
+                            site='Automatic', lst=lst, autodec=10.0))
     # Return observation definition
     return obsdef
@@ -222,6 +290,16 @@ def set_extgal(separation=3.0, lst=True):
 def set_gc(lst=True):
     """
     Setup Galactic centre KSP.
+
+    Parameters
+    ----------
+    lst : bool, optional
+        Use LSTs
+
+    Returns
+    -------
+    obsdef : list of dict
+        List of pointing definitions
     """
     # Initialise observation definition
     obsdef = []
@@ -229,12 +307,12 @@ def set_gc(lst=True):
     # Central wobble
     obsdef.extend(set_patch(lmin=-0.5, lmax=0.5, bmin=-0.5, bmax=0.5,
                             separation=0.1, hours=525,
-                            site="South", lst=lst))
+                            site='South', lst=lst))
     # Extended region
     obsdef.extend(set_patch(lmin=-10.0, lmax=10.0, bmin=-10.0, bmax=10.0,
                             separation=1.5, hours=300,
-                            site="South", lst=lst))
+                            site='South', lst=lst))
     # Return observation definition
     return obsdef
@@ -246,6 +324,18 @@ def set_gc(lst=True):
 def set_lmc(hours=250.0, lst=True):
     """
     Setup LMC KSP.
+
+    Parameters
+    ----------
+    hours : float, optional
+        Total observation duration (h)
+    lst : bool, optional
+        Use LSTs
+
+    Returns
+    -------
+    obsdef : list of dict
+        List of pointing definitions
     """
     # Initialise observation definition
     obsdef = []
@@ -275,11 +365,11 @@ def set_lmc(hours=250.0, lst=True):
         obs = {'lon': lon, 'lat': lat, 'duration': duration}
         # Add IRF
-        caldb = "prod2"
+        caldb = 'prod2'
         if lst:
-            irf = "South_50h"
+            irf = 'South_50h'
         else:
-            irf = "South_50h"
+            irf = 'South_50h'
         obs['caldb'] = caldb
         obs['irf']   = irf
@@ -287,7 +377,7 @@ def set_lmc(hours=250.0, lst=True):
         obsdef.append(obs)
     # Dump statistics
-    print("Number of pointings: %d (%.2f sec)" % (n_pnt,duration))
+    print('Number of pointings: %d (%.2f s)' % (n_pnt,duration))
     # Return observation definition
     return obsdef
@@ -298,31 +388,38 @@ def set_lmc(hours=250.0, lst=True):
 # ======================================= #
 def write_obsdef(filename, obsdef):
     """
-    Write observation definition file.
+    Write observation definition file
+
+    Parameters
+    ----------
+    filename : str
+        Observation definition file name
+    obsdef : list of dict
+        List of pointing definitions
     """
     # Open file
     file = open(filename, 'w')
     # Write header
-    file.write("ra,dec,duration,caldb,irf\n")
+    file.write('ra,dec,duration,caldb,irf\n')
     # Loop over pointings
     for obs in obsdef:
         # If we have lon,lat then convert into RA,Dec
-        if "lon" in obs and "lat" in obs:
-            lon = obs["lon"]
-            lat = obs["lat"]
+        if 'lon' in obs and 'lat' in obs:
+            lon = obs['lon']
+            lat = obs['lat']
             dir = gammalib.GSkyDir()
             dir.lb_deg(lon,lat)
             ra  = dir.ra_deg()
             dec = dir.dec_deg()
         else:
-            ra  = obs["ra"]
-            dec = obs["dec"]
+            ra  = obs['ra']
+            dec = obs['dec']
         # Write information
-        file.write("%8.4f,%8.4f,%.4f,%s,%s\n" %
+        file.write('%8.4f,%8.4f,%.4f,%s,%s\n' %
                    (ra, dec, obs['duration'], obs['caldb'], obs['irf']))
     # Close file
@@ -342,7 +439,7 @@ if __name__ == &#39;__main__&#39;:
     # Test for command line arguments
     if (len(sys.argv) > 1):
-        if sys.argv[1] == "-h":
+        if sys.argv[1] == '-h':
             need_help = True
         else:
             obsname = sys.argv[1]
@@ -351,37 +448,37 @@ if __name__ == &#39;__main__&#39;:
     # Print help if needed and exit
     if need_help:
-        print("Usage: make_pointing.py [OPTIONS]")
-        print("     -h       Display this usage message")
-        print("     gps      Galactic plane survey (2 row scheme)")
-        print("     gps3     Galactic plane survey (3 row scheme)")
-        print("     extgal   Extragalactic survey")
-        print("     gc       Galactic centre survey")
-        print("     lmc      LMC survey")
+        print('Usage: make_pointing.py [OPTIONS]')
+        print('     -h       Display this usage message')
+        print('     gps      Galactic plane survey (2 row scheme)')
+        print('     gps3     Galactic plane survey (3 row scheme)')
+        print('     extgal   Extragalactic survey')
+        print('     gc       Galactic centre survey')
+        print('     lmc      LMC survey')
         sys.exit()
     # Galactic plane survey
-    if obsname == "gps":
+    if obsname == 'gps':
         obsdef = set_gps(lst=True)
-        write_obsdef("gps.dat", obsdef)
-    elif obsname == "gps3":
+        write_obsdef('gps.dat', obsdef)
+    elif obsname == 'gps3':
         obsdef = set_gps(separation=1.5, lst=True)
-        write_obsdef("gps3.dat", obsdef)
+        write_obsdef('gps3.dat', obsdef)
     # Extragalactic survey
-    elif obsname == "extgal":
+    elif obsname == 'extgal':
         obsdef = set_extgal(lst=True)
-        write_obsdef("extgal.dat", obsdef)
+        write_obsdef('extgal.dat', obsdef)
     # Galactic centre
-    elif obsname == "gc":
+    elif obsname == 'gc':
         obsdef = set_gc(lst=True)
-        write_obsdef("gc.dat", obsdef)
+        write_obsdef('gc.dat', obsdef)
     # LMC
-    elif obsname == "lmc":
+    elif obsname == 'lmc':
         obsdef = set_lmc(lst=True)
-        write_obsdef("lmc.dat", obsdef)
+        write_obsdef('lmc.dat', obsdef)
     # Invalid pattern
     else:
@@ -44,21 +44,33 @@ class Test(gammalib.GPythonTestSuite):
         return
     # Execute Python script
-    def _execute_python(self, name):
+    def _execute_python(self, name, args='', log=''):
         """
         Execute Python script
         Parameters
         ----------
         name : str
-            Python script name with .py extension.
+            Python script name with .py extension
+        args : str, optional
+            String with arguments to be passed to script
+        log : str, optional
+            String with name to be used for logging
         """
         # Setup command
         cmd = '../examples/' + name + '.py'
+        if len(args) > 0:
+            cmd += ' ' + args
+
+        # Set logname
+        if len(log) > 0:
+            logname = log
+        else:
+            logname = name
         # Execute Python script, make sure we catch any exception
         try:
-            rc = os.system(cmd+' > example_'+name+'.log 2>&1')
+            rc = os.system(cmd+' > example_'+logname+'.log 2>&1')
         except:
             pass
@@ -78,6 +90,7 @@ class Test(gammalib.GPythonTestSuite):
         # Append tests
         self.append(self.test_generate_prod3_irfs, 'Test generate_prod3_irfs')
+        self.append(self.test_make_pointings, 'Test make_pointings')
         self.append(self.test_pipeline_binned_disk, 'Test pipeline_binned_disk')
         self.append(self.test_pipeline_binned_mem, 'Test pipeline_binned_mem')
         self.append(self.test_pipeline_stacked_disk, 'Test pipeline_stacked_disk')
@@ -88,6 +101,27 @@ class Test(gammalib.GPythonTestSuite):
         # Return
         return
+    # Test make_pointings
+    def test_make_pointings(self):
+        """
+        Test make_pointings
+        """
+        # Set script name
+        script = 'make_pointings'
+
+        # Execute script
+        self._execute_python(script, args='-h', log='make_pointings_help')
+        self._execute_python(script, args='gps', log='make_pointings_gps')
+        self._execute_python(script, args='gps3', log='make_pointings_gps3')
+        self._execute_python(script, args='extgal', log='make_pointings_extgal')
+        self._execute_python(script, args='gc', log='make_pointings_gc')
+        self._execute_python(script, args='lmc', log='make_pointings_lmc')
+
+        #TODO: Do any testing
+
+        # Return
+        return
+
     # Test generate_prod3_irfs
     def test_generate_prod3_irfs(self):
         """