GTools.cpp

Go to the documentation of this file.

/***************************************************************************
 *                       GTools.cpp - GammaLib tools                       *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2008-2026 by 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/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GTools.cpp
 * @brief Gammalib tools implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <unistd.h>        // geteuid(), access() functions
#include <fcntl.h>         // fcntl() functions
#include <sys/stat.h>
#include <sys/time.h>      // timeval
#include <sys/select.h>    // select() function
#include <sys/types.h>
#include <sys/socket.h>    // recv() function
#include <netinet/in.h>    // struct sockaddr_in, struct sockaddr
#include <netdb.h>         // struct hostent, gethostbyname
#include <pwd.h>           // getpwuid()
#include <cmath>
#include <cfloat>
#include <cctype>
#include <ctime>
#include <cstdlib>         // std::getenv() function
#include <cstring>         // std::strlen() and std::memcpy functions
#include <iostream>
#include <sstream>
#include <algorithm>
#include <cstdio>          // sprintf
#include "GTools.hpp"
#include "GException.hpp"
#include "GFits.hpp"
#include "GEnergy.hpp"
#include "GFilename.hpp"
#include "GMath.hpp"
#include "GXmlElement.hpp"
#include "GXmlNode.hpp"
#include "GXmlText.hpp"
 
/* __ OpenMP section _____________________________________________________ */
#ifdef _OPENMP
#include <omp.h>
#endif
 
/* __ Includes for memory usage determination ____________________________ */
#if defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
#include <unistd.h>
#include <sys/resource.h>
#if defined(__APPLE__) && defined(__MACH__)
#include <mach/mach.h>
#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
#include <fcntl.h>
#include <procfs.h>
#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
#include <stdio.h>
#endif
#endif
 
/* __ Compile options ____________________________________________________ */
//#define G_CHECK_FOR_NAN
//#define G_USE_CURL      //!< Use curl in host_country()
 
/* __ Function name definitions __________________________________________ */
#define G_XML2STR                           "gammalib::xml2str(std::string&)"
#define G_HTTP_QUERY       "gammalib::http_query(str::string&, std::string&)"
 
/* __ Coding definitions _________________________________________________ */
#define G_PARFORMAT_LENGTH 29
#define G_CHAR_BUFFER      256
 
 
/***********************************************************************//**
 * @brief Strip leading and trailing whitespace from string
 *
 * @param[in] arg String from which whitespace should be stripped.
 * @return String with stripped whitespace.
 ***************************************************************************/
std::string gammalib::strip_whitespace(const std::string& arg)
{
    // Return result
    return (strip_chars(arg, " "));
}
 
 
/***********************************************************************//**
 * @brief Strip leading and trailing character from string
 *
 * @param[in] arg String from which character should be stripped.
 * @param[in] chars Character(s) to be stripped.
 * @return String with stripped characters.
 ***************************************************************************/
