diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index eb0e6b8..fb22b0e 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -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 \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index fe35728..152bb3d 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -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
diff --git a/cscripts/csviscube.par b/cscripts/csviscube.par
new file mode 100644
index 0000000..ce09b9a
--- /dev/null
+++ b/cscripts/csviscube.par
@@ -0,0 +1,40 @@
+#############################################################################
+# #
+# 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"
diff --git a/cscripts/csviscube.py b/cscripts/csviscube.py
new file mode 100755
index 0000000..c7b3ad3
--- /dev/null
+++ b/cscripts/csviscube.py
@@ -0,0 +1,289 @@
+#! /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 .
+#
+# ==========================================================================
+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()
diff --git a/test/test_csviscube.py b/test/test_csviscube.py
new file mode 100755
index 0000000..6c6bc87
--- /dev/null
+++ b/test/test_csviscube.py
@@ -0,0 +1,104 @@
+#! /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 .
+#
+# ==========================================================================
+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
diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py
index 82efb66..e457a87 100755
--- a/test/test_python_cscripts.py
+++ b/test/test_python_cscripts.py
@@ -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