CECoordinates.cpp

Go to the documentation of this file.

/***************************************************************************
 *  CECoordinates.cpp: CppEphem                                            *
 * ----------------------------------------------------------------------- *
 *  Copyright © 2016-2019 JCardenzana                                      *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  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/>.  *
 *                                                                         *
 ***************************************************************************/
 
#include "CECoordinates.h"
#include "CEObserver.h"
 
/**********************************************************************/
CECoordinates::CECoordinates()
{
    init_members();
}
 
 
/**********************************************************************/
CECoordinates::CECoordinates(const CEAngle& xcoord, 
                             const CEAngle& ycoord,
                             const CECoordinateType& coord_type)
{
    init_members();
    CECoordinates::SetCoordinates(xcoord, ycoord, coord_type);
}
 
 
/**********************************************************************/
CECoordinates::CECoordinates(const std::vector<double>& xcoord,
                             const std::vector<double>& ycoord,
                             const CECoordinateType& coord_type)
{
    init_members();
    CEAngle x;
    if ((coord_type == CECoordinateType::GALACTIC) || 
        (coord_type == CECoordinateType::OBSERVED)) {
        x = CEAngle::Dms(xcoord);
    } else {
        x = CEAngle::Hms(xcoord);
    }
    CEAngle y = CEAngle::Dms(ycoord);
    CECoordinates::SetCoordinates(x, y, coord_type);
}
 
 
/**********************************************************************/
CECoordinates::CECoordinates(const CECoordinateType& coord_type)
{
    init_members();
    coord_type_ = coord_type;
}
 
/**********************************************************************/
CECoordinates::CECoordinates(const CECoordinates& other)
{
    free_members();
    init_members();
    copy_members(other);
}
 
/**********************************************************************/
CECoordinates::~CECoordinates()
{}
 
 
/**********************************************************************/
CECoordinates& CECoordinates::operator=(const CECoordinates& other)
{
    if (this != &other) {
        free_members();
        init_members();
        copy_members(other);
    }
    return *this;
}
 
 
/**********************************************************************/
void CECoordinates::CIRS2ICRS(double input_ra, double input_dec,
                              double *return_ra, double *return_dec,
                              const CEDate& date, const CEAngleType& angle_type)
{
    // Convert to radians if that's what is passed
    if (angle_type==CEAngleType::DEGREES) {
        input_ra  *= DD2R ;
        input_dec *= DD2R ;
    }
    
    // Use the sofa library to convert these coordinates
    double eo ; // Equation of the origins
    double tdb1(0.0);
    double tdb2(0.0);
    CEDate::UTC2TDB(date.MJD(), &tdb1, &tdb2);
    iauAtic13(input_ra, input_dec, tdb1, tdb2, return_ra, return_dec, &eo) ;
 
    // Subtract the eo from RA if J2000 coordinates are desired
    //*return_ra -= eo ;
    
    // Return the coordinates in the requested units
    if (angle_type == CEAngleType::DEGREES) {
        *return_ra  *= DR2D ;
        *return_dec *= DR2D ;
    }
}
 
/**********************************************************************/
void CECoordinates::CIRS2Galactic(double input_ra, double input_dec, double *glon, double *glat,
                                  const CEDate& date, const CEAngleType& angle_type)
{
    // Convert to radians if that's what is passed
    if (angle_type==CEAngleType::DEGREES) {
        input_ra  *= DD2R ;
        input_dec *= DD2R ;
    }
    
    // In order to do this with the sofa package, we must first
    // convert from CIRS -> ICRS
    double ICRS_ra(0.0), ICRS_dec(0.0) ;
    CIRS2ICRS(input_ra, input_dec, &ICRS_ra, &ICRS_dec, date, CEAngleType::RADIANS);
    
    // Now we can convert to galactic
    ICRS2Galactic(ICRS_ra, ICRS_dec, glon, glat, CEAngleType::RADIANS) ;
    
    // Convert to the desired units
    if (angle_type == CEAngleType::DEGREES) {
        *glon *= DR2D ;
        *glat *= DR2D ;
    }
}
 
/**********************************************************************/
int CECoordinates::CIRS2Observed(double ra, double dec,
                                 double *az, double *zen,
                                 const CEDate&     date,
                                 const CEObserver& observer,
                                 const CEAngleType& angle_type,
                                 double wavelength,
                                 double *observed_ra,
                                 double *observed_dec,
                                 double *hour_angle)
{
    // If we've passed in angles that are in degrees we need to convert to radians
    if (angle_type == CEAngleType::DEGREES) {
        ra  *= DD2R ;
        dec *= DD2R ;
    }
    
    // Setup the observed RA, Dec and hour_angle variables
    double temp_ra, temp_dec, temp_hour_angle ;
    // If values were passed, point these variables at the passed ones
    if (observed_ra  == nullptr) observed_ra  = &temp_ra;
    if (observed_dec == nullptr) observed_dec = &temp_dec;
    if (hour_angle   == nullptr) hour_angle   = &temp_hour_angle;
    
    int err_code = CIRS2Observed(ra, dec,
                                 az, zen,
                                 date.JD(),
                                 observer.Longitude_Rad(),
                                 observer.Latitude_Rad(),
                                 observer.Elevation_m(),
                                 observer.Pressure_hPa(),
                                 observer.Temperature_C(),
                                 observer.RelativeHumidity(),
                                 date.dut1(),
                                 date.xpolar(),
                                 date.ypolar(),
                                 wavelength,
                                 observed_ra, 
                                 observed_dec, 
                                 hour_angle) ;
    
    // Now convert back to degrees if that's what we were passed
    if (angle_type == CEAngleType::DEGREES) {
        *az              *= DR2D;
        *zen             *= DR2D;
        *observed_ra     *= DR2D;
        *observed_dec    *= DR2D;
        *hour_angle      *= DR2D;
    }
    
    return err_code ;
}
 
