Merge branch 'devel' into 1044-implement-onoff-fitting

# Please enter a commit message to explain why this merge is necessary, # especially if it merges an updated upstream into a topic branch. # # Lines starting with '#' will be ignored, and an empty message aborts # the commit.

Merge branch 'devel' into 1044-implement-onoff-fitting
# Please enter a commit message to explain why this merge is necessary, # especially if it merges an updated upstream into a topic branch. # # Lines starting with '#' will be ignored, and an empty message aborts # the commit.
Pierrick Martin
2 parents 0e5bb9e3 85d55d08
Showing 13 changed files with 334 additions and 98 deletions
AUTHORS
ChangeLog
NEWS
README
configure.ac
dev/release_gammalib.sh
doc/source/user_manual/modules/model.rst
include/GModelSpectralGauss.hpp
inst/com/test/test_COM.cpp
inst/cta/test/test_CTA.cpp
src/fits/GFits.cpp
src/fits/GFitsHDU.cpp
src/sky/GWcs.cpp
@@ -6,3 +6,4 @@ Pierrick Martin (pierrick.martin@irap.omp.eu)
 Michael Mayer (michael.mayer@desy.de)
 Ellis Owen (Ellis.Owen@mpi-hd.mpg.de)
 Anneli Schulz (anneli.schulz@desy.de)
+
-2014-01-17  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-01-18  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
  
         * Version 0.8.0 released
         ========================
@@ -23,6 +23,8 @@
         Put likelihood computation in GObservation class
         Add broken power law spectral model
         Add Gaussian spectral model
+        Update parameter files if necessary
+        Implement exact sky map pixel solid angle computation
         Rename the following classes:
         - GWcs => GSkyProjection
         - GWcslib => GWcs
-New Features and Important Changes in GammaLib 00-08-00
+vNew Features and Important Changes in GammaLib 00-08-00
  
 Author(s) : Juergen Knoedlseder
  
-17 January 2014
+18 January 2014
  
  
 1. Introduction
@@ -375,6 +375,11 @@ installed. This allows providing absolute paths for unit testing, a need
 for the special valgrind rule in case that Sonar is used for quality
 checking.
  
+The configuration script now updates any parameter files in the users
+pfiles directory that is older than the system parameter files. This
+corrects for problems that arose when parameter files in the pfiles
+directory did not include all parameters necessary for the latest version.
+
  
 4. Application module
 ---------------------
@@ -560,7 +565,12 @@ The GSkyPixel class now supports 1D and 2D pixelisations, allowing
 a more generic usage of the class. Integer and double precision
 values can be assigned from and to GSkyPixel for 1D pixelisations.
  
-The GWcsAIT class has been added for Aitoff-Hammer projections.
+The GWcsAIT class has been added for Aitoff projections.
+
+The solid angle of sky map pixels are now computed using the spherical
+excess. This provides accurate solid angles even in case that the
+pixels are highly distorted, for example at the poles of the Aitoff
+projections.
  
  
 13. Support module
 GammaLib information
 ====================
-Version: 00-08-00 (January 31st, 2013)
+Version: 00-08-00 (18 January 2014)
 Author:  Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
  
  
