Retrun to document top page
nemo.h
/* nemo.h: Nemo library of C functions, header file. */
#ifndef NEMO_H_INCLUDED
#define NEMO_H_INCLUDED
#include <stdlib.h>
#include <stddef.h>
#include <stdint.h>
#include <memory.h>
#include <math.h>
/* ==========================================================================
Library release is the publication day in YYYY.DDD format, cf. date +%Y.%j */
#define NEMO_LIBRARY_DATE 2023.215 /* (so it compares like a decimal number) */
/* ========================================================================== */
/* This header file defines Nemo Library C language implementation API
(constants, macros, types and function prototypes) in seven sections:
(0) Preliminaries
(1) Generic programming elements
(2) Trigonometry and numerical paraphernalia
(3) Spherical surface computations
(4) Concise Direction Cosines global location encoding
(5) Ellipsoid of rotation computations
(6) Spherical cartography
Naming typography convention used by the library
(and suggested for applications using it):
------------------------------------------
NEMO_MACRO_NAME manifest constant
NEMO_MacroNameXyz() macro with arguments
nemo_FunctionName() library function
nemoSomeType defined object type
Abbreviations, data types and function name mnemonics:
------------------------------------------------------
Pln Point coordinates (x,y) in plane
Cc3 As above, (x,y,z) in 3-dimensional Cartesian coordinate system
Ncs Near-conformal sphere (i,j,k) location as unit vector
Cfs As above, conformal sphere
U64 Location on Ncs in form of an 64-bit (8 byte) unsigned integer
Elr Ellipsoid of rotation
Ell Ellipsoid coordinates as φ (latitude) and λ (longitude)
Enr Ellipsoid point defined as surface normal (i,j,k) unit vector
V2 generic 2-component vector
V3 generic 3-component vector
Library documentation constants and mathematical symbols using Greek
alphabet letters and/or other frequently used UTF-8 graphemes, present
in nearly every font set and tolerated in comments and string litarals
by most current C compilers:
α Angle or azimuth
β Commonly used complement of the above or second angle of a triangle
γ Third angle of a triangle
π Circle circumference to diameter ratio
φ Latitude
λ Longitude
τ Golden ratio (sectio aurea)
Δ Difference (typically between exact value and an approximation)
σ Standard deviation
ε Spherical "excess" or tolerance or magnitude of numerical resolution
× Multiplication operator
÷ Division operator
² Square oerator
√ Square root oerator
° Degree sign
± Plus-or-minus (+/-) sign
Programmer: Hrvoje Lukatela, <www.lukatela.com/hrvoje>
Distributed under 'BSD type' license, the full text of which can be
found in the "license.text" file, co-resident with this source code file.
To contact the author, see: <www.lukatela.com/hrvoje>
*/
/* =======================
(Sect.00) Preliminaries
=======================
*/
/* Trigonometry and numerical paraphernalia
---------------------------------------- */
#define NEMO_SQRT2 1.41421356237309504880
#define NEMO_SQRT3 1.73205080756887719318
#define NEMO_PI 3.14159265358979323846
#define NEMO_TWOPI 6.28318530717958647692
#define NEMO_PIHALF 1.57079632679489661923
#define NEMO_PIQUART 0.78539816339744830962
#define NEMO_DEG2RAD (NEMO_PI / 180.0) /* convert decimal degrees to radians */
#define NEMO_RAD2DEG (180.0 / NEMO_PI) /* ...and vice versa */
#define NEMO_SECTIO_AUREA 1.61803398874989490253
/* C standard (as well as most compilers) have no direct support for NaN
("not-a-number"). This "flag value" is safe with IEEE floating point
format - used by the majority of contemporary hardware architectures: */
#define NEMO_DOUBLE_UNDEF -1.79769313486231e+308
/* number larger than any that will appear in the application context: */
#define NEMO_DOUBLE_HUGE 1.79769313486230e+308
/* loss of numerical significance; corresponding to 0.010 (millimeters) on... */
#define NEMO_FUZZ 1.57e-12 /* ...Earth's... */
#define NEMO_FUZZ_SQ 2.46e-24 /* ...surface */
#define NEMO_FUZZ_SQDBL (NEMO_FUZZ_SQ + NEMO_FUZZ_SQ) /* vertex concidence */
/* Nemo point coordinates (on plane, sphere or ellipsdois) are short arrays,
that represent coordinate components. Code may be easier to read if the
following notation is used for array element indices:
*/
#define NEMO_X 0 /* x, y in plane */
#define NEMO_Y 1
#define NEMO_Z 2 /* z (in some solid geometry productions) */
#define NEMO_EST 0 /* East, North in plane */
#define NEMO_NRT 1
#define NEMO_LAT 0 /* Latitude (φ), longitude (λ) on sphere or ellipsoid */
#define NEMO_LNG 1
#define NEMO_I 0 /* Direction cosine vector components */
#define NEMO_J 1
#define NEMO_K 2
/* Most recent hardware provides efficient simultaneous evaluation of
sine and cosine, but C language standard math library does not offer
a method to use it. With __USE_GNU possible alternative is gcc's
"void sincos (double x, double *sinx, double *cosx)" whose signature
this macro mimics. (The macro could be extended for non-gcc compilers).
*/
#if defined __GNUC__
#define NEMO_SinCos(ang, pSin, pCos) \
asm("fsincos" : "=t" (*pCos), "=u" (*pSin) : "0" (ang))
#else
#define NEMO_SinCos(ang, pSin, pCos) { *pSin = sin(ang); *pCos = cos(ang); }
#endif
/* used in testing and simulations: real (C double) random number
in the range rMin to rMax: */
#define NEMO_RandomInRange(rMin, rMax) \
((rMin) + ((rMax) - (rMin)) * ((double)(rand()) / (double)RAND_MAX))
/* The fundamental point coordinates; in plane and in 3-dimensional space: */
typedef struct {
double r[2]; /* "r": planar rectangular coordinates... */
} nemoPtPln; /* ... x, y (or east, north) in a plane */
typedef struct {
double s[3]; /* "s": Cartesian coordinates in a 3D system... */
} nemoPtCc3; /* ... x, y, z, usually geocentric space */
/* Azimuth in tangential or cartographic plane, as two direction cosines */
typedef struct {
double dc[2]; /* "dc": line direction cosines relative to... */
} nemoDxPln; /* ... x, y (or east, north) coordinate axes */
/* Generic direction in in 3-dimensional space, as three direction cosines */
typedef struct {
double dc[3]; /* "dc": line direction cosines relative to... */
} nemoDxCc3; /* ... x, y, z coordinate axes */
/* Two frequently used macros: square of chord between two points...: */
/* ...on unit circle */
#define NEMO_ChordSq2(a, b) \
(((b)[0] - (a)[0]) * ((b)[0] - (a)[0]) + \
((b)[1] - (a)[1]) * ((b)[1] - (a)[1]))
/* ...on unit sphere */
#define NEMO_ChordSq3(a, b) \
(((b)[0] - (a)[0]) * ((b)[0] - (a)[0]) + \
((b)[1] - (a)[1]) * ((b)[1] - (a)[1]) + \
((b)[2] - (a)[2]) * ((b)[2] - (a)[2]))
/* ======================================================================
(Sect.01) Elementary programming functions (not in C standard library)
======================================================================
*/
/* Sub-section A: Workspace memory */
#define NEMO_WORKMEM_ALIGN 8 /* work memory allocation granularity */
/* adjust (size_t, ssize_t) integer variable upwards to memory granularity */
#define NEMO_SizeAlign(s) { while (s % NEMO_WORKMEM_ALIGN) s++; }
typedef void *nemoWorkMem; /* workspace memory instance */
int nemo_WorkMemInit(size_t, nemoWorkMem);
void *nemo_WorkMemGet(size_t, nemoWorkMem);
void *nemo_WorkMemCheck(size_t, const nemoWorkMem);
void nemo_WorkMemRevert(void *, nemoWorkMem);
size_t nemo_WorkMemAvailable(const nemoWorkMem);
size_t nemo_WorkMemUsedMax(const nemoWorkMem);
/* The following macro simplifies the return from a function that encountered
an out-of-work-memory condition: reset memory control block (wm) next-free
pointer to the value it had (rp) when receiving the control, and return to
the invoker with the error indicator (ir, typically set to 9):
*/
#define NEMO_WorkMemBail(rp, wm, ir) { nemo_WorkMemRevert(rp, wm); return(ir); }
/* Sub-section B: Stack (LIFO list) */
/* Stack and queue "anchors" are "opaque" to the invoker */
typedef unsigned char nemoStack[3 * sizeof(void *)]; /* stack anchor */
typedef struct nemoStackNode {
struct nemoStackNode *previous; /* pointer to previous node */
unsigned char data[];
} nemoStackNode; /* item on stack (stack node) */
void nemo_StackInit(nemoStack);
int nemo_Push(nemoStack, unsigned char[], int, nemoWorkMem);
nemoStackNode *nemo_Pop(nemoStack);
nemoStackNode *nemo_PeekTop(const nemoStack);
nemoStackNode *nemo_PeekBottom(const nemoStack);
void nemo_StackInvert(nemoStack);
int nemo_StackVerify(const nemoStack);
size_t nemo_StackSize(const nemoStack);
/* Sub-section C: Queue (FIFO list) */
typedef unsigned char nemoQueue[3 * sizeof(void *)]; /* queue anchor */
typedef struct nemoQueueNode {
struct nemoQueueNode *next; /* pointer to next node */
unsigned char data[];
} nemoQueueNode; /* item in queue (queue node) */
void nemo_QueueInit(nemoQueue);
int nemo_EnQueue(nemoQueue, unsigned char[], int, nemoWorkMem);
nemoQueueNode *nemo_DeQueue(nemoQueue);
nemoQueueNode *nemo_QueuePeekFirst(const nemoQueue);
nemoQueueNode *nemo_QueuePeekLast(const nemoQueue);
int nemo_QueueVerify(const nemoQueue);
size_t nemo_QueueSize(const nemoQueue);
/* ==================================================
(Sect.02) Computational geometry utility functions
==================================================
*/
void nemo_NormalizeV2(double[2]);
void nemo_NormalizeV3(double[3]);
void nemo_MidV3(const double[3], const double[3], double[3]);
double nemo_ArcV3(const double[3], const double[3]);
double nemo_ChordSqToArcApprox(double);
double nemo_ArcToChordApprox(double);
void nemo_LatLongToDcos3(const double[2], double[3]);
void nemo_Dcos3ToLatLong(const double[3], double[2]);
double nemo_DirectionToAzimuth(const double[2]);
/* ======================================================================
(Sect.03) Near-conformal sphere ("NCS") - or any unit (r = 1.0) sphere
======================================================================
*/
typedef struct {
double dc[3]; /* "dc": i, j, k unit vector direction cosines */
} nemoPtNcs; /* point as a unit vector on NCS */
/* For possible comparison between NCS and (rigorous) ConFormal Sphere (CFS)*/
typedef struct {
double dc[3]; /* "dc": i, j, k unit vector direction cosines */
} nemoPtCfs; /* point as a unit vector on CFS */
#define NEMO_GNOMONIC_PCNT 9 /* doubly defined, as used in local random */
nemoPtNcs *nemo_SphereRandomLocation(nemoPtNcs *);
void nemo_SphereRandomLine(const nemoPtNcs *, double, nemoPtNcs *, nemoPtNcs *);
nemoPtNcs *nemo_RandomLocal(const nemoPtNcs *, double,
double[NEMO_GNOMONIC_PCNT], nemoPtNcs *);
int nemo_SphereCircumcenter(const nemoPtNcs *, const nemoPtNcs *,
const nemoPtNcs *, nemoPtNcs *);
int nemo_SphereChordCoords(const nemoPtNcs *, const nemoPtNcs *,
const nemoPtNcs *,
nemoPtNcs *, double *, double *, double *);
void nemo_SphereTangentialPlane(const nemoPtNcs *, nemoDxCc3 *, nemoDxCc3 *);
int nemo_SphereChordDirect(const nemoPtNcs *, const nemoDxPln *, double,
nemoPtNcs *);
double nemo_SphereChordInverse(const nemoPtNcs *, const nemoPtNcs *,
nemoDxPln *, nemoDxPln *);
int nemo_SphereSegsXsct(const nemoPtNcs *, const nemoPtNcs *,
const nemoPtNcs *, const nemoPtNcs *,
double, nemoPtNcs *, int *);
/* ===========================================================================
(Sect.04) Concise Direction Cosines (CDC): point as unsigned 8-byte integer
===========================================================================
*/
typedef uint64_t nemoPtU64; /* "Concise coord's": global loc. as 8-byte uint */
#define NEMO_U64_UNDEF UINT64_MAX /* undefined/unknown pt8 location */
/* CDC function prototypes */
nemoPtU64 nemo_NcsToU64(const nemoPtNcs *);
void nemo_U64ToNcs(nemoPtU64, nemoPtNcs *);
/* =============================================================
(Sect.05) Ellipsoid of rotation (sometimes called "spheroid")
=============================================================
*/
#define NEMO_WGS84_A 6378137.0 /* wgs84 semi-major axis */
/* wgs84 ellipsoid is defined by semi-major axis and flattening ((a - b) / a).
Consequently, the second line, derived from flattening, is a better
numerical representation than the first one, commonly used.
(1/7-th of a mm is however of no practical significance). */
/* #def NEMO_WGS84_B 6356752.31414 commonly published wgs84 minor axis */
#define NEMO_WGS84_B 6356752.3142451793 /* better, derived from flattening */
/* in some very approximate computations Earth can be assumed to be a sphere: */
#define NEMO_EARTH_RADIUS ((NEMO_WGS84_A + NEMO_WGS84_A + NEMO_WGS84_B) / 3.0)
/* as above, for an arbitrary ellipsoidal planet (watch explicit indices): */
#define NEMO_PlanetRadius(pPln) ((pPln->p[0] + pPln->p[0] + pPln->p[1]) / 3.0)
#define NEMO_ELR_PCNT 8
typedef struct {
double p[NEMO_ELR_PCNT]; /* parameters frequently used in computations */
} nemoElRot; /* planet, an ellipsoid of rotation */
typedef struct {
double a[2]; /* "a": angle, φ and λ respectively, in radians */
} nemoPtEll; /* point on ellipsoid as φ and λ */
typedef struct {
double dc[3]; /* "dc": direction cosines of ellipsoid normal */
} nemoPtEnr; /* point on ellipsoid defined by its normal unit vector */
/* If the ellipsoid is Earth, the first of the following two coefficients can
be used as a multiplier of geodesic length to yield the arc, and from it
the chord squared, which can be used as rejection criterion for points
guaranteed to be further away from a given point than is the length of
that geodesic, for all geodesics not exceeding the quadrant and a quarter;
the second can be used in a similar manner to include points guaranteed to
be closer than the given geodesic length. In other words, only points with
squared chord between two arc values need to be subjected to the expensive
ellipsoid coordinate transformation and geodesic length computation.
*/
#define NEMO_GEOARC_MAX 1.00225
#define NEMO_GEOARC_MIN 0.99775
void nemo_ElrInit(double, double, nemoElRot *);
const nemoElRot *nemo_ElrWgs84(void);
double nemo_MeridianArcLength(const nemoElRot *, const nemoPtEnr *);
double nemo_GeodesicVincenty(const nemoElRot *, const nemoPtEnr *,
const nemoPtEnr *, int *);
/* Ellipsoid function prototypes: initallization and NCS transformations */
void nemo_EllToNcs(const nemoElRot *, const nemoPtEll *, nemoPtNcs *);
void nemo_NcsToEll(const nemoElRot *, const nemoPtNcs *, nemoPtEll *);
double nemo_NcsElrScale(const nemoElRot *, const nemoPtNcs *);
double nemo_NcsElrRadius(const nemoElRot *, const nemoPtNcs *);
/* Less frequently used 'ellipsoid normal as i,j,k vector' transformations */
void nemo_NcsToEnr(const nemoElRot *, const nemoPtNcs *, nemoPtEnr *);
void nemo_EnrToNcs(const nemoElRot *, const nemoPtEnr *, nemoPtNcs *);
/* Basic geometry of the ellipsoid of rotation */
double nemo_EnrToCc3(const nemoElRot *, const nemoPtEnr *, double,
nemoPtCc3 *, nemoPtCc3 *);
double nemo_Cc3ToEnr(const nemoElRot *, const nemoPtCc3 *, nemoPtEnr *,
double, int *);
void nemo_EllipsoidPrincipalRadii(const nemoElRot *, const nemoPtEnr *,
double *, double *);
int nemo_EllipsoidChordDirect(const nemoElRot *,
const nemoPtEnr *, const nemoDxPln *, double,
nemoPtEnr *, double, int *);
double nemo_EllipsoidChordInverse(const nemoElRot *,
const nemoPtEnr *, const nemoPtEnr *,
nemoDxPln *, nemoDxPln *);
double nemo_EllipsoidRadius(const nemoElRot *, const nemoDxPln *,
double, double);
/* Ellipsoid of rotation to/from conformal sphere transformations */
void nemo_EnrToCfs(const nemoElRot *, const nemoPtEnr *, nemoPtCfs *);
void nemo_CfsToEnr(const nemoElRot *, const nemoPtCfs *, nemoPtEnr *);
/* ====================================
(Sect.06) Nemo spherical cartography
====================================
2 "must-have" + 3 "nice-to-have" cartographic projections. Other projections
may be useful, but stereographic and "stylindrical" are essential for any
application that depicts large volume, transient geographical data on the
computer display. Ortographic is useful for depicting orbital geometry,
gnomonic is often used in sphere/plane conformal approximations and
conical projection is traditionally used in mid-latitudes for regions
extending mostly in longitude direction.
*/
/* Manifest constants, projection indicators */
#define NEMO_PROJECTION_UNDEF 0
#define NEMO_GNOMONIC 1
#define NEMO_STEREOGRAPHIC 2
#define NEMO_ORTHOGRAPHIC 3
#define NEMO_CONICAL 4
#define NEMO_STYLINDRICAL 5
/* Manifest constants, projection parameter array sizes */
#define NEMO_GNOMONIC_PCNT 9
#define NEMO_STEREOGRAPHIC_PCNT 9
#define NEMO_ORTHOGRAPHIC_PCNT 12
#define NEMO_CONICAL_PCNT 8
#define NEMO_STYLINDRICAL_PCNT 8
/* Stereographic projection - tangency antipodes immediate vicinity */
#define NEMO_ANTIPODAL_FUZZ 5.0e-7
/* Maximum number of doubles required for projection parameters */
#define NEMO_MAX_PCNT 14
int nemo_GnomonicInit(const nemoPtNcs *, nemoPtU64,
double[NEMO_GNOMONIC_PCNT]);
int nemo_GnomonicDirect(const double[NEMO_GNOMONIC_PCNT],
const nemoPtNcs *, nemoPtU64, nemoPtPln *);
int nemo_GnomonicInverse(const double[NEMO_GNOMONIC_PCNT],
const nemoPtPln *, nemoPtNcs *, nemoPtU64 *);
int nemo_StereographicInit(const nemoPtNcs *, nemoPtU64,
double[NEMO_STEREOGRAPHIC_PCNT]);
int nemo_StereographicDirect(const double[NEMO_STEREOGRAPHIC_PCNT],
const nemoPtNcs *, nemoPtU64, nemoPtPln *);
int nemo_StereographicInverse(const double[NEMO_STEREOGRAPHIC_PCNT],
const nemoPtPln *, nemoPtNcs *, nemoPtU64 *);
int nemo_OrthographicInit(const nemoPtNcs *, nemoPtU64,
const nemoPtNcs *, nemoPtU64, double,
double[NEMO_ORTHOGRAPHIC_PCNT]);
int nemo_OrthographicDirect(const double[NEMO_ORTHOGRAPHIC_PCNT],
const nemoPtNcs *, nemoPtU64, nemoPtPln *);
int nemo_OrthographicInverse(const double[NEMO_ORTHOGRAPHIC_PCNT],
const nemoPtPln *, nemoPtNcs *, nemoPtU64 *);
int nemo_ConicalInit(const nemoPtNcs *, nemoPtU64, double[NEMO_CONICAL_PCNT]);
int nemo_ConicalDirect(const double[NEMO_CONICAL_PCNT],
const nemoPtNcs *, nemoPtU64, nemoPtPln *);
int nemo_ConicalInverse(const double[NEMO_CONICAL_PCNT],
const nemoPtPln *, nemoPtNcs *, nemoPtU64 *);
int nemo_StylindricalInit(const nemoPtNcs *, nemoPtU64,
double[NEMO_STYLINDRICAL_PCNT]);
int nemo_StylindricalDirect(const double[NEMO_STYLINDRICAL_PCNT],
const nemoPtNcs *, nemoPtU64, nemoPtPln *);
int nemo_StylindricalInverse(const double[NEMO_STYLINDRICAL_PCNT],
const nemoPtPln *, nemoPtNcs *, nemoPtU64 *);
#endif