Add initial csviscube version

Jürgen Knödlseder
1 parent b815cccb
Showing 6 changed files with 441 additions and 0 deletions
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/csviscube.par
cscripts/csviscube.py
test/test_csviscube.py
test/test_python_cscripts.py
@@ -37,6 +37,7 @@ pythonscripts = $(builddir)/__init__.py \
                 $(srcdir)/cstsdist.py \
                 $(srcdir)/cstsmapmerge.py \
                 $(srcdir)/cstsmapsplit.py \
+                $(srcdir)/csviscube.py \
                 $(srcdir)/csworkflow.py \
                 $(srcdir)/csfindobs.py \
                 $(srcdir)/csiactcopy.py \
@@ -63,6 +64,7 @@ parfiles = $(srcdir)/cscaldb.par \
            $(srcdir)/cstsdist.par \
            $(srcdir)/cstsmapmerge.par \
            $(srcdir)/cstsmapsplit.par \
+           $(srcdir)/csviscube.par \
            $(srcdir)/csworkflow.par \
            $(srcdir)/csfindobs.par \
            $(srcdir)/csiactcopy.par \
@@ -86,6 +88,7 @@ linkscripts = cscaldb \
               cstsdist \
               cstsmapmerge \
               cstsmapsplit \
+              csviscube \
               csworkflow \
               csfindobs \
               csiactcopy \
@@ -111,6 +114,7 @@ testscripts = $(srcdir)/tests/__init__.py \
               $(top_srcdir)/test/test_cstsdist.py \
               $(top_srcdir)/test/test_cstsmapmerge.py \
               $(top_srcdir)/test/test_cstsmapsplit.py \
+              $(top_srcdir)/test/test_csviscube.py \
               $(top_srcdir)/test/test_csworkflow.py \
               $(top_srcdir)/test/test_csfindobs.py \
               $(top_srcdir)/test/test_csiactcopy.py \