/**********************************************************************/
void CECoordinates::ICRS2CIRS(double input_ra, double input_dec,
                              double *return_ra, double *return_dec,
                              const CEDate& date, const CEAngleType& angle_type)
{
    // Convert to radians if necessary
    if (angle_type==CEAngleType::DEGREES) {
        input_ra  *= DD2R ;
        input_dec *= DD2R ;
    }
    
    // Store the equation of the origins
    double eo ; // Equation of the origins
    
    // Use the sofa library to convert these coordinates
    try {
        double tdb1(0.0);
        double tdb2(0.0);
        CEDate::UTC2TDB(date.MJD(), &tdb1, &tdb2);
        iauAtci13(input_ra, input_dec, 
                  0.0, 0.0, 0.0, 0.0, 
                  tdb1, tdb2, 
                  return_ra, return_dec, &eo) ;
        
        // Subtract the equation of the origins if J2000 coordinates are desired
        //*return_ra -= eo ;
    } catch (std::exception &e) {
        throw CEException::sofa_exception("CECoordinates::ICRS2CIRS",
                                          "iauAtci13", 
                                          "Exception thrown");
    }
    
    // Convert the returned coordinates to the correct angle type
    if (angle_type == CEAngleType::DEGREES) {
        *return_ra *= DR2D ;
        *return_dec *= DR2D ;
    }
    
    return ;
}
 
/**********************************************************************/
void CECoordinates::ICRS2Galactic(double input_ra, double input_dec, 
                                  double *glon, double *glat,
                                  const CEAngleType& angle_type)
{
    // Convert to radians if necessary
    if (angle_type == CEAngleType::DEGREES) {
        input_ra *= DD2R ;
        input_dec *= DD2R ;
    }
    
    // Use the sofa method to convert the coordinates
    iauIcrs2g(input_ra, input_dec, glon, glat) ;
    
    // Convert back to degrees if necessary
    if (angle_type == CEAngleType::DEGREES) {
        *glon *= DR2D ;
        *glat *= DR2D ;
    }
}
 
/**********************************************************************/
int CECoordinates::ICRS2Observed(double ra, double dec,
                                 double *az, double *zen,
                                 const CEDate&     date,
                                 const CEObserver& observer,
                                 const CEAngleType& angle_type,
                                 double wavelength,
                                 double *observed_ra,
                                 double *observed_dec,
                                 double *hour_angle)
{
    // Setup the observed RA, Dec and hour_angle variables
    double temp_ra(0.0);
    double temp_dec(0.0);
    double temp_hour_angle(0.0);
 
    // If values were not passed, point them at the temporary ones
    if (observed_ra  == nullptr) observed_ra  = &temp_ra;
    if (observed_dec == nullptr) observed_dec = &temp_dec;
    if (hour_angle   == nullptr) hour_angle   = &temp_hour_angle;
 
    // If we've passed in angles that are in degrees we need to convert to radians
    if (angle_type == CEAngleType::DEGREES) {
        ra  *= DD2R ;
        dec *= DD2R ;
    }
 
    // Call the necessary sofa method
    int err_code = ICRS2Observed(ra, dec,
                                 az, zen,
                                 date.JD(),
                                 observer.Longitude_Rad(),
                                 observer.Latitude_Rad(),
                                 observer.Elevation_m(),
                                 observer.Pressure_hPa(),
                                 observer.Temperature_C(),
                                 observer.RelativeHumidity(),
                                 date.dut1(),
                                 date.xpolar(),
                                 date.ypolar(),
                                 wavelength,
                                 observed_ra, observed_dec, hour_angle) ;
 
    // Now convert back to degrees if that's what we were passed
    if (angle_type == CEAngleType::DEGREES) {
        *az           *= DR2D ;
        *zen          *= DR2D ;
        *observed_ra  *= DR2D ;
        *observed_dec *= DR2D ;
        *hour_angle   *= DR2D ;
    }
    
    return err_code ;
}
 
 
/*********************************************************
 * Convert GALACTIC Coordinates
 *********************************************************/
 