@@ -92,7 +92,7 @@ The file INSTALL details more about using configure. Also try
      % ./configure --help.
  
 The `make check' command will run an extensive unit test to verify that
-GammaLib was correctly built. Make sure that all tests were successful. 
+GammaLib was correctly built.  Make sure that all tests were successful. 
  
  
 Macintosh OS X Installation
@@ -145,14 +145,14 @@ will generate a universal build which contains 32-bit and 64-bit code.
  
 BSD Installation
 ================
-GammaLib has been tested on FreeBSD successfully. Follow the Linux
+GammaLib has been tested on FreeBSD successfully.  Follow the Linux
 installation instructions above.
  
  
 Solaris Installation
 ====================
 GammaLib compiles on Solaris and openSolaris using the GNU compiler
-collection. Follow the Linux installation instructions above.
+collection.  Follow the Linux installation instructions above.
 GammaLib also compiles using the Sun Studio compiler, but due to
 some unexpected handling of global variables, it is not fully
 working (see "Known problems" below).
@@ -175,11 +175,25 @@ extensive unit test by typing:
 If everything works successfully you should see
  
     ===================
-    All 18 tests passed
+    All 20 tests passed
     ===================
  
-at the end of the test. If no Python support was compiled in, the
-number of tests performed will be reduced by one.
+or
+
+============================================================================
+Testsuite summary for gammalib 0.8.0
+============================================================================
+# TOTAL: 20
+# PASS:  20
+# SKIP:  0
+# XFAIL: 0
+# FAIL:  0
+# XPASS: 0
+# ERROR: 0
+============================================================================
+
+at the end of the test (depending on your automake version).  If no Python
+support was compiled in, the number of tests performed will be reduced by one.
  
  
 Setting up your environment
@@ -227,43 +241,50 @@ Documentation
 =============
 The doc directory (usually at /usr/local/gamma/share/doc/gammalib)
 contains the most recent set of updated documentation for this release.
-A detailed documentation can be created using Doxygen by typing:
+A detailed documentation can be created by typing:
+
+    % make doc
  
-    % make doxygen
+before installing the library. Two types of documentation exist: 
+* code documentation
+* user documentation
  
-before installing the library. Doxygen can be obtained from
-http://www.stack.nl/~dimitri/doxygen/
+Code documentation is created using Doxygen.  You need Doxygen on your
+system to generate the code documentation.  This includes man pages.
+Doxygen can be obtained from http://www.stack.nl/~dimitri/doxygen/.
  
-This is a development release and the documentation is largely, but 
-not entirely up to date. We are working on it. Please be patient or
-volunteer to help.
+User documentation is created using Sphinx.  You need Sphinx on your
+system to generate the user documentation.  Sphinx can be obtained
+from http://sphinx-doc.org/install.html.
  
  
 Bug reports
 ===========
 To report or search for bugs, please use the GammaLib Bug Tracker at
-https://sourceforge.net/tracker/?group_id=282143&atid=1196795.
+https://cta-redmine.irap.omp.eu/projects/gammalib.  Before using the
+tracker, please read 
+https://cta-redmine.irap.omp.eu/projects/gammalib/wiki/Submission_guidelines
  
  
 Known problems
 ==============
 Python support
   GammaLib comes with Python wrappers so that all classes can be
-  directly used from Python. To compile-in Python support, GammaLib
+  directly used from Python.  To compile-in Python support, GammaLib
   needs the Python.h header file, which on many distributions is not
   installed by default.
   To make Python.h available, install the Python developer package
-  in your distribution. Otherwise you will not be able to use GammaLib
+  in your distribution.  Otherwise you will not be able to use GammaLib
   from Python.
  
 Readline support
   Many distributions do not have the readline header files and symbolic
-  link to the shared library set up by default. The same is true for the
-  ncurses library that is needed by readline. To properly compile-in
+  link to the shared library set up by default.  The same is true for the
+  ncurses library that is needed by readline.  To properly compile-in
   readline support in GammaLib, make sure that the redline header files
-  exist and that the symbolic links are set. In many distributions,
+  exist and that the symbolic links are set.  In many distributions,
   the shared redline and ncurses are located in /lib, while the symbolic
-  links are in /usr/lib. Make sure that you have symbolic links with
+  links are in /usr/lib.  Make sure that you have symbolic links with
   names
  
      libreadline.so
@@ -279,23 +300,23 @@ Readline support
 Solaris
   Although GammaLib builds on Solaris using the Sun compiler, there are
   problems with global symbols in the shared library that prevent the
-  model registry to work correctly. Furthermore, GammaLib is not able
+  model registry to work correctly.  Furthermore, GammaLib is not able
   to catch its own exceptions, which prevents the FITS interface to work
   correctly.
   GammaLib has however been built and tested successfully using the GNU
   compiler, and this is the only build method that is currently supported.
   Problems have also been encountered when compiling cfitsio versions
-  more recent than 3.250. The problems have been reported to the cfitsio
-  developer team, and are likely to be solved in the future. For the
+  more recent than 3.250.  The problems have been reported to the cfitsio
+  developer team, and are likely to be solved in the future.  For the
   time being, it is recommended to use cfitsio version 3.250 on Solaris.
  
 OpenSolaris
   On OpenSolaris, the same problems concerning the SunStudio compiler
   occur as for Solaris, and also here, the GNU compiler is the recommended
-  tool to build GammaLib. Also here, cfitsio version 3.250 is the
+  tool to build GammaLib.  Also here, cfitsio version 3.250 is the
   recommended library as more recent version feature relocation
   problems. GammaLib has been tested using gcc 4.3.2 on 
-  OpenSolaris 2009.06. Make sure to create the symbolic links
+  OpenSolaris 2009.06.  Make sure to create the symbolic links
  
      ln -s /usr/bin/gcc4.3.2 /usr/bin/gcc
      ln -s /usr/bin/g++4.3.2 /usr/bin/g++
@@ -303,12 +324,12 @@ OpenSolaris
   which are not there by default to avoid excess warnings during
   compilation. 
   For OpenSolaris 2009.06, no readline package is available, hence 
-  readline support is not readily available. Readline support can however
+  readline support is not readily available.  Readline support can however
   be enabled by directly installing the GNU ncurses and readline packages
-  from source (http://ftp.gnu.org/gnu). This has been done successfully
-  on OpenSolaris 2009.06 for ncurses-5.9 and readline-6.2. Make sure
+  from source (http://ftp.gnu.org/gnu).  This has been done successfully
+  on OpenSolaris 2009.06 for ncurses-5.9 and readline-6.2.  Make sure
   to specify the `--with-shared' option when configuring ncurses, so