std::string gammalib::strip_chars(const std::string& arg,
                                  const std::string& chars)
{
    // Initialise empty result string
    std::string result;
 
    // Continue only if argument is not empty
    if (!arg.empty()) {
 
        // Get start and stop
        std::string::size_type start = arg.find_first_not_of(chars);
        std::string::size_type stop  = arg.find_last_not_of(chars);
 
        // Continue only if start and stop are valid
        if (start != std::string::npos && stop != std::string::npos) {
 
            // Continue only if stop is larger then or equal to  start
            if (start <= stop) {
                result = arg.substr(start, stop-start+1);
            }
 
        } // endif: start and stop were valid
 
    } // endif: argument was not empty
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Strip trailing character from string
 *
 * @param[in] arg String from which character should be stripped.
 * @param[in] chars Character(s) to be stripped.
 * @return String with stripped characters.
 ***************************************************************************/
std::string gammalib::rstrip_chars(const std::string& arg,
                                   const std::string& chars)
{
    // Initialise empty result string
    std::string result;
 
    // Continue only if argument is not empty
    if (!arg.empty()) {
 
        // Get stop
        std::string::size_type stop = arg.find_last_not_of(chars);
 
        // Continue only if stop is valid
        if (stop != std::string::npos) {
            result = arg.substr(0, stop+1);
        }
 
    } // endif: argument was not empty
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Replace string segment in string
 *
 * @param[in] arg String in which character segements are to be replaced
 * @param[in] segment String segment to be replaced.
 * @param[in] replacement Replacement string.
 * @return String with replaced segements.
 *
 * Replaces string segments by a replacement string in a given string.
 *
 * If the input string @p arg is "Wonderful", the @p segment is "onder" and
 * the @p replacement is "ish" the method will return "Wishful".
 ***************************************************************************/
std::string gammalib::replace_segment(const std::string& arg,
                                      const std::string& segment,
                                      const std::string& replacement)
{
    // Initialise result string by argument
    std::string result = arg;
 
    // Initialise character pointer
    std::string::size_type pos   = 0;
    std::string::size_type start = 0;
 
    // Loop over string
    while (start != std::string::npos) {
        start = result.find(segment, pos);
        if (start != std::string::npos) {
            result = result.replace(start, segment.length(), replacement);
            pos    = start + replacement.length();
        }
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Expand environment variables in string
 *
 * @param[in] arg String.
 * @return String with expected environment variables.
 *
 * Expands any environment variable that is found in a string. Valid
 * delimiters for environment variables are $ENV{<name>}, $ENV(<name>),
 * ${<name>}, $(<name>) and $<name> (in the last case the terminating
 * delimiter is either a / or a blank character or the end of the string).
 * Environment variables occuring within single quotes (') are ignored.
 * Environment variables that do not exist will be kept as specified.
 *
 * The method also replaces ~ or ~user by the user's home directory, ~+
 * by the value of the PWD environment variable and ~- by the value of the
 * OLDPWD variable. If the user or the PWD or OLDPWD variable are not found,
 * no replacement is done.
 *
 * This function has been inspired by the function ape_util_expand_env_var
 * from ape_util.c in the ape software developed at HEASARC.
 ***************************************************************************/
std::string gammalib::expand_env(const std::string& arg)
{
    // Set environment variable delimiters
    static const char* begin_delim[] = { "$ENV{", "$ENV(", "${", "$(", "$", "~" };
    static const char* end_delim[]   = { "}", ")", "}", ")", "/", "/" };
    static const int   num_delim     = 6;
 
    // Initialise result with argument
    std::string result = arg;
 
    // Initialise parse parameters
    size_t index    = 0;
    bool   in_quote = false;
 
    // Loop over string
    while (index < result.length()) {
 
        // If we have an escaped character then skip the current character
        // and the next one
        if (result[index] == '\\') {
            index += 2;
            continue;
        }
 
        // If we have a single quote then toggle the quote state
        if (result[index] == '\'') {
            in_quote = !in_quote;
            index++;
            continue;
        }
 
        // Don't expand environment variables inside single quotes. Note
        // that double quotes are ok.
        if (in_quote) {
            index++;
            continue;
        }
 
        // Find delimiter which indicates beginning of an environment variable
        size_t begin_length = 0;
        int    delim_idx    = 0;
        for (; delim_idx < num_delim; ++delim_idx) {
            size_t len = std::strlen(begin_delim[delim_idx]);
            if (result.compare(index, len, begin_delim[delim_idx]) == 0) {
                begin_length = len;
                break;
            }
        }
 
        // If we found a delimiter then process the environment variable
        if (begin_length > 0) {
 
            // Search for the termination delimiter of the environment
            // variable. There is a special case for delimiters 4 and 5:
            // It has always an end_length of zero as the / is not a real
            // delimiter, but just an indicator that the environment variable
            // ends. Another indicator is a blank. If the end of the
            // string has been reached this is also acceptable.
            size_t i_start    = index + begin_length;
            size_t i_end      = i_start;
            size_t end_length = 0;
            if (delim_idx == 4 || delim_idx == 5) {
                while (i_end < result.length() &&
                       result.compare(i_end, 1, "/") != 0 &&
                       result.compare(i_end, 1, " ") != 0) {
                    i_end++;
                }
            }
            else {
                end_length = std::strlen(end_delim[delim_idx]);
                while (i_end < result.length() &&
                       result.compare(i_end, end_length, end_delim[delim_idx]) != 0) {
                    i_end++;
                }
            }
 
            // If termination delimiter has been found then expand the
            // environment variable
            if (i_end < result.length() || delim_idx == 4) {
 
                // Extract environment variable name
                std::string name        = result.substr(i_start, i_end-i_start);
                size_t      name_length = name.length();
 
                // Initialise pointer on environment variable
                const char* env = NULL;
 
                // Handle ~
                if (delim_idx == 5) {
                    if (name_length == 0) {
                        if ((env = std::getenv("HOME")) == NULL) {
                            struct passwd *pw = getpwuid(getuid());
                            if (pw != NULL) {
                                env = pw->pw_dir;
                            }
                        }
                    }
                    else {
                        if (name == "+") {
                            env = std::getenv("PWD");
                        }
                        else if (name == "-") {
                            env = std::getenv("OLDPWD");
                        }
                        else {
                            struct passwd *pw = getpwnam(name.c_str());
                            if (pw != NULL) {
                                env = pw->pw_dir;
                            }
                        }
                    }
                }
 
                // ... otherwise get the environment variable
                else {
                    env = std::getenv(name.c_str());
                }
 
                // If the environment variable has been found then replace
                // it by its value
                if (env != NULL) {
 
                    // Erase delimiters and environment variable
                    result.erase(index, begin_length+name_length+end_length);
 
                    // Set replacement string and its length
                    std::string replace(env);
                    size_t      replace_length = replace.length();
 
                    // Insert replacement string
                    result.insert(index, replace);
 
                    // Advance pointer
                    index += replace_length;
 
                } // endif: environment variable has been found
 
                // If no environment variable has been found then set
                // index=i_end+1
                else {
                    index = i_end + 1;
                }
 
            } // endif: termination delimiter found
 
            // If no environment variable has been found then set index=i_end+1
            else {
                index = i_end + 1;
            }
 
        } // endif: we found an environment variable delimiter
 
        // ... otherwise advance to next character
        else {
            index++;
        }
 
    } // endwhile: looped over string
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Build file path from path name and file name
 *
 * @param[in] pathname Path name.
 * @param[in] filename File name.
 * @return File path.
 *
 * Builds a file path by combining the @p pathname and the @p filename
 * following
 *
 *      filepath = pathname/filename
 *
 * If @p pathname is an empty string, the method simply returns the
 * @p filename.
 ***************************************************************************/
std::string gammalib::filepath(const std::string& pathname,
                               const std::string& filename)
{
    // Initialise filepath
    std::string filepath;
 
    // If path name is empty, simply return the file name
    if (pathname.empty()) {
        filepath = filename;
    }
 
    // ... otherwise combine both
    else {
        filepath = pathname + "/" + filename;
    }
 
    // Return the file path
    return filepath;
}
 
 
/***********************************************************************//**
 * @brief Return temporary file name
 *
 * @return Temporary file name.
 *
 * Returns a temporary file name.
 ***************************************************************************/
std::string gammalib::tmpnam(void)
{
    // Set default temporary directory
    char default_tmpdir[] = "/tmp";
 
    // Get pointer to name of temporary directory by searching various
    // possible environment variables
    char *tmpdir = NULL;
    if ((tmpdir = std::getenv("TEMP")) == NULL) {
        if ((tmpdir = std::getenv("TMP")) == NULL) {
            if ((tmpdir = std::getenv("TMPDIR")) == NULL) {
                tmpdir = default_tmpdir;
            }
        }
    }
 
    // Allocate empty filename
    char filename[256];
    filename[0] = '\0';
 
    // Combine temporary directory with temporary filename
    strcat(filename, tmpdir);
    strcat(filename, "/gammalibXXXXXX");
 
    // Create temporary file
    int         fd = mkstemp(filename);
    std::string tmpname(filename);
    close(fd);
    unlink(filename);
 
    // Return temporary file name
    return tmpname;
}
 
 
/***********************************************************************//**
 * @brief Return value of environment variable
 *
 * @param[in] arg Environment variable
 * @return Value of environment variable.
 *
 * Returns the value of an environment variable @p arg. If the environment
 * variable is not found an empty string will be returned.
 ***************************************************************************/
std::string gammalib::getenv(const std::string& arg)
{
    // Initialise environment variable
    std::string env;
 
    // Get pointer to environment variable string
    const char* ptr = std::getenv(arg.c_str());
 
    // If pointer is not NULL then extract value of environment variable
    if (ptr != NULL) {
        env = std::string(ptr);
    }
 
    // Return environment variable
    return env;
}
 
 
/***********************************************************************//**
 * @brief Convert unsigned short integer value into string
 *
 * @param[in] value Unsigned short integer to be converted into string.
 * @return String with unsigned short integer value.
 ***************************************************************************/
std::string gammalib::str(const unsigned short int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert unsigned integer value into string
 *
 * @param[in] value Unsigned integer to be converted into string.
 * @return String with unsigned integer value.
 ***************************************************************************/
std::string gammalib::str(const unsigned int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert unsigned long integer value into string
 *
 * @param[in] value Unsigned long integer to be converted into string.
 * @return String with unsigned long integer value.
 ***************************************************************************/
std::string gammalib::str(const unsigned long int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert unsigned long long integer value into string
 *
 * @param[in] value Unsigned long long integer to be converted into string.
 * @return String with unsigned long long integer value.
 ***************************************************************************/
std::string gammalib::str(const unsigned long long int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert short integer value into string
 *
 * @param[in] value Short integer to be converted into string.
 * @return String with short integer value.
 ***************************************************************************/
std::string gammalib::str(const short int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert integer value into string
 *
 * @param[in] value Integer to be converted into string.
 * @param[in] fmt Format string.
 * @return String with integer value.
 ***************************************************************************/
std::string gammalib::str(const int& value, const std::string& fmt)
{
    // Allocate character buffer
    char buffer[G_CHAR_BUFFER];
 
    // Put integer into buffer (not C-98 compliant)
    gammalib_snprintf(buffer, G_CHAR_BUFFER, fmt.c_str(), value);
 
    // Convert buffer into string
    std::string str_buffer(buffer);
 
    // Return string buffer
    return str_buffer;
}
 
 
/***********************************************************************//**
 * @brief Convert long integer value into string
 *
 * @param[in] value Long integer to be converted into string.
 * @return String with long integer value.
 ***************************************************************************/
std::string gammalib::str(const long int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert long long integer value into string
 *
 * @param[in] value Long long integer to be converted into string.
 * @return String with long long integer value.
 ***************************************************************************/
std::string gammalib::str(const long long int& value)
{
    std::ostringstream s_value;
    s_value << value;
    return s_value.str();
}
 
 
/***********************************************************************//**
 * @brief Convert single precision value into string
 *
 * @param[in] value Single precision value to be converted into string.
 * @param[in] precision Floating point precision.
 * @return String with single precision value.
 *
 * Converts a single precision value into a string. Any positive
 * @p precision argument specifies the exact number of digits after the
 * comma.
 ***************************************************************************/
std::string gammalib::str(const float& value, const int& precision)
{
    // Allocate output stream
    std::ostringstream s_value;
 
    // If specified then set the requested fixed point precision. Otherwise
    // use a precision that should be sufficient for floating point values.
    if (precision > 0) {
        s_value.precision(precision);
        s_value.setf(std::ios::fixed, std::ios::floatfield);
    }
    else {
        s_value.precision(7);
    }
 
    // Put floating point value in stream
    s_value << value;
 
    // Convert to a string
    std::string result = s_value.str();
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert double precision value into string
 *
 * @param[in] value Double precision value to be converted into string.
 * @param[in] precision Floating point precision.
 * @return String with double precision value.
 *
 * Converts a double precision value into a string. Any positive
 * @p precision argument specifies the exact number of digits after the
 * comma.
 ***************************************************************************/
std::string gammalib::str(const double& value, const int& precision)
{
    // Allocate output stream
    std::ostringstream s_value;
 
    // If specified then set the requested fixed point precision. Otherwise
    // use a precision that should be sufficient for floating point values.
    if (precision > 0) {
        s_value.precision(precision);
        s_value.setf(std::ios::fixed, std::ios::floatfield);
    }
    else {
        s_value.precision(15);
    }
 
    // Put double precision floating point value in stream
    s_value << value;
 
    // Convert to a string
    std::string result = s_value.str();
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert complex value into string
 *
 * @param[in] value Complex value to be converted into string.
 * @param[in] precision Floating point precision.
 * @return String with complex value.
 *
 * Converts a complex value into a string. Any positive @p precision argument
 * specifies the exact number of digits after the comma.
 ***************************************************************************/
std::string gammalib::str(const std::complex<double>& value,
                          const int&                  precision)
{
    // Allocate output stream
    std::ostringstream s_value;
 
    // If specified then set the requested fixed point precision. Otherwise
    // use a precision that should be sufficient for floating point values.
    if (precision > 0) {
        s_value.precision(precision);
        s_value.setf(std::ios::fixed, std::ios::floatfield);
    }
    else {
        s_value.precision(15);
    }
 
    // Put double precision floating point value in stream
    s_value << value.real();
    if (value.imag() < 0.0) {
        s_value << "-";
    }
    else {
        s_value << "+";
    }
    s_value << std::abs(value.imag()) << "j";
 
    // Convert to a string
    std::string result = s_value.str();
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Return current date
 *
 * Returns the current date as string in the format yyyy-mm-ddThh:mm:ss.
 ***************************************************************************/
std::string gammalib::strdate(void)
{
    // Allocate variables
    struct std::tm timeStruct;
    std::time_t    now;
    char           buffer[100];
 
    // Get time
    now = std::time(NULL);
    #ifdef HAVE_GMTIME_R
    std::gmtime_r(&now, &timeStruct);
    #else
    std::memcpy(&timeStruct, gmtime(&now), sizeof(struct tm));
    #endif
 
    // Write message type, time and task name to buffer
    std::snprintf(buffer, 100, "%04d-%02d-%02dT%02d:%02d:%02d",
                  timeStruct.tm_year + 1900,
                  timeStruct.tm_mon + 1,
                  timeStruct.tm_mday,
                  timeStruct.tm_hour,
                  timeStruct.tm_min,
                  timeStruct.tm_sec);
 
    // Build string from buffer
    std::string date = buffer;
 
    // Return date
    return date;
}
 
 
/***********************************************************************//**
 * @brief Convert string into short value
 *
 * @param[in] arg String to be converted.
 * @return Short value.
 ***************************************************************************/
short gammalib::toshort(const std::string& arg)
{
    std::istringstream iss(arg);
    short              result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into unsigned short value
 *
 * @param[in] arg String to be converted.
 * @return Unsigned short value.
 ***************************************************************************/
unsigned short gammalib::toushort(const std::string& arg)
{
    std::istringstream iss(arg);
    unsigned short     result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into integer value
 *
 * @param[in] arg String to be converted.
 * @return Integer value.
 ***************************************************************************/
int gammalib::toint(const std::string& arg)
{
    std::istringstream iss(arg);
    int                result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into unsigned integer value
 *
 * @param[in] arg String to be converted.
 * @return Unsigned integer value.
 ***************************************************************************/
unsigned int gammalib::touint(const std::string& arg)
{
    std::istringstream iss(arg);
    unsigned int       result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into long value
 *
 * @param[in] arg String to be converted.
 * @return Long value.
 ***************************************************************************/
long gammalib::tolong(const std::string& arg)
{
    std::istringstream iss(arg);
    long               result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into unsigned long value
 *
 * @param[in] arg String to be converted.
 * @return Unsigned long value.
 ***************************************************************************/
unsigned long gammalib::toulong(const std::string& arg)
{
    std::istringstream iss(arg);
    unsigned long      result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into long long value
 *
 * @param[in] arg String to be converted.
 * @return Long long value.
 ***************************************************************************/
long long gammalib::tolonglong(const std::string& arg)
{
    std::istringstream iss(arg);
    long long          result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into unsigned long long value
 *
 * @param[in] arg String to be converted.
 * @return Unsigned long long value.
 ***************************************************************************/
unsigned long long gammalib::toulonglong(const std::string& arg)
{
    std::istringstream iss(arg);
    unsigned long long result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into single precision value
 *
 * @param[in] arg String to be converted.
 * @return Single precision value.
 ***************************************************************************/
float gammalib::tofloat(const std::string& arg)
{
    std::istringstream iss(arg);
    float              result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string into double precision value
 *
 * @param[in] arg String to be converted.
 * @return Double precision value.
 ***************************************************************************/
double gammalib::todouble(const std::string& arg)
{
    std::istringstream iss(arg);
    double             result;
    iss >> std::dec >> result;
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string to upper case
 *
 * @param[in] arg String to be converted to upper case.
 * @return String converted to upper case.
 ***************************************************************************/
std::string gammalib::toupper(const std::string& arg)
{
    std::string s = arg;
    std::transform(s.begin(), s.end(), s.begin(), ::toupper);
    return s;
}
 
 
/***********************************************************************//**
 * @brief Convert string to lower case
 *
 * @param[in] arg String to be converted to upper case.
 * @return String converted to lower case.
 ***************************************************************************/
std::string gammalib::tolower(const std::string& arg)
{
    std::string s = arg;
    std::transform(s.begin(), s.end(), s.begin(), ::tolower);
    return s;
}
 
 
/***********************************************************************//**
 * @brief Split string
 *
 * @param[in] s String to be splitted.
 * @param[in] sep Separator(s).
 * @return Vector of split strings.
 *
 * Splits a string on the basis of one or multiple separator characters. The
 * separator characters are provided by the @p sep argument. Subsequent
 * separator characters that are not seperated by some other characters will
 * lead to an empty string element, except for a blank separator where
 * subsequent blanks are takens as a single separator. Below a few examples
 * that illustrate how the function will split a given string.
 *
 *     "Name;RA;DEC" => ["Name","RA","DEC"] (sep=";")
 *     "My house is    red" => ["My","house","is","red"] (sep=" ")
 *     "IRF::FRONT" => ["IRF","","FRONT"] (sep=":")
 *     "Fields;RA,DEC,Flux" => ["Fields","RA","DEC","Flux"] (sep=";,")
 *     "Last;Field;" => ["Last","Field",""] (sep=";")
 ***************************************************************************/
std::vector<std::string> gammalib::split(const std::string& s,
                                         const std::string& sep)
{
    // Allocate result string vector
    std::vector<std::string> result;
 
    // Initialise counters
    std::size_t pos = 0;
    std::size_t len = s.length();
 
    // Loop over string
    while (pos < len && pos != std::string::npos) {
 
        // Get index of first separator occurence and preset the length
        // of the substring to the end of the string
        std::size_t index = s.find_first_of(sep, pos);
        std::size_t n     = std::string::npos;
 
        // If we did not reach the end then compute now the length of the
        // substring
        if (index != std::string::npos) {
            n = index-pos;
        }
 
        // If we have no whitespace separator and the length of the
        // substring is zero then push back an empty string. If the
        // length of the substring is positive then push back the
        // substring.
        if (sep != " " && n == 0) {
            result.push_back("");
        }
        else if (n > 0) {
            result.push_back(s.substr(pos, n));
        }
 
        // Go to the string position after the last separator
        pos = (index != std::string::npos) ? index + 1 : std::string::npos;
 
        // If the position is pointing right beyong the last string
        // character we terminated with a separator, hence we need to push
        // back one more empty string before we leave
        if (sep != " " && pos == len) {
            result.push_back("");
        }
 
    } // endwhile: there were still characters in the string
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Fill string with n strings of same type
 *
 * @param[in] s String to be filled.
 * @param[in] n Number of fillings.
 * @return Filled strings.
 *
 * Replicates a given string n time.
 ***************************************************************************/
std::string gammalib::fill(const std::string& s, const int& n)
{
    // Initialise result
    std::string result = "";
 
    // Replicate string
    for (int i = 0; i < n; ++i) {
        result.append(s);
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Left justify string to achieve a length of n characters
 *
 * @param[in] s String to be left justified.
 * @param[in] n Requested total width.
 * @param[in] c Fill character.
 * @return Left justified string.
 *
 * Left justify string by adding @p c to the right to achieve a length of
 * @p n characters.
 ***************************************************************************/
std::string gammalib::left(const std::string& s, const int& n, const char& c)
{
    // Compute number of characters to fill right
    int n_right  = n - s.length();
 
    // Set result
    std::string result = s + fill(std::string(1,c), n_right);
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Right justify string to achieve a length of n characters
 *
 * @param[in] s String to be right justified.
 * @param[in] n Requested total width.
 * @param[in] c Fill character.
 * @return Right justified string.
 *
 * Right justify string by adding @p c to the left to achieve a length of
 * @p n characters.
 ***************************************************************************/
std::string gammalib::right(const std::string& s, const int& n, const char& c)
{
    // Compute number of characters to fill right
    int n_left  = n - s.length();
 
    // Set result
    std::string result = fill(std::string(1,c), n_left) + s;
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Centre string to achieve a length of n characters
 *
 * @param[in] s String to be centred.
 * @param[in] n Requested total width.
 * @param[in] c Fill character.
 * @return Centred string.
 *
 * Centre string by adding @p c to the left and the right to achieve a length
 * of @p n characters.
 ***************************************************************************/
std::string gammalib::centre(const std::string& s, const int& n, const char& c)
{
    // Compute number of characters to fill left and right
    int n_right = (n-s.length()) / 2;
    int n_left  = n - s.length() - n_right;
 
    // Set result
    std::string result = fill(std::string(1,c), n_left) + s + 
                         fill(std::string(1,c), n_right);
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert string in parameter format
 *
 * @param[in] s String to be converted.
 * @param[in] indent Indentation of parameter.
 * @return Parameter string.
 *
 * Converts and string into the parameter format of type "s ......: " with a
 * total length of G_PARFORMAT_LENGTH.
 ***************************************************************************/
std::string gammalib::parformat(const std::string& s, const int& indent)
{
    // Compute number of characters to fill right. Do not clip the string if
    // it is too long since we do not want to loose information.
    int n_right = G_PARFORMAT_LENGTH - s.length() - 3 - indent;
 
    // Set result
    std::string result = " " + s + " " + fill(".", n_right) + ": ";
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert singular noun into number noun
 *
 * @param[in] noun Singular noun.
 * @param[in] number Number of instance of noun.
 * @return Converted noun.
 *
 * Converts a singular noun into a number noun by appending a "s" to the
 * noun if the @p number of the instances of the noun is not one.
 ***************************************************************************/
std::string gammalib::number(const std::string& noun, const int& number)
{
    // Copy input noun
    std::string result(noun);
 
    // Append "s" if number if not one
    if (number != 1) {
        result += "s";
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert number into order noun
 *
 * @param[in] number Number.
 * @return Order string derived from number.
 *
 * Converts a number into an order string, such as "0th", "1st", "2nd",
 * "3rd", "4th", etc. This function only works correctly on non-negative
 * numbers. No checking is done on whether the number is non-negative.
 ***************************************************************************/
std::string gammalib::order(const int& number)
{
    // Initialise result
    std::string result;
 
    // Set results
    if (number == 0) {
        result = "0th";
    }
    else if (number == 1) {
        result = "1st";
    }
    else if (number == 2) {
        result = "2nd";
    }
    else if (number == 3) {
        result = "3rd";
    }
    else {
        result = gammalib::str(number)+"th";
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Compute photon flux between two energies for a power law
 *
 * @param[in] emin Minimum energy.
 * @param[in] emax Maximum energy.
 * @param[in] epivot Pivot energy.
 * @param[in] gamma Spectral index.
 * @return Photon flux under power law.
 *
 * Analytically computes
 * \f[\int_{E_{\rm min}}^{E_{\rm max}} 
 *                             \left( E/E_{\rm pivot} \right)^{\gamma} dE\f]
 * where
 * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
 * energy, respectively, and
 * \f$E_{\rm pivot}\f$ is the pivot energy, and
 * \f$\gamma\f$ is the spectral index.
 ***************************************************************************/
double gammalib::plaw_photon_flux(const double& emin, const double& emax,
                                  const double& epivot, const double& gamma)
{
    // Initialise flux
    double flux = 0.0;
 
    // Continue only if emax > emin
    if (emax > emin) {
 
        // Compute photon flux. Computations dependend on the exponent. We
        // add here a kluge to assure numerical accuracy.
        double exponent = gamma + 1.0;
        if (std::abs(exponent) > 1.0e-11) {
            double xmin = emin   / epivot;
            double xmax = emax   / epivot;
            flux        = epivot / exponent * (std::pow(xmax, exponent) -
                                               std::pow(xmin, exponent));
        }
        else {
            double ratio = emax / emin;
            flux         = epivot * std::log(ratio);
        }
 
    } // endif: emax > emin
 
    // Return result
    return flux;
}
 
 
/***********************************************************************//**
 * @brief Compute energy flux between two energies for a power law
 *
 * @param[in] emin Minimum energy.
 * @param[in] emax Maximum energy.
 * @param[in] epivot Pivot energy.
 * @param[in] gamma Spectral index.
 * @return Energy flux under power law.
 *
 * Analytically computes
 * \f[\int_{E_{\rm min}}^{E_{\rm max}} 
 *                           \left( E/E_{\rm pivot} \right)^{\gamma} E dE\f]
 * where
 * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
 * energy, respectively, and
 * \f$E_{\rm pivot}\f$ is the pivot energy, and
 * \f$\gamma\f$ is the spectral index.
 ***************************************************************************/
double gammalib::plaw_energy_flux(const double& emin, const double& emax,
                                  const double& epivot, const double& gamma)
{
    // Initialise flux
    double flux = 0.0;
 
    // Continue only if emax > emin
    if (emax > emin) {
 
        // Compute energy flux. Computations dependend on the exponent. We
        // add here a kluge to assure numerical accuracy.
        double exponent = gamma + 2.0;
        double epivot2  = epivot * epivot;
        if (std::abs(exponent) > 1.0e-11) {
            double xmin = emin    / epivot;
            double xmax = emax    / epivot;
            flux        = epivot2 / exponent * (std::pow(xmax, exponent) -
                                                std::pow(xmin, exponent));
        }
        else {
            double ratio = emax / emin;
            flux         = epivot2 * (std::log(ratio));
        }
 
    } // endif: emax > emin
 
    // Return result
    return flux;
}
 
 
/***********************************************************************//**
 * @brief Computes log mean energy
 *
 * @param[in] a First energy.
 * @param[in] b Second energy.
 * @return Log mean energy.
 *
 * Computes the logarithmic mean energy
 * \f$10^{0.5 * (\log E_{\rm a} + \log E_{\rm b})}\f$
 * for two energies.
 ***************************************************************************/
GEnergy gammalib::elogmean(const GEnergy& a, const GEnergy& b)
{
    // Compute logarithmic mean energy
    GEnergy elogmean;
    double  eloga = a.log10MeV();
    double  elogb = b.log10MeV();
    elogmean.MeV(std::pow(10.0, 0.5 * (eloga + elogb)));
 
    // Return
    return elogmean;
}
 
 
/***********************************************************************//**
 * @brief Checks if directory exists
 *
 * @param[in] dirname Directory name.
 * @return True if directory exists, false otherwise.
 *
 * Checks if a directory exists. The function expands any environment
 * variable prior to checking.
 ***************************************************************************/
bool gammalib::dir_exists(const std::string& dirname)
{
    // Initialise result
    bool result = false;
 
    // Allocate file information structure
    struct stat info;
 
    // Get file information structure
    int ret = stat(gammalib::expand_env(dirname).c_str(), &info);
 
    // Check if we have a directory
    if (ret == 0 && S_ISDIR(info.st_mode)) {
        result = true;
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Checks if a substring is in a string
 *
 * @param[in] str String you want to search in.
 * @param[in] substring String you are looking for in @p str.
 * @return True if a string contains a sub string.
 *
 * Checks if substring is contained in str
 ***************************************************************************/
bool gammalib::contains(const std::string& str, const std::string& substring)
{
    // Initialise result
    bool result = false;
 
    // checks if substring is in str
    if (str.find(substring) != std::string::npos) {
        result = true;
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Checks if a string is contained in a vector of strings
 *
 * @param[in] strings Vector of strings you want to search in.
 * @param[in] string string you are looking for in strings.
 * @return True if a string is contained a vector.
 *
 * Checks if a string is contained in a vector of strings
 ***************************************************************************/
bool gammalib::contains(const std::vector<std::string>& strings,
                        const std::string&              string)
{
    // Compute result
    bool result = std::find(strings.begin(), strings.end(), string) !=
                  strings.end();
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Emits warning
 *
 * @param[in] origin Name of method that emits warning.
 * @param[in] message Warning message.
 *
 * Writes a warning to the console.
 ***************************************************************************/
void gammalib::warning(const std::string& origin,
                       const std::string& message)
{
    // Compile option: enable/disable warnings
    #if defined(G_WARNINGS)
 
    // Set warning message
    std::string warning = "+++ WARNING in " + origin + ": " + message;
 
    // Writes warning to the console
    std::cout << warning << std::endl;
 
    // End of compile option
    #endif
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Convert XML character references in string to characters
 *
 * @param[in] arg String containing XML character references.
 * @return String with character reference replaced by respective characters.
 *
 * @exception GException::invalid_argument
 *            Could not extract number from numerical character reference.
 *
 * Converts all character references found in a string in their respective
 * characters. For more information about XML character references read
 * http://en.wikipedia.org/wiki/List_of_XML_and_HTML_character_entity_references
 ***************************************************************************/
std::string gammalib::xml2str(const std::string& arg)
{
    // Initialise string
    std::string result;
 
    // Initialise position, lengths and flags
    size_t length = arg.length();
    size_t pos    = 0;
    size_t start  = 0;
    size_t stop   = 0;
    size_t len    = 0;      // Length of string preceeding char. reference
    bool   found  = false;
 
    // Loop over string
    while (pos < length) {
 
        // If we have not yet found a character reference then search for
        // the next one. If we did not find one we break here.
        if (!found) {
            start = arg.find("&", pos);
            if (start != std::string::npos) {
                len   = start - pos;
                pos   = start;
                found = true;
            }
            else {
                break;
            }
        }
 
        // If we have found a character reference then attach the string
        // that is preceeding the character reference to the result and
        // search for the end character of the character reference.
        if (found) {
 
            // Initialse successful decoding flag
            bool decoded = false;
 
            // Attach preceeding string to the result
            if (len > 0) {
                result += arg.substr(start-len, len);
            }
 
            // Search for the end character
            stop = arg.find(";", pos);
 
            // If the end character was found then try decoding it
            if (stop != std::string::npos) {
            
                // Extract the sequence between start and end character as
                // the character reference
                std::string cref = arg.substr(start+1, stop-start-1);
                len              = cref.length();
 
                // Check whether we have numerical character reference
                //TODO: Check that there are only valid characters in
                //      numerical field
                if (len >= 2 && cref[0] == '#') {
                    int number = -1;
                    if (cref[1] == 'x') {
                        number = (int)std::strtol(cref.substr(2,len-2).c_str(), NULL, 16);
                    }
                    else {
                        number = toint(cref.substr(1,len-1));
                    }
                    if (number != -1) {
                        result.push_back((char)number);
                        decoded = true;
                    }
                    else {
                        std::string msg = "Could not extract number from "
                                          "numerical character reference &"+
                                          cref+";";
                        throw GException::invalid_argument(G_XML2STR, msg);
                    }
                }
 
                // ... otherwise check for a character entity reference.
                // If we found out then push back the corresponding character
                // to the result string
                else {
                    if (cref == "quot") {
                        result.push_back((char)34);
                        decoded = true;
                    }
                    else if (cref == "amp") {
                        result.push_back((char)38);
                        decoded = true;
                    }
                    else if (cref == "apos") {
                        result.push_back((char)39);
                        decoded = true;
                    }
                    else if (cref == "lt") {
                        result.push_back((char)60);
                        decoded = true;
                    }
                    else if (cref == "gt") {
                        result.push_back((char)62);
                        decoded = true;
                    }
                } // endelse: checked for character entity reference
 
                // If decoding was successful then skip the character
                // reference, otherwise add ampersand to result and
                // continue searching for an ampersand
                if (decoded) {
                    found = false;
                    pos   = stop + 1;
                }
                else {
                    found   = false;
                    result += "&";
                    pos++;
                }
 
            } // endif: end of character reference found
 
            // ... otherwise add ampersand to result and continue
            // searching for an ampersand
            else {
                found   = false;
                result += "&";
                pos++;
            }
 
        } // endelse: search for end of character reference
 
    } // endwhile
 
    // Append any pending string to the result
    if (pos < length) {
        len     = length - pos;
        result += arg.substr(pos, len);
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Convert special characters in string to XML character references
 *
 * @param[in] arg String.
 * @return String with special characters replaced by character references.
 *
 * Converts all special characters found in a string into character
 * references. For more information about XML character references read
 * http://en.wikipedia.org/wiki/List_of_XML_and_HTML_character_entity_references
 ***************************************************************************/
std::string gammalib::str2xml(const std::string& arg)
{
    // Initialise string
    std::string result;
 
    // Loop over string
    for (int i = 0; i < arg.length(); ++i) {
 
        // Initialize string to add
        std::string character = arg.substr(i, 1);
 
        // Replace special characters
        if (character == "\"") {
            character = "&quot;";
        }
        else if (character == "&") {
            character = "&amp;";
        }
        else if (character == "'") {
            character = "&apos;";
        }
        else if (character == "<") {
            character = "&lt;";
        }
        else if (character == ">") {
            character = "&gt;";
        }
 
        // Append character
        result += character;
 
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Checks if parameter with given name in XML element exists
 *
 * @param[in] xml XML element.
 * @param[in] name Parameter name.
 * @return True if parameter exists, false otherwise.
 *
 * Checks whether a parameter with given @p name exists in XML element.
 ***************************************************************************/
bool gammalib::xml_has_par(const GXmlElement& xml, const std::string& name)
{
    // Initialize flag
    bool found = false;
 
    // Get number of elements
    int n = xml.elements("parameter");
 
    // Search for parameter with given name
    for (int i = 0; i < n; ++i) {
        const GXmlElement* element = xml.element("parameter", i);
        if (element->attribute("name") == name) {
            found = true;
            break;
        }
    }
 
    // Return
    return found;
}
 
 
/***********************************************************************//**
 * @brief Return pointer to parameter with given name in XML element
 *
 * @param[in] origin Method requesting parameter.
 * @param[in] xml XML element.
 * @param[in] name Parameter name.
 * @return Pointer to parameter XML element.
 *
 * @exception GException::invalid_value
 *            Invalid XML format encountered.
 *
 * Returns pointer to parameter with given @p name in XML element. If the
 * @p name is not found, a parameter with the given @p name is added. In
 * that respect the function differs from xml_get_par which does not add a
 * parameter element.
 *
 * The function checks for multiple occurences of a parameter and throws an
 * exception in case that more than one parameter with a given name is found.
 ***************************************************************************/
GXmlElement* gammalib::xml_need_par(const std::string& origin,
                                    GXmlElement&       xml,
                                    const std::string& name)
{
    // Initialize XML element pointer
    GXmlElement* par = NULL;
 
    // Number of elements
    int number = 0;
 
    // Get number of elements in XML element
    int n = xml.elements("parameter");
 
    // Search for parameter with given name
    for (int i = 0; i < n; ++i) {
        GXmlElement* element = xml.element("parameter", i);
        if (element->attribute("name") == name) {
            par = element;
            number++;
        }
    }
 
    // Create parameter if none was found
    if (number == 0) {
        par = static_cast<GXmlElement*>(xml.append(GXmlElement("parameter name=\""+name+"\"")));
        number++;
    }
 
    // Check that there are no multiple parameters
    gammalib::xml_check_par(origin, name, number);
 
    // Return
    return par;
}
 
 
/***********************************************************************//**
 * @brief Return pointer to parameter with given name in XML element
 *
 * @param[in] origin Method requesting parameter.
 * @param[in] xml XML element.
 * @param[in] name Parameter name.
 * @return Pointer to parameter XML element.
 *
 * @exception GException::invalid_value
 *            Invalid XML format encountered.
 *
 * Returns pointer to parameter with given @p name in XML element. The
 * function checks whether the parameter has been found and throws an
 * exception if no parameter or multiple occurences of a parameter with given
 * @p name are found.
 ***************************************************************************/
const GXmlElement* gammalib::xml_get_par(const std::string& origin,
                                         const GXmlElement& xml,
                                         const std::string& name)
{
    // Initialize XML element pointer
    const GXmlElement* par = NULL;
 
    // Number of elements
    int number = 0;
 
    // Get number of elements in XML element
    int n = xml.elements("parameter");
 
    // Search for parameter with given name
    for (int i = 0; i < n; ++i) {
        const GXmlElement* element = xml.element("parameter", i);
        if (element->attribute("name") == name) {
            par = element;
            number++;
        }
    }
 
    // Check that there is at least one parameter and that there are no
    // multiple parameters
    gammalib::xml_check_par(origin, name, number);
 
    // Return
    return par;
}
 
 
/***********************************************************************//**
 * @brief Return attribute value for a given parameter in XML element
 *
 * @param[in] origin Method requesting parameter.
 * @param[in] xml XML element.
 * @param[in] name Parameter name.
 * @param[in] attribute Attribute name.
 * @return Value of attribute.
 *
 * @exception GException::invalid_value
 *            Attribute not found.
 *
 * Returns the value of @p attribute of parameter @p name in XML element.
 * The function checks whether the parameter has been found and throws an
 * exception if no parameter or multiple occurences of a parameter with given
 * @p name are found. The function furthermore checks whether the attribute
 * exists. 
 ***************************************************************************/
std::string gammalib::xml_get_attr(const std::string& origin,
                                   const GXmlElement& xml,
                                   const std::string& name,
                                   const std::string& attribute)
{
    // Initialise attribute value
    std::string value = "";
 
    // Get parameter
    const GXmlElement* par = gammalib::xml_get_par(origin, xml, name);
 
    // Throw an exception if a parameter has not the requested attribute
    if (!par->has_attribute(attribute)) {
        std::string msg = "Attribute \""+attribute+"\" not found in XML "
                          "parameter \""+name+"\". Please verify the XML "
                          "format.";
        throw GException::invalid_value(origin, msg);
    }
 
    // Extract attribute
    value = par->attribute(attribute);
 
    // Return attribute value
    return value;
}
 
 
/***********************************************************************//**
 * @brief Checks number of parameters
 *
 * @param[in] origin Method performing the check.
 * @param[in] xml XML element.
 * @param[in] number Expected number of parameters.
 *
 * @exception GException::invalid_value
 *            Invalid XML format encountered.
 *
 * Checks the number of parameter in an XML element.
 ***************************************************************************/
void gammalib::xml_check_parnum(const std::string& origin,
                                const GXmlElement& xml,
                                const int&         number)
{
    // Throw exception if number of elements does not correspond to the
    // expected number
    if (xml.elements() != number) {
        std::string msg = "Number of "+gammalib::str(xml.elements())+
                          " child elements in XML file does not correspond "
                          "to expected number of "+gammalib::str(number)+
                          " elements. Please verify the XML format.";
        throw GException::invalid_value(origin, msg);
    }
 
    // Throw exception if number of "parameter" elements does not correspond
    // to the xpected number
    int npars = xml.elements("parameter");
    if (npars != number) {
        std::string msg = "Number of "+gammalib::str(npars)+" \"parameter\" "
                          "child elements in XML file does not correspond to "
                          "expected number of "+gammalib::str(number)+
                          " elements. Please verify the XML format.";
        throw GException::invalid_value(origin, msg);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Checks the model type
 *
 * @param[in] origin Method performing the check.
 * @param[in,out] xml XML element.
 * @param[in] type Expected model typeN.
 *
 * @exception GException::invalid_value
 *            Invalid XML format encountered.
 *
 * Checks the number of parameter in an XML element.
 ***************************************************************************/
void gammalib::xml_check_type(const std::string& origin,
                              GXmlElement&       xml,
                              const std::string& type)
{
    // If XML element has no model type then set model type
    if (xml.attribute("type") == "") {
        xml.attribute("type", type);
    }
 
    // Throw an exception if the model type is not the expected one
    if (xml.attribute("type") != type) {
        std::string msg = "Model type \""+xml.attribute("type")+"\" is not "
                          "the expected type \""+type+"\". Please verify "
                          "the XML format.";
        throw GException::invalid_value(origin, msg);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Checks whether a parameter has occured once
 *
 * @param[in] origin Method performing the check.
 * @param[in] name Parameter name.
 * @param[in] number Number of occurences of parameter.
 *
 * @exception GException::invalid_value
 *            Invalid XML format encountered.
 *
 * Throws an exception if a given parameter has not exactly occured once.
 * The exception text is adapted to the case that none or multiple parameters
 * have been found.
 ***************************************************************************/
void gammalib::xml_check_par(const std::string& origin,
                             const std::string& name,
                             const int&         number)
{
    // Throw case dependent exception
    if (number < 1) {
        std::string msg = "Parameter \""+name+"\" not found in XML element."
                          " Please verify the XML format.";
        throw GException::invalid_value(origin, msg);
    }
    else if (number > 1) {
        std::string msg = "Parameter \""+name+"\" found "+
                          gammalib::str(number)+" times in XML element."
                          " Please verify the XML format.";
        throw GException::invalid_value(origin, msg);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Expand file name provided as XML attribute for loading
 *
 * @param[in] xml XML element.
 * @param[in] filename File name.
 * @return Expanded file name.
 *
 * Expands file name provided as XML attribute for loading. The XML file
 * access path will be prepended to the file name access path if the file
 * name has not an absolute path and if the file name access path does not
 * start with the XML file access path.
 ***************************************************************************/
GFilename gammalib::xml_file_expand(const GXmlElement& xml,
                                    const std::string& filename)
{
    // Set file name
    GFilename fname(filename);
 
    // If the file name is not empty and is not an absolute path then we
    // assume that the filename is a relative file name with respect to the
    // XML file access path and we therefore prepend the XML file access path
    // to the file name
    if (!fname.is_empty()) {
        size_t access_length = fname.path().length();
        size_t xml_length    = xml.filename().path().length();
 
        // If filename has no access path and the XML file has an access path
        // then preprend the XML file access path
        if (access_length == 0) {
            if (xml_length > 0) {
                fname = xml.filename().path() + fname;
            }
        }
 
        // If filename assess path is not absolute then continue ...
        else if (fname.path()[0] != '/') {
 
            // If XML file access path is not contained in file name access
            // path then prepend the XML file access path
            if (access_length >= xml_length) {
                if (fname.path().compare(0, xml_length, xml.filename().path()) != 0) {
                    fname = xml.filename().path() + fname;
                }
            }
 
            // If file name access pass is shorter than XML file access path
            // then prepend the XML file access path
            else {
                fname = xml.filename().path() + fname;
            }
        }
    }
 
    // Return file name
    return fname;
}
 
 
/***********************************************************************//**
 * @brief Reduce file name provided for writing as XML attribute
 *
 * @param[in] xml XML element.
 * @param[in] filename File name.
 * @return Reduced file name.
 *
 * Reduces file name provided for writing as XML attribute. If the file name
 * is not empty and its access path starts has the same access path as the
 * XML file the XML file access path is stripped from the file name.
 ***************************************************************************/
GFilename gammalib::xml_file_reduce(const GXmlElement& xml,
                                    const std::string& filename)
{
    // Set file name
    GFilename fname(filename);
 
    // If the file name is not empty and the access path starts with the
    // same characters as the XML file access path then we assume that the
    // access path has been added and hence we strip it now
    if (!fname.is_empty()) {
        size_t access_length = fname.path().length();
        size_t xml_length    = xml.filename().path().length();
        if (xml_length > 0) {
            if (access_length == xml_length) {
                if (fname.path() == xml.filename().path()) {
                    fname = fname.file();
                }
            }
            else if (access_length > xml_length) {
                if (fname.path().compare(0, xml_length, xml.filename().path()) == 0) {
                    std::string relpath = fname.path().substr(xml_length, access_length-xml_length);
                    fname               = GFilename(relpath + fname.file());
                }
            }
        }
    }
 
    // Return file name
    return fname;
}
 
 
/***********************************************************************//**
 * @brief Extract name / value pair from XML node
 *
 * @param[in] node Pointer to XML node.
 * @param[out] name Name string.
 * @param[out] value Value string.
 *
 * Extracts a name / value pair from a XML node. If the XML node pointer is
 * NULL, the name and value strings will be empty.
 ***************************************************************************/
void gammalib::xml_get_name_value_pair(const GXmlNode* node,
                                       std::string&    name,
                                       std::string&    value)
{
    // Clear name and value strings
    name.clear();
    value.clear();
 
    // Continue only if node is valid
    if (node != NULL) {
 
        // Get name node and extract text content
        const GXmlNode* ptr = node->element("name", 0);
        if (ptr != NULL) {
            const GXmlText* text = static_cast<const GXmlText*>((*ptr)[0]);
            if (text != NULL) {
                name = text->text();
            }
        }
 
        // Get value node and extract text content
        ptr = node->element("value", 0);
        if (ptr != NULL) {
            const GXmlText* text = static_cast<const GXmlText*>((*ptr)[0]);
            if (text != NULL) {
                value = text->text();
            }
        }
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Return XML element for a given tag and name/value attribute
 *
 * @param[in] element Pointer to XML element in which to search for element.
 * @param[in] tag Tag of element to search for.
 * @param[in] name Attribute name of element to search for.
 * @param[in] value Attribute value of element to search for.
 * @return Pointer to XML element.
 *
 * Searches the XML structure from the specified @p element on for a given
 * @p tag that contains an attribute with @p name and @p value. This function
 * works recursevely down the entire XML structure and stops once the first
 * element with a given name/value attribute was found.
 ***************************************************************************/
const GXmlElement* gammalib::xml_get_element_by_attribute(const GXmlElement* element,
                                                          const std::string& tag,
                                                          const std::string& name,
                                                          const std::string& value)
{
    // Initialise result element
    const GXmlElement* result = NULL;
 
    // First check if element is the one searched for
    if (element->name() == tag) {
        if (element->has_attribute(name)) {
            if (element->attribute(name) == value) {
                result = element;
            }
        }
    }
 
    // If we have no result then search all childern
    if (result == NULL) {
 
        // Get number of elements in element
        int elements = element->elements();
 
        // Loop over elements
        for (int i = 0; i < elements; ++i) {
 
            // Get child element
            const GXmlElement* child = element->element(i);
 
            // Search in child element
            result = xml_get_element_by_attribute(child, tag, name, value);
 
            // If element was found then break
            if (result != NULL) {
                break;
            }
 
        } // endfor: looped over elements
 
    } // endif: searched all children
 
    // Return result element
    return result;
}
 
 
/***********************************************************************//**
 * @brief Checks whether a parameter has occured once
 *
 * @param[in] fd Socket file descriptor.
 * @param[out] buffer Buffer to hold data.
 * @param[in] len Maximum number of bytes to recv().
 * @param[in] flags Flags (as the fourth param to recv() ).
 * @param[in] timeout Timeout in milliseconds.
 * @return recv() error code, -2 == timeout
 *
 * This function implements the recv() function with a timeout. The timeout
 * is specified in milliseconds.
 ***************************************************************************/
int gammalib::recv(int fd, char *buffer, int len, int flags, int timeout)
{
    // Initialise error code with time out
    int error = -2;
 
    // Initialize the set
    fd_set readset;
    FD_ZERO(&readset);
    FD_SET(fd, &readset);
 
    // Initialize time out struct
    struct timeval tv;
    if (timeout >= 1000) {
        tv.tv_sec  = timeout/1000;
        tv.tv_usec = 0;
    }
    else {
        tv.tv_sec  = 0;
        tv.tv_usec = timeout*1000;
    }
 
    // select()
    int result = select(fd+1, &readset, NULL, NULL, &tv);
 
    // Check status
    if (result < 0) {
        error = -1;
    }
    else if (result > 0 && FD_ISSET(fd, &readset)) {
 
        // Initialise flags
        int iof = -1;
 
        // Set non-blocking mode
        if ((iof = fcntl(fd, F_GETFL, 0)) != -1) {
            fcntl(fd, F_SETFL, iof | O_NONBLOCK);
        }
 
        // Receive data
        result = ::recv(fd, buffer, len, flags);
 
        // Set as before
        if (iof != -1) {
            fcntl(fd, F_SETFL, iof);
        }
 
        // Set error
        error = result;
    }
 
    // Return error
    return error;
}
 
 
/***********************************************************************//**
 * @brief Returns length of circular arc within circular ROI
 *
 * @param[in] rad Circle radius in radians (<pi).
 * @param[in] dist Circle centre distance to ROI centre (<pi).
 * @param[in] cosdist Cosine of circle centre distance to ROI centre.
 * @param[in] sindist Sinus of circle centre distance to ROI centre.
 * @param[in] roi Radius of ROI in radians.
 * @param[in] cosroi Cosine of ROI radius.
 *
 * This method returns the arclength in radians of a circle of radius 'rad'
 * with a centre that is offset by 'dist' from the ROI centre, where the ROI
 * radius is given by 'roi'. To speed-up computations, the cosines and sinus
 * of 'roi' and 'psf' should be calculated by the client and be passed to
 * the method.
 ***************************************************************************/
double gammalib::roi_arclength(const double& rad,     const double& dist,
                               const double& cosdist, const double& sindist,
                               const double& roi,     const double& cosroi)
{
    // Declare arclength
    double arclength;
 
    // Handle special case that circle centre matches ROI centre
    if (dist == 0.0) {
        if (rad > roi) arclength = 0.0;             // Circle outside ROI
        else           arclength = gammalib::twopi; // Circle inside ROI
    }
 
    // ... otherwise circle and ROI centres are not identical
    else {
 
        // Handle special case that we evaluate exactly at the circle
        // centre. In this case we have in fact a point, and if this point
        // falls within the ROI it has a formal arclength of 2pi.
        //
        if (rad == 0.0) {
            if (dist > roi) arclength = 0.0;             // Circle centre outside ROI
            else            arclength = gammalib::twopi; // Circle centre inside ROI
        }
 
        // ... otherwise we have to handle the general case
        else {
            double d = roi - dist;
            if (-rad >= d) {
                arclength = 0.0;
            }
            else if (rad <= d) {
                arclength = gammalib::twopi;
            }
            else {
                double cosrad = std::cos(rad);
                double sinrad = std::sin(rad);
                double cosang = (cosroi - cosdist*cosrad) / (sindist*sinrad);
                arclength     = 2.0 * gammalib::acos(cosang);
                #if defined(G_CHECK_FOR_NAN)
                if (gammalib::is_infinite(arclength) ||
                    gammalib::is_notanumber(arclength)) {
                    std::cout << "roi_arclength: NaN occured";
                    std::cout << " rad=" << rad;
                    std::cout << " sinrad=" << sinrad;
                    std::cout << " cosrad=" << cosrad;
                    std::cout << " sindist=" << sindist;
                    std::cout << " cosdist=" << cosdist;
                    std::cout << " cosang=" << cosang;
                    std::cout << std::endl;
                }
                #endif
            }
        }
 
    } // endelse: Circle and ROI centres were not identical
 
    // Return arclength
    return arclength;
}
 
 
/***********************************************************************//**
 * @brief Compare two floating point values with tolerance
 *
 * @param[in] a First floating point value.
 * @param[in] b Second floating point value.
 * @param[in] tol Relative tolerance.
 * @return True if both floating point values are identical within
 *         the tolerance.
 *
 * This method tests whether @p a is identical to @p b within a given
 * relative tolerance @p tol. The test verifies whether
 *
 * \f[|a - b| \le |a \times tol|\f]
 *
 * for \f$a \ne 0\f$ or whether
 *
 * \f[|a - b| \le |b \times tol|\f]
 *
 * for \f$b \ne 0\f$. If both \f$a = 0\f$ and \f$b = 0\f$ the method
 * returns true.
 ***************************************************************************/
bool gammalib::compare(const double& a, const double& b, const double& tol)
{
    // Initialise identity
    bool identity = false;
 
    // Case A: a is not zero
    if (a != 0.0) {
        identity = std::abs(a - b) <= std::abs(a * tol);
    }
 
    // Case B: b is not zero
    else if (b != 0.0) {
        identity = std::abs(a - b) <= std::abs(b * tol);
    }
 
    // Case C: a and b are zero
    else {
        identity = true;
    }
 
    // Return identity
    return identity;
}
 
 
/***********************************************************************//**
 * @brief Return response to a HTTP query
 *
 * @param[in] host Host address.
 * @param[in] query Query string.
 * @return Response to a HTTP query.
 *
 * @exception GException::runtime_error
 *            Unable to open or to connect to socket.
 * @exception GException::invalid_argument
 *            Host not found.
 *
 * Returns response to a HTTP query. Be aware that this method will not
 * work for https servers.
 ***************************************************************************/
std::string gammalib::http_query(const std::string& host, const std::string& query)
{
    // Set constants
    const int portno = 80;
 
    // Initialise buffers
    char message[1024];
    char response[4096];
    memset(message,  0, sizeof(message));
    memset(response, 0, sizeof(response));
 
    // Build message string (C++98 compliant)
    #if ((defined __STDC_VERSION__ && __STDC_VERSION__ >= 199901L) || \
     (defined __cplusplus && __cplusplus >= 201103L))
    std::snprintf(message, sizeof(message), "GET /%s HTTP/1.1\r\nHost: %s\r\n\r\n",
                  query.c_str(), host.c_str());
    #else
    std::sprintf(message, "GET /%s HTTP/1.1\r\nHost: %s\r\n\r\n",
                 query.c_str(), host.c_str());
    #endif
 
    // Create the socket
    int sockfd = socket(AF_INET, SOCK_STREAM, 0);
    if (sockfd < 0) {
        std::string msg = "Unable to open socket.";
        throw GException::runtime_error(G_HTTP_QUERY, msg);
    }
 
    // Lookup IP address
    struct hostent* server = gethostbyname(host.c_str());
    if (server == NULL) {
        std::string msg = "Unable to find host \""+host+"\".";
        throw GException::invalid_argument(G_HTTP_QUERY, msg);
    }
 
    // Fill in structure
    struct sockaddr_in serveraddr;
    memset(&serveraddr, 0, sizeof(serveraddr));
    serveraddr.sin_family = AF_INET;
    memcpy(&serveraddr.sin_addr.s_addr, server->h_addr, server->h_length);
    serveraddr.sin_port   = htons(portno);
 
    // Connect the socket
    if (connect(sockfd, (struct sockaddr*)&serveraddr, sizeof(serveraddr)) < 0) {
        std::string msg = "Unable to connect to socket for host \""+host+"\".";
        throw GException::runtime_error(G_HTTP_QUERY, msg);
    }
 
    // Send request
    write(sockfd, message, strlen(message));
 
    // Receive response
    read(sockfd, response, sizeof(response)-1);
 
    // Close socket
    close(sockfd);
 
    // Return response
    return std::string(response);
}
 
 
#if defined(G_USE_CURL)
/***********************************************************************//**
 * @brief Return two-digit host country code
 *
 * @param[in] force_query Force query of host country code?
 * @return Host two-digit host country code.
 *
 * Returns two-digit host country code, either by reading the code from the
 * $HOME/.gamma/host-country file, or if the file does not exist, by using
 * the http query
 * http://ip-api.com/line/?fields=countryCode.
 *
 * The result is saved in a static variable, hence once the country code is
 * retrieved no queries will be executed anymore. The http query is issued
 * using the curl tool with a timeout of one second. If the query fails,
 * or the curl tool is not available, the method returns an empty string.
 *
 * If no host-country file exists and the method retrieved a two-digit
 * country code it will write the results into a host-country file.
 ***************************************************************************/
#else
/***********************************************************************//**
 * @brief Return two-digit host country code
 *
 * @param[in] force_query Force query of host country code?
 * @return Host two-digit host country code.
 *
 * Returns two-digit host country code, either by reading the code from the
 * $HOME/.gamma/host-country file, or if the file does not exist, by using
 * the http query
 * http://ip-api.com/line/?fields=countryCode.
 *
 * The result is saved in a static variable, hence once the country code is
 * retrieved no queries will be executed anymore. The http query is issued
 * using the http_query() method. If the query fails the method returns an
 * empty string.
 *
 * If no host-country file exists and the method retrieved a two-digit
 * country code it will write the results into a host-country file.
 ***************************************************************************/
#endif
std::string gammalib::host_country(const bool& force_query)
{
    // Initialise country as static variable
    static std::string country;
 
    // Continue only if host country is empty
    if (country.empty()) {
 
        // Attempt to get host country from $HOME/.gamma/host-country file
        GFilename filename = gammalib::gamma_filename("host-country");
 
        // Continue only if file exists
        if ((!force_query) && (access(filename.url().c_str(), R_OK) == 0)) {
 
            // OpenMP critical zone to read host-country file
            #pragma omp critical(GammaLib_host_country)
            {
 
                // Get file lock. Continue only in case of success.
                struct flock lock;
                lock.l_type   = F_RDLCK;  // Want a read lock
                lock.l_whence = SEEK_SET; // Want beginning of file
                lock.l_start  = 0;        // No offset, lock entire file ...
                lock.l_len    = 0;        // ... to the end
                lock.l_pid    = getpid(); // Current process ID
                int fd        = open(filename.url().c_str(), O_RDONLY);
                if (fd != -1) {
 
                    // Lock file
                    fcntl(fd, F_SETLKW, &lock);
 
                    // Open host-country file. Continue only if opening was
                    // successful.
                    FILE* fptr = fopen(filename.url().c_str(), "r");
                    if (fptr != NULL) {
 
                        // Allocate line buffer
                        const int n = 1000;
                        char      line[n];
 
                        // Read line
                        fgets(line, n-1, fptr);
 
                        // Close file
                        fclose(fptr);
 
                        // Extract country string
                        country = gammalib::strip_chars(std::string(line), "\n");
 
                    } // endif: host-country file opened successfully
 
                    // Unlock file
                    lock.l_type = F_UNLCK;
                    fcntl(fd, F_SETLK, &lock);
                    close(fd);
 
                } // endif: file locking successful
 
            } // end of OMP critial zone
 
        } // endif: file existed
 
        // ... otherwise there is no host-country file then get the
        // information from a geolocalisation web site
        else {
 
            // If curl should be used, the curl tool is available and the
            // country code is not set then determine the country code
            // from a geolocalisation web site
            #if defined(G_USE_CURL)
            #if defined(G_HAS_CURL)
            // If curl is available then use curl to infer the country code
            char command[256];
 
            // Setup curl command with timeout if possible to avoid lock up
            #if defined(G_HAS_PERL)
            sprintf(command, "a=$(perl -e 'alarm shift; exec \"curl --silent"
                             " http://ip-api.com/line/?fields=countryCode"
                             " 2>/dev/null\"' \"1\";);echo $a");
            #elif defined(G_HAS_TIMEOUT)
            sprintf(command, "timeout 1s curl --silent http://ip-api.com/line/"
                             "?fields=countryCode 2>/dev/null");
            #else
            sprintf(command, "curl --silent http://ip-api.com/line/?fields=countryCode"
                             " 2>/dev/null");
            #endif
 
            // Open the process with given 'command' for reading
            FILE* file = popen(command, "r");
 
            // Retrieve curl output
            char output[1024];
            if (fgets(output, sizeof(output)-1, file) != NULL) {
                country = strip_chars(std::string(output), "\n");
            }
 
            // Close process
            pclose(file);
            #endif
 
            // ... otherwise use the http_query method. Catch any exceptions.
            #else
            try {
                std::string response = http_query("ip-api.com", "line/?fields=countryCode");
                std::vector<std::string> lines = split(response, "\n");
                bool body = false;
                for (int i = 0; i < lines.size(); ++i) {
                    if (body) {
                        country = lines[i];
                        break;
                    }
                    if ((lines[i] == "\n") || (lines[i] == "\r") || (lines[i] == " ")) {
                        body = true;
                    }
                }
            } // endtry
            catch (const std::exception& e) {
                //
            }
            #endif
 
            // If we have a two-digit country code then write the code into
            // host-country file
            if (country.length() == 2) {
 
                // OpenMP critical zone to write host-country file
                #pragma omp critical(GammaLib_host_country)
                {
 
                    // Open host-country file, and in case of success, write
                    // country
                    FILE* fptr = fopen(filename.url().c_str(), "w");
                    if (fptr != NULL) {
                        fprintf(fptr, "%s\n", country.c_str());
                        fclose(fptr);
                    }
 
                } // end of OMP critial zone
 
            } // endif: had no two digit country code
 
        } // endelse: there was no host country file
 
    } // endif: host country was empty
 
    // Return country
    return country;
}
 
 
/***********************************************************************//**
 * @brief Returns filename in .gamma directory
 *
 * @param[in] name Name of file in .gamma directory.
 * @return Filename in .gamma directory.
 *
 * Returns the filename of the name of a file in the .gamma directory.
 ***************************************************************************/
GFilename gammalib::gamma_filename(const std::string& name)
{
    // Initialise lock filename
    GFilename filename;
 
    // First try accessing the user ID to get the home directory
    uid_t uid         = geteuid();
    struct passwd* pw = getpwuid(uid);
    if (pw != NULL) {
        filename = std::string(pw->pw_dir) + "/.gamma/" + name;
    }
 
    // ... otherwise try using the $HOME environment variable
    else {
 
        // First attempt fetching $HOME environment variable
        char* home_ptr = std::getenv("HOME");
        if (home_ptr != NULL) {
            filename = std::string(home_ptr) + "/.gamma/" + name;
        }
 
        // ... otherwise use ~ as home directory
        else {
            filename = "~/.gamma/" + name;
        }
    }
 
    // Return filename
    return filename;
}
 
 
/***********************************************************************//**
 * @brief Get current resident set size (physical memory use) in Bytes
 *
 * @return Physical memory use in Bytes.
 ***************************************************************************/
size_t gammalib::get_current_rss(void)
{
    // Initialize resident set size
    size_t rss = 0;
 
    // Determine resident set size (architecture dependent)
    // OSX
    #if defined(__APPLE__) && defined(__MACH__)
    #ifdef MACH_TASK_BASIC_INFO
    struct mach_task_basic_info info;
    mach_msg_type_number_t      infoCount = MACH_TASK_BASIC_INFO_COUNT;
    if (task_info(mach_task_self( ), MACH_TASK_BASIC_INFO,
        (task_info_t)&info, &infoCount) == KERN_SUCCESS) {
        rss = (size_t)info.resident_size;
    }
    #else
    struct task_basic_info info;
    mach_msg_type_number_t info_count = TASK_BASIC_INFO_COUNT;
    if (task_info(mach_task_self(), TASK_BASIC_INFO,
        (task_info_t)&info, &info_count) == KERN_SUCCESS) {
        rss = (size_t)info.resident_size;
    }
    #endif
    // Linux
    #elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
    FILE* fp = NULL;
    if ((fp = fopen( "/proc/self/statm", "r" )) != NULL) {
        if (fscanf( fp, "%*s%ld", &rss ) == 1) {
            rss *= (size_t)sysconf(_SC_PAGESIZE);
        }
        fclose(fp);
    }
    // AIX, BSD, Solaris, and Unknown OS
    #else
    rss = 0;
    #endif
 
    // Return resident set size
    return rss;
}
 
 
/***********************************************************************//**
 * @brief Get current clock in seconds
 *
 * @return Current clock in seconds.
 ***************************************************************************/
double gammalib::get_current_clock(void)
{
    // Get current clock in seconds
    #ifdef _OPENMP
    double time = omp_get_wtime();
    #else
    double time = double(clock()) / (double)CLOCKS_PER_SEC;
    #endif
 
    // Return time
    return time;
}