/**********************************************************************/
void CECoordinates::Galactic2CIRS(double glon, double glat, double *ra, double *dec,
                                  const CEDate& date, const CEAngleType& angle_type)
{
    // Check for degrees
    if (angle_type == CEAngleType::DEGREES) {
        glon *= DD2R ;
        glat *= DD2R ;
    }
    
    // Do the Galactic -> ICRS converstion
    Galactic2ICRS(glon, glat, ra, dec, CEAngleType::RADIANS) ;
    
    // Now convert ICRS -> CIRS
    double tmp_ra(*ra);
    double tmp_dec(*dec);
    ICRS2CIRS(tmp_ra, tmp_dec, ra, dec, date, CEAngleType::RADIANS) ;
    
    // Now make sure to return the coordinates in the correct format
    if (angle_type == CEAngleType::DEGREES) {
        *ra  *= DR2D ;
        *dec *= DR2D ;
    }
}
 
/**********************************************************************/
void CECoordinates::Galactic2ICRS(double glon, double glat, 
                                  double *ra, double *dec,
                                  const CEAngleType& angle_type)
{
    // Check for degrees
    if (angle_type == CEAngleType::DEGREES) {
        glon *= DD2R ;
        glat *= DD2R ;
    }
    
    // Do the Galactic -> ICRS converstion
    iauG2icrs(glon, glat, ra, dec) ;
    
    // Now make sure to return the coordinates in the correct format
    if (angle_type == CEAngleType::DEGREES) {
        *ra  *= DR2D ;
        *dec *= DR2D ;
    }
}
 
/**********************************************************************/
int CECoordinates::Galactic2Observed(double glon, double glat,
                                     double *az, double *zen,
                                     const CEDate&     date,
                                     const CEObserver& observer,
                                     const CEAngleType& angle_type,
                                     double wavelength,
                                     double *observed_glon,
                                     double *observed_glat,
                                     double *hour_angle)
{
    // Convert angles to radians if necessary
    if (angle_type == CEAngleType::DEGREES) {
        glon *= DD2R ;
        glat *= DD2R ;
    }
    
    // Do the conversion
    int error_code = Galactic2Observed(glon, glat,
                                       az, zen,
                                       date.JD(),
                                       observer.Longitude_Rad(),
                                       observer.Latitude_Rad(),
                                       observer.Elevation_m(),
                                       observer.Pressure_hPa(),
                                       observer.Temperature_C(),
                                       observer.RelativeHumidity(),
                                       date.dut1(),
                                       date.xpolar(),
                                       date.ypolar(),
                                       wavelength,
                                       observed_glon, observed_glat, hour_angle) ;
    
    // Convert back to degrees if necessary
    if (angle_type == CEAngleType::DEGREES) {
        *az  *= DR2D ;
        *zen *= DR2D ;
        if (observed_glon != nullptr) *observed_glon *= DR2D ;
        if (observed_glat != nullptr) *observed_glat *= DR2D ;
        if (hour_angle    != nullptr) *hour_angle    *= DR2D ;
    }
    
    return error_code;
}
 
/*********************************************************
 * Convert OBSERVED Coordinates
 *********************************************************/
 
 
