Functions setting up the ellipsoid parameters, performing the geometry
   computations on its surface and mapping a point from ellipsoid to
   NCS ("near-conformal-sphere") and vice versa.

Ellipsoid of rotation (x²·b²+y²·b²+z²·a²=a²·b²) meridian section

Description:
   A short list of important geometry elements of the meridian ellipse,
   and of an arbitrary point on the ellipsoid surface:
   ·
   P:         Point on planetary surface
   O:         x/O/z coordinate system origin
   F, F'      Ellipse focal points
   ·
   F↔B = O↔A: Ellipse semi major axis
   OB:        Ellipse semi minor axis
   ·
   n:         Ellipse normal in P
   t:         Ellipse tangent in P
   N:         Normal "toe point"
   ·
   P↔T:      Tangent free term
   O↔T:      Normal free term
   T:        Paracenter of P
   ·
   Ellipsoid coordinates are normally encoded as angular φ snd λ measured
   in radians, while NCS coordinates are encoded as an i, j, k vector,
   i.e., the direction cosines of radius-vector with X, Y and Z axis of
   a unit sphere.
   .
   Unlike on the sphere, the normal to the tangential plane of the ellipsoid
   surface in some point in general case does not pass through the coordinate
   origin; it does so only for two Poles and any point on the Equator. The
   point on the normal closest to the coordinate origin is in the documentation
   of functions of this Library called "paracenter".
   ·
   The longitude of NCS coordinates is identical to ellipsoid longitude,
   while the latitude is numerically very close to the spherical latitude
   that would be obtained by rigorous ellipsoid to sphere conformal mapping.
   The transformation between the spherical and ellipsoidal coordinate
   domains is much faster for NCS than it would be for rigorously computed
   ellipsoid-to-sphere (and inverse) conformal mapping.
   ·
   For additional discussion of the Near-Conformal Sphere, see:
   www.lukatela.com/hrvoje/papers/ncsphere.html

Synopsis:
   #include <nemo.h>
   void nemo_ElrInit(double a,
                     double b,
                     nemoElRot *pElr);
Description:
   Given ellipsoid semi-major and semi-minor axes, populate a structure
   that contains a short set of numerical parameters frequently used
   in evaluation of spatial relationships on the surface of ellipsoid.
Arguments:
   a, b:
      Doubles, given semi-major and semi-minor ellipsoid axes.
   pElr:
      Pointer to a receiving structure of nemoElRot type, returned
      ellipsoid parameters. Once computed by this function, the
      parameters should not be changed by the application code.
See Also:
   nemo_ElrWgs84()

Synopsis:
   #include <nemo.h>
   nemoElRot *nemo_ElrWgs84();
