vectors.h docs + more conversion functions

Former-commit-id: 28a2b617273f4fd9dea0c78f93edab55c29ca791

qpms/vectors.h

@@ -1,3 +1,6 @@
+/*! \file vectors.h
+ * \brief Coordinate transforms and vector arithmetics.
+ */
 #ifndef VECTORS_H
 #define VECTORS_H
 #include <math.h>
@@ -5,27 +8,32 @@
 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
 #endif
 #include "qpms_types.h"
+#include "qpms_error.h"
 
 //static inline double vectors_h_sq(double x) {return x*x;}
 
 static const cart2_t CART2_ZERO = {0, 0};
 static const cart3_t CART3_ZERO = {0, 0, 0};
 
+/// 2D vector addition.
 static inline cart2_t cart2_add(const cart2_t a, const cart2_t b) {
 	cart2_t res = {a.x+b.x, a.y+b.y};
 	return res;
 }
 
+/// 2D vector substraction.
 static inline cart2_t cart2_substract(const cart2_t a, const cart2_t b) {
 	cart2_t res = {a.x-b.x, a.y-b.y};
 	return res;
 }
 
+/// 2D vector scaling.
 static inline cart2_t cart2_scale(const double c, const cart2_t v) {
 	cart2_t res = {c * v.x, c * v.y};
 	return res;
 }
 
+/// 2D vector dot product.
 static inline double cart2_dot(const cart2_t a, const cart2_t b) {
 	return a.x * b.x + a.y * b.y;
 }
@@ -34,10 +42,12 @@ static inline double cart2_normsq(const cart2_t a) {
 	return cart2_dot(a, a);
 }
 
+/// 2D vector euclidian norm.
 static inline double cart2norm(const cart2_t v) {
 	return hypot(v.x, v.y); //sqrt(v.x*v.x + v.y*v.y);
 }
 