/**********************************************************************/
int CECoordinates::Observed2CIRS(double az, double zen,
                            double *ra, double *dec,
                            const CEDate&     date,
                            const CEObserver& observer,
                            const CEAngleType& angle_type)
{
    // Convert to radians if necessary
    if (angle_type == CEAngleType::DEGREES) {
        az *= DD2R ;
        zen *= DD2R ;
    }
    
    // Call the raw method
    int err_code = Observed2CIRS(az, zen, ra, dec, date.JD(),
                  observer.Longitude_Rad(), observer.Latitude_Rad(),
                  observer.Elevation_m(),
                  observer.Pressure_hPa(), observer.Temperature_C(),
                  observer.RelativeHumidity(),
                  date.dut1(),
                  date.xpolar(), 
                  date.ypolar(),
                  observer.Wavelength_um()) ;
    
    // Convert back to degrees if necessary
    if (angle_type == CEAngleType::DEGREES) {
        *ra *= DR2D ;
        *dec *= DR2D ;
    }
    
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::Observed2ICRS(double az, double zen,
                                 double *ra, double *dec,
                                 const CEDate&     date,
                                 const CEObserver& observer,
                                 const CEAngleType& angle_type)
{
    // Convert to radians if necessary
    if (angle_type == CEAngleType::DEGREES) {
        az *= DD2R ;
        zen *= DD2R ;
    }
    
    // Call the raw method
    int err_code = Observed2ICRS(az, zen, ra, dec, 
                  date.JD(),
                  observer.Longitude_Rad(), 
                  observer.Latitude_Rad(),
                  observer.Elevation_m(),
                  observer.Pressure_hPa(), 
                  observer.Temperature_C(),
                  observer.RelativeHumidity(),
                  date.dut1(),
                  date.xpolar(), 
                  date.ypolar(),
                  observer.Wavelength_um()) ;
    
    // Convert back to degrees if necessary
    if (angle_type == CEAngleType::DEGREES) {
        *ra *= DR2D ;
        *dec *= DR2D ;
    }
    
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::Observed2Galactic(double az, double zen,
                                 double *glon, double *glat,
                                 const CEDate&     date,
                                 const CEObserver& observer,
                                 const CEAngleType& angle_type)
{
    // Convert to radians if necessary
    if (angle_type == CEAngleType::DEGREES) {
        az *= DD2R ;
        zen *= DD2R ;
    }
    
    // Call the raw method
    int err_code = Observed2Galactic(az, zen, glon, glat, date.JD(),
                  observer.Longitude_Rad(), observer.Latitude_Rad(),
                  observer.Elevation_m(),
                  observer.Pressure_hPa(), observer.Temperature_C(),
                  observer.RelativeHumidity(),
                  date.dut1(),
                  date.xpolar(), date.ypolar(),
                  observer.Wavelength_um()) ;
    
    // Convert back to degrees if necessary
    if (angle_type == CEAngleType::DEGREES) {
        *glon *= DR2D ;
        *glat *= DR2D ;
    }
    
    return err_code ;
}
 
 
/*********************************************************
 * Main routines for converting TO OBSERVED Coordinates
 *********************************************************/
 
/**********************************************************************/
int CECoordinates::CIRS2Observed(double ra, double dec,
                                 double *az, double *zen,
                                 double julian_date,
                                 double longitude,
                                 double latitude,
                                 double elevation_m,
                                 double pressure_hPa,
                                 double temperature_celsius,
                                 double relative_humidity,
                                 double dut1,
                                 double xp, double yp,
                                 double wavelength_um,
                                 double *observed_ra,
                                 double *observed_dec,
                                 double *hour_angle)
{
    // Setup the observed RA, Dec and hour_angle variables
    double temp_ra(0.0);
    double temp_dec(0.0);
    double temp_hour_angle(0.0);
 
    // If values were not passed, point them at the temporary ones
    if (observed_ra  == nullptr) observed_ra  = &temp_ra;
    if (observed_dec == nullptr) observed_dec = &temp_dec;
    if (hour_angle   == nullptr) hour_angle   = &temp_hour_angle;
 
    // Call the necessary sofa method
    int err_code = 0;
    try {
        CEDate date(julian_date, CEDateType::JD);
        err_code = iauAtio13(ra, dec,
                             CEDate::GetMJD2JDFactor(), date.MJD(),
                             dut1,
                             longitude,
                             latitude,
                             elevation_m,
                             xp, yp,
                             pressure_hPa,
                             temperature_celsius,
                             relative_humidity,
                             wavelength_um,
                             az, zen,
                             hour_angle, 
                             observed_dec, 
                             observed_ra) ;
    } catch (std::exception &e) {
        throw CEException::sofa_exception("CECoordinates::CIRS2Observed",
                                          "iauAtio13",
                                          e.what());
    }
 
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::Observed2CIRS(double az, double zen,
                                 double *ra, double *dec,
                                 double julian_date,
                                 double longitude,
                                 double latitude,
                                 double elevation_m,
                                 double pressure_hPa,
                                 double temperature_celsius,
                                 double relative_humidity,
                                 double dut1,
                                 double xp, double yp,
                                 double wavelength_um)
{
    CEDate date(julian_date, CEDateType::JD);
    int err_code = iauAtoi13("A", az, zen,
                             CEDate::GetMJD2JDFactor(), date.MJD(),
                             dut1,
                             longitude, latitude,
                             elevation_m,
                             xp, yp,
                             pressure_hPa,
                             temperature_celsius,
                             relative_humidity,
                             wavelength_um,
                             ra, dec) ;
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::ICRS2Observed(double ra, double dec,
                                 double *az, double *zen,
                                 double julian_date,
                                 double longitude,
                                 double latitude,
                                 double elevation_m,
                                 double pressure_hPa,
                                 double temperature_celsius,
                                 double relative_humidity,
                                 double dut1,
                                 double xp, double yp,
                                 double wavelength_um,
                                 double *observed_ra,
                                 double *observed_dec,
                                 double *hour_angle)
{
    // Setup the observed RA, Dec and hour_angle variables
    double temp_ra(0.0);
    double temp_dec(0.0);
    double temp_hour_angle(0.0);
 
    // If values were passed, point these variables at the passed ones
    if (observed_ra  == nullptr) observed_ra  = &temp_ra;
    if (observed_dec == nullptr) observed_dec = &temp_dec;
    if (hour_angle   == nullptr) hour_angle   = &temp_hour_angle;
    
    // First convert the ICRS coordinates to CIRS coordinates
    CEDate date(julian_date, CEDateType::JD) ;
    ICRS2CIRS(ra, dec, &ra, &dec, date) ;
    
    // Call the necessary sofa method
    int err_code = CIRS2Observed(ra, dec,
                                 az, zen,
                                 julian_date,
                                 longitude,
                                 latitude,
                                 elevation_m,
                                 pressure_hPa,
                                 temperature_celsius,
                                 relative_humidity,
                                 dut1, xp, yp,
                                 wavelength_um,
                                 observed_ra,
                                 observed_dec,
                                 hour_angle) ;
    
    // Convert the apparent CIRS RA,Dec to ICRS RA,Dec
    CIRS2ICRS(*observed_ra, *observed_dec, observed_ra, observed_dec, 
              date, CEAngleType::RADIANS) ;
    
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::Observed2ICRS(double az, double zen,
                                 double *ra, double *dec,
                                 double julian_date,
                                 double longitude,
                                 double latitude,
                                 double elevation_m,
                                 double pressure_hPa,
                                 double temperature_celsius,
                                 double relative_humidity,
                                 double dut1,
                                 double xp, double yp,
                                 double wavelength_um)
{
    // Convert from Observed -> CIRS
    int err_code = Observed2CIRS(az, zen,
                                 ra, dec,
                                 julian_date,
                                 longitude, latitude,
                                 elevation_m, pressure_hPa,
                                 temperature_celsius, relative_humidity,
                                 dut1, xp, yp, wavelength_um) ;
    if (err_code == 0) {
        // Convert from CIRS -> ICRS
        CIRS2ICRS(*ra, *dec, ra, dec,
                  CEDate(julian_date,CEDateType::JD),
                  CEAngleType::RADIANS) ;
    }
    
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::Galactic2Observed(double glon, double glat,
                                     double *az, double *zen,
                                     double julian_date,
                                     double longitude,
                                     double latitude,
                                     double elevation_m,
                                     double pressure_hPa,
                                     double temperature_celsius,
                                     double relative_humidity,
                                     double dut1,
                                     double xp, double yp,
                                     double wavelength,
                                     double *observed_glon,
                                     double *observed_glat,
                                     double *hour_angle)
{
    // Setup the observed RA, Dec and hour_angle variables
    double temp_glon(0.0);
    double temp_glat(0.0);
    double temp_ra(0.0);
    double temp_dec(0.0);
    double temp_hour_angle(0.0);
 
    // If values were passed, point these variables at the passed ones
    if (observed_glon == nullptr) observed_glon = &temp_glon;
    if (observed_glat == nullptr) observed_glat = &temp_glat;
    if (hour_angle    == nullptr) hour_angle    = &temp_hour_angle;
    
    // Convert GALACTIC to ICRS
    double ra(0.0);
    double dec(0.0);
    CEDate date(julian_date, CEDateType::JD) ;
    Galactic2CIRS(glon, glat, &ra, &dec, date) ;
 
    // Call the necessary sofa method
    int err_code = CIRS2Observed(ra, dec,
                                 az, zen,
                                 julian_date,
                                 longitude,
                                 latitude,
                                 elevation_m,
                                 pressure_hPa,
                                 temperature_celsius,
                                 relative_humidity,
                                 dut1, xp, yp,
                                 wavelength,
                                 &temp_ra,
                                 &temp_dec,
                                 &temp_hour_angle) ;
 
    // Convert the apparent RA,Dec to galactic longitude,latitude
    CIRS2Galactic(temp_ra, temp_dec, observed_glon, observed_glat, date) ;
    
    return err_code ;
}
 
/**********************************************************************/
int CECoordinates::Observed2Galactic(double az, double zen,
                                 double *glon, double *glat,
                                 double julian_date,
                                 double longitude,
                                 double latitude,
                                 double elevation_m,
                                 double pressure_hPa,
                                 double temperature_celsius,
                                 double relative_humidity,
                                 double dut1,
                                 double xp, double yp,
                                 double wavelength_um)
{
    // Convert from Observed -> ICRS
    int err_code = Observed2ICRS(az, zen,
                                 glon, glat,
                                 julian_date,
                                 longitude, latitude,
                                 elevation_m, pressure_hPa,
                                 temperature_celsius, relative_humidity,
                                 dut1, xp, yp, wavelength_um) ;
    if (err_code == 0) {
        // Convert from ICRS -> Galactic
        ICRS2Galactic(*glon, *glat, glon, glat,
                      CEAngleType::RADIANS) ;
    }
    
    return err_code ;
}
 
/**********************************************************************/
CECoordinates CECoordinates::GetObservedCoords(const double& julian_date,
                                const double& longitude,
                                const double& latitude,
                                const double& elevation_m,
                                const double& pressure_hPa,
                                const double& temperature_celsius,
                                const double& relative_humidity,
                                const double& dut1,
                                const double& xp, const double& yp,
                                const double& wavelength) const
{
    // Preliminary variables
    double azimuth(0);
    double zenith(0);
    double observed1(0);
    double observed2(0);
    double observed3(0);
    
    if (coord_type_ == CECoordinateType::CIRS) {
        // Convert CIRS to Observed
        CIRS2Observed(xcoord_, ycoord_, &azimuth, &zenith,
                      julian_date, longitude, latitude,
                      elevation_m, pressure_hPa, temperature_celsius,
                      relative_humidity, dut1, xp, yp, wavelength,
                      &observed1, &observed2, &observed3) ;
    } else if (coord_type_ == CECoordinateType::ICRS) {
        // Convert ICRS to Observed
        ICRS2Observed(xcoord_, ycoord_, &azimuth, &zenith,
                      julian_date, longitude, latitude,
                      elevation_m, pressure_hPa, temperature_celsius,
                      relative_humidity, dut1, xp, yp, wavelength,
                      &observed1, &observed2, &observed3) ;
    } else if (coord_type_ == CECoordinateType::GALACTIC) {
        // Convert Galactic to Observed
        Galactic2Observed(xcoord_, ycoord_, &azimuth, &zenith,
                          julian_date, longitude, latitude,
                          elevation_m, pressure_hPa, temperature_celsius,
                          relative_humidity, dut1, xp, yp, wavelength,
                          &observed1, &observed2, &observed3) ;
    } else if (coord_type_ == CECoordinateType::OBSERVED) {
        azimuth = xcoord_;
        zenith  = ycoord_; 
    }
    
    // Create the CECoordinates object to be returned
    return CECoordinates(azimuth, zenith, CECoordinateType::OBSERVED) ;
}
 
/**********************************************************************/
CECoordinates CECoordinates::GetObservedCoords(const CEDate& date,
                                               const CEObserver& observer) const
{
    return GetObservedCoords(date,
                             observer.Longitude_Rad(),
                             observer.Latitude_Rad(),
                             observer.Elevation_m(),
                             observer.Pressure_hPa(),
                             observer.Temperature_C(),
                             observer.RelativeHumidity(),
                             date.dut1(), 
                             date.xpolar(), 
                             date.ypolar(), 
                             observer.Wavelength_um());
}
 
 
/**********************************************************************/
CEAngle CECoordinates::AngularSeparation(const CECoordinates& coords) const
{
    return AngularSeparation(*this, coords);
}
 
/**********************************************************************/
CEAngle CECoordinates::AngularSeparation(const CECoordinates& coords1, 
                                         const CECoordinates& coords2)
{
    // Make sure the coordinates are in the same frame
    if (coords1.GetCoordSystem() != coords2.GetCoordSystem()) {
        throw CEException::invalid_value("CECoordinates::AngularSeparation(CECoords&, CECoords&)",
                                         "Supplied coordinates are in different frames");
    }
 
    // Get the appropriate Y-Coordinates
    double y1 = coords1.YCoordinate_Rad();
    double y2 = coords2.YCoordinate_Rad();
    if (coords1.GetCoordSystem() == CECoordinateType::OBSERVED) {
        y1 = M_PI_2 - y1;
        y2 = M_PI_2 - y2;
    }
 
    // Convert the second coordinates to be the same type as the first set of coordinates
    CEAngle angsep = AngularSeparation(coords1.XCoord(), y1,
                                       coords2.XCoord(), y2) ;
    
    return angsep ;
}
 
/**********************************************************************/
CEAngle CECoordinates::AngularSeparation(const CEAngle& xcoord_first, 
                                         const CEAngle& ycoord_first,
                                         const CEAngle& xcoord_second, 
                                         const CEAngle& ycoord_second)
{
    // Call the 'iauSeps' SOFA algorithm
    CEAngle angsep(iauSeps(xcoord_first, ycoord_first,
                           xcoord_second, ycoord_second)) ;
    
    return angsep ;
}
 
/**********************************************************************/
CECoordinates CECoordinates::ConvertTo(CECoordinateType output_coord_type,
                                       const CEObserver& observer, 
                                       const CEDate& date)
{
    CEDate tmpdate(date);
    return ConvertTo(output_coord_type,
                     tmpdate,
                     observer.Longitude_Rad(),
                     observer.Latitude_Rad(),
                     observer.Elevation_m(),
                     observer.Pressure_hPa(),
                     observer.Temperature_C(),
                     observer.RelativeHumidity(),
                     tmpdate.dut1(),
                     tmpdate.xpolar(),
                     tmpdate.ypolar(),
                     observer.Wavelength_um()) ;
}
 
 
/**********************************************************************/
CECoordinates CECoordinates::ConvertTo(CECoordinateType output_coord_type,
                                       double jd,
                                       double longitude,
                                       double latitude,
                                       double elevation_m,
                                       double pressure_hPa,
                                       double temperature_celsius,
                                       double relative_humidity,
                                       double dut1,
                                       double xp, double yp,
                                       double wavelength_um)
{
    // Create return coordinate object
    CECoordinates out_coord;
 
    // Do conversion to CIRS
    if (output_coord_type == CECoordinateType::CIRS) {
        out_coord = ConvertToCIRS(jd, longitude, latitude,
                                  elevation_m, pressure_hPa,
                                  temperature_celsius, relative_humidity,
                                  dut1, xp, yp, wavelength_um) ;
    } else if (output_coord_type == CECoordinateType::ICRS) {
        out_coord = ConvertToICRS(jd, longitude, latitude,
                                  elevation_m, pressure_hPa,
                                  temperature_celsius, relative_humidity,
                                  dut1, xp, yp, wavelength_um) ;
    } else if (output_coord_type == CECoordinateType::GALACTIC) {
        out_coord = ConvertToGalactic(jd, longitude, latitude,
                                      elevation_m, pressure_hPa,
                                      temperature_celsius, relative_humidity,
                                      dut1, xp, yp, wavelength_um) ;
    } else if (output_coord_type == CECoordinateType::OBSERVED) {
        out_coord = ConvertToObserved(jd, longitude, latitude,
                                      elevation_m, pressure_hPa,
                                      temperature_celsius, relative_humidity,
                                      dut1, xp, yp, wavelength_um) ;
    }
 
    return out_coord;
}
 
/**********************************************************************/
CECoordinates CECoordinates::ConvertToCIRS(double jd,
                            double longitude,
                            double latitude,
                            double elevation_m,
                            double pressure_hPa,
                            double temperature_celsius,
                            double relative_humidity,
                            double dut1,
                            double xp, double yp,
                            double wavelength_um)
{
    double xcoord_new(0.0);
    double ycoord_new(0.0) ;
    if (coord_type_ == CECoordinateType::CIRS) {
        // CIRS -> CIRS
        xcoord_new = XCoordinate_Rad() ;
        ycoord_new = YCoordinate_Rad() ;
    } else if (coord_type_ == CECoordinateType::ICRS) {
        // ICRS -> CIRS
        ICRS2CIRS(XCoordinate_Rad(), YCoordinate_Rad(),
                  &xcoord_new, &ycoord_new,
                  jd, CEAngleType::RADIANS) ;
    } else if (coord_type_ == CECoordinateType::GALACTIC) {
        // Galactic -> CIRS
        Galactic2CIRS(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new,
                      jd, CEAngleType::RADIANS) ;
    } else if (coord_type_ == CECoordinateType::OBSERVED) {
        // Observed -> CIRS
        Observed2CIRS(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new,
                      jd, longitude, latitude,
                      elevation_m, pressure_hPa,
                      temperature_celsius, relative_humidity,
                      dut1, xp, yp, wavelength_um) ;
    }
    
    return CECoordinates(xcoord_new, 
                         ycoord_new,
                         CECoordinateType::CIRS) ;
}
 
/**********************************************************************/
CECoordinates CECoordinates::ConvertToICRS(double jd,
                            double longitude,
                            double latitude,
                            double elevation_m,
                            double pressure_hPa,
                            double temperature_celsius,
                            double relative_humidity,
                            double dut1,
                            double xp, double yp,
                            double wavelength_um)
{
    double xcoord_new(0.0);
    double ycoord_new(0.0);
    if (coord_type_ == CECoordinateType::CIRS) {
        CIRS2ICRS(XCoordinate_Rad(), YCoordinate_Rad(),
                  &xcoord_new, &ycoord_new, jd, CEAngleType::RADIANS) ;
    } else if (coord_type_ == CECoordinateType::ICRS) {
        xcoord_new = XCoordinate_Rad() ;
        ycoord_new = YCoordinate_Rad() ;
    } else if (coord_type_ == CECoordinateType::GALACTIC) {
        Galactic2ICRS(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new, CEAngleType::RADIANS) ;
    } else if (coord_type_ == CECoordinateType::OBSERVED) {
        Observed2ICRS(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new,
                      jd, longitude, latitude,
                      elevation_m, pressure_hPa,
                      temperature_celsius, relative_humidity,
                      dut1, xp, yp, wavelength_um) ;
    }
    
    return CECoordinates(xcoord_new, ycoord_new,
                         CECoordinateType::ICRS) ;
}
 
/**********************************************************************/
CECoordinates CECoordinates::ConvertToGalactic(double jd,
                                double longitude,
                                double latitude,
                                double elevation_m,
                                double pressure_hPa,
                                double temperature_celsius,
                                double relative_humidity,
                                double dut1,
                                double xp, double yp,
                                double wavelength_um)
{
    double xcoord_new(0.0);
    double ycoord_new(0.0);
    if (coord_type_ == CECoordinateType::CIRS) {
        CIRS2Galactic(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new, jd, CEAngleType::RADIANS) ;
    } else if (coord_type_ == CECoordinateType::ICRS) {
        ICRS2Galactic(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new, CEAngleType::RADIANS) ;
    } else if (coord_type_ == CECoordinateType::GALACTIC) {
        xcoord_new = XCoordinate_Rad() ;
        ycoord_new = YCoordinate_Rad() ;
    } else if (coord_type_ == CECoordinateType::OBSERVED) {
        Observed2Galactic(XCoordinate_Rad(), YCoordinate_Rad(),
                          &xcoord_new, &ycoord_new,
                          jd, longitude, latitude,
                          elevation_m, pressure_hPa,
                          temperature_celsius, relative_humidity,
                          dut1, xp, yp, wavelength_um) ;
    }
    
    return CECoordinates(xcoord_new, ycoord_new, CECoordinateType::GALACTIC) ;
}
 
/**********************************************************************/
CECoordinates CECoordinates::ConvertToObserved(double jd,
                                double longitude,
                                double latitude,
                                double elevation_m,
                                double pressure_hPa,
                                double temperature_celsius,
                                double relative_humidity,
                                double dut1,
                                double xp, double yp,
                                double wavelength_um)
{
    double xcoord_new(0.0);
    double ycoord_new(0.0);
    double apparent_x, apparent_y, apparent_hourangle;
    if (coord_type_ == CECoordinateType::CIRS) {
        CIRS2Observed(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new,
                      jd,
                      longitude, latitude,
                      elevation_m, pressure_hPa,
                      temperature_celsius,
                      relative_humidity,
                      dut1, xp, yp,
                      wavelength_um,
                      &apparent_x, &apparent_y, &apparent_hourangle) ;
    } else if (coord_type_ == CECoordinateType::ICRS) {
        ICRS2Observed(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new,
                      jd,
                      longitude, latitude,
                      elevation_m, pressure_hPa,
                      temperature_celsius,
                      relative_humidity,
                      dut1, xp, yp,
                      wavelength_um,
                      &apparent_x, &apparent_y, &apparent_hourangle) ;
    } else if (coord_type_ == CECoordinateType::GALACTIC) {
        Galactic2Observed(XCoordinate_Rad(), YCoordinate_Rad(),
                      &xcoord_new, &ycoord_new,
                      jd,
                      longitude, latitude,
                      elevation_m, pressure_hPa,
                      temperature_celsius,
                      relative_humidity,
                      dut1, xp, yp,
                      wavelength_um,
                      &apparent_x, &apparent_y, &apparent_hourangle) ;
    } else {
        // Throw an exception
        std::string msg = "[ERROR] ";
        if (coord_type_ == CECoordinateType::OBSERVED) {
            // Cant convert from observed to observed without additional
            // observer information
            msg += "Unable to convert to OBSERVED from OBSERVED";
        }
        throw CEException::invalid_value("CECoordinates::ConvertToObserved(long form)",
                                         msg);
    }
 
    return CECoordinates(xcoord_new, ycoord_new, CECoordinateType::OBSERVED) ;
}
 
/**********************************************************************/
std::vector<double> CECoordinates::GetDMS(const double& angle,
                                          const CEAngleType& angle_type)
{
    CEAngle ang;
    ang.SetAngle(angle, angle_type);
    return ang.DmsVect() ;
}
 
/**********************************************************************/
std::vector<double> CECoordinates::GetHMS(const double& angle,
                                          const CEAngleType& angle_type)
{
    // Convert to degrees if passed radians
    CEAngle ang;
    ang.SetAngle(angle, angle_type);
    return ang.HmsVect();
}
 
 
/**********************************************************************/
double CECoordinates::HMSToAngle(const std::vector<double>& angle,
                                 const CEAngleType& return_type)
{
    CEAngle ret_angle = CEAngle::Hms(angle);
 
    if (return_type == CEAngleType::DEGREES) {
        return ret_angle.Deg();
    } else {
        return ret_angle.Rad();
    }
}
 
 
/**********************************************************************/
double CECoordinates::DMSToAngle(const std::vector<double>& angle,
                                 const CEAngleType& return_type)
{
    CEAngle ret_angle = CEAngle::Dms(angle);
    
    if (return_type == CEAngleType::DEGREES) {
        return ret_angle.Deg();
    } else {
        return ret_angle.Rad();
    }
}
 
 
/**********************************************************************/
void CECoordinates::SetCoordinates(const CEAngle& xcoord, 
                                   const CEAngle& ycoord,
                                   const CECoordinateType& coord_type)
{
    xcoord_     = xcoord ;
    ycoord_     = ycoord ;
    coord_type_ = coord_type ;
}
 
/**********************************************************************/
void CECoordinates::SetCoordinates(const CECoordinates& coords)
{
    copy_members(coords);
}
 
 
/**********************************************************************/
std::string CECoordinates::print(void) const
{
    std::string msg = "Coordinates:\n";
    msg += "   - System : " + std::to_string(int(coord_type_)) + "\n";
    msg += "   - X-coord: " + std::to_string(XCoordinate_Deg()) + " deg\n";
    msg += "   - Y-coord: " + std::to_string(YCoordinate_Deg()) + " deg\n";
    return msg;
}
 
 
/**********************************************************************/
void CECoordinates::copy_members(const CECoordinates& other)
{
    coord_type_ = other.coord_type_;
    xcoord_     = other.xcoord_;
    ycoord_     = other.ycoord_;
}
 
 
/**********************************************************************/
void CECoordinates::free_members(void)
{
    // Nothing to do here
}
 
 
/**********************************************************************/
void CECoordinates::init_members(void)
{
    coord_type_ = CECoordinateType::ICRS;
    xcoord_     = 0.0;
    ycoord_     = 0.0;
}
 
 
/**********************************************************************/
bool operator==(const CECoordinates& lhs, const CECoordinates& rhs)
{
    bool are_equal = true;
    
    // Check the coordinate systems are equivalent
    if (lhs.GetCoordSystem() != rhs.GetCoordSystem()) {
        are_equal = false;
    }
    // Check that the x-coordinate and the y-coordinate are the same
    else {
        // Check how far appart the coordinates are from each other
        CEAngle angsep = CECoordinates::AngularSeparation(lhs, rhs);
        // Currently require separation < 0.03 arcsec
        double marcsec_rad = 4.848e-6;
        if (angsep > 3.0*marcsec_rad) {
            are_equal = false;
        }
    }
 
    return are_equal;
}
 
 
/**********************************************************************/
bool operator!=(const CECoordinates& lhs, const CECoordinates& rhs)
{
    return !(lhs == rhs);
}