@@ -40,6 +40,7 @@ __all__ = [
     "cstsdist",
     "cstsmapmerge",
     "cstsmapsplit",
+    "csviscube",
     "csworkflow",
     "obsutils",
     "modutils",
@@ -69,6 +70,7 @@ from cscripts.csspec       import csspec
 from cscripts.cstsdist     import cstsdist
 from cscripts.cstsmapmerge import cstsmapmerge
 from cscripts.cstsmapsplit import cstsmapsplit
+from cscripts.csviscube    import csviscube
 from cscripts.csworkflow   import csworkflow
 from cscripts              import obsutils
 from cscripts              import modutils
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      csviscube.par                                                 #
+#  Version:   1.0.0                                                         #
+#                                                                           #
+#  Author:    Juergen Knoedlseder                                           #
+#             Jurgen.Knodlseder@irap.omp.eu                                 #
+#             IRAP                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the csviscube script                       #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+tmin,       r, a,   0.0,,, "Start time (CTA MET in seconds)"
+tmax,       r, a,  31557600.0,,, "End time (CTA MET in seconds)"
+geolon,     r, a,  79.4041,,, "Geographic longitude of array (deg)"
+geolat,     r, a, -24.6272,,, "Geographic latitude of array (deg)"
+sunzenith,  r, h, 105.0,,, "Minimum Sun zenith angle for observations (deg)"
+moonzenith, r, h, 90.0,,, "Minimum Moon zenith angle for observations (deg)"
+outfile,    f, a, viscube.fits,,, "Visibility cube output file"
+
+#
+# Visibility cube parameters
+#===========================
+binsz, r, h, 1.0,,, "Image scale (in degrees/pixel)"
+
+#
+# Standard parameters
+#====================
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug,   b, h, no,,, "Debugging mode activated"
+mode,    s, h, ql,,, "Mode of automatic parameters"
+logfile, f, h, csviscube.log,,, "Log filename"
+#! /usr/bin/env python
+# ==========================================================================
+# Compute a visibility cube
+#
+# Copyright (C) 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
+# 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import sys
+import glob
+import os
+import math
+import gammalib
+import ctools
+
+
+# =============== #
+# csviscube class #
+# =============== #
+class csviscube(ctools.cscript):
+    """
+    Compute a visibility cube
+    """
+
+    # Constructor
+    def __init__(self, *argv):
+        """
+        Constructor.
+        
+        Parameters
+        ----------
+        argv : list of str
+            List of IRAF command line parameter strings of the form
+            ``parameter=3``.
+
+        Raises
+        ------
+        TypeError
+            An invalid number of command line arguments was provided.
+        """
+        # Set name
+        self._name    = 'csviscube'
+        self._version = '1.0.0'
+
+        # Initialise application by calling the appropriate class
+        # constructor.
+        self._init_cscript(argv)
+
+        # Return
+        return
+
+
+    # Private methods
+    def _get_parameters(self):
+        """
+        Get all parameters
+        """
+        # Query parameters
+        self['tmin'].real()
+        self['tmax'].real()
+        self['geolon'].real()
+        self['geolat'].real()
+        self['sunzenith'].real()
+        self['moonzenith'].real()
+        self['outfile'].filename()
+        self['binsz'].real()
+
+        #  Write input parameters into logger
+        self._log_parameters(gammalib.TERSE)
+
+        # Return
+        return
+
+    def _hour_angle_weight(self):
+        """
+        Compute hour angle weight vector
+        
+        Compute during which time the array was observing due to hour angle
+        constraints. We take here a very simple model where the night lasts
+        between 6 and 10 hours (or 90-150 degrees), with a sinusoidal variation
+        along the year. Every day the hour angle of the Sun set is advancing by
+        4 minutes.
+        """
+        # Set number of hour angle bins
+        nsteps = 1000
+        d_min  =  6.0/24.0*nsteps
+        d_max  = 10.0/24.0*nsteps
+    
+        # Initialise hour angle array and weight
+        hours  = [0.0 for i in range(nsteps)]
+        weight = 24.0/float(nsteps)  # Weight per hour angle step
+
+        # Loop over all days in one year
+        for day in range(nsteps):
+            h_set  = day
+            h_rise = day + \
+                     int(d_min + (d_max-d_min)*math.sin(day/float(nsteps)*gammalib.pi))
+            for i in range(h_set,h_rise):
+                if i >= nsteps:
+                    hours[i-nsteps] += weight
+                else:
+                    hours[i] += weight
+
+        # Return
+        return hours
+
+    def _visibility_cube(self):
+        """
+        Compute visibility cube
+        
+        Compute visibility cube by displacing the zenith angle map for all
+        hour angles. The visibility cube contains the number of hours a given
+        celestial position is visible under a given zenith angle interval.
+        Summing over all zenith angle intervals specifies for how long a given
+        celestial position will be visible.
+        """
+        # Compute zenith angle map
+        map = self._zenith_angle_map()
+        
+        # Set zenith angle range and bin size
+        dz = 1.0
+        nz = int(60.0/dz)
+    
+        # Initialise visibility cube
+        self._cube = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180,nz)
+
+        # Setup hour angle weights. They specify for how many hours a given
+        # hour angle will be observed during one year.
+        hours = self._hour_angle_weight()
+
+        # Compute normalisation factor
+        norm = 365.0/float(len(hours))
+
+        # Loop over all hour angle weights
+        for h, weight in enumerate(hours):
+
+            # Compute temporal shift
+            shift = 360.0/float(len(hours)) * float(h)
+
+            # Initialise shifted zenith angle map
+            map_shift = gammalib.GSkyMap('CAR','CEL',shift,0.0,-1.0,1.0,360,180)
+
+            # Merge zenith angle map into shifted map
+            map_shift += map
+
+            # Loop over all pixels of shifted map
+            for i, zenith in enumerate(map_shift):
+                iz     = int(zenith/dz)
+                if iz < nz:
+                    self._cube[i,iz] += weight * norm
+
+        # Return
+        return
+
+    def _zenith_angle_map(self):
+        """
+        Compute zenith angle map
+
+        The zenith angle of a position (ra,dec) depends on the declination and
+        the hour angle h and is given by
+
+        zenith(h,dec) = arccos( sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(h) )
+
+        The hour angle h (or local hour angle, LHA) is defined as the difference
+        between local siderial time (LST) and the Right Ascension
+    
+        h = LST - ra
+
+        LST is the difference between the Greenwich siderial time (GST) and the
+        geographic longitude of the observer
+    
+        LST = GST - lon
+    
+        where lon is counted positive west from the meridian. Hence
+    
+        h = GST - lon - ra
+
+        The map is computed for h=-ra which is equivalent to GST=lon (or LST=0).
+        In other words, the map corresponds to the time when ra=0 passes the
+        local meridian.
+        """
+        # Initialise zenith angle map
+        map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180)
+
+        # Set hour angle and declination vectors
+        hours = [float(i) for i in range(360)]
+        decs  = [float(i)-89.5 for i in range(180)]
+
+        # Get array geographic longitude and latitude
+        geolon = self['geolon'].real()
+        geolat = self['geolat'].real()
+
+        # Precompute latitude terms
+        cos_lat = math.cos(geolat*gammalib.deg2rad)
+        sin_lat = math.sin(geolat*gammalib.deg2rad)
+
+        # Loop over all declination and hour angles and compute the zenith
+        # angle. Store the zenith angle in the map
+        index = 0
+        for dec in decs:
+            cos_dec = math.cos(dec*gammalib.deg2rad)
+            sin_dec = math.sin(dec*gammalib.deg2rad)
+            for h in hours:
+                cos_h  = math.cos(h*gammalib.deg2rad)
+                zenith = math.acos(sin_lat*sin_dec +
+                                   cos_lat*cos_dec*cos_h)*gammalib.rad2deg
+                map[index] = zenith
+                index     += 1
+
+        # Return zenith angle map
+        return map
+
+
+    # Public methods
+    def run(self):
+        """
+        Run the script
+        """
+        # Switch screen logging on in debug mode
+        if self._logDebug():
+            self._log.cout(True)
+
+        # Get parameters
+        self._get_parameters()
+
+        # Compute visibility cube
+        self._visibility_cube()
+
+        # Return
+        return
+
+    def save(self):
+        """
+        Save the visibility cube
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save visibility cube')
+
+        # Get outfile parameter
+        outfile = self['outfile'].filename()
+        
+        # Log file name
+        self._log_value(gammalib.TERSE, 'Visibility cube file', outfile.url())
+
+        # Save the visibility cube
+        self._cube.save(outfile, self['clobber'].boolean())
+
+        # Return
+        return
+
+    def execute(self):
+        """
+        Execute the script
+        """
+        # Open logfile
+        self.logFileOpen()
+
+        # Run the script
+        self.run()
+
+        # Save the cube
+        self.save()
+
+        # Return
+        return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+    # Create instance of application
+    app = csviscube(sys.argv)
+
+    # Execute application
+    app.execute()
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the csviscube script
+#
+# Copyright (C) 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
+# 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import cscripts
+from testing import test
+
+
+# =============================== #
+# Test class for csviscube script #
+# =============================== #
+class Test(test):
+    """
+    Test class for csviscube script
+    """
+
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        test.__init__(self)
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name('csviscube')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test csviscube on command line')
+        self.append(self._test_python, 'Test csviscube from Python')
+
+        # Return
+        return
+
+    # Test csviscube on command line
+    def _test_cmd(self):
+        """
+        Test csviscube on the command line.
+        """
+        # Set script name
+        csviscube = self._script('csviscube')
+
+        # Setup csviscube command
+        cmd = csviscube+' tmin=0.0 tmax=31557600.0 geolon=79.4041'+\
+                        ' geolat=-24.6272 outfile=csviscube_cmd1.fits'+\
+                        ' logfile="csviscube_cmd1.log" chatter=1'
+
+        # Check if execution of wrong command fails
+        self.test_assert(self._execute('command_that_does_not_exist') != 0,
+             'Self test of test script')
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Return
+        return
+
+    # Test csviscube from Python
+    def _test_python(self):
+        """
+        Test csviscube from Python.
+        """
+        # Set-up csviscube
+        viscube = cscripts.csviscube()
+        viscube['tmin']    =  0.0
+        viscube['tmax']    =  31557600.0
+        viscube['geolon']  =  79.4041
+        viscube['geolat']  = -24.6272
+        viscube['outfile'] = 'csviscube_py1.fits'
+        viscube['logfile'] = 'csviscube_py1.log'
+        viscube['chatter'] = 2
+
+        # Run script
+        viscube.logFileOpen()   # Make sure we get a log file
+        viscube.run()
+
+        # Return
+        return
@@ -41,6 +41,7 @@ import test_csspec
 import test_cstsdist
 import test_cstsmapmerge
 import test_cstsmapsplit
+import test_csviscube
 import test_csworkflow
 import test_csiactcopy
 import test_csiactdata
@@ -124,6 +125,7 @@ def test(installed=False, debug=False):
              test_cstsdist.Test(),
              test_cstsmapmerge.Test(),
              test_cstsmapsplit.Test(),
+             test_csviscube.Test(),
              test_csworkflow.Test()]
     # Append tests for Python 2.6+ (the IACT cscripts depend on the json