EarthShape.h

Go to the documentation of this file.

//- ****************************************************************************
//-
//- Copyright 2009 Sandia Corporation. Under the terms of Contract
//- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
//- retains certain rights in this software.
//-
//- BSD Open Source License.
//- All rights reserved.
//-
//- Redistribution and use in source and binary forms, with or without
//- modification, are permitted provided that the following conditions are met:
//-
//-    * Redistributions of source code must retain the above copyright notice,
//-      this list of conditions and the following disclaimer.
//-    * Redistributions in binary form must reproduce the above copyright
//-      notice, this list of conditions and the following disclaimer in the
//-      documentation and/or other materials provided with the distribution.
//-    * Neither the name of Sandia National Laboratories nor the names of its
//-      contributors may be used to endorse or promote products derived from
//-      this software without specific prior written permission.
//-
//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//- POSSIBILITY OF SUCH DAMAGE.
//-
//- ****************************************************************************
 
#ifndef EARTH_SHAPE_H
#define EARTH_SHAPE_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <iostream>
#include <string>
 
// use standard library objects
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "CPPUtils.h"
#include "GeoTessException.h"
 
// **** _BEGIN GEOTESS NAMESPACE_ **********************************************
 
namespace geotess
{
 
// **** _FORWARD REFERENCES_ ***************************************************
 
// **** _CLASS DEFINITION_ *****************************************************
 
class GEOTESS_EXP_IMP EarthShape
{
private:
 
    string shapeName;
 
    double equatorialRadius;
 
    double inverseFlattening;
 
    double eccentricitySqr;
 
    double e1;
 
    double e2;
 
    bool constantRadius;
 
    bool sphere;
 
public:
 
    EarthShape(const string& earthShape = "WGS84") { setEarthShape(earthShape); }
 
    void setEarthShape(const string& earthShape)
    {
        shapeName = earthShape;
 
        if (earthShape == "SPHERE")
        {
            equatorialRadius = 6371;
            inverseFlattening = 1e99;
            constantRadius = true;
        }
        else if (earthShape == "GRS80")
        {
            equatorialRadius = 6378.137;
            inverseFlattening = 298.257222101;
            constantRadius = false;
        }
        else if (earthShape == "GRS80_RCONST")
        {
            equatorialRadius = 6371.;
            inverseFlattening = 298.257222101;
            constantRadius = true;
        }
        else if (earthShape == "WGS84")
        {
            equatorialRadius = 6378.137;;
            inverseFlattening = 298.257223563;
            constantRadius = false;
        }
        else if (earthShape == "WGS84_RCONST")
        {
            equatorialRadius = 6371.;
            inverseFlattening = 298.257223563;
            constantRadius = true;
        }
        else if (earthShape == "IERS2003")
        {
            equatorialRadius = 6378.1366;
            inverseFlattening = 298.25642;
            constantRadius = false;
        }
        else if (earthShape == "IERS2003_RCONST")
        {
            equatorialRadius = 6371.;
            inverseFlattening = 298.25642;
            constantRadius = true;
        }
        else
        {
            ostringstream os;
            os << endl << "ERROR in EarthShape::setEarthShape" << endl
                    << earthShape << " is not a recognized EarthShape" << endl
                    << "Valid EarthShapes include SPHERE, GRS80, GRS80_RCONST, WGS84, WGS84_RCONST, IERS2003 and IERS2003_RCONST" << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 9001);
        }
 
        sphere = shapeName == "SPHERE";
        eccentricitySqr = (2.- 1./inverseFlattening)/inverseFlattening;
        e1 = 1.-eccentricitySqr;
        e2 = eccentricitySqr/(1.-eccentricitySqr);
    }
 
    virtual ~EarthShape()
    {
    }
 
    static string class_name()
    { return "EarthShape"; }
 
    virtual int class_size() const
    { return (int) sizeof(EarthShape); }
 
    double getEarthRadius(const double* const v)
    { return constantRadius ? equatorialRadius : equatorialRadius / sqrt(1 + e2 * v[2] * v[2]); }
 
    double getGeocentricLat(const double& lat)
    { return sphere ? lat : atan(tan(lat) * e1); }
 
    double getGeographicLat(const double& lat)
    { return sphere ? lat : atan(tan(lat) / e1); }
 
    double getGeocentricLatDegrees(const double& lat)
    { return sphere ? lat : CPPUtils::toDegrees(atan(tan(CPPUtils::toRadians(lat)) * e1)); }
 
    double getGeographicLatDegrees(const double& lat)
    { return sphere ? lat : CPPUtils::toDegrees(atan(tan(CPPUtils::toRadians(lat)) / e1)); }
 
    double getLat(const double* const v)
    { return getGeographicLat(asin(v[2])); }
 
    double getLatDegrees(const double* const v)
    { return CPPUtils::toDegrees(getLat(v)); }
 
    double getLon(const double* const v)
    { return atan2(v[1], v[0]); }
 
    double getLonDegrees(const double* const v)
    { return CPPUtils::toDegrees(atan2(v[1], v[0])); }
 
    void getVectorDegrees(const double& lat, const double& lon, double* v)
    { getVector(CPPUtils::toRadians(lat), CPPUtils::toRadians(lon), v); }
 
    double* getVectorDegrees(const double& lat, const double& lon)
    { return getVector(CPPUtils::toRadians(lat), CPPUtils::toRadians(lon)); }
 
    double* getVector(const double& lat, const double& lon)
    {
        double* v = new double[3];
        getVector(lat, lon, v);
        return v;
    }
 
    double* getVector(const double& lat, const double& lon, double* v)
    {
        // convert lat from geographic to geocentric latitude.
        double temp = getGeocentricLat(lat);
 
        // z component of v is sin of geocentric latitude.
        v[2] = sin(temp);
 
        // set lat = to cos of geocentric latitude
        temp = cos(temp);
 
        // compute x and y components of v
        v[0] = temp * cos(lon);
        v[1] = temp * sin(lon);
 
        return v;
    }
 
    string getLatLonString(const double* const v)
    {
      char s[300];
      string frmt("%9.5f %10.5f");
      sprintf(s, frmt.c_str(), getLatDegrees(v), getLonDegrees(v));
      return s;
    }
 
    bool isConstantRadius() const { return constantRadius; }
 
    double getEccentricitySqr() const { return eccentricitySqr; }
 
    double getEquatorialRadius() const { return equatorialRadius; }
 
    double getInverseFlattening() const { return inverseFlattening; }
 
    bool isSphere() const { return sphere; }
 
    const string& getShapeName() const { return shapeName; }
 
};
// end class EarthShape
 
} // end namespace geotess
 
#endif  // EARTH_SHAPE_H