+/// 2D cartesian to polar coordinates conversion. See @ref coord_conversions.
 static inline pol_t cart2pol(const cart2_t cart) {
 	pol_t pol;
 	pol.r = cart2norm(cart);
@@ -45,6 +55,7 @@ static inline pol_t cart2pol(const cart2_t cart) {
 	return pol;
 }
 
+/// Polar to spherical coordinates conversion. See @ref coord_conversions.
 static inline sph_t pol2sph_equator(const pol_t pol) {
 	sph_t sph;
 	sph.r = pol.r;
@@ -53,6 +64,7 @@ static inline sph_t pol2sph_equator(const pol_t pol) {
 	return sph;
 }
 
+/// 2D cartesian to spherical coordinates conversion. See @ref coord_conversions.
 static inline sph_t cart22sph(const cart2_t cart) {
 	sph_t sph;
 	sph.r = cart2norm(cart);
@@ -61,6 +73,19 @@ static inline sph_t cart22sph(const cart2_t cart) {
 	return sph;
 }
 
+/// 1D cartesian to spherical coordinates conversion. See @ref coord_conversions.
+static inline sph_t cart12sph_zaxis(double z) {
+	sph_t sph = {fabs(z), z < 0 ? M_PI : 0, 0};
+	return sph;
+}
+
+/// 1D to 3D cartesian coordinates conversion. See @ref coord_conversions.
+static inline cart3_t cart12cart3z(double z) {
+	cart3_t c = {0, 0, z};
+	return c;
+}
+
+/// 2D to 3D cartesian coordinates conversion. See @ref coord_conversions.
 static inline cart3_t cart22cart3xy(const cart2_t a) {
 	cart3_t c;
 	c.x = a.x;
@@ -74,11 +99,13 @@ static inline cart2_t cart3xy2cart2(const cart3_t a) {
 	return c;
 }
 
+/// 3D vector euclidian norm.
 static inline double cart3norm(const cart3_t v) {
 	return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
 }
 
 
+/// 3D cartesian to spherical coordinates conversion. See @ref coord_conversions.
 static inline sph_t cart2sph(const cart3_t cart) {
 	sph_t sph;
 	sph.r = cart3norm(cart);
@@ -87,6 +114,7 @@ static inline sph_t cart2sph(const cart3_t cart) {
 	return sph;
 }
 
+/// Spherical to 3D cartesian coordinates conversion. See @ref coord_conversions.
 static inline cart3_t sph2cart(const sph_t sph) {
 	cart3_t cart;
 	double sin_th = 
@@ -100,6 +128,7 @@ static inline cart3_t sph2cart(const sph_t sph) {
 	return cart;
 }
 
+/// Polar to 2D cartesian coordinates conversion. See @ref coord_conversions.
 static inline cart2_t pol2cart(const pol_t pol) {
 	cart2_t cart;
 	cart.x = pol.r * cos(pol.phi);
@@ -107,21 +136,32 @@ static inline cart2_t pol2cart(const pol_t pol) {
 	return cart;
 }
 
+/// Polar to 3D cartesian coordinates conversion. See @ref coord_conversions.
+static inline cart3_t pol2cart3_equator(const pol_t pol) {
+	cart2_t c = pol2cart(pol);
+	cart3_t cart3 = {c.x, c.y, 0};
+	return cart3;
+}
+
+/// 3D vector addition.
 static inline cart3_t cart3_add(const cart3_t a, const cart3_t b) {
 	cart3_t res = {a.x+b.x, a.y+b.y, a.z+b.z};
 	return res;
 }
 
+/// 3D vector substraction.
 static inline cart3_t cart3_substract(const cart3_t a, const cart3_t b) {
 	cart3_t res = {a.x-b.x, a.y-b.y, a.z-b.z};
 	return res;
 }
 
+/// 3D vector scaling
 static inline cart3_t cart3_scale(const double c, const cart3_t v) {
 	cart3_t res = {c * v.x, c * v.y, c * v.z};
 	return res;
 }
 
+/// Euclidian distance between two 3D points.
 static inline double cart3_dist(const cart3_t a, const cart3_t b) {
 	return cart3norm(cart3_substract(a,b));
 }
@@ -130,52 +170,62 @@ static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol,
 	return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5;
 }
 
+/// Complex 3D vector scaling.
 static inline ccart3_t ccart3_scale(const complex  double c, const ccart3_t v) {
 	ccart3_t res = {c * v.x, c * v.y, c * v.z};
 	return res;
 }
 
+/// Complex 3D vector adition.
 static inline ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b) {
 	ccart3_t res = {a.x+b.x, a.y+b.y, a.z+b.z};
 	return res;
 }
 
+/// Complex 3D vector substraction.
 static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) {
 	ccart3_t res = {a.x-b.x, a.y-b.y, a.z-b.z};
 	return res;
 }
 
+/// Complex 3D vector (geographic coordinates) addition.
 static inline csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) {
 	csphvec_t res = {a.rc + b.rc, a.thetac + b.thetac, a.phic + b.phic};
 	return res;
 }
 
+/// Complex 3D vector (geographic coordinates) substraction.
 static inline csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b) {
 	csphvec_t res = {a.rc - b.rc, a.thetac - b.thetac, a.phic - b.phic};
 	return res;
 }
 
+/// Complex 3D vector (geographic coordinates) scaling.
 static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) {
 	csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic};
 	return res;
 }
 
+/// Complex 3D vector (geographic coordinates) "dot product" without conjugation.
 static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
 	//N.B. no complex conjugation done here
 	return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic;
 }
 
