* converts from polar to cartesian based on Airy 1830

* converts from cartesian to polar based on Airy 1830
* need to handle different ellipsoids

Author: Nathaniel McHugh

geospatial.c

@@ -41,13 +41,26 @@ ZEND_BEGIN_ARG_INFO_EX(helmert_args, 0, 0, 3)
     ZEND_ARG_INFO(0, z)
 ZEND_END_ARG_INFO()
 
+ZEND_BEGIN_ARG_INFO_EX(polar_to_cartesian_args, 0, 0, 2)
+    ZEND_ARG_INFO(0, latitude)
+    ZEND_ARG_INFO(0, longitude)
+ZEND_END_ARG_INFO()
+
+ZEND_BEGIN_ARG_INFO_EX(cartesian_to_polar_args, 0, 0, 3)
+    ZEND_ARG_INFO(0, x)
+    ZEND_ARG_INFO(0, y)
+    ZEND_ARG_INFO(0, z)
+ZEND_END_ARG_INFO()
+
 /* {{{ geospatial_functions[]
  *
  * Every user visible function must have an entry in geospatial_functions[].
  */
 const zend_function_entry geospatial_functions[] = {
     PHP_FE(haversine, haversine_args)
 	PHP_FE(helmert, helmert_args)
+    PHP_FE(polar_to_cartesian, polar_to_cartesian_args)
+    PHP_FE(cartesian_to_polar, cartesian_to_polar_args)
     /* End of functions */
 	{ NULL, NULL, NULL }
 };
@@ -145,6 +158,61 @@ PHP_FUNCTION(helmert)
     add_next_index_double(return_value, zOut);
 }
 
+PHP_FUNCTION(polar_to_cartesian)
+{
+    double latitude, longitude;
+    double x, y, z;
+    if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &latitude, &longitude) == FAILURE) {
+        return;
+    }
+
+    array_init(return_value);
+    double phi = latitude * GEO_DEG_TO_RAD;
+    double lambda = longitude * GEO_DEG_TO_RAD;
+    double eSq = ((AIRY_1830_A * AIRY_1830_A)  - (AIRY_1830_B * AIRY_1830_B))  /  (AIRY_1830_A * AIRY_1830_A);
+    double nu = AIRY_1830_A / sqrt(1 - (eSq * sin(phi) * sin(phi)));
+    x = nu + HEIGHT;
+    x *= cos(phi) * cos(lambda);
+    y = nu + HEIGHT;
+    y *= cos(phi) * sin(lambda);
+    z = ((1 - eSq) * nu) + HEIGHT;
+    z*= sin(phi);
+    add_next_index_double(return_value, x);
+    add_next_index_double(return_value, y);
+    add_next_index_double(return_value, z);
+
+}
+
+PHP_FUNCTION(cartesian_to_polar)
+{
+    double latitude, longitude;
+    double x, y, z;
+    double nu, lambda, h;
+
+    if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &y, &z) == FAILURE) {
+        return;
+    }
+    //aiming for 1m accuracy
+    double precision = 0.1 / AIRY_1830_A;
+    array_init(return_value);
+    double eSq = ((AIRY_1830_A * AIRY_1830_A)  - (AIRY_1830_B * AIRY_1830_B))  /  (AIRY_1830_A * AIRY_1830_A);
+    double p = sqrt(x * x + y * y);
+    double phi = atan2(z, p * (1 - eSq));
+    double phiP = 2 * M_PI;
+    while (abs(phi - phiP) > precision) {
+        nu = AIRY_1830_A / sqrt(1 - eSq * sin(phi) * sin(phi));
+        phiP = phi;
+        phi = atan2(z + eSq * nu * sin(phi), p);
+    }
+    lambda = atan2(y ,x);
+    h = p / cos(phi) - nu;
+
+    add_assoc_double(return_value, "lat", phi / GEO_DEG_TO_RAD);
+    add_assoc_double(return_value, "long", lambda / GEO_DEG_TO_RAD);
+    add_assoc_double(return_value, "height", h);
+
+}
+
 /* {{{ proto haversine(double fromLat, double fromLong, double toLat, double toLong [, double radius ])
  * Calculates the greater circle distance between the two lattitude/longitude pairs */
 PHP_FUNCTION(haversine)

php_geospatial.h

@@ -56,12 +56,16 @@ extern zend_module_entry geospatial_module_entry;
 #define ROTATION_Y -0.2470;
 #define ROTATION_Z -0.8421;
 
+#define HEIGHT 24.7;
+
 
 PHP_MINIT_FUNCTION(geospatial);
 PHP_MINFO_FUNCTION(geospatial);
 
 PHP_FUNCTION(haversine);
 PHP_FUNCTION(helmert);
+PHP_FUNCTION(polar_to_cartesian);
+PHP_FUNCTION(cartesian_to_polar);
 
 #endif	/* PHP_GEOSPATIAL_H */
 

tests/cartesian_to_polar.phpt

@@ -0,0 +1,17 @@
+--TEST--
+Test polar to cartesian
+--FILE--
+<?php
+$x = 3810891.6734396;
+$y = 97591.624686311;
+$z = 5095766.3939034;
+
+$polar = cartesian_to_polar($x, $y, $z);
+echo round($polar['lat'] ,6),PHP_EOL;
+echo round($polar['long'] ,6),PHP_EOL;
+echo round($polar['height'] ,3),PHP_EOL;
+?>
+--EXPECT--
+53.383611
+1.466944
+24.7

tests/helmert.phpt

@@ -11,9 +11,9 @@ var_dump(helmert($x, $y, $z));
 --EXPECT--
 array(3) {
   [0]=>
-  float(3909460.068)
+  float(3909460.0676711)
   [1]=>
-  float(-146987.30)
+  float(-146987.30138174)
   [2]=>
-  float(5019888.07)
+  float(5019888.0705933)
 }

tests/polar_to_cartesian.phpt

@@ -0,0 +1,18 @@
+--TEST--
+Test polar to cartesian
+--FILE--
+<?php
+$lat = 53.38361111111;
+$long = 1.4669444444;
+
+var_dump(polar_to_cartesian($lat, $long));
+?>
+--EXPECT--
+array(3) {
+  [0]=>
+  float(3810891.6734396)
+  [1]=>
+  float(97591.624686311)
+  [2]=>
+  float(5095766.3939034)
+}