-  that the shared library will be built and installed. Otherwise,
+  that the shared library will be built and installed.  Otherwise,
   relocation errors may occur during compilation of GammaLib.
  
  
@@ -2,7 +2,7 @@
 # Configuration for GammaLib - a versatile toolbox for high-level analysis  #
 #                              of astronomical gamma-ray data               #
 # ------------------------------------------------------------------------- #
-# Copyright (C) 2006-2013 Juergen Knoedlseder                               #
+# Copyright (C) 2006-2014 Juergen Knoedlseder                               #
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  This program is free software: you can redistribute it and/or modify     #
@@ -65,20 +65,6 @@ fi
 echo "Create LaTeX documentation"
 cd $sourcedir/doc
 #
-cd um/latex
-latex gammalib_um.tex
-dvips gammalib_um.dvi -o gammalib_um.ps
-ps2pdf gammalib_um.ps
-cd ../..
-mv um/latex/gammalib_um.pdf .
-#
-cd dev/coding
-latex gammalib_coding.tex
-dvips gammalib_coding.dvi -o gammalib_coding.ps
-ps2pdf gammalib_coding.ps
-cd ../..
-mv dev/coding/gammalib_coding.pdf .
-#
 cd dev/inst
 latex gammalib_inst.tex
 dvips gammalib_inst.dvi -o gammalib_inst.ps
@@ -89,11 +89,11 @@ Power law
 The ``GModelSpectralPlaw`` class implements the power law function
  
 .. math::
-    \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma}
+    \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
  
 where the parameters in the XML definition have the following mappings:
  
-* :math:`N_0` = ``Prefactor``
+* :math:`k_0` = ``Prefactor``
 * :math:`\gamma` = ``Index``
 * :math:`E_0` = ``Scale``
  
@@ -145,12 +145,12 @@ Exponentially cut-off power law
 The ``GModelSpectralExpPlaw`` class implements the exponentially cut-off power law function
  
 .. math::
-    \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma}
+    \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
                     \exp \left( \frac{-E}{\tt E_{\rm cut}} \right)
  
 where the parameters in the XML definition have the following mappings:
  
-* :math:`N_0` = ``Prefactor``
+* :math:`k_0` = ``Prefactor``
 * :math:`\gamma` = ``Index``
 * :math:`E_0` = ``Scale``
 * :math:`E_{\rm cut}` = ``Cutoff``
@@ -174,7 +174,7 @@ The ``GModelSpectralBrokenPlaw`` class implements the broken power law function
  
 .. math::
  
