CassiniSoldner.hpp 9.39 KB
Newer Older
Valentin Platzgummer's avatar
Valentin Platzgummer committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
/**
 * \file CassiniSoldner.hpp
 * \brief Header for GeographicLib::CassiniSoldner class
 *
 * Copyright (c) Charles Karney (2009-2019) <charles@karney.com> and licensed
 * under the MIT/X11 License.  For more information, see
 * https://geographiclib.sourceforge.io/
 **********************************************************************/

#if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP)
#define GEOGRAPHICLIB_CASSINISOLDNER_HPP 1

#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicLine.hpp>
#include <GeographicLib/Constants.hpp>

namespace GeographicLib {

  /**
   * \brief Cassini-Soldner projection
   *
   * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e
   * lon0, on the ellipsoid.  This projection is a transverse cylindrical
   * equidistant projection.  The projection from (\e lat, \e lon) to easting
   * and northing (\e x, \e y) is defined by geodesics as follows.  Go north
   * along a geodesic a distance \e y from the central point; then turn
   * clockwise 90&deg; and go a distance \e x along a geodesic.
   * (Although the initial heading is north, this changes to south if the pole
   * is crossed.)  This procedure uniquely defines the reverse projection.  The
   * forward projection is constructed as follows.  Find the point (\e lat1, \e
   * lon1) on the meridian closest to (\e lat, \e lon).  Here we consider the
   * full meridian so that \e lon1 may be either \e lon0 or \e lon0 +
   * 180&deg;.  \e x is the geodesic distance from (\e lat1, \e lon1) to
   * (\e lat, \e lon), appropriately signed according to which side of the
   * central meridian (\e lat, \e lon) lies.  \e y is the shortest distance
   * along the meridian from (\e lat0, \e lon0) to (\e lat1, \e lon1), again,
   * appropriately signed according to the initial heading.  [Note that, in the
   * case of prolate ellipsoids, the shortest meridional path from (\e lat0, \e
   * lon0) to (\e lat1, \e lon1) may not be the shortest path.]  This procedure
   * uniquely defines the forward projection except for a small class of points
   * for which there may be two equally short routes for either leg of the
   * path.
   *
   * Because of the properties of geodesics, the (\e x, \e y) grid is
   * orthogonal.  The scale in the easting direction is unity.  The scale, \e
   * k, in the northing direction is unity on the central meridian and
   * increases away from the central meridian.  The projection routines return
   * \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the
   * reciprocal of the scale in the northing direction.
   *
   * The conversions all take place using a Geodesic object (by default
   * Geodesic::WGS84()).  For more information on geodesics see \ref geodesic.
   * The determination of (\e lat1, \e lon1) in the forward projection is by
   * solving the inverse geodesic problem for (\e lat, \e lon) and its twin
   * obtained by reflection in the meridional plane.  The scale is found by
   * determining where two neighboring geodesics intersecting the central
   * meridian at \e lat1 and \e lat1 + \e dlat1 intersect and taking the ratio
   * of the reduced lengths for the two geodesics between that point and,
   * respectively, (\e lat1, \e lon1) and (\e lat, \e lon).
   *
   * Example of use:
   * \include example-CassiniSoldner.cpp
   *
   * <a href="GeodesicProj.1.html">GeodesicProj</a> is a command-line utility
   * providing access to the functionality of AzimuthalEquidistant, Gnomonic,
   * and CassiniSoldner.
   **********************************************************************/

  class GEOGRAPHICLIB_EXPORT CassiniSoldner {
  private:
    typedef Math::real real;
    Geodesic _earth;
    GeodesicLine _meridian;
    real _sbet0, _cbet0;
    static const unsigned maxit_ = 10;

  public:

    /**
     * Constructor for CassiniSoldner.
     *
     * @param[in] earth the Geodesic object to use for geodesic calculations.
     *   By default this uses the WGS84 ellipsoid.
     *
     * This constructor makes an "uninitialized" object.  Call Reset to set the
     * central latitude and longitude, prior to calling Forward and Reverse.
     **********************************************************************/
    explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84());