+/// 3D vector dot product.
 static inline double cart3_dot(const cart3_t a, const cart3_t b) {
 	return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
-static inline csph_t sph_cscale(complex double c, const sph_t s) {
+/// Spherical coordinate system scaling.
-	csph_t res = {c * s.r, s.theta, s.phi};
+static inline sph_t sph_scale(double c, const sph_t s) {
+	sph_t res = {c * s.r, s.theta, s.phi};
 	return res;
 }
 
-static inline sph_t sph_scale(double c, const sph_t s) {
+/// "Complex spherical" coordinate system scaling.
-	sph_t res = {c * s.r, s.theta, s.phi};
+static inline csph_t sph_cscale(complex double c, const sph_t s) {
+	csph_t res = {c * s.r, s.theta, s.phi};
 	return res;
 }
 
@@ -244,13 +294,16 @@ void print_sph(sph_t);
 
 // kahan sums for various types... TODO make generic code using macros
 
+/// Kanan sum initialisation for ccart3_t.
 static inline void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation) {
 	sum->x = sum->y = sum->z = compensation->x = compensation->y = compensation->z = 0;
 }
+/// Kanan sum initialisation for csphvec_t.
 static inline void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation) {
 	sum->rc = sum->thetac = sum->phic = compensation->rc = compensation->thetac = compensation->phic = 0;
 }
 
+/// Add element to Kahan sum (ccart3_t).
 static inline void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input) {
 	ccart3_t comped_input = ccart3_substract(input, *compensation);
 	ccart3_t nsum = ccart3_add(*sum, comped_input);
@@ -258,6 +311,7 @@ static inline void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const
 	*sum = nsum;
 }
 
+/// Add element to Kahan sum (csphvec_t).
 static inline void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input) {
 	csphvec_t comped_input = csphvec_substract(input, *compensation);
 	csphvec_t nsum = csphvec_add(*sum, comped_input);
@@ -265,6 +319,7 @@ static inline void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, con
 	*sum = nsum;
 }
 
+/// Euclidian norm of a vector in geographic coordinates.
 static inline double csphvec_norm(const csphvec_t a) {
 	return  sqrt(creal(a.rc * conj(a.rc) + a.thetac * conj(a.thetac) + a.phic * conj(a.phic)));
 }
@@ -280,6 +335,151 @@ static inline double csphvec_reldiff(const csphvec_t a, const csphvec_t b) {
 	return csphvec_reldiff_abstol(a, b, 0);
 }
 