-    \frac{dN}{dE} = N_0 \times \left \{
+    \frac{dN}{dE} = k_0 \times \left \{
     \begin{eqnarray}
       \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\
       \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise}
@@ -183,7 +183,7 @@ The ``GModelSpectralBrokenPlaw`` class implements the broken power law function
  
 where the parameters in the XML definition have the following mappings:
  
-* :math:`N_0` = ``Prefactor``
+* :math:`k_0` = ``Prefactor``
 * :math:`\gamma_1` = ``Index1``
 * :math:`\gamma_2` = ``Index2``
 * :math:`E_b` = ``BreakValue``
@@ -199,6 +199,33 @@ The XML format for specifying a broken power law is:
     <parameter name="Index2"     scale="-1"    value="2.70" min="0.01"  max="1000.0" free="1"/>
    </spectrum>
  
+
+Gaussian
+^^^^^^^^
+
+The ``GModelSpectralGauss`` class implements the gaussian function
+
+.. math::
+    \frac{dN}{dE} = \frac{N_0}{\sqrt{2\pi}\sigma}
+                    \exp \left( \frac{-(E-\bar{E})^2}{2 \sigma^2} \right)
+
+where the parameters in the XML definition have the following mappings:
+
+* :math:`N_0` = ``Normalization``
+* :math:`\bar{E}` = ``Mean``
+* :math:`\sigma` = ``Sigma``
+
+The XML format for specifying a Gaussian is:
+
+.. code-block:: xml
+
+   <spectrum type="Gaussian">
+    <parameter name="Normalization" scale="1e-10" value="1.0"  min="1e-07" max="1000.0" free="1"/>
+    <parameter name="Mean"          scale="1e6"   value="5.0"  min="0.01"  max="100.0"  free="1"/>
+    <parameter name="Sigma"         scale="1e6"   value="1.0"  min="0.01"  max="100.0"  free="1"/>
+   </spectrum>
+
+
 Node function
 ^^^^^^^^^^^^^
  
@@ -231,11 +258,11 @@ Log parabola
 The ``GModelSpectralLogParabola`` class implements the log parabola function
  
 .. math::
-    \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma+\eta \ln(E/E_0)}
+    \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma+\eta \ln(E/E_0)}
  
 where the parameters in the XML definition have the following mappings:
  
-* :math:`N_0` = ``Prefactor``
+* :math:`k_0` = ``Prefactor``
 * :math:`\gamma` = ``Index``
 * :math:`\eta` = ``Curvature``
 * :math:`E_0` = ``Scale``
@@ -182,7 +182,7 @@ GEnergy GModelSpectralGauss::mean(void) const
 /***********************************************************************//**
  * @brief Set mean energy
  *
- * @param[in] Mean energy.
+ * @param[in] mean Mean energy.
  ***************************************************************************/
 inline
 void GModelSpectralGauss::mean(const GEnergy& mean)
