aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/pj_obs_api_test.c33
-rw-r--r--examples/s45b.pol402
-rwxr-xr-xnad/testvarious4
-rw-r--r--src/Makefile.am2
-rw-r--r--src/PJ_calcofi.c4
-rw-r--r--src/PJ_cart.c110
-rw-r--r--src/PJ_helmert.c492
-rw-r--r--src/PJ_horner.c953
-rw-r--r--src/PJ_krovak.c2
-rw-r--r--src/PJ_pipeline.c54
-rw-r--r--src/lib_proj.cmake2
-rw-r--r--src/makefile.vc3
-rw-r--r--src/pj_fwd.c85
-rw-r--r--src/pj_fwd3d.c96
-rw-r--r--src/pj_init.c76
-rw-r--r--src/pj_inv.c60
-rw-r--r--src/pj_inv3d.c92
-rw-r--r--src/pj_list.h2
-rw-r--r--src/pj_obs_api.c81
-rw-r--r--src/pj_param.c12
-rw-r--r--src/pj_strerrno.c1
-rw-r--r--src/proj.def41
-rw-r--r--src/proj.h19
23 files changed, 2307 insertions, 319 deletions
diff --git a/examples/pj_obs_api_test.c b/examples/pj_obs_api_test.c
index c75b4a21..bb792c62 100644
--- a/examples/pj_obs_api_test.c
+++ b/examples/pj_obs_api_test.c
@@ -44,8 +44,6 @@
int pj_pipeline_selftest (void);
-
-
int main (void) {
PJ *P;
PJ_OBS a, b;
@@ -71,7 +69,7 @@ int main (void) {
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude) to (12, 55) */
- a = pj_obs_null;
+ a = pj_obs_null ();
/* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */
a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(55);
@@ -100,13 +98,12 @@ int main (void) {
dist = pj_xyz_dist (a.coo.xyz, b.coo.xyz);
printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);
-
/* Invalid projection */
a = pj_trans (P, 42, a);
if (a.coo.lpz.lam!=DBL_MAX)
printf ("%15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);
- err = pj_error (P);
- printf ("pj_error: %d\n", err);
+ err = pj_err_level (P, PJ_ERR_TELL);
+ printf ("pj_err_level: %d\n", err);
/* Clean up */
pj_free (P);
@@ -117,7 +114,7 @@ int main (void) {
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
- a = b = pj_obs_null;
+ a = b = pj_obs_null ();
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
a.coo.lpz.z = 100;
@@ -152,17 +149,18 @@ int main (void) {
P I P E L I N E T E S T S
- ****************************************************************************
-
-
***************************************************************************/
/* forward-reverse geo->utm->geo */
- P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +step +proj=utm +inv");
+ P = pj_create (
+ "+proj=pipeline +ellps=GRS80 +zone=32 +step "
+ "+proj=utm +step "
+ "+proj=utm +inv"
+ );
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
- a = b = pj_obs_null;
+ a = b = pj_obs_null ();
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
a.coo.lpz.z = 100;
@@ -180,12 +178,15 @@ int main (void) {
/* And now the back-to-back situation utm->geo->utm */
- P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +inv +step +proj=utm");
+ P = pj_create (
+ "+proj=pipeline +ellps=GRS80 +zone=32 +step "
+ "+proj=utm +inv +step "
+ "+proj=utm");
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (easting, northing) to utm(12, 55) */
- a = b = pj_obs_null;
+ a = b = pj_obs_null ();
a.coo.xy = cph_utm32;
printf ("PRE: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);
@@ -207,7 +208,7 @@ int main (void) {
return puts ("Oops"), 0;
/* zero initialize everything, then set (easting, northing) to utm(12, 55) */
- a = b = pj_obs_null;
+ a = b = pj_obs_null ();
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
printf ("PRE: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
@@ -221,7 +222,6 @@ int main (void) {
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
-
/* Geodesic distance between two points with angular 2D coordinates */
a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(60);
@@ -237,7 +237,6 @@ int main (void) {
dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
printf ("1 deg at Equator: %15.9f\n", dist);
-
pj_free (P);
pj_pipeline_selftest ();
return 0;
diff --git a/examples/s45b.pol b/examples/s45b.pol
new file mode 100644
index 00000000..9b70fbf7
--- /dev/null
+++ b/examples/s45b.pol
@@ -0,0 +1,402 @@
+################################################################################
+#
+# Horner polynomial evaluation demo data
+#
+################################################################################
+#
+# These are the polynomial coefficients used for transforming to and from
+# the Danish legacy system s45b "System 45 Bornholm"
+#
+################################################################################
+
+<s45b>
+proj=pipeline
+
+step init=./s45b.pol:s45b_tc32
+# step init=./s45b.pol:tc32_utm32
+
+step proj=utm inv ellps=intl zone=32
+step proj=cart ellps=intl
+
+step proj=helmert ellps=GRS80
+ x=-81.0703 y=-89.3603 z=-115.7526
+ rx=-484.88 ry=-24.36 rz=-413.21 s=-0.540645
+
+step proj=cart inv ellps=GRS80
+step proj=utm ellps=GRS80 zone=32
+
+
+
+
+################################################################################
+#
+# S 4 5 B -> T C 3 2
+#
+################################################################################
+<s45b_tc32>
+
+proj=horner
+ellps=intl
+range=500000
+
+fwd_origin=47022.563745,51779.260103
+inv_origin=878354.943082,6125305.175366
+deg=6
+
+
+# static double C_ttb[]
+# tc32_ed50 -> s45b
+# m_lim_gen: 0.153 red = 0 OBS = 1074
+# m = 1.51 cm my_loss = +3 y_enp = +8.4
+# m = 1.53 cm mx_loss = +4 x_enp = +8.4
+
+# mht C_ttb er
+# fwd-inv ombyttet ifht original Poder/Engsager-kode
+# For at opnå at to fwd transform fører fra s45b->tc32->utm32 (->ETRS89)
+
+inv_v=\
+# Poly NORTH :: e-degree = 0 : n-degree = 6
+ 5.1779004699e+04,\
+ 9.9508320295e-01,\
+ -2.9453823207e-10,\
+ 1.9995084102e-14,\
+ -1.4895751366e-18,\
+ -9.9734812211e-23,\
+ 1.1194218845e-26,\
+
+# Poly NORTH :: e-degree = 1 : n-degree = 5
+ -8.4285679515e-02,\
+ -7.9623049286e-09,\
+ -3.7190046062e-14,\
+ -2.3324127411e-18,\
+ -1.1150449763e-22,\
+ 2.8703154270e-27,\
+
+# Poly NORTH :: e-degree = 2 : n-degree = 4
+ 8.7160434140e-10,\
+ -3.3634602927e-14,\
+ -5.5718245313e-18,\
+ 6.2611750909e-23,\
+ -2.1011243838e-26,\
+
+# Poly NORTH :: e-degree = 3 : n-degree = 3
+ 1.0905463989e-14,\
+ -4.3960034360e-18,\
+ 3.6121595001e-22,\
+ -1.3493066011e-27,\
+
+# Poly NORTH :: e-degree = 4 : n-degree = 2
+ -1.3360171462e-18,\
+ 1.0780850646e-22,\
+ 4.5118286607e-26,\
+
+# Poly NORTH :: e-degree = 5 : n-degree = 1
+ -1.3718883973e-22,\
+ 1.6263920750e-26,\
+
+# Poly NORTH :: e-degree = 6 : n-degree = 0
+ -5.1004217526e-27
+
+# tcy 6125305.175366
+
+inv_u=\
+# Poly EAST :: n-degree = 0 : e-degree = 6
+ 4.7022495967e+04,\
+ -9.9508282498e-01,\
+ 3.2436283039e-09,\
+ -2.6276394334e-15,\
+ 8.6318533291e-18,\
+ -3.8327518550e-23,\
+ -2.5704924282e-26,\
+
+# Poly EAST :: n-degree = 1 : e-degree = 5
+ -8.4285975934e-02,\
+ 5.7098765263e-10,\
+ -6.0863955939e-14,\
+ 2.3608788740e-18,\
+ 6.8899581969e-24,\
+ -1.1429511179e-26,\
+
+# Poly EAST :: n-degree = 2 : e-degree = 4
+ -4.6079778412e-09,\
+ 1.5072604543e-14,\
+ 5.4063862378e-18,\
+ 1.2591327827e-22,\
+ 7.9336388691e-27,\
+
+# Poly EAST :: n-degree = 3 : e-degree = 3
+ -2.9479268638e-14,\
+ 1.7090049434e-18,\
+ 2.8413337985e-22,\
+ -3.3577391552e-27,\
+
+# Poly EAST :: n-degree = 4 : e-degree = 2
+ 3.0434879273e-18,\
+ -1.8081673510e-22,\
+ -2.3651419850e-26,\
+
+# Poly EAST :: n-degree = 5 : e-degree = 1
+ 9.2060044804e-23,\
+ 3.7807953325e-27,\
+
+# Poly EAST :: n-degree = 6 : e-degree = 0
+ -4.9415665221e-27
+
+# tcx 878354.943082
+
+
+# static double C_btt[]
+# s45b -> tc32_ed50
+# m_lim_gen: 0.154 red = 0 OBS = 1074
+# m = 1.50 cm my_loss = +3 y_enp = +8.5
+# m = 1.54 cm mx_loss = +4 x_enp = +8.3
+
+fwd_v=\
+# Poly NORTH :: e-degree = 0 : n-degree = 6
+ 6.1253054245e+06,\
+ 9.9778251908e-01,\
+ -7.7346152025e-10,\
+ -2.5359789369e-14,\
+ 1.5614918228e-18,\
+ 9.8091134295e-23,\
+ -1.1092581145e-26,\
+
+# Poly NORTH :: e-degree = 1 : n-degree = 5
+ -8.4514352088e-02,\
+ -7.9847579284e-09,\
+ -2.6865560962e-14,\
+ -2.0731372756e-18,\
+ -1.3660341123e-22,\
+ 1.1244836340e-26,\
+
+# Poly NORTH :: e-degree = 2 : n-degree = 4
+ 8.0551988135e-11,\
+ 3.6661500679e-14,\
+ 5.4247705403e-18,\
+ 8.4494604807e-23,\
+ 1.3334858516e-26,\
+
+# Poly NORTH :: e-degree = 3 : n-degree = 3
+ 8.3889821184e-15,\
+ -4.8124202237e-18,\
+ 2.9088188830e-22,\
+ -2.0129874264e-26,\
+
+# Poly NORTH :: e-degree = 4 : n-degree = 2
+ 2.4716463766e-18,\
+ -2.1717177513e-22,\
+ -3.2828537638e-26,\
+
+# Poly NORTH :: e-degree = 5 : n-degree = 1
+ -1.2080655753e-22,\
+ 2.5050435391e-26,\
+
+# Poly NORTH :: e-degree = 6 : n-degree = 0
+ 1.1383483826e-27
+
+# tcy 51779.260103,
+
+fwd_u=\
+# Poly EAST :: n-degree = 0 : e-degree = 6
+ 8.7835485387e+05,\
+ -9.9778289691e-01,\
+ 3.2537215213e-09,\
+ 6.9217640616e-15,\
+ 8.6268883840e-18,\
+ 4.6748156909e-23,\
+ -2.6492402009e-26,\
+
+# Poly EAST :: n-degree = 1 : e-degree = 5
+ -8.4514648771e-02,\
+ 1.4399520180e-09,\
+ -6.0423329711e-14,\
+ 6.9816167332e-20,\
+ 6.7729233542e-23,\
+ -5.3308251880e-27,\
+
+# Poly EAST :: n-degree = 2 : e-degree = 4
+ -4.5697800099e-09,\
+ -1.5194038814e-14,\
+ 5.1112653016e-18,\
+ -2.0307532869e-22,\
+ 1.0374125432e-26,\
+
+# Poly EAST :: n-degree = 3 : e-degree = 3
+ -2.8983003841e-14,\
+ -1.6414425785e-18,\
+ 1.7874983379e-22,\
+ 1.5492164174e-26,\
+
+# Poly EAST :: n-degree = 4 : e-degree = 2
+ 2.7919197366e-18,\
+ 1.9218613279e-22,\
+ -2.1007264634e-26,\
+
+# Poly EAST :: n-degree = 5 : e-degree = 1
+ 1.0032412389e-22,\
+ -5.9007997846e-27,\
+
+# Poly EAST :: n-degree = 6 : e-degree = 0
+ -4.4410970979e-27
+
+# tcx 47022.563745
+
+
+
+
+
+################################################################################
+#
+# T C 3 2 -> U T M 3 2
+#
+################################################################################
+
+<tc32_utm32>
+
+proj=horner
+ellps=intl
+range=500000
+
+
+fwd_origin=877605.269066,6125810.306769
+inv_origin=877605.760036,6125811.281773
+
+
+# tc32_ed50 -> utm32_ed50 : Bornholm
+
+deg=4
+
+# ttu_n and ttu_e are based on static double C_ttu_b[]
+# m_lim_gen: 0.086 red = 0 OBS = 852
+# m = 1.38 cm my_loss = +2 y_enp = +10.5
+# m = 1.44 cm mx_loss = +2 x_enp = +10.4
+# static double ttu_n[]
+
+fwd_v=\
+# Poly NORTH :: e-degree = 0 : n-degree = 0..4
+ 6.1258112678e+06,\
+ 9.9999971567e-01,\
+ 1.5372750011e-10,\
+ 5.9300860915e-15,\
+ 2.2609497633e-19,\
+
+# Poly NORTH :: e-degree = 1 : n-degree = 0..3
+ 4.3188227445e-05,\
+ 2.8225130416e-10,\
+ 7.8740007114e-16,\
+ -1.7453997279e-19,\
+
+# Poly NORTH :: e-degree = 2 : n-degree = 0..2
+ 1.6877465415e-10,\
+ -1.1234649773e-14,\
+ -1.7042333358e-18,\
+
+# Poly NORTH :: e-degree = 3 : n-degree = 0..1
+ -7.9303467953e-15,\
+ -5.2906832535e-19,\
+
+# Poly NORTH :: e-degree = 4 : n-degree = 0
+ 3.9984284847e-19
+
+# tcy 6125810.306769
+
+
+
+# static double ttu_e[]
+fwd_u=\
+# Poly EAST :: n-degree = 0 : e-degree = 0..4
+ 8.7760574982e+05,\
+ 9.9999752475e-01,\
+ 2.8817299305e-10,\
+ 5.5641310680e-15,\
+ -1.5544700949e-18,\
+
+# Poly EAST :: n-degree = 1 : e-degree = 3
+ -4.1357045890e-05,\
+ 4.2106213519e-11,\
+ 2.8525551629e-14,\
+ -1.9107771273e-18,\
+
+# Poly EAST :: n-degree = 2 : e-degree = 2
+ 3.3615590093e-10,\
+ 2.4380247154e-14,\
+ -2.0241230315e-18,\
+
+# Poly EAST :: n-degree = 3 : e-degree = 1
+ 1.2429019719e-15,\
+ 5.3886155968e-19,\
+
+# Poly EAST :: n-degree = 4 : e-degree = 0
+ -1.0167505000e-18
+
+# tcx 877605.760036
+
+
+# utt_n and utt_e are based on static double C_utt_b[]
+# utm32_ed50 -> tc32_ed50 : Bornholm
+# m_lim_gen: 0.086 red = 0 OBS = 852
+# m = 1.38 cm my_loss = +2 y_enp = +10.8
+# m = 1.44 cm mx_loss = +2 x_enp = +10.7
+# static double utt_n[]
+
+inv_v=\
+# Poly NORTH :: e-degree = 0 : n-degree = 4
+ 6.1258103208e+06,\
+ 1.0000002826e+00,\
+ -1.5372762184e-10,\
+ -5.9304261011e-15,\
+ -2.2612705361e-19,\
+
+# Poly NORTH :: e-degree = 1 : n-degree = 3
+ -4.3188331419e-05,\
+ -2.8225549995e-10,\
+ -7.8529116371e-16,\
+ 1.7476576773e-19,\
+
+# Poly NORTH :: e-degree = 2 : n-degree = 2
+ -1.6875687989e-10,\
+ 1.1236475299e-14,\
+ 1.7042518057e-18,\
+
+# Poly NORTH :: e-degree = 3 : n-degree = 1
+ 7.9300735257e-15,\
+ 5.2881862699e-19,\
+
+# Poly NORTH :: e-degree = 4 : n-degree = 0
+ -3.9990736798e-19
+
+# tcy 6125811.281773
+
+
+
+
+# static double utt_e[]
+inv_u=\
+# Poly EAST :: n-degree = 0 : e-degree = 0..4
+ 8.7760527928e+05,\
+ 1.0000024735e+00,\
+ -2.8817540032e-10,\
+ -5.5627059451e-15,\
+ 1.5543637570e-18,\
+
+# Poly EAST :: n-degree = 1 : e-degree = 0..3
+ 4.1357152105e-05,\
+ -4.2114813612e-11,\
+ -2.8523713454e-14,\
+ 1.9109017837e-18,\
+
+# Poly EAST :: n-degree = 2 : e-degree = 0..2
+ -3.3616407783e-10,\
+ -2.4382678126e-14,\
+ 2.0245020199e-18,\
+
+# Poly EAST :: n-degree = 3 : e-degree = 0..1
+ -1.2441377565e-15,\
+ -5.3885232238e-19,\
+
+# Poly EAST :: n-degree = 4 : e-degree = 0
+ 1.0167203661e-18
+
+# tcx 877605.760036
+
+<>
diff --git a/nad/testvarious b/nad/testvarious
index 58770347..f11bc806 100755
--- a/nad/testvarious
+++ b/nad/testvarious
@@ -283,8 +283,8 @@ $EXE +proj=hammer +datum=WGS84 \
EOF
echo "##############################################################" >> ${OUT}
echo "Test healpix forward projection on sphere" >> ${OUT}
-$EXE +proj=latlong +a=1 +lon_0=0 \
- +to +proj=healpix +a=1 +lon_0=0 -f '%.'5'f' \
+$EXE +proj=latlong +R=1 +lon_0=0 \
+ +to +proj=healpix +R=1 +lon_0=0 -f '%.'5'f' \
-E >>${OUT} <<EOF
0 41.81031
-90 0
diff --git a/src/Makefile.am b/src/Makefile.am
index 69e0e44e..a2508fe8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -74,7 +74,7 @@ libproj_la_SOURCES = \
jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \
pj_strtod.c \
\
- pj_obs_api.c PJ_cart.c PJ_pipeline.c
+ pj_obs_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c
install-exec-local:
rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_calcofi.c b/src/PJ_calcofi.c
index b9a8ca22..0f206155 100644
--- a/src/PJ_calcofi.c
+++ b/src/PJ_calcofi.c
@@ -174,8 +174,8 @@ int pj_calcofi_selftest (void) {
double tolerance_lp = 1e-10;
double tolerance_xy = 1e-7;
- char e_args[] = {"+proj=calcofi +ellps=GRS80 +lat_1=0.5 +lat_2=2"};
- char s_args[] = {"+proj=calcofi +a=6400000 +lat_1=0.5 +lat_2=2"};
+ char e_args[] = {"+proj=calcofi +ellps=GRS80 +lat_1=0.5 +lat_2=2 +no_defs"};
+ char s_args[] = {"+proj=calcofi +R=6400000 +lat_1=0.5 +lat_2=2 +no_defs"};
LP fwd_in[] = {
{ 2, 1},
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 4ec20ced..c585fbeb 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -88,6 +88,16 @@ PROJ_HEAD(cart, "Geodetic/cartesian conversions");
(TF, below).
+ Close to the poles, we avoid singularities by switching to an
+ approximation requiring knowledge of the geocentric radius
+ at the given latitude. For this, we use an adaptation of the
+ formula given in:
+
+ Wikipedia: Earth Radius
+ https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude
+ (Derivation and commentary at http://gis.stackexchange.com/questions/20200/how-do-you-compute-the-earths-radius-at-a-given-geodetic-latitude)
+
+ (WP2, below)
These routines are probably not as robust at those in
geocent.c, at least thay haven't been through as heavy
@@ -103,13 +113,6 @@ static void freeup (PJ *P) {
return;
}
-
-/*********************************************************************/
-static double semiminor_axis (double a, double es) {
-/*********************************************************************/
- return a * sqrt (1 - es);
-}
-
/*********************************************************************/
static double normal_radius_of_curvature (double a, double es, double phi) {
/*********************************************************************/
@@ -120,15 +123,11 @@ static double normal_radius_of_curvature (double a, double es, double phi) {
return a / sqrt (1 - es*s*s);
}
-
/*********************************************************************/
-static double second_eccentricity_squared (double a, double es) {
-/**********************************************************************
- Follows the definition, e'2 = (a2-b2)/b2, but uses the identity
- (a2-b2) = (a-b)(a+b), for improved numerical precision.
-***********************************************************************/
- double b = semiminor_axis (a, es);
- return (a - b) * (a + b) / (b*b);
+static double geocentric_radius (double a, double b, double phi) {
+/*********************************************************************/
+ /* This is from WP2, but uses hypot() for potentially better numerical robustness */
+ return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi));
}
@@ -148,21 +147,6 @@ static XYZ cartesian (LPZ geodetic, PJ *P) {
xyz.y = (N + h) * cosphi * sin(lam);
xyz.z = (N * (1 - P->es) + h) * sin(phi);
- /*********************************************************************/
- /* */
- /* For historical reasons, in proj, plane coordinates are measured */
- /* in units of the semimajor axis. Since 3D handling is grafted on */
- /* later, this is not the case for heights. And even though this */
- /* coordinate really is 3D cartesian, the z-part looks like a height */
- /* to proj. Hence, we have the somewhat unusual situation of having */
- /* a point coordinate with differing units between dimensions. */
- /* */
- /* The scaling and descaling is handled by the pj_fwd/inv functions. */
- /* */
- /*********************************************************************/
- xyz.x /= P->a;
- xyz.y /= P->a;
-
return xyz;
}
@@ -173,15 +157,12 @@ static LPZ geodetic (XYZ cartesian, PJ *P) {
double N, b, p, theta, c, s, e2s;
LPZ lpz;
- cartesian.x *= P->a;
- cartesian.y *= P->a;
-
/* Perpendicular distance from point to Z-axis (HM eq. 5-28) */
p = hypot (cartesian.x, cartesian.y);
/* Ancillary ellipsoidal parameters */
- b = semiminor_axis (P->a, P->es);
- e2s = second_eccentricity_squared (P->a, P->es);
+ b = P->b;
+ e2s = P->e2s;
/* HM eq. (5-37) */
theta = atan2 (cartesian.z * P->a, p * b);
@@ -196,11 +177,12 @@ static LPZ geodetic (XYZ cartesian, PJ *P) {
c = cos(lpz.phi);
if (fabs(c) < 1e-6) {
- /* poleward of 89.99994 deg, we avoid division by zero */
- /* by computing the height as the cartesian z value */
- /* minus the semiminor axis length */
- /* (potential improvement: approx. by 2nd order poly) */
- lpz.z = cartesian.z - (cartesian.z > 0? b: -b);
+ /* poleward of 89.99994 deg, we avoid division by zero */
+ /* by computing the height as the cartesian z value */
+ /* minus the geocentric radius of the Earth at the given */
+ /* latitude */
+ double r = geocentric_radius (P->a, b, lpz.phi);
+ lpz.z = fabs (cartesian.z) - r;
}
else
lpz.z = p / c - N;
@@ -239,6 +221,8 @@ PJ *PROJECTION(cart) {
P->inv3d = geodetic;
P->fwd = cart_forward;
P->inv = cart_reverse;
+ P->left = PJ_IO_UNITS_RADIANS;
+ P->right = PJ_IO_UNITS_METERS;
return P;
}
@@ -255,9 +239,6 @@ int pj_cart_selftest (void) {
char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
- /* Log everything libproj offers to log for you */
- pj_log_level (0, PJ_LOG_TRACE);
-
/* An utm projection on the GRS80 ellipsoid */
P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==P)
@@ -269,7 +250,7 @@ int pj_cart_selftest (void) {
/* Same projection, now using argc/argv style initialization */
P = pj_create_argv (3, args);
if (0==P)
- return puts ("Oops"), 0;
+ return 2;
/* Turn off logging */
pj_log_level (0, PJ_LOG_NONE);
@@ -294,18 +275,18 @@ int pj_cart_selftest (void) {
dist = pj_xy_dist (a.coo.xy, b.coo.xy);
if (dist > 2e-9)
- return 2;
+ return 3;
/* Invalid projection */
a = pj_trans (P, 42, a);
- if (a.coo.lpz.lam!=DBL_MAX)
- return 3;
- err = pj_error (P);
- if (0==err)
+ if (a.coo.lpz.lam!=HUGE_VAL)
return 4;
+ err = pj_err_level (P, PJ_ERR_TELL);
+ if (0==err)
+ return 5;
/* Clear error */
- pj_error_set (P, 0);
+ pj_err_level (P, 0);
/* Clean up */
pj_free (P);
@@ -313,7 +294,7 @@ int pj_cart_selftest (void) {
/* Now do some 3D transformations */
P = pj_create ("+proj=cart +ellps=GRS80");
if (0==P)
- return 5;
+ return 6;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
a = b = pj_obs_null;
@@ -328,9 +309,32 @@ int pj_cart_selftest (void) {
dist = pj_roundtrip (P, PJ_FWD, 10000, a);
dist = pj_roundtrip (P, PJ_INV, 10000, b);
if (dist > 2e-9)
- return 6;
+ return 7;
+
+
+ /* Test at the North Pole */
+ a = b = pj_obs_null;
+ a.coo.lpz.lam = TORAD(0);
+ a.coo.lpz.phi = TORAD(90);
+ a.coo.lpz.z = 100;
+
+ /* Forward projection: Ellipsoidal-to-3D-Cartesian */
+ dist = pj_roundtrip (P, PJ_FWD, 1, a);
+ if (dist > 1e-12)
+ return 8;
+
+ /* Test at the South Pole */
+ a = b = pj_obs_null;
+ a.coo.lpz.lam = TORAD(0);
+ a.coo.lpz.phi = TORAD(-90);
+ a.coo.lpz.z = 100;
+
+ /* Forward projection: Ellipsoidal-to-3D-Cartesian */
+ dist = pj_roundtrip (P, PJ_FWD, 1, a);
+ if (dist > 1e-12)
+ return 9;
- /* Inverse projection: Ellipsoidal-to-3D-Cartesian */
+ /* Inverse projection: 3D-Cartesian-to-Ellipsoidal */
b = pj_trans (P, PJ_INV, b);
/* Move p to another context */
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
new file mode 100644
index 00000000..72cefe90
--- /dev/null
+++ b/src/PJ_helmert.c
@@ -0,0 +1,492 @@
+/***********************************************************************
+
+ Helmert transformations, 3- and 7-parameter
+
+ Thomas Knudsen, 2016-05-24
+
+************************************************************************
+
+ Implements 3- and 7-parameter Helmert transformations for 3D data.
+
+ Primarily useful for implementation of datum shifts in transformation
+ pipelines.
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-24/06-05
+Last update: 2016-11-24
+
+************************************************************************
+* Copyright (c) 2016, Thomas Knudsen / SDFE
+*
+* Permission is hereby granted, free of charge, to any person obtaining a
+* copy of this software and associated documentation files (the "Software"),
+* to deal in the Software without restriction, including without limitation
+* the rights to use, copy, modify, merge, publish, distribute, sublicense,
+* and/or sell copies of the Software, and to permit persons to whom the
+* Software is furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included
+* in all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+* DEALINGS IN THE SOFTWARE.
+*
+***********************************************************************/
+
+#define PJ_LIB__
+#include <proj.h>
+#include <projects.h>
+#include <geocent.h>
+#include <assert.h>
+#include <stddef.h>
+#include <errno.h>
+PROJ_HEAD(helmert, "3- and 7-parameter Helmert shift");
+
+static XYZ helmert_forward_3d (LPZ lpz, PJ *P);
+static LPZ helmert_reverse_3d (XYZ xyz, PJ *P);
+
+
+
+static void *freeup_msg (PJ *P, int errlev) { /* Destructor */
+ if (0==P)
+ return 0;
+
+ if (0!=P->ctx)
+ pj_ctx_set_errno (P->ctx, errlev);
+
+ if (0==P->opaque)
+ return pj_dealloc (P);
+
+ pj_dealloc (P->opaque);
+ return pj_dealloc(P);
+}
+
+
+/* Adapts pipeline_freeup to the format defined for the PJ object */
+static void freeup (PJ *P) {
+ freeup_msg (P, 0);
+ return;
+}
+
+
+/***********************************************************************/
+struct pj_opaque_helmert {
+/************************************************************************
+ Projection specific elements for the "helmert" PJ object
+************************************************************************/
+ XYZ xyz;
+ PJ_OPK opk;
+ double scale;
+ double R[3][3];
+ int no_rotation;
+};
+
+
+/* Make the maths of the rotation operations somewhat more readable and textbook like */
+#define R00 (Q->R[0][0])
+#define R01 (Q->R[0][1])
+#define R02 (Q->R[0][2])
+
+#define R10 (Q->R[1][0])
+#define R11 (Q->R[1][1])
+#define R12 (Q->R[1][2])
+
+#define R20 (Q->R[2][0])
+#define R21 (Q->R[2][1])
+#define R22 (Q->R[2][2])
+
+
+
+
+
+
+/***********************************************************************/
+static XY helmert_forward (LP lp, PJ *P) {
+/************************************************************************
+
+ The 2D implementation is a slightly silly stop-gap, returning a
+ 2D sub-space of the 3D transformation.
+
+ In order to follow the "principle of least astonishment", the 2D
+ interface should provide an implementation of the 4 parameter
+ 2D Helmert shift.
+
+************************************************************************/
+ PJ_TRIPLET point;
+ point.lp = lp;
+ point.lpz.z = 0;
+ point.xyz = helmert_forward_3d (point.lpz, P);
+ return point.xy;
+}
+
+
+/***********************************************************************/
+static LP helmert_reverse (XY xy, PJ *P) {
+/***********************************************************************/
+ PJ_TRIPLET point;
+ point.xy = xy;
+ point.xyz.z = 0;
+ point.lpz = helmert_reverse_3d (point.xyz, P);
+ return point.lp;
+}
+
+
+/***********************************************************************/
+static XYZ helmert_forward_3d (LPZ lpz, PJ *P) {
+/***********************************************************************/
+ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
+ PJ_TRIPLET point;
+ double X, Y, Z, scale;
+ if (Q->no_rotation) {
+ point.xyz.x = lpz.lam + Q->xyz.x;
+ point.xyz.y = lpz.phi + Q->xyz.y;
+ point.xyz.z = lpz.z + Q->xyz.z;
+ return point.xyz;
+ }
+
+ scale = 1 + Q->scale * 1e-6;
+
+ X = lpz.lam;
+ Y = lpz.phi;
+ Z = lpz.z;
+
+ point.lpz = lpz;
+
+ point.xyz.x = scale * ( R00 * X + R01 * Y + R02 * Z);
+ point.xyz.y = scale * ( R10 * X + R11 * Y + R12 * Z);
+ point.xyz.z = scale * ( R20 * X + R21 * Y + R22 * Z);
+
+ point.xyz.x += Q->xyz.x;
+ point.xyz.y += Q->xyz.y;
+ point.xyz.z += Q->xyz.z;
+
+ return point.xyz;
+}
+
+
+/***********************************************************************/
+static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) {
+/***********************************************************************/
+ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
+ PJ_TRIPLET point;
+ double X, Y, Z, scale;
+
+ point.xyz = xyz;
+
+ if (Q->no_rotation) {
+ point.xyz.x = xyz.x - Q->xyz.x;
+ point.xyz.y = xyz.y - Q->xyz.y;
+ point.xyz.z = xyz.z - Q->xyz.z;
+ return point.lpz;
+ }
+
+ scale = 1 + Q->scale * 1e-6;
+
+ /* Unscale and deoffset */
+ X = (xyz.x - Q->xyz.x) / scale;
+ Y = (xyz.y - Q->xyz.y) / scale;
+ Z = (xyz.z - Q->xyz.z) / scale;
+
+ /* Inverse rotation through transpose multiplication */
+ point.xyz.x = ( R00 * X + R10 * Y + R20 * Z);
+ point.xyz.y = ( R01 * X + R11 * Y + R21 * Z);
+ point.xyz.z = ( R02 * X + R12 * Y + R22 * Z);
+
+ return point.lpz;
+}
+
+
+static PJ_OBS helmert_forward_obs (PJ_OBS point, PJ *P) {
+ point.coo.xyz = helmert_forward_3d (point.coo.lpz, P);
+ return point;
+}
+
+
+static PJ_OBS helmert_reverse_obs (PJ_OBS point, PJ *P) {
+ point.coo.lpz = helmert_reverse_3d (point.coo.xyz, P);
+ return point;
+}
+
+
+/* Milliarcsecond to radians */
+#define MAS_TO_RAD (DEG_TO_RAD / 3600000.0)
+
+
+/***********************************************************************/
+PJ *PROJECTION(helmert) {
+/***********************************************************************/
+ XYZ xyz = {0, 0, 0};
+ PJ_OPK opk = {0, 0, 0};
+ double scale = 0;
+ int approximate = 0, transpose = 0;
+
+ struct pj_opaque_helmert *Q = pj_calloc (1, sizeof (struct pj_opaque_helmert));
+ if (0==Q)
+ return freeup_msg (P, ENOMEM);
+ P->opaque = (void *) Q;
+
+ P->fwdobs = helmert_forward_obs;
+ P->invobs = helmert_reverse_obs;
+ P->fwdobs = 0;
+ P->invobs = 0;
+ P->fwd3d = helmert_forward_3d;
+ P->inv3d = helmert_reverse_3d;
+ P->fwd = helmert_forward;
+ P->inv = helmert_reverse;
+
+ P->left = PJ_IO_UNITS_METERS;
+ P->right = PJ_IO_UNITS_METERS;
+
+ /* Translations */
+ if (pj_param (P->ctx, P->params, "tx").i)
+ xyz.x = pj_param (P->ctx, P->params, "dx").f;
+
+ if (pj_param (P->ctx, P->params, "ty").i)
+ xyz.y = pj_param (P->ctx, P->params, "dy").f;
+
+ if (pj_param (P->ctx, P->params, "tz").i)
+ xyz.z = pj_param (P->ctx, P->params, "dz").f;
+
+ /* Rotations */
+ if (pj_param (P->ctx, P->params, "trx").i)
+ opk.o = pj_param (P->ctx, P->params, "drx").f * MAS_TO_RAD;
+
+ if (pj_param (P->ctx, P->params, "try").i)
+ opk.p = pj_param (P->ctx, P->params, "dry").f * MAS_TO_RAD;
+
+ if (pj_param (P->ctx, P->params, "trz").i)
+ opk.k = pj_param (P->ctx, P->params, "drz").f * MAS_TO_RAD;
+
+
+ /* Scale */
+ if (pj_param (P->ctx, P->params, "ts").i)
+ scale = pj_param (P->ctx, P->params, "ds").f;
+
+ /* Use small angle approximations? */
+ if (pj_param (P->ctx, P->params, "bapprox").i)
+ approximate = 1;
+
+ /* Use "other" rotation sign convention? */
+ if (pj_param (P->ctx, P->params, "ttranspose").i)
+ transpose = 1;
+
+ Q->xyz = xyz;
+ Q->opk = opk;
+ Q->scale = scale;
+
+ if ((opk.o==0) && (opk.p==0) && (opk.k==0) && (scale==0)) {
+ Q->no_rotation = 1;
+ return P;
+ }
+
+
+ /**************************************************************************
+
+ Build rotation matrix.
+ ----------------------
+
+ Here we rename rotation indices from omega, phi, kappa (opk), to
+ fi (i.e. phi), theta, psi (ftp), in order to reduce the mental agility
+ needed to implement the expression for the rotation matrix derived over
+ at https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions
+
+ If the option "approximate" is set, small angle approximations are used:
+ The matrix elements are approximated by expanding the trigonometric
+ functions to linear order (i.e. cos(x) = 1, sin(x) = x), and discarding
+ products of second order.
+
+ This was a useful hack when calculating by hand was the only option,
+ but in general, today, should be avoided because:
+
+ 1. It does not save much computation time, as the rotation matrix
+ is built only once, and probably used many times.
+
+ 2. The error induced may be too large for ultra high accuracy
+ applications: the Earth is huge and the linear error is
+ approximately the angular error multiplied by the Earth radius.
+
+ However, in many cases the approximation is necessary, since it has
+ been used historically: Rotation angles from older published datum
+ shifts may actually be a least squares fit to the linearized rotation
+ approximation, hence not being strictly valid for deriving the full
+ rotation matrix.
+
+ So in order to fit historically derived coordinates, the access to
+ the approximate rotation matrix is necessary - at least in principle.
+
+ Also, when using any published datum transformation information, one
+ should always check which convention (exact or approximate rotation
+ matrix) is expected, and whether the induced error for selecting
+ the opposite convention is acceptable (which it often is).
+
+
+ Sign conventions
+ ----------------
+
+ Take care: Two different sign conventions exist.
+
+ Conceptually they relate to whether we rotate the coordinate system
+ or the "position vector" (the vector going from the coordinate system
+ origin to the point being transformed, i.e. the point coordinates
+ interpreted as vector coordinates).
+
+ Switching between the "position vector" and "coordinate system"
+ conventions is simply a matter of switching the sign of the rotation
+ angles, which algebraically also translates into a transposition of
+ the rotation matrix.
+
+ Hence, as geodetic constants should preferably be referred to exactly
+ as published, the "transpose" option provides the ability to switch
+ between the conventions.
+
+ ***************************************************************************/
+ do {
+ double f, t, p; /* phi/fi , theta, psi */
+ double cf, ct, cp; /* cos (fi, theta, psi) */
+ double sf, st, sp; /* sin (fi, theta, psi) */
+
+ /* rename (omega, phi, kappa) to (fi, theta, psi) */
+ f = opk.o;
+ t = opk.p;
+ p = opk.k;
+
+ if (approximate) {
+ R00 = 1;
+ R01 = p;
+ R02 = -t;
+
+ R10 = -p;
+ R11 = 1;
+ R12 = f;
+
+ R20 = t;
+ R21 = -f;
+ R22 = 1;
+ break;
+ }
+
+ cf = cos(f);
+ sf = sin(f);
+ ct = cos(t);
+ st = sin(t);
+ cp = cos(p);
+ sp = sin(p);
+
+
+ R00 = ct*cp;
+ R01 = cf*sp + sf*st*cp;
+ R02 = sf*sp - cf*st*cp;
+
+ R10 = -ct*sp;
+ R11 = cf*cp - sf*st*sp;
+ R12 = sf*cp + cf*st*sp;
+
+ R20 = st;
+ R21 = -sf*ct;
+ R22 = cf*ct;
+
+ /*
+ For comparison: Description from Engsager/Poder implementation
+ in set_dtm_1.c (trlib)
+
+ DATUM SHIFT:
+ TO = scale * ROTZ * ROTY * ROTX * FROM + TRANSLA
+
+ ( cz sz 0) (cy 0 -sy) (1 0 0)
+ ROTZ=(-sz cz 0), ROTY=(0 1 0), ROTX=(0 cx sx)
+ ( 0 0 1) (sy 0 cy) (0 -sx cx)
+
+ trp->r11 = cos_ry*cos_rz;
+ trp->r12 = cos_rx*sin_rz + sin_rx*sin_ry*cos_rz;
+ trp->r13 = sin_rx*sin_rz - cos_rx*sin_ry*cos_rz;
+
+ trp->r21 = -cos_ry*sin_rz;
+ trp->r22 = cos_rx*cos_rz - sin_rx*sin_ry*sin_rz;
+ trp->r23 = sin_rx*cos_rz + cos_rx*sin_ry*sin_rz;
+
+ trp->r31 = sin_ry;
+ trp->r32 = -sin_rx*cos_ry;
+ trp->r33 = cos_rx*cos_ry;
+
+ trp->scale = 1.0 + scale;
+ */
+
+ } while (0);
+
+
+ if (transpose) {
+ double r;
+ r = R01; R01 = R10; R10 = r;
+ r = R02; R02 = R20; R20 = r;
+ r = R12; R12 = R21; R21 = r;
+ }
+
+ return P;
+}
+
+
+
+
+#ifndef PJ_SELFTEST
+
+int pj_helmert_selftest (void) {return 0;} /* self test stub */
+
+#else
+
+
+static int test (char *args, PJ_TRIPLET in, PJ_TRIPLET expect) {
+ PJ_TRIPLET out;
+ PJ *P = pj_init_plus (args);
+
+ if (0==P)
+ return 5;
+
+ out.xyz = pj_fwd3d (in.lpz, P);
+ if (pj_xyz_dist (out.xyz, expect.xyz) > 1e-4)
+ return 1;
+
+ out.lpz = pj_inv3d (out.xyz, P);
+ if (pj_xyz_dist (out.xyz, in.xyz) > 1e-4)
+ return 2;
+
+ return 0;
+}
+
+
+
+int pj_helmert_selftest (void) {
+
+ /* This example is from
+ Lotti Jivall:
+ Simplified transformations from ITRF2008/IGS08 to ETRS89 for maritime applications */
+ PJ_TRIPLET in1 = {{3565285.0000, 855949.0000, 5201383.0000}};
+ PJ_TRIPLET expect1 = {{3565285.41342351, 855948.67986759, 5201382.72939791}};
+ char args1[] = {
+ " +proj=helmert +ellps=GRS80"
+ " +x=0.67678 +y=0.65495 +z=-0.52827"
+ " +rx=-22.742 +ry=12.667 +rz=22.704 +s=-0.01070"
+ };
+
+
+ /* This example is a random point, transformed from ED50 to ETRS89 using KMStrans2 */
+ PJ_TRIPLET in2 = {{3494994.3012, 1056601.9725, 5212382.1666}};
+ PJ_TRIPLET expect2 = {{3494909.84026368, 1056506.78938633, 5212265.66699761}};
+ char args2[] = {
+ " +proj=helmert +ellps=GRS80"
+ " +x=-81.0703 +y=-89.3603 +z=-115.7526"
+ " +rx=-484.88 +ry=-24.36 +rz=-413.21 +s=-0.540645"
+ };
+ int ret;
+
+ ret = test (args1, in1, expect1); if (ret) return ret;
+ ret = test (args2, in2, expect2); if (ret) return ret + 10;
+ return 0;
+}
+
+#endif
diff --git a/src/PJ_horner.c b/src/PJ_horner.c
new file mode 100644
index 00000000..f39d2c2c
--- /dev/null
+++ b/src/PJ_horner.c
@@ -0,0 +1,953 @@
+#define PJ_LIB__
+#include <proj.h>
+#include <projects.h>
+#include <assert.h>
+#include <stddef.h>
+#include <math.h>
+#include <errno.h>
+PROJ_HEAD(horner, "Horner polynomial evaluation");
+#define HORNER_SILENCE
+
+/* The next few hundred lines comprises a direct cut-and-paste from the horner.h header library */
+
+/***********************************************************************
+
+ Interfacing to a classic piece of geodetic software
+
+************************************************************************
+
+ gen_pol is a highly efficient, classic implementation of a generic
+ 2D Horner's Scheme polynomial evaluation routine by Knud Poder and
+ Karsten Engsager, originating in the vivid geodetic environment at
+ what was then (1960-ish) the Danish Geodetic Institute.
+
+ The original Poder/Engsager gen_pol implementation (where
+ the polynomial degree and two sets of polynomial coefficients
+ are packed together in one compound array, handled via a simple
+ double pointer) is compelling and "true to the code history":
+
+ It has a beautiful classical 1960s ring to it, not unlike the
+ original fft implementations, which revolutionized spectral
+ analysis in twenty lines of code.
+
+ The Poder coding sound, as classic 1960s as Phil Spector's Wall
+ of Sound, is beautiful and inimitable.
+
+ On the other hand: For the uninitiated, the gen_pol code is hard
+ to follow, despite being compact.
+
+ Also, since adding metadata and improving maintainability
+ of the code are among the implied goals of a current SDFE/DTU Space
+ project, the material in this file introduces a version with a more
+ more modern (or at least 1990s) look, introducing a "double 2D
+ polynomial" data type, HORNER.
+
+ Despite introducing a new data type for handling the polynomial
+ coefficients, great care has been taken to keep the coefficient
+ array organization identical to that of gen_pol.
+
+ Hence, on one hand, the HORNER data type helps improving the
+ long term maintainability of the code by making the data
+ organization more mentally accessible.
+
+ On the other hand, it allows us to preserve the business end of
+ the original gen_pol implementation - although not including the
+ famous "Poder dual autocheck" in all its enigmatic elegance.
+
+ The original code has, however, been included in the conditionally
+ compiled TEST-section.
+
+ This is partially for validation of the revised version, partially
+ to enable more generic enjoyment of an interesting piece of
+ ingenious geodetic code - simplistic and enigmatic at the same time.
+
+ **********************************************************************
+
+ The material included here was written by Knud Poder, starting
+ around 1960, and Karsten Engsager, starting around 1970. It was
+ originally written in Algol 60, later (1980s) reimplemented in C.
+
+ The HORNER data type interface, and the organization as a header
+ library was implemented by Thomas Knudsen, starting around 2015.
+
+ ***********************************************************************
+ *
+ * The gen_pol routine comes with this legal statement (ISC/OpenBSD):
+ *
+ * Copyright (c) 2011, National Survey and Cadastre, Denmark
+ * (Kort- og Matrikelstyrelsen), kms@kms.dk
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ *
+ *
+ *******************************************************************************
+ *
+ * The remaining parts are...
+ *
+ * Copyright (c) 2016, Thomas Knudsen / Karsten Engsager / SDFE http://www.sdfe.dk
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *
+ *****************************************************************************/
+
+
+ #ifndef HORNER_H
+ #define HORNER_H
+ #ifdef __cplusplus
+ extern "C" {
+ #endif
+
+
+
+#if defined(PROJ_H) || defined(PROJECTS_H)
+#define horner_dealloc(x) pj_dealloc(x)
+#define horner_calloc(n,x) pj_calloc(n,x)
+#else
+#define horner_dealloc(x) free(x)
+#define horner_calloc(n,x) calloc(n,x)
+typedef struct {double u,v;} UV;
+#endif
+
+struct horner;
+typedef struct horner HORNER;
+static UV horner (const HORNER *transformation, int direction, UV position);
+static HORNER *horner_alloc (size_t order);
+static void horner_free (HORNER *h);
+
+struct horner {
+ int order; /* maximum degree of polynomium */
+ int coefs; /* number of coefficients for each polynomium */
+ double range; /* radius of the region of validity */
+
+ double *fwd_u; /* coefficients for the forward transformations */
+ double *fwd_v; /* i.e. latitude/longitude to northing/easting */
+
+ double *inv_u; /* coefficients for the inverse transformations */
+ double *inv_v; /* i.e. northing/easting to latitude/longitude */
+
+ UV *fwd_origin; /* False longitude/latitude */
+ UV *inv_origin; /* False easting/northing */
+};
+
+/* e.g. degree = 2: a + bx + cy + dxx + eyy + fxy, i.e. 6 coefficients */
+#define horner_number_of_coefficients(order) \
+ (((order + 1)*(order + 2)/2))
+
+static int horner_degree_u (int order, int index);
+static int horner_degree_v (int order, int index);
+static int horner_index (int order, int degree_u, int degree_v);
+
+#ifndef HORNER_HEADER_ONLY
+
+/***************************************************************************/
+static int horner_index (int order, int degree_1, int degree_2) {
+/****************************************************************************
+
+ Returns the index of the polynomial coefficient, C, for the element
+
+ C * pow (c_1, degree_2) * pow (c_2, degree_2),
+
+ given that degree_1 > -1, degree_2 > -1, degree_1 + degree_2 <= order.
+
+ Otherwise returns -1 and sets errno to EDOM.
+
+ The range of the index is [0 : (order + 1) * (order + 2) / 2 - 1].
+
+ A very important thing to note is that the order of the coordinates
+ c_1 and c_2 depend on the polynomium:
+
+ For the fwd and inv polynomia for the "u" coordinate,
+ u is first (degree_1), v is second (degree_2).
+ For the fwd and inv polynomia for the "v" coordinate,
+ v is first (degree_1), u is second (degree_2).
+
+****************************************************************************/
+ if ( (degree_1 < 0) || (degree_2 < 0) || (degree_1 + degree_2 > order) ) {
+ errno = EDOM;
+ return -1;
+ }
+
+ return ( horner_number_of_coefficients(order) - 1
+ - (order - degree_1)*(order - degree_1 + 1)/2
+ - (order - degree_1 - degree_2));
+}
+
+#define index_u(h, u, v) horner_index (h->order, u, v)
+#define index_v(h, u, v) horner_index (h->order, v, u)
+
+
+static int horner_degree_u (int order, int index) {
+ int n = horner_number_of_coefficients(order);
+ int i, j;
+ if ((order < 0) || (index >= n)) {
+ errno = EDOM;
+ return -1;
+ }
+ for (i = 0; i <= order; i++)
+ for (j = 0; j <= order - i; j++)
+ if (index == horner_index (order, i, j))
+ return i;
+ return -1;
+}
+
+
+static int horner_degree_v (int order, int index) {
+ int n = horner_number_of_coefficients(order);
+ int i, j;
+ if ((order < 0) || (index >= n)) {
+ errno = EDOM;
+ return -1;
+ }
+ for (i = 0; i <= order; i++)
+ for (j = 0; j <= order - i; j++)
+ if (index == horner_index (order, i, j))
+ return j;
+ return -1;
+}
+
+
+
+
+static void horner_free (HORNER *h) {
+ horner_dealloc (h->inv_v);
+ horner_dealloc (h->inv_u);
+ horner_dealloc (h->fwd_v);
+ horner_dealloc (h->fwd_u);
+ horner_dealloc (h);
+}
+
+
+static HORNER *horner_alloc (size_t order) {
+ /* size_t is unsigned, so we need not check for order > 0 */
+ int n = horner_number_of_coefficients(order);
+
+ HORNER *h = horner_calloc (1, sizeof (HORNER));
+
+ if (0==h)
+ return 0;
+
+ h->order = order;
+ h->coefs = n;
+
+ h->fwd_u = horner_calloc (n, sizeof(double));
+ h->fwd_v = horner_calloc (n, sizeof(double));
+ h->inv_u = horner_calloc (n, sizeof(double));
+ h->inv_v = horner_calloc (n, sizeof(double));
+
+ h->fwd_origin = horner_calloc (1, sizeof(UV));
+ h->inv_origin = horner_calloc (1, sizeof(UV));
+
+ if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v && h->fwd_origin && h->inv_origin)
+ return h;
+
+ /* safe, since all pointers are null-initialized (by calloc) */
+ horner_free (h);
+ return 0;
+}
+
+
+
+
+/**********************************************************************/
+static UV horner (const HORNER *transformation, int direction, UV position) {
+/***********************************************************************
+
+A reimplementation of the classic Engsager/Poder 2D Horner polynomial
+evaluation engine "gen_pol".
+
+This version omits the inimitable Poder "dual autocheck"-machinery,
+which here is intended to be implemented at a higher level of the
+library: We separate the polynomial evaluation from the quality
+control (which, given the limited MTBF for "computing machinery",
+typical when Knud Poder invented the dual autocheck method,
+was not defensible at that time).
+
+Another difference from the original version is that we return the
+result on the stack, rather than accepting pointers to result variables
+as input. This results in code that is easy to read:
+
+ projected = horner (s34j, 1, geographic);
+ geographic = horner (s34j, -1, projected );
+
+and experiments have shown that on contemporary architectures, the time
+taken for returning even comparatively large objects on the stack (and
+the UV is not that large - typically only 16 bytes) is negligibly
+different from passing two pointers (i.e. typically also 16 bytes) the
+other way.
+
+The polynomium has the form:
+
+P = sum (i = [0 : order])
+ sum (j = [0 : order - i])
+ pow(par_1, i) * pow(par_2, j) * coef(index(order, i, j))
+
+For numerical stability, the summation is carried out backwards,
+summing the tiny high order elements first.
+
+***********************************************************************/
+
+ /* These variable names follow the Engsager/Poder implementation */
+ int sz; /* Number of coefficients per polynomial */
+ double *tcx, *tcy; /* Coefficient pointers */
+ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */
+ double n, e;
+ UV uv_error;
+ uv_error.u = uv_error.v = HUGE_VAL;
+
+ if (0==transformation)
+ return uv_error;
+
+ /* Check for valid value of direction (-1, 0, 1) */
+ switch (direction) {
+ case 0: /* no-op */
+ return position;
+ case 1: /* forward */
+ case -1: /* inverse */
+ break;
+ default: /* invalid */
+ errno = EINVAL;
+ return uv_error;
+ }
+
+ /* Prepare for double Horner */
+ sz = horner_number_of_coefficients(transformation->order);
+ range = transformation->range;
+
+
+ if (direction==1) { /* forward */
+ tcx = transformation->fwd_u + sz;
+ tcy = transformation->fwd_v + sz;
+ e = position.u - transformation->fwd_origin->u;
+ n = position.v - transformation->fwd_origin->v;
+ } else { /* inverse */
+ tcx = transformation->inv_u + sz;
+ tcy = transformation->inv_v + sz;
+ e = position.u - transformation->inv_origin->u;
+ n = position.v - transformation->inv_origin->v;
+ }
+
+ if ((fabs(n) > range) || (fabs(e) > range)) {
+ errno = EDOM;
+ return uv_error;
+ }
+
+ /* The melody of this block is straight out of the great Engsager/Poder songbook */
+ else {
+ int g = transformation->order;
+ int r = g, c;
+ double u, v, N, E;
+
+ /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */
+ for (N = *--tcy, E = *--tcx; r > 0; r--) {
+ for (c = g, u = *--tcy, v = *--tcx; c >= r; c--) {
+ u = n*u + *--tcy;
+ v = e*v + *--tcx;
+ }
+ N = e*N + u;
+ E = n*E + v;
+ }
+
+ position.u = E;
+ position.v = N;
+ }
+
+ return position;
+}
+
+
+#ifdef HORNER_SILENCE
+/**********************************************************************/
+static int horner_silence (int i) {
+/***********************************************************************
+ useless function that silences coompiler warnings about unused stuff
+***********************************************************************/
+ HORNER *Q;
+ UV uv_error;
+ if (i==0)
+ return i;
+ uv_error.u = uv_error.v = HUGE_VAL;
+ horner(0, 1, uv_error);
+ Q = horner_alloc (2);
+ if (Q)
+ horner_free (Q);
+ if (horner_degree_u (2,1))
+ return horner_degree_v (2,1);
+ return i;
+}
+#endif /* def HORNER_SILENCE */
+
+
+
+
+#ifdef HORNER_TEST_ORIGINAL_GEN_POL_CODE
+
+double fwd_u[] = {1,2,3};
+double fwd_v[] = {4,5,6};
+double inv_u[] = {4,6,5};
+double inv_v[] = {1,3,2};
+UV uv_origin_fwd = {0, 0};
+UV uv_origin_inv = {0, 0};
+HORNER uv = {1, 3, 500000.0, fwd_u, fwd_v, inv_u, inv_v, &uv_origin_fwd, &uv_origin_inv};
+
+void show_test (void);
+void tuut_b_test (void);
+
+
+int main (int argc, char **argv) {
+ int i, j, order = 10;
+ UV res = {1,1};
+
+ if (argc==2)
+ order = atoi (argv[1]);
+
+ printf ("Testing %d combinations\n", horner_number_of_coefficients(order));
+
+ for (i = 0; i <= order; i++)
+ for (j = 0; j <= order - i; j++) {
+ int hh = horner_index (order, i, j);
+ int ii = horner_degree_u (order, hh);
+ int jj = horner_degree_v (order, hh);
+ assert (i==ii);
+ assert (j==jj);
+ printf ("%2.2d %1d%1du %1d%1dv\n", hh, i, ii, j, jj);
+ }
+
+ tuut_b_test ();
+
+ puts ("Forward...");
+
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, 1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v);
+ assert ( 6==res.u); /* fwd_u: a + bu + cv = 1 + 2*1 + 3*1 = 6 */
+ assert (15==res.v); /* fwd_v: a + bu + cv = 4 + 5*1 + 6*1 = 15 */
+
+ res.u = 2;
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, 1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v); /* u = 2, v = 15 */
+ assert (50==res.u); /* fwd_u: a + bu + cv = 1 + 2*(u=2) + 3*(v=15) = 50 */
+ assert (91==res.v); /* fwd_v: a + bu + cv = 4 + 5*(v=15) + 6*(u=2) = 91 */
+
+ puts ("Inverse...");
+
+ res.u = 1;
+ res.v = 1;
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, -1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v); /* u = 1, v = 1 */
+ assert (15==res.u); /* inv_u: a + bu + cv = 4 + 5*1 + 6*1 = 15 */
+ assert ( 6==res.v); /* inv_v: a + bu + cv = 1 + 2*1 + 3*1 = 6 */
+
+ res.u = 2;
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, -1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v); /* u = 2, v = 6 */
+ assert (46==res.u); /* inv_u: a + bu + cv = 4 + 6*(u=2) + 5*(v=6) = 46 */
+ assert (23==res.v); /* inv_v: a + bu + cv = 1 + 2*(u=2) + 3*(v=6) = 23 */
+
+/* show_test ();*/
+ return 0;
+}
+
+
+
+
+
+
+/* ttu_n and ttu_e are based on static double C_ttu_b[] */
+static double ttu_n[] = {
+ /* tc32_ed50 -> utm32_ed50 : Bornholm */
+ /* m_lim_gen: 0.086 red = 0 OBS = 852 */
+ /* m = 1.38 cm my_loss = +2 y_enp = +10.5 */
+ /* m = 1.44 cm mx_loss = +2 x_enp = +10.4 */
+
+ /*deg 4.0,*/
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258112678e+06, 9.9999971567e-01, 1.5372750011e-10,
+ /* 3*/ 5.9300860915e-15, 2.2609497633e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ 4.3188227445e-05, 2.8225130416e-10, 7.8740007114e-16,
+ /* 8*/ -1.7453997279e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ 1.6877465415e-10, -1.1234649773e-14, -1.7042333358e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ -7.9303467953e-15, -5.2906832535e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ 3.9984284847e-19,
+
+ /*tcy 6125810.306769, */
+};
+
+static double ttu_e[] = {
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760574982e+05, 9.9999752475e-01, 2.8817299305e-10,
+ /* 3*/ 5.5641310680e-15, -1.5544700949e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ -4.1357045890e-05, 4.2106213519e-11, 2.8525551629e-14,
+ /* 8*/ -1.9107771273e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ 3.3615590093e-10, 2.4380247154e-14, -2.0241230315e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ 1.2429019719e-15, 5.3886155968e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ -1.0167505000e-18,
+
+ /* tcx 877605.269066 */
+ };
+
+
+/* utt_n and utt_e are based on static double C_utt_b[] */
+static double utt_n[] = {
+ /* utm32_ed50 -> tc32_ed50 : Bornholm */
+ /* m_lim_gen: 0.086 red = 0 OBS = 852 */
+ /* m = 1.38 cm my_loss = +2 y_enp = +10.8 */
+ /* m = 1.44 cm mx_loss = +2 x_enp = +10.7 */
+
+ /*deg 4.0,*/
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258103208e+06, 1.0000002826e+00, -1.5372762184e-10,
+ /* 3*/ -5.9304261011e-15, -2.2612705361e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ -4.3188331419e-05, -2.8225549995e-10, -7.8529116371e-16,
+ /* 8*/ 1.7476576773e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ -1.6875687989e-10, 1.1236475299e-14, 1.7042518057e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ 7.9300735257e-15, 5.2881862699e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ -3.9990736798e-19,
+
+ /*tcy 6125811.281773,*/
+};
+
+static double utt_e[] = {
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760527928e+05, 1.0000024735e+00, -2.8817540032e-10,
+ /* 3*/ -5.5627059451e-15, 1.5543637570e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ 4.1357152105e-05, -4.2114813612e-11, -2.8523713454e-14,
+ /* 8*/ 1.9109017837e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ -3.3616407783e-10, -2.4382678126e-14, 2.0245020199e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ -1.2441377565e-15, -5.3885232238e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ 1.0167203661e-18,
+
+ /*tcx 877605.760036 */
+};
+
+UV tuut_b_origin_fwd = {877605.269066, 6125810.306769};
+UV tuut_b_origin_inv = {877605.760036, 6125811.281773};
+HORNER tuut_b = {4, 15, 500000.0, ttu_e, ttu_n, utt_e, utt_n, &tuut_b_origin_fwd, &tuut_b_origin_inv};
+
+
+/* Prototype and forward declarations of the material needed for cross-check with original implementation */
+int gen_pol(double *C_f, double *C_r, double N_in, double E_in, double *Nout, double *Eout);
+static double C_ttu_b[];
+static double C_utt_b[];
+
+
+void gen_pol_roundtrip (double *C_ttu_b, double *C_utt_b, UV fwd) {
+ UV res, hrn;
+ int ret;
+
+ ret = gen_pol(C_ttu_b, C_utt_b, fwd.v, fwd.u, &res.v, &res.u);
+ printf ("\n------\n\n");
+ if (0!=ret) printf ("ret: %d\n", ret);
+ printf ("inp: %11.3f %11.3f\n", fwd.u, fwd.v);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ hrn = horner (&tuut_b, 1, fwd);
+ printf ("hrn: %11.3f %11.3f\n", hrn.u, hrn.v);
+ assert (hrn.u==res.u);
+ assert (hrn.v==res.v);
+
+ ret = gen_pol(C_utt_b, C_ttu_b, res.v, res.u, &res.v, &res.u);
+ hrn = horner (&tuut_b, -1, hrn);
+ printf ("hrn: %11.3f %11.3f\n", hrn.u, hrn.v);
+ if (0!=ret) printf ("ret: %d\n", ret);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ assert (hrn.u==res.u);
+ assert (hrn.v==res.v);
+ printf ("inp: %11.3f %11.3f (%.3g mm)\n", fwd.u, fwd.v, 1e3*hypot(fwd.u-res.u, fwd.v-res.v));
+}
+
+void tuut_b_test (void) {
+ UV fwd = tuut_b_origin_fwd;
+ UV res;
+
+ puts ("Bornholm");
+ printf ("fwd: %11.3f %11.3f\n", fwd.u, fwd.v);
+ res = horner (&tuut_b, 1, fwd);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ res = horner (&tuut_b, -1, res);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ printf ("fwd: %11.3f %11.3f\n", fwd.u, fwd.v);
+
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+ fwd.u = 877000;
+ fwd.v = 6125000;
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+ fwd.u = 800000;
+ fwd.v = 6100000;
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+ fwd.u = 850000;
+ fwd.v = 6200000;
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+
+}
+
+
+
+
+
+#define FLOATLIMIT 5.0e5
+#define TRF_AREA_ EDOM
+
+/********************************************************************************************/
+int gen_pol(double *C_f, double *C_r, double N_in, double E_in, double *Nout, double *Eout) {
+/*********************************************************************************************
+
+ This is the original Poder/Engsager implementation of gen_pol.
+ It is included here for test-comparison with the horner() routine.
+
+*********************************************************************************************/
+ double N, E, n, e;
+ double *Cp, *tcy, *tcx;
+ double tol = 1.0e-4;
+ int i;
+ int g;
+ int sz;
+ int r, c;
+ int res = 0;
+
+ /* Preserve input for reverse check */
+ N = N_in;
+ E = E_in;
+ Cp = C_f;
+
+ /* Transformation loop */
+ for (i = -1; i <= 1 && res == 0; ++i)
+ if (i) {
+
+ /* Prepare for double Horner */
+ g = (int) *Cp;
+ sz = (g + 1)*(g + 2)/2 + 1;
+ tcy = Cp + sz;
+ tcx = tcy + sz;
+ /* Double Horner's scheme */
+ /* N = n*Cy*e -> yout, E = e*Cx*n -> xout */
+ n = N - *tcy;
+ e = E - *tcx;
+ if ((fabs(n) < FLOATLIMIT) && (fabs(e) < FLOATLIMIT)) {
+
+ for ( r = g, N = *--tcy, E = *--tcx; r > 0; r--) {
+ double u, v;
+ for (c = g, u = *--tcy, v = *--tcx; c >= r; c--) {
+ u = n*u + *--tcy;
+ v = e*v + *--tcx;
+ }
+ N = e*N + u;
+ E = n*E + v;
+ }
+ } else res = TRF_AREA_;
+ }
+ else { /* collect output coord,switch to reverse checking */
+ *Nout = N;
+ *Eout = E;
+ Cp = C_r;
+ }
+
+ /* tol-check of results*/
+ if (res == 0 && (fabs(N - N_in) < tol && fabs(E - E_in) < tol))
+ return (0);
+ else if (res == 0) res = TRF_AREA_;
+ return(res);
+
+#undef FLOATLIMIT
+
+}
+
+
+/* s45b polynomia */
+ static double C_ttu_b[] = {
+ /* tc32_ed50 -> utm32_ed50 : Bornholm */
+ /* m_lim_gen: 0.086 red = 0 OBS = 852 */
+ /* m = 1.38 cm my_loss = +2 y_enp = +10.5 */
+ /* m = 1.44 cm mx_loss = +2 x_enp = +10.4 */
+
+ /*deg*/ 4.0,
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258112678e+06, 9.9999971567e-01, 1.5372750011e-10,
+ /* 3*/ 5.9300860915e-15, 2.2609497633e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ 4.3188227445e-05, 2.8225130416e-10, 7.8740007114e-16,
+ /* 8*/ -1.7453997279e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ 1.6877465415e-10, -1.1234649773e-14, -1.7042333358e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ -7.9303467953e-15, -5.2906832535e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ 3.9984284847e-19,
+
+ /*tcy*/ 6125810.306769,
+
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760574982e+05, 9.9999752475e-01, 2.8817299305e-10,
+ /* 3*/ 5.5641310680e-15, -1.5544700949e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ -4.1357045890e-05, 4.2106213519e-11, 2.8525551629e-14,
+ /* 8*/ -1.9107771273e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ 3.3615590093e-10, 2.4380247154e-14, -2.0241230315e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ 1.2429019719e-15, 5.3886155968e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ -1.0167505000e-18,
+
+ /*tcx*/ 877605.269066
+ };
+
+ static double C_utt_b[] = {
+ /* utm32_ed50 -> tc32_ed50 : Bornholm */
+ /* m_lim_gen: 0.086 red = 0 OBS = 852 */
+ /* m = 1.38 cm my_loss = +2 y_enp = +10.8 */
+ /* m = 1.44 cm mx_loss = +2 x_enp = +10.7 */
+
+ /*deg*/ 4.0,
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258103208e+06, 1.0000002826e+00, -1.5372762184e-10,
+ /* 3*/ -5.9304261011e-15, -2.2612705361e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ -4.3188331419e-05, -2.8225549995e-10, -7.8529116371e-16,
+ /* 8*/ 1.7476576773e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ -1.6875687989e-10, 1.1236475299e-14, 1.7042518057e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ 7.9300735257e-15, 5.2881862699e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ -3.9990736798e-19,
+
+ /*tcy*/ 6125811.281773,
+
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760527928e+05, 1.0000024735e+00, -2.8817540032e-10,
+ /* 3*/ -5.5627059451e-15, 1.5543637570e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ 4.1357152105e-05, -4.2114813612e-11, -2.8523713454e-14,
+ /* 8*/ 1.9109017837e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ -3.3616407783e-10, -2.4382678126e-14, 2.0245020199e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ -1.2441377565e-15, -5.3885232238e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ 1.0167203661e-18,
+
+ /*tcx*/ 877605.760036
+ };
+
+#endif /* HORNER_TEST_ORIGINAL_GEN_POL_CODE */
+
+
+#endif /* ndef HORNER_HEADER_ONLY */
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* ndef HORNER_H */
+
+
+
+
+
+
+
+static PJ_OBS horner_forward_obs (PJ_OBS point, PJ *P) {
+ point.coo.uv = horner ((HORNER *) P->opaque, 1, point.coo.uv);
+ return point;
+}
+
+static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) {
+ point.coo.uv = horner ((HORNER *) P->opaque, -1, point.coo.uv);
+ return point;
+}
+
+
+
+
+
+
+static void *horner_freeup (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ if (0==P->opaque)
+ return pj_dealloc (P);
+ horner_free ((HORNER *) P->opaque);
+ return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+ horner_freeup (P);
+ return;
+}
+
+
+static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
+ char buf[20], *init, *next;
+ int i;
+ sprintf (buf, "t%s", param);
+ if (0==pj_param (P->ctx, P->params, buf).i)
+ return 0;
+ sprintf (buf, "s%s", param);
+ init = pj_param(P->ctx, P->params, buf).s;
+
+ for (i = 0; i < ncoefs; i++) {
+ if (i > 0) {
+ if (','!=*next) {
+ pj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs);
+ return 0;
+ }
+ init = ++next;
+ }
+ coefs[i] = pj_strtod (init, &next);
+ }
+ return 1;
+}
+
+
+/*********************************************************************/
+PJ *PROJECTION(horner) {
+/*********************************************************************/
+ int degree = 0, n;
+ HORNER *Q;
+ P->fwdobs = horner_forward_obs;
+ P->invobs = horner_reverse_obs;
+ P->fwd3d = 0;
+ P->inv3d = 0;
+ P->fwd = 0;
+ P->inv = 0;
+ P->left = P->right = PJ_IO_UNITS_METERS;
+ /* silence a few compiler warnings */
+ horner_silence (0);
+
+ /* Polynomial degree specified? */
+ if (pj_param (P->ctx, P->params, "tdeg").i) /* degree specified? */
+ degree = pj_param(P->ctx, P->params, "ideg").i;
+ else {
+ pj_log_debug (P, "Horner: Need to specify polynomial degree, (+deg=n)");
+ return horner_freeup (P);
+ }
+
+ Q = horner_alloc (degree);
+ P->opaque = (void *) Q;
+
+ n = horner_number_of_coefficients (degree);
+
+ if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, Q->inv_u, "inv_u", n))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, Q->inv_v, "inv_v", n))
+ return horner_freeup (P);
+
+ if (0==parse_coefs (P, (double *)(Q->fwd_origin), "fwd_origin", 2))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, &Q->range, "range", 1))
+ Q->range = 500000;
+
+ return P;
+}
+
+#ifndef PJ_SELFTEST
+/* selftest stub */
+int pj_horner_selftest (void) {return 0;}
+#else
+
+char tc32_utm32[] = {
+ " +proj=horner"
+ " +ellps=intl"
+ " +range=500000"
+ " +fwd_origin=877605.269066,6125810.306769"
+ " +inv_origin=877605.760036,6125811.281773"
+ " +deg=4"
+ " +fwd_v=6.1258112678e+06,9.9999971567e-01,1.5372750011e-10,5.9300860915e-15,2.2609497633e-19,4.3188227445e-05,2.8225130416e-10,7.8740007114e-16,-1.7453997279e-19,1.6877465415e-10,-1.1234649773e-14,-1.7042333358e-18,-7.9303467953e-15,-5.2906832535e-19,3.9984284847e-19"
+ " +fwd_u=8.7760574982e+05,9.9999752475e-01,2.8817299305e-10,5.5641310680e-15,-1.5544700949e-18,-4.1357045890e-05,4.2106213519e-11,2.8525551629e-14,-1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1.2429019719e-15,5.3886155968e-19,-1.0167505000e-18"
+ " +inv_v=6.1258103208e+06,1.0000002826e+00,-1.5372762184e-10,-5.9304261011e-15,-2.2612705361e-19,-4.3188331419e-05,-2.8225549995e-10,-7.8529116371e-16,1.7476576773e-19,-1.6875687989e-10,1.1236475299e-14,1.7042518057e-18,7.9300735257e-15,5.2881862699e-19,-3.9990736798e-19"
+ " +inv_u=8.7760527928e+05,1.0000024735e+00,-2.8817540032e-10,-5.5627059451e-15,1.5543637570e-18,4.1357152105e-05,-4.2114813612e-11,-2.8523713454e-14,1.9109017837e-18,-3.3616407783e-10,-2.4382678126e-14,2.0245020199e-18,-1.2441377565e-15,-5.3885232238e-19,1.0167203661e-18"
+};
+
+
+int pj_horner_selftest (void) {
+ PJ *P;
+ PJ_OBS a, b;
+ double dist;
+
+ /* The TC32 for "System 45 Bornholm" */
+ /* TC32 -> UTM32" */
+ P = pj_create (tc32_utm32);
+ if (0==P)
+ return 10;
+
+ a = b = pj_obs_null;
+ a.coo.uv.v = 6125305.4245;
+ a.coo.uv.u = 878354.8539;
+
+ /* Forward projection */
+ b = pj_trans (P, PJ_FWD, a);
+ b = pj_trans (P, PJ_INV, b);
+
+ /* Check roundtrip precision for 1 iteration each way */
+ dist = pj_roundtrip (P, PJ_FWD, 1, a);
+ if (dist > 0.01)
+ return 1;
+ return 0;
+}
+#endif
diff --git a/src/PJ_krovak.c b/src/PJ_krovak.c
index 136978a0..c469346f 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -239,7 +239,7 @@ int pj_krovak_selftest (void) {
double tolerance_lp = 1e-10;
double tolerance_xy = 1e-7;
- char e_args[] = {"+proj=krovak +ellps=GRS80"};
+ char e_args[] = {"+proj=krovak +ellps=GRS80 +no_defs"};
LP fwd_in[] = {
{ 2, 1},
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index ef86ccd3..27fe6ec9 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -114,6 +114,8 @@ struct pj_opaque {
int *reverse_step;
int *omit_forward;
int *omit_inverse;
+ char **argv;
+ char **current_argv;
PJ_OBS stack[PIPELINE_STACK_SIZE];
PJ **pipeline;
};
@@ -175,13 +177,17 @@ static LP pipeline_reverse (XY xyz, PJ *P);
static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
- int i, first_step, last_step, incr;
+ int i, first_step, last_step;
first_step = 1;
last_step = P->opaque->steps + 1;
- incr = 1;
- for (i = first_step; i != last_step; i += incr) {
+ for (i = first_step; i != last_step; i++) {
+ pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s",
+ i-first_step, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
+ P->opaque->pipeline[i]->descr
+ );
+
if (P->opaque->omit_forward[i])
continue;
if (P->opaque->reverse_step[i])
@@ -191,19 +197,23 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
if (P->opaque->depth < PIPELINE_STACK_SIZE)
P->opaque->stack[P->opaque->depth++] = point;
}
+ pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
P->opaque->depth = 0; /* Clear the stack */
return point;
}
+
static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
- int i, first_step, last_step, incr;
+ int i, first_step, last_step;
first_step = P->opaque->steps;
last_step = 0;
- incr = -1;
-
- for (i = first_step; i != last_step; i += incr) {
+ for (i = first_step; i != last_step; i--) {
+ pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f",
+ i - 1, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
+ P->opaque->pipeline[i]->a, P->opaque->pipeline[i]->rf
+ );
if (P->opaque->omit_inverse[i])
continue;
if (P->opaque->reverse_step[i])
@@ -213,6 +223,7 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
if (P->opaque->depth < PIPELINE_STACK_SIZE)
P->opaque->stack[P->opaque->depth++] = point;
}
+ pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
P->opaque->depth = 0; /* Clear the stack */
return point;
@@ -264,7 +275,7 @@ static void *pipeline_freeup (PJ *P, int errlev) { /* Destructor */
if (0==P)
return 0;
- pj_error_set (P, errlev);
+ pj_err_level (P, errlev);
if (0==P->opaque)
return pj_dealloc (P);
@@ -275,6 +286,8 @@ static void *pipeline_freeup (PJ *P, int errlev) { /* Destructor */
pj_dealloc (P->opaque->reverse_step);
pj_dealloc (P->opaque->omit_forward);
pj_dealloc (P->opaque->omit_inverse);
+ pj_dealloc (P->opaque->argv);
+ pj_dealloc (P->opaque->current_argv);
pj_dealloc (P->opaque->pipeline);
pj_dealloc (P->opaque);
@@ -313,8 +326,9 @@ static PJ *pj_create_pipeline (PJ *P, size_t steps) {
return P;
}
+
/* count the number of args in pipeline definition */
-size_t argc_params (paralist *params) {
+static size_t argc_params (paralist *params) {
size_t argc = 0;
for (; params != 0; params = params->next)
argc++;
@@ -322,18 +336,18 @@ size_t argc_params (paralist *params) {
}
/* Sentinel for argument list */
-static char argv_sentinel[5] = "step";
+static char argv_sentinel[] = "step";
-/* turn paralist int argc/argv style argument list */
-char **argv_params (paralist *params) {
+/* turn paralist into argc/argv style argument list */
+static char **argv_params (paralist *params, size_t argc) {
char **argv;
- size_t argc = 0;
- argv = pj_calloc (argc_params (params), sizeof (char *));
+ size_t i = 0;
+ argv = pj_calloc (argc, sizeof (char *));
if (0==argv)
return 0;
for (; params != 0; params = params->next)
- argv[argc++] = params->param;
- argv[argc++] = argv_sentinel;
+ argv[i++] = params->param;
+ argv[i++] = argv_sentinel;
return argv;
}
@@ -356,17 +370,15 @@ PJ *PROJECTION(pipeline) {
return 0;
argc = argc_params (P->params);
- argv = argv_params (P->params);
+ P->opaque->argv = argv = argv_params (P->params, argc);
if (0==argv)
return pipeline_freeup (P, ENOMEM);
- /* The elements of current_argv are not used - we just use argv_params */
- /* as allocator for a "large enough" container needed later */
- current_argv = argv_params (P->params);
+ P->opaque->current_argv = current_argv = pj_calloc (argc, sizeof (char *));
if (0==current_argv)
return pipeline_freeup (P, ENOMEM);
- /* Do some syntactic sanity checking */
+ /* Do some syntactical sanity checking */
for (i = 0; i < argc; i++) {
if (0==strcmp ("step", argv[i])) {
if (-1==i_pipeline) {
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 9e1ac4f7..2fcc4e31 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -77,6 +77,8 @@ SET(SRC_LIBPROJ_PJ
PJ_gstmerc.c
PJ_hammer.c
PJ_hatano.c
+ PJ_helmert.c
+ PJ_horner.c
PJ_igh.c
PJ_isea.c
PJ_imw_p.c
diff --git a/src/makefile.vc b/src/makefile.vc
index a2869217..e1eacf04 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -59,7 +59,7 @@ support = \
pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj
pipeline = \
- pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj
+ pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj
geodesic = geodesic.obj
@@ -147,6 +147,7 @@ install: all
copy *.exe $(INSTDIR)\bin
copy *.dll $(INSTDIR)\bin
copy *.lib $(INSTDIR)\lib
+ copy proj.h $(INSTDIR)\include
copy proj_api.h $(INSTDIR)\include
copy projects.h $(INSTDIR)\include
copy geodesic.h $(INSTDIR)\include
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 469e9a47..598c6dd2 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -6,41 +6,52 @@
XY /* forward projection entry */
pj_fwd(LP lp, PJ *P) {
XY xy;
- double t;
-
- /* check for forward and latitude or longitude overange */
- if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
- xy.x = xy.y = HUGE_VAL;
- pj_ctx_set_errno( P->ctx, -14);
- } else { /* proceed with projection */
- P->ctx->last_errno = 0;
- pj_errno = 0;
- errno = 0;
-
- if (fabs(t) <= EPS)
- lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
- else if (P->geoc)
- lp.phi = atan(P->rone_es * tan(lp.phi));
- lp.lam -= P->lam0; /* compute del lp.lam */
- if (!P->over)
- lp.lam = adjlon(lp.lam); /* adjust del longitude */
-
- /* Check for NULL pointer */
- if (P->fwd != NULL)
- {
- xy = (*P->fwd)(lp, P); /* project */
- if ( P->ctx->last_errno )
- xy.x = xy.y = HUGE_VAL;
- /* adjust for major axis and easting/northings */
- else {
- xy.x = P->fr_meter * (P->a * xy.x + P->x0);
- xy.y = P->fr_meter * (P->a * xy.y + P->y0);
- }
- }
- else
- {
- xy.x = xy.y = HUGE_VAL;
- }
- }
- return xy;
+ XY err;
+ double t;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.x = err.y = HUGE_VAL;
+
+ if (0==P->fwd)
+ return err;
+
+ P->ctx->last_errno = pj_errno = errno = 0;
+
+ /* Check input coordinates if angular */
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+
+ /* check for forward and latitude or longitude overange */
+ t = fabs(lp.phi)-M_HALFPI;
+ if (t > EPS || fabs(lp.lam) > 10.) {
+ pj_ctx_set_errno( P->ctx, -14);
+ return err;
+ }
+
+ /* Clamp latitude to -90..90 degree range */
+ if (fabs(t) <= EPS)
+ lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
+ else if (P->geoc) /* Maybe redundant and never used. */
+ lp.phi = atan(P->rone_es * tan(lp.phi));
+ lp.lam -= P->lam0; /* compute del lp.lam */
+ if (!P->over)
+ lp.lam = adjlon(lp.lam); /* adjust del longitude */
+ }
+
+ /* Do the transformation */
+ xy = (*P->fwd)(lp, P);
+ if ( P->ctx->last_errno )
+ return err;
+
+ /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xy.x *= P->a;
+ xy.y *= P->a;
+ }
+
+ /* Handle false eastings/northings and non-metric linear units */
+ xy.x = P->fr_meter * (xy.x + P->x0);
+ xy.y = P->fr_meter * (xy.y + P->y0);
+ /* z is not scaled since this is handled by vto_meter outside */
+
+ return xy;
}
diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c
index 63fc529e..0fb2f633 100644
--- a/src/pj_fwd3d.c
+++ b/src/pj_fwd3d.c
@@ -1,47 +1,59 @@
-/* general forward projection */
#define PJ_LIB__
#include <projects.h>
#include <errno.h>
# define EPS 1.0e-12
- XYZ /* forward projection entry */
-pj_fwd3d(LPZ lpz, PJ *P) {
- XYZ xyz;
- double t;
-
- /* check for forward and latitude or longitude overange */
- if ((t = fabs(lpz.phi)-M_HALFPI) > EPS || fabs(lpz.lam) > 10.) {
- xyz.x = xyz.y = xyz.z = HUGE_VAL;
- pj_ctx_set_errno( P->ctx, -14);
- } else { /* proceed with projection */
- P->ctx->last_errno = 0;
- pj_errno = 0;
- errno = 0;
-
- if (fabs(t) <= EPS)
- lpz.phi = lpz.phi < 0. ? -M_HALFPI : M_HALFPI;
- else if (P->geoc) /* Maybe redundant and never used. */
- lpz.phi = atan(P->rone_es * tan(lpz.phi));
- lpz.lam -= P->lam0; /* compute del lp.lam */
- if (!P->over)
- lpz.lam = adjlon(lpz.lam); /* adjust del longitude */
-
- /* Check for NULL pointer */
- if (P->fwd3d != NULL)
- {
- xyz = (*P->fwd3d)(lpz, P); /* project */
- if ( P->ctx->last_errno )
- xyz.x = xyz.y = xyz.z = HUGE_VAL;
- /* adjust for major axis and easting/northings */
- else {
- xyz.x = P->fr_meter * (P->a * xyz.x + P->x0);
- xyz.y = P->fr_meter * (P->a * xyz.y + P->y0);
- /* z is not scaled since this handled by vto_meter outside */
- }
- }
- else
- {
- xyz.x = xyz.y = xyz.z = HUGE_VAL;
- }
- }
- return xyz;
+
+/* 3D Forward transformation */
+
+XYZ pj_fwd3d(LPZ lpz, PJ *P) {
+ XYZ xyz;
+ XYZ err;
+ double t;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.x = err.y = err.z = HUGE_VAL;
+
+ if (0==P->fwd3d)
+ return err;
+
+
+ P->ctx->last_errno = pj_errno = errno = 0;
+
+ /* Check input coordinates if angular */
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+
+ /* check for forward and latitude or longitude overange */
+ t = fabs(lpz.phi)-M_HALFPI;
+ if (t > EPS || fabs(lpz.lam) > 10.) {
+ pj_ctx_set_errno( P->ctx, -14);
+ return err;
+ }
+
+ /* Clamp latitude to -90..90 degree range */
+ if (fabs(t) <= EPS)
+ lpz.phi = lpz.phi < 0. ? -M_HALFPI : M_HALFPI;
+ else if (P->geoc) /* Maybe redundant and never used. */
+ lpz.phi = atan(P->rone_es * tan(lpz.phi));
+ lpz.lam -= P->lam0; /* compute del lp.lam */
+ if (!P->over)
+ lpz.lam = adjlon(lpz.lam); /* adjust del longitude */
+ }
+
+ /* Do the transformation */
+ xyz = (*P->fwd3d)(lpz, P);
+ if ( P->ctx->last_errno )
+ return err;
+
+ /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xyz.x *= P->a;
+ xyz.y *= P->a;
+ }
+
+ /* Handle false eastings/northings and non-metric linear units */
+ xyz.x = P->fr_meter * (xyz.x + P->x0);
+ xyz.y = P->fr_meter * (xyz.y + P->y0);
+ /* z is not scaled since this is handled by vto_meter outside */
+
+ return xyz;
}
diff --git a/src/pj_init.c b/src/pj_init.c
index c468f1d7..b6be52d9 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -35,10 +35,12 @@
#include <errno.h>
#include <ctype.h>
+/* Maximum size of files using the "escape carriage return" feature */
+#define MAX_CR_ESCAPE 65537
typedef struct {
projCtx ctx;
PAFile fid;
- char buffer[8193];
+ char buffer[MAX_CR_ESCAPE];
int buffer_filled;
int at_eof;
} pj_read_state;
@@ -51,6 +53,7 @@ static const char *fill_buffer(pj_read_state *state, const char *last_char)
{
size_t bytes_read;
size_t char_remaining, char_requested;
+ char *r, *w;
/* -------------------------------------------------------------------- */
/* Don't bother trying to read more if we are at eof, or if the */
@@ -85,7 +88,36 @@ static const char *fill_buffer(pj_read_state *state, const char *last_char)
state->buffer[state->buffer_filled + bytes_read] = '\0';
}
- state->buffer_filled += bytes_read;
+/* -------------------------------------------------------------------- */
+/* Line continuations: skip whitespace after escaped newlines */
+/* -------------------------------------------------------------------- */
+ r = state->buffer;
+ w = state->buffer;
+ while (*r) {
+ /* Escaped newline? */
+ while ((r[0]=='\\') && ((r[1]=='\n') || (r[1]=='\r'))) {
+ r += 2;
+ while (isspace (*r))
+ r++;
+ /* we also skip comments immediately after an escaped newline */
+ while (*r=='#') {
+ while( *r && (*r != '\n') )
+ r++;
+ while (isspace (*r))
+ r++;
+ /* Reaching end of buffer while skipping continuation comment is currently an error */
+ if (0==*r) {
+ pj_ctx_set_errno (state->ctx, -2);
+ pj_log (state->ctx, PJ_LOG_ERROR, "init file too big");
+ return 0;
+ }
+ }
+ }
+ *w++ = *r++;
+ }
+ *w = 0;
+ state->buffer_filled += (bytes_read - (r-w));
+
return last_char;
}
@@ -96,11 +128,11 @@ static paralist *
get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
int *found_def) {
pj_read_state *state = (pj_read_state*) calloc(1,sizeof(pj_read_state));
- char sword[302];
+ char sword[MAX_CR_ESCAPE];
+ char *pipeline;
int len;
int in_target = 0;
const char *next_char = NULL;
-
state->fid = fid;
state->ctx = ctx;
next_char = fill_buffer(state, NULL);
@@ -110,6 +142,9 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
len = strlen(name);
*sword = 't';
+ if (0==next_char)
+ return 0;
+
/* loop till we find our target keyword */
while (*next_char)
{
@@ -120,6 +155,8 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
next_char++;
next_char = fill_buffer(state, next_char);
+ if (0==next_char)
+ return 0;
/* for comments, skip past end of line. */
if( *next_char == '#' )
@@ -128,6 +165,8 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
next_char++;
next_char = fill_buffer(state, next_char);
+ if (0==next_char)
+ return 0;
if (*next_char == '\n')
next_char++;
if (*next_char == '\r')
@@ -171,7 +210,7 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
}
/* capture parameter */
- while( *next_char && !isspace(*next_char) )
+ while ( *next_char && !isspace(*next_char) )
{
next_char++;
word_len++;
@@ -180,8 +219,11 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
strncpy(sword+1, start_of_word, word_len);
sword[word_len+1] = '\0';
- /* do not override existing parameter value of same name */
- if (!pj_param(ctx, *start, sword).i) {
+ /* do not override existing parameter value of same name - unless in pipeline definition */
+ pipeline = pj_param(ctx, *start, "sproj").s;
+ if (pipeline && strcmp (pipeline, "pipeline"))
+ pipeline = 0;
+ if (pipeline || !pj_param(ctx, *start, sword).i) {
/* don't default ellipse if datum, ellps or any earth model
information is set. */
if( strncmp(sword+1,"ellps=",6) != 0
@@ -194,8 +236,9 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
{
next = next->next = pj_mkparam(sword+1);
}
+ else
+ next = next->next = pj_mkparam(sword+1);
}
-
}
else
{
@@ -388,6 +431,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
PJ *(*proj)(PJ *);
paralist *curr;
int i;
+ int defer_init_expansion = 0;
PJ *PIN = 0;
ctx->last_errno = 0;
@@ -403,9 +447,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
if (ctx->last_errno) goto bum_call;
/* check if +init present */
- if (pj_param(ctx, start, "tinit").i) {
+ if (pj_param(ctx, start, "tinit").i && ! defer_init_expansion) {
int found_def = 0;
-
+ /* avoid expanding additional inits (as could happen in a pipeline) */
+ defer_init_expansion = 1;
if (!(curr = get_init(ctx,&start, curr,
pj_param(ctx, start, "sinit").s,
&found_def)))
@@ -444,7 +489,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
if (pj_datum_set(ctx, start, PIN)) goto bum_call;
/* set ellipsoid/sphere parameters */
- if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) goto bum_call;
+ if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) {
+ pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
+ goto bum_call;
+ }
PIN->a_orig = PIN->a;
PIN->es_orig = PIN->es;
@@ -464,11 +512,11 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
/* flattening */
PIN->f = 1 - cos (PIN->alpha); /* = 1 - sqrt (1 - PIN->es); */
- PIN->rf = PIN->f? 1/PIN->f: HUGE_VAL;
+ PIN->rf = PIN->f? 1.0/PIN->f: HUGE_VAL;
/* second flattening */
- PIN->f2 = 1/cos (PIN->alpha) - 1;
- PIN->rf = PIN->f2? 1/PIN->f2: HUGE_VAL;
+ PIN->f2 = 1/cos (PIN->alpha) - 1;
+ PIN->rf2 = PIN->f2? 1/PIN->f2: HUGE_VAL;
/* third flattening */
PIN->n = pow (tan (PIN->alpha/2), 2);
diff --git a/src/pj_inv.c b/src/pj_inv.c
index ed05618a..d5393261 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -7,35 +7,49 @@
/* inverse projection entry */
LP pj_inv(XY xy, PJ *P) {
LP lp;
+ LP err;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.lam = err.phi = HUGE_VAL;
+
+ if (0==P->inv)
+ return err;
/* can't do as much preliminary checking as with forward */
if (xy.x == HUGE_VAL || xy.y == HUGE_VAL) {
- lp.lam = lp.phi = HUGE_VAL;
pj_ctx_set_errno( P->ctx, -15);
- return lp;
+ return err;
+ }
+ P->ctx->last_errno = errno = pj_errno = 0;
+
+ /* de-scale and de-offset */
+ xy.x = (xy.x * P->to_meter - P->x0);
+ xy.y = (xy.y * P->to_meter - P->y0);
+
+ /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
+ /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */
+ /* stomps on a and hence depends on this */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xy.x *= P->ra;
+ xy.y *= P->ra;
}
- errno = pj_errno = 0;
- P->ctx->last_errno = 0;
-
- /* descale and de-offset */
- xy.x = (xy.x * P->to_meter - P->x0) * P->ra;
- xy.y = (xy.y * P->to_meter - P->y0) * P->ra;
-
- /* Check for NULL pointer */
- if (P->inv != NULL) {
- lp = (*P->inv)(xy, P); /* inverse project */
- if (P->ctx->last_errno ) {
- lp.lam = lp.phi = HUGE_VAL;
- } else {
- lp.lam += P->lam0; /* reduce from del lp.lam */
- if (!P->over)
- lp.lam = adjlon(lp.lam); /* adjust longitude to CM */
- if (P->geoc && fabs(fabs(lp.phi)-M_HALFPI) > EPS)
- lp.phi = atan(P->one_es * tan(lp.phi));
- }
- } else {
- lp.lam = lp.phi = HUGE_VAL;
+ /* Do inverse transformation */
+ lp = (*P->inv) (xy, P);
+ if (P->ctx->last_errno)
+ return err;
+
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+ /* reduce from del lp.lam */
+ lp.lam += P->lam0;
+
+ /* adjust longitude to central meridian */
+ if (!P->over)
+ lp.lam = adjlon(lp.lam);
+
+ /* This may be redundant and never used */
+ if (P->geoc && fabs(fabs(lp.phi)-M_HALFPI) > EPS)
+ lp.phi = atan(P->one_es * tan(lp.phi));
}
return lp;
diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c
index bf6a9324..c6052ad4 100644
--- a/src/pj_inv3d.c
+++ b/src/pj_inv3d.c
@@ -1,45 +1,57 @@
-/* general inverse projection */
#define PJ_LIB__
#include <projects.h>
#include <errno.h>
# define EPS 1.0e-12
- LPZ /* inverse projection entry */
-pj_inv3d(XYZ xyz, PJ *P) {
- LPZ lpz;
-
- /* can't do as much preliminary checking as with forward */
- if (xyz.x == HUGE_VAL || xyz.y == HUGE_VAL || xyz.z == HUGE_VAL ) {
- lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
- pj_ctx_set_errno( P->ctx, -15);
- return lpz;
- }
-
- errno = pj_errno = 0;
- P->ctx->last_errno = 0;
-
- xyz.x = (xyz.x * P->to_meter - P->x0) * P->ra; /* descale and de-offset */
- xyz.y = (xyz.y * P->to_meter - P->y0) * P->ra;
- /* z is not scaled since that is handled by vto_meter before we get here */
-
- /* Check for NULL pointer */
- if (P->inv3d != NULL)
- {
- lpz = (*P->inv3d)(xyz, P); /* inverse project */
- if (P->ctx->last_errno )
- lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
- else {
- lpz.lam += P->lam0; /* reduce from del lp.lam */
- if (!P->over)
- lpz.lam = adjlon(lpz.lam); /* adjust longitude to CM */
-
- /* This may be redundant and never used */
- if (P->geoc && fabs(fabs(lpz.phi)-M_HALFPI) > EPS)
- lpz.phi = atan(P->one_es * tan(lpz.phi));
- }
- }
- else
- {
- lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
- }
- return lpz;
+
+/* 3D inverse transformation */
+
+LPZ pj_inv3d (XYZ xyz, PJ *P) {
+ LPZ lpz;
+ LPZ err;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.lam = err.phi = err.z = HUGE_VAL;
+
+ if (0==P->inv3d)
+ return err;
+
+ /* can't do as much preliminary checking as with forward */
+ if (xyz.x == HUGE_VAL || xyz.y == HUGE_VAL || xyz.z == HUGE_VAL ) {
+ pj_ctx_set_errno( P->ctx, -15);
+ return err;
+ }
+
+ P->ctx->last_errno = errno = pj_errno = 0;
+
+ /* de-scale and de-offset */
+ /* z is not de-scaled since that is handled by vto_meter before we get here */
+ xyz.x = (xyz.x * P->to_meter - P->x0);
+ xyz.y = (xyz.y * P->to_meter - P->y0);
+ /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
+ /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */
+ /* stomps on a and hence depends on this */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xyz.x *= P->ra;
+ xyz.y *= P->ra;
+ }
+
+ /* Do inverse transformation */
+ lpz = (*P->inv3d) (xyz, P);
+ if (P->ctx->last_errno)
+ return err;
+
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+ /* reduce from del lp.lam */
+ lpz.lam += P->lam0;
+
+ /* adjust longitude to central meridian */
+ if (!P->over)
+ lpz.lam = adjlon(lpz.lam);
+
+ /* This may be redundant and never used */
+ if (P->geoc && fabs(fabs(lpz.phi)-M_HALFPI) > EPS)
+ lpz.phi = atan(P->one_es * tan(lpz.phi));
+ }
+
+ return lpz;
}
diff --git a/src/pj_list.h b/src/pj_list.h
index f49fd574..96db13d7 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -52,6 +52,8 @@ PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff")
PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area")
PROJ_HEAD(healpix, "HEALPix")
PROJ_HEAD(rhealpix, "rHEALPix")
+PROJ_HEAD(helmert, "3- and 7-parameter Helmert shift")
+PROJ_HEAD(horner, "Horner polynomial evaluation")
PROJ_HEAD(igh, "Interrupted Goode Homolosine")
PROJ_HEAD(imw_p, "International Map of the World Polyconic")
PROJ_HEAD(isea, "Icosahedral Snyder Equal Area")
diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
index 716c6095..4a6e98a5 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -32,7 +32,7 @@
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
*****************************************************************************/
-#define PJ_OBS_C
+#define PJ_OBS_API_C
#include <proj.h>
#include <projects.h>
#include <geodesic.h>
@@ -40,16 +40,9 @@
#include <math.h>
-/* Used as return value in case of errors */
-const PJ_OBS pj_obs_error = {
- /* Cannot use HUGE_VAL here: MSVC misimplements HUGE_VAL as something that is not compile time constant */
- {{DBL_MAX,DBL_MAX,DBL_MAX,DBL_MAX}},
- {{DBL_MAX,DBL_MAX,DBL_MAX}},
- 0, 0
-};
-
/* Used for zero-initializing new objects */
-const PJ_OBS pj_obs_null = {
+const PJ_COORD pj_coo_null = {{0, 0, 0, 0}};
+const PJ_OBS pj_obs_null = {
{{0, 0, 0, 0}},
{{0, 0, 0}},
0, 0
@@ -77,6 +70,23 @@ double pj_xyz_dist (XYZ a, XYZ b) {
}
+/* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */
+PJ_COORD pj_coo_error (void) {
+ PJ_COORD c;
+ c.v[0] = c.v[1] = c.v[2] = c.v[3] = HUGE_VAL;
+ return c;
+}
+
+PJ_OBS pj_obs_error (void) {
+ PJ_OBS obs;
+ obs.coo = pj_coo_error ();
+ obs.anc.v[0] = obs.anc.v[1] = obs.anc.v[2] = HUGE_VAL;
+ obs.id = obs.flags = 0;
+ return obs;
+}
+
+
+
PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
if (0!=P->fwdobs) {
obs = P->fwdobs (obs, P);
@@ -90,8 +100,8 @@ PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
obs.coo.xy = pj_fwd (obs.coo.lp, P);
return obs;
}
- pj_error_set (P, EINVAL);
- return pj_obs_error;
+ pj_err_level (P, EINVAL);
+ return pj_obs_error ();
}
@@ -108,8 +118,8 @@ PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) {
obs.coo.lp = pj_inv (obs.coo.xy, P);
return obs;
}
- pj_error_set (P, EINVAL);
- return pj_obs_error;
+ pj_err_level (P, EINVAL);
+ return pj_obs_error ();
}
@@ -129,8 +139,8 @@ PJ_OBS pj_trans (PJ *P, enum pj_direction direction, PJ_OBS obs) {
break;
}
- pj_error_set (P, EINVAL);
- return pj_obs_error;
+ pj_err_level (P, EINVAL);
+ return pj_obs_error ();
}
@@ -143,7 +153,7 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
return HUGE_VAL;
if (n < 1) {
- pj_error_set (P, EINVAL);
+ pj_err_level (P, EINVAL);
return HUGE_VAL;
}
@@ -152,15 +162,15 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
for (i = 0; i < n; i++) {
switch (direction) {
case PJ_FWD:
- u.coo.xyz = pj_fwd3d (o.coo.lpz, P);
- o.coo.lpz = pj_inv3d (u.coo.xyz, P);
+ u = pj_fwdobs (o, P);
+ o = pj_invobs (u, P);
break;
case PJ_INV:
- u.coo.lpz = pj_inv3d (o.coo.xyz, P);
- o.coo.xyz = pj_fwd3d (u.coo.lpz, P);
+ u = pj_invobs (o, P);
+ o = pj_fwdobs (u, P);
break;
default:
- pj_error_set (P, EINVAL);
+ pj_err_level (P, EINVAL);
return HUGE_VAL;
}
}
@@ -179,18 +189,24 @@ PJ *pj_create_argv (int argc, char **argv) {
}
-/* From here: Minimum viable support for contexts. The first four functions */
+/* Below: Minimum viable support for contexts. The first four functions */
/* relate to error reporting, debugging, and logging, hence being generically */
/* useful. The remaining is a compact implementation of the more low level */
/* proj_api.h thread contexts, which may or may not be useful */
-int pj_error (PJ *P) {
- return pj_ctx_get_errno (pj_get_ctx(P));
-}
-
-
-void pj_error_set (PJ *P, int err) {
- pj_ctx_set_errno (pj_get_ctx(P), err);
+/* Set error level 0-3, or query current level */
+int pj_err_level (PJ *P, int err_level) {
+ int previous;
+ PJ_CONTEXT *ctx;
+ if (0==P)
+ ctx = pj_get_default_ctx();
+ else
+ ctx = pj_get_ctx (P);
+ previous = pj_ctx_get_errno (ctx);
+ if (PJ_ERR_TELL==err_level)
+ return previous;
+ pj_ctx_set_errno (ctx, err_level);
+ return previous;
}
@@ -211,7 +227,7 @@ enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level) {
/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */
-void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *)) {
+void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *)) {
PJ_CONTEXT *ctx = pj_get_ctx (P);
ctx->app_data = app_data;
if (0!=log)
@@ -223,7 +239,7 @@ void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *))
int pj_context_renew (PJ *P) {
PJ_CONTEXT *ctx = pj_ctx_alloc ();
if (0==ctx) {
- pj_error_set (P, ENOMEM);
+ pj_err_level (P, ENOMEM);
return 1;
}
@@ -231,6 +247,7 @@ int pj_context_renew (PJ *P) {
return 0;
}
+
/* Move daughter to mother's context */
void pj_context_inherit (PJ *mother, PJ *daughter) {
pj_set_ctx (daughter, pj_get_ctx (mother));
diff --git a/src/pj_param.c b/src/pj_param.c
index 434d4535..24c2354e 100644
--- a/src/pj_param.c
+++ b/src/pj_param.c
@@ -23,13 +23,13 @@ pj_mkparam(char *str) {
/* character in `opt' is a parameter type which can take the */
/* values: */
/* */
-/* `t' - test for presence, return TRUE/FALSE in PROJVALUE.i */
-/* `i' - integer value returned in PROJVALUE.i */
-/* `d' - simple valued real input returned in PROJVALUE.f */
+/* `t' - test for presence, return TRUE/FALSE in PROJVALUE.i */
+/* `i' - integer value returned in PROJVALUE.i */
+/* `d' - simple valued real input returned in PROJVALUE.f */
/* `r' - degrees (DMS translation applied), returned as */
-/* radians in PROJVALUE.f */
-/* `s' - string returned in PROJVALUE.s */
-/* `b' - test for t/T/f/F, return in PROJVALUE.i */
+/* radians in PROJVALUE.f */
+/* `s' - string returned in PROJVALUE.s */
+/* `b' - test for t/T/f/F, return in PROJVALUE.i */
/* */
/************************************************************************/
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index 8a2a9c4b..bef42882 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -82,4 +82,5 @@ char *pj_strerrno(int err) {
return note;
}
}
+ return 0;
}
diff --git a/src/proj.def b/src/proj.def
index 3a5e3de9..b01b846c 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -89,23 +89,24 @@ EXPORTS
geod_polygon_testpoint @87
geod_polygon_clear @88
pj_run_selftests @89
- pj_error @90
- pj_create @91
- pj_create_argv @92
- pj_trans @93
- pj_fwdobs @94
- pj_invobs @95
- pj_roundtrip @96
- pj_log_level @97
- pj_error_set @98
- pj_log_set @99
- pj_context_renew @100
- pj_context_inherit @101
- pj_context_free @102
- pj_fileapi_set @103
- pj_log_error @104
- pj_log_debug @105
- pj_log_trace @106
- pj_lp_dist @107
- pj_xy_dist @108
- pj_xyz_dist @109
+ pj_create @90
+ pj_create_argv @91
+ pj_trans @92
+ pj_fwdobs @93
+ pj_invobs @94
+ pj_roundtrip @95
+ pj_log_level @96
+ pj_err_level @97
+ pj_log_func @98
+ pj_context_renew @99
+ pj_context_inherit @100
+ pj_context_free @101
+ pj_fileapi_set @102
+ pj_log_error @103
+ pj_log_debug @104
+ pj_log_trace @105
+ pj_lp_dist @106
+ pj_xy_dist @107
+ pj_xyz_dist @108
+ pj_coo_error @109
+ pj_obs_error @110
diff --git a/src/proj.h b/src/proj.h
index 6c33a916..f6b3bd2a 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -271,10 +271,13 @@ double pj_xy_dist (XY a, XY b);
double pj_xyz_dist (XYZ a, XYZ b);
-#ifndef PJ_OBS_C
-extern const PJ_OBS pj_obs_error;
-extern const PJ_OBS pj_obs_null;
-extern const PJ *pj_shutdown;
+#ifndef PJ_OBS_API_C
+/* Part of MSVC workaround: Make pj_*_null look function-like for symmetry with pj_*_error */
+#define pj_coo_null(x) pj_coo_null
+#define pj_obs_null(x) pj_obs_null
+extern const PJ_COORD pj_coo_null;
+extern const PJ_OBS pj_obs_null;
+extern const PJ *pj_shutdown;
#endif
#ifndef TODEG
@@ -301,17 +304,19 @@ enum pj_log_level {
PJ_LOG_DEBUG_MINOR = 3 /* for proj_api.h compatibility */
};
+/* Set or read error level */
+#define PJ_ERR_TELL -56789
+int pj_err_level (PJ *P, int err_level);
+
/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level);
void pj_log_error (PJ *P, const char *fmt, ...);
void pj_log_debug (PJ *P, const char *fmt, ...);
void pj_log_trace (PJ *P, const char *fmt, ...);
+void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));
-void pj_error_set (PJ *P, int err);
-void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *));
-
/* Lower level functionality for handling thread contexts */
int pj_context_renew (PJ *P);
void pj_context_inherit (PJ *mother, PJ *daughter);
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 12:13:26 +0000