SpatialUtils.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-

/* 
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 * 
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * 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 LSST License Statement and 
 * the GNU General Public License along with this program.  If not, 
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#include "lsst/ap/utils/SpatialUtils.h"

#include <cfloat>

#include "lsst/pex/exceptions.h"
#include "lsst/ap/constants.h"

using lsst::pex::exceptions::InvalidParameterError;
using lsst::pex::exceptions::RuntimeError;

using lsst::afw::geom::HALFPI;
using lsst::afw::geom::PI;
using lsst::afw::geom::TWOPI;
using lsst::afw::geom::Angle;
using lsst::afw::geom::radians;
using lsst::afw::coord::IcrsCoord;


namespace lsst { namespace ap { namespace utils {

void thetaRangeReduce(Angle &min, Angle &max) {
    if (min > max) {
        if (max < 0.0 || min >= TWOPI) {
            throw LSST_EXCEPT(InvalidParameterError,
                              "Invalid longitude angle interval");
        }
    } else if (max - min >= TWOPI) {
        min = 0.0 * radians;
        max = TWOPI * radians;
    } else if (min < 0.0 || max >= TWOPI) {
        // range reduce
        min.wrap();
        max.wrap();
    }
}

Angle const maxAlpha(
    Angle radius,   
    Angle centerPhi 
) {
    static const double POLE_EPSILON = 1e-7;

    if (radius < 0.0 || radius > HALFPI) {
        throw LSST_EXCEPT(InvalidParameterError,
                          "radius must be in range [0, PI/2] deg");
    }
    if (radius == 0.0) {
        return 0.0 * radians;
    }
    centerPhi = clampPhi(centerPhi);
    if (std::fabs(centerPhi.asRadians()) + radius.asRadians() > HALFPI - POLE_EPSILON) {
        return PI*(1 + 2.0*DBL_EPSILON) * radians;
    }
    double y = std::sin(radius.asRadians());
    double x = std::sqrt(std::fabs(std::cos(centerPhi.asRadians() - radius.asRadians()) *
                                   std::cos(centerPhi.asRadians() + radius.asRadians())));
    return std::fabs(std::atan(y / x)) * radians;
}

void positionAndVelocity(
    Eigen::Vector3d &p, 
    Eigen::Vector3d &v, 
    Angle ra,           
    Angle decl,         
    double muRa,        
    double muDecl,      
    double vRadial,     
    Angle parallax      
) {
    // distance (AU)
    double r = 1.0 / parallax.asRadians();

    // convert to position and velocity vector (AU, AU/day)
    double sinRa = sin(ra.asRadians());
    double cosRa = cos(ra.asRadians());
    double sinDecl = sin(decl.asRadians());
    double cosDecl = cos(decl.asRadians());
    double s = r*cosDecl;
    double t = r*muDecl*sinDecl;
    double u = cosDecl*vRadial;
    p = Eigen::Vector3d(s*cosRa, s*sinRa, r*sinDecl);
    v = Eigen::Vector3d(cosRa*u - p.y()*muRa - cosRa*t,
                        sinRa*u + p.x()*muRa - sinRa*t,
                        sinDecl*vRadial + s*muDecl);
    if (v.squaredNorm() > 0.25*C_AU_PER_DAY*C_AU_PER_DAY) {
        throw LSST_EXCEPT(RuntimeError,
                          "star velocity vector magnitude exceeds half "
                          "the speed of light");
    }
}

Eigen::Vector2d const cartesianToSpherical(Eigen::Vector3d const &v) {
    double d2 = v.x()*v.x() + v.y()*v.y();
    double theta = (d2 == 0.0) ? 0.0 : std::atan2(v.y(), v.x());
    double phi = (v.z() == 0.0) ? 0.0 : std::atan2(v.z(), std::sqrt(d2));
    return Eigen::Vector2d(theta, phi);
}

}}} // namespace lsst::ap::utils