@@ -208,7 +208,7 @@ GEnergy GModelSpectralGauss::sigma(void) const
 /***********************************************************************//**
  * @brief Set energy width
  *
- * @param[in] Energy width
+ * @param[in] sigma Energy width.
  ***************************************************************************/
 inline
 void GModelSpectralGauss::sigma(const GEnergy& sigma)
 /***************************************************************************
  *                       test_COM.cpp - test COM classes                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -580,7 +580,7 @@ void TestGCOMOptimize::test_binned_optimizer(void)
                             47, 0,
                             29.9617, 0.377506,
                             49, 0,
-                            27.9993, 0.372603};
+                            27.9992, 0.372603};
     test_try("Perform LM optimization");
     try {
         GOptimizerLM opt;
 /***************************************************************************
  *                       test_CTA.cpp - Test CTA classes                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -335,7 +335,8 @@ void TestGCTAResponse::test_response_irf_diffuse(void)
 {
     // Set reference value
     //double ref = 13803.800313356;
-    const double ref = 13803.5186374;
+    //const double ref = 13803.5186374;
+    const double ref = 13803.6932774; // After GWcs::solidangle improvement
  
     // Set parameters
     double src_ra  = 201.3651;
@@ -412,7 +413,8 @@ void TestGCTAResponse::test_response_npred_diffuse(void)
 {
     // Set reference value
     //double ref = 11212.26274; // npred_spec precision of 1e-6
-    const double ref = 11212.437464; // npred_spec precision of 1e-5
+    //const double ref = 11212.437464; // npred_spec precision of 1e-5
+    const double ref = 11212.4370702; // After GWcs::solidangle improvement
  
     // Set parameters
     double src_ra  = 201.3651;
@@ -766,10 +768,10 @@ void TestGCTAOptimize::test_binned_optimizer(void)
     double fit_results[] = {83.6331, 0,
                             22.0145, 0,
                             5.616410411e-16, 1.904730785e-17,
-                            -2.481781246, -0.02580905077,
+                            -2.48178, -0.02580905077,
                             300000, 0,
                             1, 0,
-                            2.933677595, 0.06639644824,
+                            2.93368, 0.06639644824,
                             6.550723074e-05, 1.945714239e-06,
                             -1.833781187, -0.0160819,
                             1000000, 0,
@@ -1213,19 +1213,32 @@ void GFits::saveto(const std::string&amp; filename, const bool&amp; clobber)
     // Expand environment variables
     std::string fname = gammalib::expand_env(filename);
  
+    // Create gzipped file name version
+    std::string gzfname = fname + ".gz";
+
     // Debug header
     #if defined(G_DEBUG)
     std::cout << "GFits::saveto(\"" << fname << "\", " << clobber << ")"
               << " (size=" << size() << ") -->" << std::endl;
     #endif
  
-    // If overwriting has been specified then remove any existing file ...
+    // If overwriting has been specified then remove any existing file. We
+    // consider here also the possibility that the file is gzipped, but
+    // the file name does not contain the .gz extension.
     if (clobber) {
-        std::remove(fname.c_str());
+        if (gammalib::file_exists(fname)) {
+            std::remove(fname.c_str());
+        }
+        else {
+            if (gammalib::file_exists(gzfname)) {
+                std::remove(gzfname.c_str());
+            }
+        }
     }
  
-    // ... otherwise, if file exists then throw an exception
-    else if (gammalib::file_exists(fname)) {
+    // Otherwise, if file exists then throw an exception
+    else if (gammalib::file_exists(fname) ||
+             gammalib::file_exists(gzfname)) {
         throw GException::fits_file_exist(G_SAVETO, fname);
     }
  
@@ -1393,7 +1406,7 @@ void GFits::free_members(void)
     // If FITS file has been opened then close it now
     if (m_fitsfile != NULL) {
  
-        // If file has been created but not yet save then delete the file
+        // If file has been created but not yet saved then delete the file
         // now. We do not worry about the status in this case.
         if (m_created) {
             int status = 0;
 /***************************************************************************
  *                   GFitsHDU.cpp  - FITS HDU handling class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2008-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2008-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -29,6 +29,7 @@
 #include <config.h>
 #endif
 #include <iostream>
+#include <cstdlib>
 #include "GException.hpp"
 #include "GTools.hpp"
 #include "GFitsCfitsio.hpp"
@@ -511,10 +512,3 @@ void GFitsHDU::free_members(void)
     // Return
     return;
 }
-
-
-/*==========================================================================
- =                                                                         =
- =                                 Friends                                 =
- =                                                                         =
- ==========================================================================*/
@@ -56,6 +56,7 @@
  
 /* __ Coding definitions _________________________________________________ */
 //#define G_LIN_MATINV_FORCE_PC                             // Force PC usage
+#define G_SOLIDANGLE_OPTION 2                // 1=sin/cos, 2=tan, 3=vectors
  
 /* __ Debug definitions __________________________________________________ */
 //#define G_DIR2XY_DEBUG                                      // Debug dir2xy
@@ -425,34 +426,213 @@ void GWcs::write(GFitsHDU&amp; hdu) const
  *
  * @param[in] pixel Pixel index (x,y)
  *