Description:
   Applications that model the spatial geometry on the planet Earth,
   currenrly almost universaly represented by WGS84 ellipsoid, can use
   this function instead of instantiating the ellipsoid using semi-exes
   and nemo_ElrInit(() to populate the ellipsoid parameters structure.
Return Value:
   Pointer to an internal static structure of nemoElRot type, parameters
   representing WGS84 ellipsoid.
See Also:
   nemo_ElrInit()

Synopsis:
   #include <nemo.h>
   double nemo_GeodesicVincenty(const nemoElRot *pElr,
                                const nemoPtEnr *ptA,
                                const nemoPtEnr *ptB,
                                int *iterCount);
Description:
   Given two points on the surface of an ellipsoid of rotation, compute the
   length of geodesic: the shortest line (on the surface of the ellipsoid)
   between them. The geodesic has the following important property:
   ·
   Each point on geodesic and the two very close points which are also on
   geodesic, one on each side of it, define a plane which contains the normal
   to the ellipddoid surface in the mid-point ("osculating plane"). However,
   each line on ellipsoid that has this property is not necessarily the
   shortest line between its endpoints, thus, in the strict sense, not a
   geodesic. The example of such line is the Equator: is not the shortest
   distance for any two points on it.
   ·
   Computing the length of geodesic is difficult. Even more difficult is
   finding the coordinates of a point which divides a given geodesic in an
   arbitrary ratio, finding the distance between an arbtrary point on the
   ellipsoid surface and a given geodesic or finding the point of intersection
   between two geodesics. No closed formulae definding the geometry of any of
   those productions exists; thus all solutiona are either (to some degree)
   approximate, or iterative and dependent on some pre-defined "solution
   proximity" criterion.
   ·
   The geodesic length computation method used in this function is that
   proposed by Tadeusz Szpila (later in life known as Thaddeus Vincenty),
   published in "Survey Review" No. 176, April 1975.
   ·
   Vincenty's method is considered slow but very accurate. It is commonly
   used, together with that of Rudoe's (Bomford, Geodyssey) and Karney's
   (Journal of Geodesy, 87/2013), for verifying the results of ellipsoid
   geometry propositions obtained by various approximate, but faster
   computations or methods.
   ·
   The procedure is iterative; the function will return an undefined
   signal value if the solution fails to converge.
   ·
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation
      parameters, as initialized by nemo_InitEllipsoid() or as returned
      by nemo_ElrWgs84().
   ptA:
      Pointer to nemoPtEnr structure, first given point ellipsoidal
      vector coordinates (normal i,j,k).
   ptB:
      Pointer to nemoPtEnr structure, second given point coordinates.
   iterCount:
      Pointer to a natural integer variable, returned number of iterations
      required to reach the solution, ignored if NULL.
Return Value:
   Distance along the geodesic between ptA and ptB, measured in
   the same units as are used for first two pElr array elements.
   (NEMO_DOUBLE_UNDEF if the solution (unexpectedly) failed to
   converge in VINCENTY_ITERATIONS_LIMIT number of steps).

Synopsis:
   #include <nemo.h>
   double nemo_MeridianArcLength(const nemoElRot *pElr,
                                 const nemoPtEnr *pLat);
Description:
   Find the length of ellipsdoid arc from the Equator to a given latitude.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pLat:
      Pointer to nemoPtEnr structure, given ellipsoid coordinates as the unit
      vector of its ellipsoid normal. The meridian arc the length of which is
      required starts at the Equator (φ = 0) and ends at the latitude of the
      given point. (The longitude of the point is ignored; but note that the
      point must be a normalized vector.)
Return Value:
   Distance of the meridian arc in the same units as pElr.
Example:
   ...
   /* Hello, Round World: A historical curio - or -
      GPS vs. chevalier Delambre et citoyen Méchain - or -
      what was the length of the meridian that was meant to define
      The Metre, measured at length from 1792 to 1799 ?
      The two landmarks coordinates - with no accuracy pretensions - may be:
      Beffroi Saint-Éloi de Dunkerque,         51.035625886, 2.376145883;
      Torre de guaita del Castell de Montjuïc, 41.363566497, 2.166566974;
      N.B.: It is not the distance between two geodetic stations we want to
      know, it is the distance between the two points with given latitudes,
      assuming they are located due north/south from each other.
    */
   ...
   double dunkerque[2] = {NEMO_DEG2RAD * 51.035625886,
                          NEMO_DEG2RAD * 2.376145883};
   double montjuic[2] =  {NEMO_DEG2RAD * 41.363566497,
                          NEMO_DEG2RAD * 2.166566974};
   nemoPtEnr ptEnrDunk, ptEnrMont;
   double eqDunk, eqMont;
   ...
   nemo_LatLongToDcos3(dunkerque, ptEnrDunk.dc);
   nemo_LatLongToDcos3(montjuic, ptEnrMont.dc);
   eqDunk = nemo_MeridianArcLength(nemo_ElrWgs84(), &ptEnrDunk);
   eqMont = nemo_MeridianArcLength(nemo_ElrWgs84(), &ptEnrMont);
   printf("Length of Delambre/Méchain meridian: %12.3f meters\n",
          eqDunk - eqMont);
   ...

Synopsis:
   #include <nemo.h>
   void nemo_EllToNcs(const nemoElRot *pElr,
                      const nemoPtEll *ptEll,
                      nemoPtNcs *ptNcs);
Description:
   Perform direct mapping of a point from ellipsoid to NCS.
   ·
   As mentioned above, the longitudes (λ) of ellipsoid, conformal sphere and
   NCS are the same, only the latitudes (φ) are different. On Earth, the
   difference between rigorous conformal sphere and NCS is rather small,
   approximately half an arc second, which is about 15 meters on the surface.
    ·
   The difference in φ between the ellipsoid and either sphere is considerably
   more: approximately 11.5 arc minutes, which is about 21 kilometres on the
   Earth's surface. This number is of consequence when constructing computer
   display (as opposed to "hardcopy") maps of high volume geographical data,
   when using fast NCS↔plane projections, instead of considerably slower
   NCS↔ellipsoid↔plane pair of projections.
   ·
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   ptEll:
      Pointer to nemoPtEll structure, given ellipsoid φ, λ coordinetes.
   ptNcs:
      Pointer to nemoPtNcs structure, returned NCS point coordinates.
See Also:
   nemo_NcsToEll()
   nemo_EnrToNcs()
   nemo_NcsToEnr()

Synopsis:
   #include <nemo.h>
   void nemo_NcsToEll(const nemoElRot *pElr,
                      const nemoPtNcs *ptNcs,
                      nemoPtEll *ptEll);
Description:
   Perform inverse mapping of a point from NCS to ellipsoid. "Direct"
   mapping is the term commonly used for coordinate transformation from
   a mathematically more complex surface (for instance, an ellipsoid) to
   a different, mathematically simpler surface (for instance, a sphere),
   while the "inverse" mapping would be coordinate transformation in the
   opposite direction.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   ptNcs:
      Pointer to nemoPtNcs structure, given NCS point coordinates.
   ptEll:
      Pointer to nemoPtEll structure, returned ellipsoid φ, λ coordinates.
See Also:
   nemo_EllToNcs()
   nemo_EnrToNcs()
   nemo_NcsToEnr()

Synopsis:
   #include <nemo.h>
   double nemo_NcsElrScale(const nemoElRot *pElr,
                            const nemoPtNcs *ptNcs);
Description:
   Return mapping scale factor between the NCS (Near Conformal Sphere)
   and the ellipsoid. Note that the factor includes both the mapping
   transformation local scale and the diffrence of magnitude between two
   surfaces (i.e., the sphere with the radius of 1.0 and the elipsoid
   defined by planetary magnitude semi-axes).
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   ptNcs:
      Pointer to nemoPtNcs structure, given NCS point coordinates.
Return Value:
   Double, the scale factor as defined above. If spherical coordinates
   are ill-defined, the function will return NEMO_DOUBLE_UNDEF.
Example:
   ...
   /* Return geodesic approximation as the NCS great circle arc */
   double geodesicApprox(nemoPtNcs *pPtA, nemoPtNcs *pPtB) {
      nemoPtNcs ptMid;
      nemo_MidV3(pPtA->dc, pPtB->dc, ptMid.dc);
      return(nemo_NcsElrScale(nemo_ElrWgs84(), &ptMid) *
             nemo_ArcV3(pPtA->dc, pPtB->dc));
      }
   ...

Synopsis:
   #include <nemo.h>
   double nemo_NcsElrRadius(const nemoElRot *pElr,
                            const nemoPtNcs *ptNcs);
Description:
   Return mean radius of the ellipsoid for a point defined on NCS (Near
   Conformal Sphere.)
   ·
   Mean radius of curvature for a point on the surface of an ellipsoid of
   rotation is defined as the geometric mean of the radii of the
   intersections of the surface with the prime vertical and the curvature
   of meridian in the same pont.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   ptNcs:
      Pointer to nemoPtNcs structure, given NCS point coordinates.
Return Value:
   Double, the mean radius of ellipsoid in the given point. If spherical
   coordinates are ill-defined, the function will return NEMO_DOUBLE_UNDEF.
See Also:
   nemo_NcsElrScale()

Synopsis:
   #include <nemo.h>
   double nemo_EnrToCc3(const nemoElRot *pElr,
                        const nemoPtEnr *pPtP,
                        double elvn,
                        nemoPtCc3 *pPtR,
                        nemoPtCc3 *pPtT);

Ellipsoid normal - coordinates of a point on the ellipsoid

Description:
   xOz:       Meridian section of the ellipsoid at P, point on its surface.
   n:         Ellipsoid normal to the ellipsoid surface in point P.
   t:         Tangential plane in the point P.
   R:         Point above P (along its normal), at elevation elvn.
   T:         Paracenter of P: point on its normal, closest to the origon O.
   With given point P and the elevation above or below it, find the Cartesian
   ccoridates of point R
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtP:
      Pointer to given ellipsoid coordinates as the unit vector of its
      ellipsoid normal.
   elvn:
      Double, elevation above or below pPtP to the point of which the
      coordinates are requited. Measured in the same units as the ellipsoid
      semi-axes. elvn of 0.0 will return point on the ellipsoid durface.
   pPtR:
      Pointer to a receiving structure of nemoPtCc3 type, returned point
      Cartesian coordinates. Ignored if NULL.
   pPtT:
      Pointer to a receiving structure of nemoPtCc3 type, returned
      Cartesian coordinates of P's paracenter'. Ignored if NULL.
Return Value:
   Double, free term of the tangential plane to ellipsoid in pPtP, the T↔P
   distance on the illustration. If the normal in pPtP is I, J, and K, and
   the free term as returned by this function is N, then the canonical
   equation of tangential plane is: I × x + J × y + K × z ± N = 0
   (if both pPtR and pPtT are NULL, this is the only item the invoker is
   interested in).
See Also:
   nemo_Cc3ToEnr()

Synopsis:
   #include <nemo.h>
   double nemo_Cc3ToEnr(const nemoElRot *pElr,
                        const nemoPtCc3 *pPtR,
                              nemoPtEnr *pPtP,
                              double iterCriter,
                              int *iterCount);
Description:
   Given a position in dimensional Cartesian coordinate system (above or
   below the surface of an ellipsoid of rotation), find the coordinates of
   corresponding point on the surface of the ellipsoid. (In this context.
   "corresponding" point is one below or above, on the same ellipsoid normal
   as the given point).
   .
   (N.B.: geometry objects in the meridian plane are the same as those
   depicted and labelled in the illustration of nemo_EnrToCc3().
   .
   The solution is iterative, but on the ellipsoids with the small amount of
   flattening typical for the planetary bodies, the convergence is very fast.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtR:
      Pointer to a structure of nemoPtCc3 type, given point Cartesian
      coordinates. The point can be above of below the ellipsoid surface,
      but note that a point close to planetary results in a numerically
      weak definition of normal of a surface point.
   pPtP:
      Pointer to a receiving structure of nemoPtCc3 type, returned
      Cartesian coordinates of P's paracenter'.
   iterCriter:
      Double, iteration criterion, linear value in the same units as pElr.
      For instance, if the value of the argument is 0.010, the iteration
      will return the value of ellipsoid normal that did not change more
      than 10 millimers (measured on the planetary surface) from the one
      calculated in the previous iteration.
      With the default limit of iteration steps, the convergence on Earth is
      assured for sub-millimetric values of the criterion, for given points
      at Moho depth (35 km) below the surface, up to geostationary orbit
      altitude (35 000 km).
   iterCount:
      Pointer to an integer variable, returned number of iterations required
      to reach the solution, ignored if NULL.
Return Value:
   Double, elevation of the given point above (positive value) or below
   (negative) the ellipsoid surface. NEMO_DOUBLE_UNDEF if the solution did
   not converge in MAX_ITER steps - the iteration criterion (iterCriter)
   was below the minimum (the value of which depends on the size of the
   planet and the numerical significance inherent in double precision real
   number arithmetic, or if the given point (ptR) was too close to the
   ellipsoid center.
See Also:
   nemo_EnrToCc3()

Synopsis:
   #include <nemo.h>
   void nemo_NcsToEnr(const nemoElRot *pElr,
                      const nemoPtNcs *pPtNcs,
                      nemoPtEnr *pPtEnr);
Description:
   Convert Near Conformal Sphere coordinates into ellipsoid normal
   expressed not in the common φ,λ angular form, but instead in form
   of an i,j,k unit vector.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtNcs:
      Pointer to nemoPtNcs structure, given NCS point coordinates.
   pPtEnt:
      Pointer to nemoPtEnt structure, returned ellipsoid coordinates
      in vector form.
See Also:
   nemo_EllToNcs()
   nemo_NcsToEll()
   nemo_EnrToNcs()

Synopsis:
   #include <nemo.h>
   void nemo_EnrToNcs(const nemoElRot *pElr,
                      const nemoPtEnr *pPtEnr,
                      nemoPtNcs *pPtNcs);
Description:
   Convert ellipsoid coordinates in vector form to a Near Conformal sphere.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtEnr:
      Pointer to nemoPtEnr structure, given ellipsoid normal in i,j,k
      form.
   pPtNcs:
      Pointer to nemoPtNcs structure, returned NCS coordinates.
See Also:
   nemo_EllToNcs()
   nemo_NcsToEll()
   nemo_NcsToEnr()

Synopsis:
   #include <nemo.h>
   void nemo_EllipsoidPrincipalRadii(const nemoElRot *pElr,
                                     const nemoPtEnr *pPtEnr,
                                     double *pPrvRad,
                                     double *pMrdRad);
Description:
   Find two principal radii of curvature at a given point on the ellipsoid
   surface: along the meridian and along the Prime Vertical (a plane
   perpendicular to the meridian plane, that contains the ellipsoid surface
   normal).
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtEnr:
      Pointer to nemoPtEnr structure, given point ellipsoid normal in i,j,k
      form.
   pPrvRad:
      pointer to a receiving double variable, radius of curvatures along the
      Prime Vertical.
   pMrdRad:
      pointer to a receiving double variable, radius of curvatures along the
      meridian.
See Also:
   nemo_EllipsoidRadius()

Synopsis:
   #include <nemo.h>
   double nemo_EllipsoidRadius(const nemoElRot *pElr,
                               const nemoDxPln *pDir,
                               double prvRad,
                               double mrdRad);
Description:
   Find the radius of curvature along a line at arbitrary azimuth (i.e., the
   "general" ellipsoid surface curvature radius). Note that the computation
   of general curvature radius requires the knowledge of two principal radii
   - as returned by a previous invocation of nemo_EllipsoidPrincipalRadii()).
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pDir:
      pointer to nemoDxPln structure, given azimuth of the line along
      which the radius of curvature in point pPtEnr is required.
   prvRad:
      double, given radius of curvatures along the Prime Vertical.
   mrdRad:
      double, radius of curvatures along the meridian.
Return Value:
   Returned radius of curvature, as a linear measure in the same units
   as pElr.
See Also:
   nemo_EllipsoidPrincipalRadii()

Synopsis:
   #include <nemo.h>
   int nemo_EllipsoidChordDirect(const nemoElRot *pElr,
                                 const nemoPtEnr *pPtEnrA,
                                 const nemoDxPln *pDxAB,
                                 double chordElrAB,
                                 nemoPtEnr *pPtEnrB,
                                 double iterCriter,
                                 int *iterCount);
Description:
   Solve the "direct geodetic problem" on the ellipsoid surface: given are the
   coordinates of a point, find the coordinates of another point at a given
   direction and distance from the given point. In this function, "distance"
   is defined as the chord between two ellipsoid surface points. Since the  
   intersection between the chord and the spheroidal surface becomes 
   ill-defined as the angular distance between two points approaches π,
   the function is restricted to distance that are less than planetary
   quadrant. Note that "direct" versus "inverse" problem are historical 
   terms, and somewhat misleading, as the inverse problem is usually solved 
   directly, while the direct problem can be solved onlt by iteration.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   ptA:
      Pointer to nemoPtEnr structure, given ellipsoid point.
   pDxAB:
      Pointer to nemoDxPln structure, given direction from ptA to ptB in
      the tangential plane of pPtEnrA.
   chordElrAB:
      Double, chord distance between pPtEnrA and pPtEnrB.
   pPtEnrB:
      Pointer to nemoPtEnr structure, returned coordinates of the unknown
      point as defined in the decsription.
   iterCriter:
      Double, iteration criterion, linear value in the same units as pElr.
      For instance, if the value of the argument is 0.010, the iteration
      will return the position of the new point when the inverse problem
      calculation on two points returns the direction and distance deltas 
      within given value (measured on the planetary surface) from the 
      computed new point ellipsoid coordinates.
   iterCount:
      Pointer to a natural integer variable, returned number of iterations
      performed (if successful, required to reach the solution, if not, the 
      maximum allowed), ignored if NULL.
Return Value:
   Integer, success/failure indicator:
      > 2 Internal assertion error.
      = 2 Failed to converge in maximum allowed number of steps.
      = 1 Given chord distance is above the limit.
      = 0 Otherwise.
See Also:
   nemo_EllipsoidChordInverse()

Synopsis:
   #include <nemo.h>
   double nemo_EllipsoidChordInverse(const nemoElRot *pElr,
                                 const nemoPtEnr *pPtEnrA,
                                 const nemoPtEnr *pPtEnrB,
                                 nemoDxPln *pDxAB, nemoDxPln *pDxBA);
Description:
   Solve the "Inverse geodetic problem" on the ellipsoid surface: given are
   coordinates of two non-coincident points, find the length of the chord
   between them and the direction of the line to the other given point in
   each given point's tangential plane.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtEnrA:
      Pointer to nemoPtEnr structure, first given ellipsoid point.
   pPtEnrB:
      Pointer to nemoPtEnr structure, second given ellipsoid point.
   pDxAB:
      Pointer to nemoDxPln structure or NULL, returned direction from
      ptA to ptB in the tangential plane of ptA. Note that the direction
      is in the vertical intersection plane of ptA. 
   pDxBA:
      Pointer to nemoDxPln structure or NULL, returned direction from
      ptB to ptA in the tangential plane of ptB. Note that the direction
      is in the vertical intersection plane of ptB. (If either of two
      direction pointers are NULL, the invoker does not require the
      respective direction).
Return Value:
   Double, square of the chord distance between ptA and ptB.
See Also:
   nemo_EllipsoidChordDirect()

Synopsis:
   #include <nemo.h>
   void nemo_EnrToCfs(const nemoElRot *pElr,
                      const nemoPtEnr *pPtEnr,
                      nemoPtCfs *pPtCfs);
Description:
   Convert ellipsoid coordinates in vector form to a Conformal sphere.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtEnr:
      Pointer to nemoPtEnr structure, given ellipsoid coordinates.
   pPtCfs:
      Pointer to nemoPtCfs structure, returned conformal sphere coordinates.
See Also:
   nemo_CfsToEnr()

Synopsis:
   #include <nemo.h>
   void nemo_CfsToEnr(const nemoElRot *pElr,
                      const nemoPtCfs *pPtCfs,
                      nemoPtEnr *pPtEnr);
Description:
   Convert ConFormal sphere coordinates into ellipsoid coordinates.
Arguments:
   pElr:
      Pointer to nemoElRot structure, given ellipsoid of rotation parameters.
   pPtCfs:
      Pointer to nemoPtCfs structure, given conformal sphere coordinates.
   pPtEnr:
      Pointer to nemoPtEnr structure, returned ellipsoid coordinates.
See Also:
   nemo_EnrToCfs()