    /**
     * Constructor for CassiniSoldner specifying a center point.
     *
     * @param[in] lat0 latitude of center point of projection (degrees).
     * @param[in] lon0 longitude of center point of projection (degrees).
     * @param[in] earth the Geodesic object to use for geodesic calculations.
     *   By default this uses the WGS84 ellipsoid.
     *
     * \e lat0 should be in the range [&minus;90&deg;, 90&deg;].
     **********************************************************************/
    CassiniSoldner(real lat0, real lon0,
                   const Geodesic& earth = Geodesic::WGS84());

    /**
     * Set the central point of the projection
     *
     * @param[in] lat0 latitude of center point of projection (degrees).
     * @param[in] lon0 longitude of center point of projection (degrees).
     *
     * \e lat0 should be in the range [&minus;90&deg;, 90&deg;].
     **********************************************************************/
    void Reset(real lat0, real lon0);

    /**
     * Forward projection, from geographic to Cassini-Soldner.
     *
     * @param[in] lat latitude of point (degrees).
     * @param[in] lon longitude of point (degrees).
     * @param[out] x easting of point (meters).
     * @param[out] y northing of point (meters).
     * @param[out] azi azimuth of easting direction at point (degrees).
     * @param[out] rk reciprocal of azimuthal northing scale at point.
     *
     * \e lat should be in the range [&minus;90&deg;, 90&deg;].  A call to
     * Forward followed by a call to Reverse will return the original (\e lat,
     * \e lon) (to within roundoff).  The routine does nothing if the origin
     * has not been set.
     **********************************************************************/
    void Forward(real lat, real lon,
                 real& x, real& y, real& azi, real& rk) const;

    /**
     * Reverse projection, from Cassini-Soldner to geographic.
     *
     * @param[in] x easting of point (meters).
     * @param[in] y northing of point (meters).
     * @param[out] lat latitude of point (degrees).
     * @param[out] lon longitude of point (degrees).
     * @param[out] azi azimuth of easting direction at point (degrees).
     * @param[out] rk reciprocal of azimuthal northing scale at point.
     *
     * A call to Reverse followed by a call to Forward will return the original
     * (\e x, \e y) (to within roundoff), provided that \e x and \e y are
     * sufficiently small not to "wrap around" the earth.  The routine does
     * nothing if the origin has not been set.
     **********************************************************************/
    void Reverse(real x, real y,
                 real& lat, real& lon, real& azi, real& rk) const;

    /**
     * CassiniSoldner::Forward without returning the azimuth and scale.
     **********************************************************************/
    void Forward(real lat, real lon,
                 real& x, real& y) const {
      real azi, rk;
      Forward(lat, lon, x, y, azi, rk);
    }

    /**
     * CassiniSoldner::Reverse without returning the azimuth and scale.
     **********************************************************************/
    void Reverse(real x, real y,
                 real& lat, real& lon) const {
      real azi, rk;
      Reverse(x, y, lat, lon, azi, rk);
    }

    /** \name Inspector functions
     **********************************************************************/
    ///@{
    /**
     * @return true if the object has been initialized.
     **********************************************************************/
    bool Init() const { return _meridian.Init(); }

    /**
     * @return \e lat0 the latitude of origin (degrees).
     **********************************************************************/
    Math::real LatitudeOrigin() const
    { return _meridian.Latitude(); }

    /**
     * @return \e lon0 the longitude of origin (degrees).
     **********************************************************************/
    Math::real LongitudeOrigin() const
    { return _meridian.Longitude(); }

    /**
     * @return \e a the equatorial radius of the ellipsoid (meters).  This is
     *   the value inherited from the Geodesic object used in the constructor.
     **********************************************************************/
    Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); }

    /**
     * @return \e f the flattening of the ellipsoid.  This is the value
     *   inherited from the Geodesic object used in the constructor.
     **********************************************************************/
    Math::real Flattening() const { return _earth.Flattening(); }

    /**
      * \deprecated An old name for EquatorialRadius().
      **********************************************************************/
    // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
    Math::real MajorRadius() const { return EquatorialRadius(); }
    ///@}

  };

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_CASSINISOLDNER_HPP