- * Estimate solid angles of pixel by computing the coordinates in the 4 pixel
- * corners. The surface is computed using a cartesian approximation:
- *           a
- *     1-----5-----2                 a+b
- *    /      | h    \    where A = h ---
- *   4-------6-------3                2
- *           b
- * This is a brute force technique that works sufficiently well for non-
- * rotated sky maps. Something more intelligent should be implemented in
- * the future.
+ * Estimates solid angles of pixels using the Girard equation for excess
+ * area - see: http://mathworld.wolfram.com/SphericalPolygon.html
+ *
+ *
+ * Below, the definiton of the pixel cornes and sides are shown as used
+ * within the code.
+ *
+ *             a12
+ *         1---------2
+ *         |\       /|
+ *         | \a13  / |
+ *         |  \   /  |
+ *         |   \ /   |
+ *      a14|    X    |a23
+ *         |   / \   |
+ *         |  /   \  |
+ *         | /a24  \ |
+ *         |/       \|
+ *         4---------3
+ *             a34
+ *
  ***************************************************************************/
 double GWcs::solidangle(const GSkyPixel& pixel) const
 {
-    // Get the sky directions of the 6 points
+    // Get the sky directions of the 4 points
     GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
     GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
     GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
     GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
-    GSkyDir dir5 = pix2dir(GSkyPixel(pixel.x(),     pixel.y()-0.5));
-    GSkyDir dir6 = pix2dir(GSkyPixel(pixel.x(),     pixel.y()+0.5));
  
-    // Compute distances between sky directions
-    double a = dir1.dist(dir2);
-    double b = dir3.dist(dir4);
-    double h = dir5.dist(dir6);
+    // Compile option 1: Use triginometric function
+    #if G_SOLIDANGLE_OPTION == 1
+
+    // Initialise solid angle
+    double solidangle = 0.0;
+
+    // Compute angular distances between pixel corners
+    double a12 = dir1.dist(dir2);
+    double a14 = dir1.dist(dir4);
+    double a23 = dir2.dist(dir3);
+    double a34 = dir3.dist(dir4);
+
+    // Special case: a12 or a14 is zero, then pixel is a triangle composed
+    // of [2,3,4]
+    if (a12 <= 0.0 || a14 <= 0.0) {
+
+        // Compute diagonals
+        double a24 = dir2.dist(dir4);
+
+        // Compute sines
+        double sin_a23 = std::sin(a23);
+        double sin_a24 = std::sin(a24);
+        double sin_a34 = std::sin(a34);
+
+        // Compute cosines
+        double cos_a23 = std::cos(a23);
+        double cos_a24 = std::cos(a24);
+        double cos_a34 = std::cos(a34);
+
+        // Compute angles
+        double angle2 = gammalib::acos((cos_a34-cos_a24*cos_a23)/(sin_a24*sin_a23));
+        double angle3 = gammalib::acos((cos_a24-cos_a34*cos_a23)/(sin_a34*sin_a23));
+        double angle4 = gammalib::acos((cos_a23-cos_a24*cos_a34)/(sin_a24*sin_a34));
+
+        // Compute excess area to determine solid angle
+        solidangle = (angle2 + angle3 + angle4) - gammalib::pi;
+
+    }
+
+    // Special case: a23 or a 34 is zero, then pixel is a triangle composed
+    // of [1,2,4]
+    else if (a23 <= 0.0 || a34 <= 0.0) {
+
+        // Compute diagonals
+        double a24 = dir2.dist(dir4);
  
-    // Compute solid angle
-    double solidangle = 0.5*(h*(a+b));
+        // Compute sines
+        double sin_a12 = std::sin(a12);
+        double sin_a14 = std::sin(a14);
+        double sin_a24 = std::sin(a24);
+
+        // Compute cosines
+        double cos_a12 = std::cos(a12);
+        double cos_a14 = std::cos(a14);
+        double cos_a24 = std::cos(a24);
+
+        // Compute angles
+        double angle1 = gammalib::acos((cos_a24-cos_a12*cos_a14)/(sin_a12*sin_a14));
+        double angle2 = gammalib::acos((cos_a14-cos_a12*cos_a24)/(sin_a12*sin_a24));
+        double angle4 = gammalib::acos((cos_a12-cos_a14*cos_a24)/(sin_a14*sin_a24));
+
+        // Compute excess area to determine solid angle
+        solidangle = (angle1 + angle2 + angle4) - gammalib::pi;
+
+    }
+
+    // Otherwise we have a polygon
+    else {
+
+        // Compute diagonals
+        double a13 = dir1.dist(dir3);
+        double a24 = dir2.dist(dir4);
+
+        // Compute sines
+        double sin_a12 = std::sin(a12);
+        double sin_a14 = std::sin(a14);
+        double sin_a23 = std::sin(a23);
+        double sin_a34 = std::sin(a34);
+
+        // Compute cosines
+        double cos_a12 = std::cos(a12);
+        double cos_a13 = std::cos(a13);
+        double cos_a14 = std::cos(a14);
+        double cos_a23 = std::cos(a23);
+        double cos_a24 = std::cos(a24);
+        double cos_a34 = std::cos(a34);
+
+        // Compute angles
+        double angle1 = std::acos((cos_a13-cos_a34*cos_a14)/(sin_a34*sin_a14));
+        double angle2 = std::acos((cos_a24-cos_a23*cos_a34)/(sin_a23*sin_a34));
+        double angle3 = std::acos((cos_a13-cos_a12*cos_a23)/(sin_a12*sin_a23));
+        double angle4 = std::acos((cos_a24-cos_a14*cos_a12)/(sin_a14*sin_a12));
+
+        // Use Girard equation for excess area to determine solid angle
+        solidangle = (angle1 + angle2 + angle3 + angle4) - gammalib::twopi;
+
+    } // endif: we had a polynom
+
+    // Compile option 2: Use Huilier's theorem
+    #elif G_SOLIDANGLE_OPTION == 2
+
+    // Initialise solid angle
+    double solidangle = 0.0;
+
+    // Compute angular distances between pixel corners
+    double a12 = dir1.dist(dir2);
+    double a14 = dir1.dist(dir4);
+    double a23 = dir2.dist(dir3);
+    double a24 = dir2.dist(dir4);
+    double a34 = dir3.dist(dir4);
+
+    // Special case: a12 or a14 is zero, then pixel is a triangle composed
+    // of [2,3,4]
+    if (a12 <= 0.0 || a14 <= 0.0) {
+        double s   = 0.5 * (a23 + a34 + a24);
+        solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
+                                               std::tan(0.5*(s-a23)) *
+                                               std::tan(0.5*(s-a34)) *
+                                               std::tan(0.5*(s-a24))));
+    }
+
+    // Special case: a23 or a 34 is zero, then pixel is a triangle composed
+    // of [1,2,4]
+    else if (a23 <= 0.0 || a34 <= 0.0) {
+        double s   = 0.5 * (a12 + a24 + a14);
+        solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
+                                               std::tan(0.5*(s-a12)) *
+                                               std::tan(0.5*(s-a24)) *
+                                               std::tan(0.5*(s-a14))));
+    }
+
+    // Otherwise we have a polygon
+    else {
+
+        // Triangle 1 [1,2,4]
+        double s1      = 0.5 * (a12 + a24 + a14);
+        double excess1 = std::atan(std::sqrt(std::tan(0.5*s1) *
+                                             std::tan(0.5*(s1-a12)) *
+                                             std::tan(0.5*(s1-a24)) *
+                                             std::tan(0.5*(s1-a14))));
+
+        // Triangle 2 [2,3,4]
+        double s2      = 0.5 * (a23 + a34 + a24);
+        double excess2 = std::atan(std::sqrt(std::tan(0.5*s2) *
+                                             std::tan(0.5*(s2-a23)) *
+                                             std::tan(0.5*(s2-a34)) *
+                                             std::tan(0.5*(s2-a24))));
+
+        // Determine solid angle
+        solidangle = 4.0 * (excess1 + excess2);
+
+    } // endif: we had a polynom
+
+    // Compile option 3: Use vectors
+    #else
+
+    // Get vectors to pixel corners
+    GVector vec1 = dir1.celvector();
+    GVector vec2 = dir2.celvector();
+    GVector vec3 = dir3.celvector();
+    GVector vec4 = dir4.celvector();
+
+    // Compute inner angles of pixel corners
+    double angle1 = gammalib::acos(cross(vec2, (cross(vec1, vec2))) * 
+                                   cross(vec2, (cross(vec3, vec2))));
+    double angle2 = gammalib::acos(cross(vec3, (cross(vec2, vec3))) *
+                                   cross(vec3, (cross(vec4, vec3))));
+    double angle3 = gammalib::acos(cross(vec4, (cross(vec3, vec4))) *
+                                   cross(vec4, (cross(vec1, vec4))));
+    double angle4 = gammalib::acos(cross(vec1, (cross(vec4, vec1))) *
+                                   cross(vec1, (cross(vec2, vec1))));
+
+    // Use Girard equation for excess area to determine solid angle
+    double solidangle = gammalib::deg2rad * gammalib::deg2rad +
+                        (angle1 + angle2 + angle3 + angle4) -
+                        gammalib::twopi;
+    #endif
  
     // Return solid angle
     return solidangle;