+
+/*! \page coord_conversions Coordinate systems and default conversions
+ *
+ * The coordinate system transformations are defined as following:
+ * 
+ * \section coordtf_same_d Equal-dimension coordinate tranforms
+ * \subsection sph_cart3 Spherical and 3D cartesian coordinates
+ * * \f$ x = r \sin \theta \cos \phi \f$,
+ * * \f$ y = r \sin \theta \sin \phi \f$,
+ * * \f$ z = r \cos \theta \f$.
+ * \subsection pol_cart2 Polar and 2D cartesian coordinates
+ * * \f$ x = r \cos \phi \f$,
+ * * \f$ y = r \sin \phi \f$.
+ *
+ * \section coordtf_123 Lower to higher dimension conversions.
+ * * The 1D coordinate is identified with the \a z 3D cartesian coordinate.
+ * * The 2D cartesian coordinates \a x, \a y are identified with the \a x, \a y
+ *   3D cartesian coordinates.
+ * * For the sake of consistency, default conversion between
+ *   1D and 2D coordinates is not allowed and yields NAN values.
+ *
+ * \section coordtf_321 Higher to lower dimension conversions.
+ * Default conversions from higher to lower-dimensional coordinate
+ * systems are not allowed. Any projections have to be done explicitly.
+ */
+
+/// Conversion from anycoord_point_t to explicitly spherical coordinates.
+/** See @ref coord_conversions for the conversion definitions.
+ */ 
+static inline sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t) {
+	switch(t & QPMS_COORDS_BITRANGE) {
+		case QPMS_COORDS_SPH:
+			return p.sph;
+			break;
+		case QPMS_COORDS_POL:
+			return pol2sph_equator(p.pol);
+			break;
+		case QPMS_COORDS_CART3:
+			return cart2sph(p.cart3);
+			break;
+		case QPMS_COORDS_CART2:
+			return cart22sph(p.cart2);
+			break;
+		case QPMS_COORDS_CART1:
+			return cart12sph_zaxis(p.z);
+			break;
+	}
+	QPMS_WTF;
+}
+
+
+/// Conversion from anycoord_point_t to explicitly 3D cartesian coordinates.
+/** See @ref coord_conversions for the conversion definitions.
+ */ 
+static inline cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t) {
+	switch(t & QPMS_COORDS_BITRANGE) {
+		case QPMS_COORDS_SPH:
+			return sph2cart(p.sph);
+			break;
+		case QPMS_COORDS_POL:
+			return pol2cart3_equator(p.pol);
+			break;
+		case QPMS_COORDS_CART3:
+			return p.cart3;
+			break;
+		case QPMS_COORDS_CART2:
+			return cart22cart3xy(p.cart2);
+			break;
+		case QPMS_COORDS_CART1:
+			return cart12cart3z(p.z);
+			break;
+	}
+	QPMS_WTF;
+}
+
+#if 0
+// Convenience identifiers for return values.
+static const cart3_t CART3_INVALID = {NAN, NAN, NAN};
+static const cart2_t CART2_INVALID = {NAN, NAN};
+static const double CART1_INVALID = NAN;
+static const sph_t SPH_INVALID = {NAN, NAN, NAN};
+static const pol_t POL_INVALID = {NAN, NAN};
+#endif
+
+/// Conversion from anycoord_point_t to explicitly polar coordinates.
+/** See @ref coord_conversions for the conversion definitions.
+ */ 
+static inline pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t) {
+	switch(t & QPMS_COORDS_BITRANGE) {
+		case QPMS_COORDS_SPH:
+		case QPMS_COORDS_CART3:
+			QPMS_PR_ERROR("Implicit conversion from 3D to 2D"
+				       	" coordinates not allowed");
+			break;
+		case QPMS_COORDS_POL:
+			return p.pol;
+			break;
+		case QPMS_COORDS_CART2:
+			return cart2pol(p.cart2);
+			break;
+		case QPMS_COORDS_CART1:
+			QPMS_PR_ERROR("Implicit conversion from 1D to 2D"
+					" coordinates not allowed");
+			break;
+	}
+	QPMS_WTF;
+}
+
+
+/// Conversion from anycoord_point_t to explicitly 2D cartesian coordinates.
+/** See @ref coord_conversions for the conversion definitions.
+ */ 
+static inline cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t) {
+	switch(t & QPMS_COORDS_BITRANGE) {
+		case QPMS_COORDS_SPH:
+		case QPMS_COORDS_CART3:
+			QPMS_PR_ERROR("Implicit conversion from 3D to 2D"
+				       	" coordinates not allowed");
+			break;
+		case QPMS_COORDS_POL:
+			return pol2cart(p.pol);
+			break;
+		case QPMS_COORDS_CART2:
+			return p.cart2;
+			break;
+		case QPMS_COORDS_CART1:
+			QPMS_PR_ERROR("Implicit conversion from 1D to 2D"
+					" coordinates not allowed");
+			break;
+	}
+	QPMS_WTF;
+}
+
+
+/// Conversion from anycoord_point_t to explicitly 1D cartesian coordinates.
+/** See @ref coord_conversions for the conversion definitions.
+ */ 
+static inline double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) {
+	if (t & QPMS_COORDS_BITRANGE == QPMS_COORDS_CART1) 
+			return p.z;
+	else
+		QPMS_PR_ERROR("Implicit conversion from nD (n > 1)"
+				" to 1D not allowed.");
+}
+
 typedef double matrix3d[3][3];
 typedef double matrix2d[2][2];
 typedef complex double cmatrix3d[3][3];