diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/PJ_aea.c | 1 | ||||
| -rw-r--r-- | src/PJ_nell_h.c | 129 |
2 files changed, 103 insertions, 27 deletions
diff --git a/src/PJ_aea.c b/src/PJ_aea.c index 4cf1030f..e25adf13 100644 --- a/src/PJ_aea.c +++ b/src/PJ_aea.c @@ -364,7 +364,6 @@ int pj_latlong_selftest (void) {return 10000;} int pj_lonlat_selftest (void) {return 10000;} int pj_longlat_selftest (void) {return 10000;} -int pj_nell_h_selftest (void) {return 10000;} int pj_nicol_selftest (void) {return 10000;} int pj_nsper_selftest (void) {return 10000;} int pj_nzmg_selftest (void) {return 10000;} diff --git a/src/PJ_nell_h.c b/src/PJ_nell_h.c index 22061ea9..1e863945 100644 --- a/src/PJ_nell_h.c +++ b/src/PJ_nell_h.c @@ -1,32 +1,109 @@ #define PJ_LIB__ -# include <projects.h> +#include <projects.h> + PROJ_HEAD(nell_h, "Nell-Hammer") "\n\tPCyl., Sph."; + #define NITER 9 #define EPS 1e-7 -FORWARD(s_forward); /* spheroid */ - (void) P; - xy.x = 0.5 * lp.lam * (1. + cos(lp.phi)); - xy.y = 2.0 * (lp.phi - tan(0.5 *lp.phi)); - return (xy); + + +static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ + XY xy = {0.0,0.0}; + (void) P; + + xy.x = 0.5 * lp.lam * (1. + cos(lp.phi)); + xy.y = 2.0 * (lp.phi - tan(0.5 *lp.phi)); + + return xy; } -INVERSE(s_inverse); /* spheroid */ - double V, c, p; - int i; - (void) P; - - p = 0.5 * xy.y; - for (i = NITER; i ; --i) { - c = cos(0.5 * lp.phi); - lp.phi -= V = (lp.phi - tan(lp.phi/2) - p)/(1. - 0.5/(c*c)); - if (fabs(V) < EPS) - break; - } - if (!i) { - lp.phi = p < 0. ? -HALFPI : HALFPI; - lp.lam = 2. * xy.x; - } else - lp.lam = 2. * xy.x / (1. + cos(lp.phi)); - return (lp); + + +static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ + LP lp = {0.0,0.0}; + double V, c, p; + int i; + (void) P; + + p = 0.5 * xy.y; + for (i = NITER; i ; --i) { + c = cos(0.5 * lp.phi); + lp.phi -= V = (lp.phi - tan(lp.phi/2) - p)/(1. - 0.5/(c*c)); + if (fabs(V) < EPS) + break; + } + if (!i) { + lp.phi = p < 0. ? -HALFPI : HALFPI; + lp.lam = 2. * xy.x; + } else + lp.lam = 2. * xy.x / (1. + cos(lp.phi)); + + return lp; +} + + +static void *freeup_new (PJ *P) { /* Destructor */ + if (0==P) + return 0; + + return pj_dealloc(P); +} + + +static void freeup (PJ *P) { + freeup_new (P); + return; } -FREEUP; if (P) pj_dalloc(P); } -ENTRY0(nell_h) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P) + + +PJ *PROJECTION(nell_h) { + P->es = 0.; + P->inv = s_inverse; + P->fwd = s_forward; + + return P; +} + + +#ifdef PJ_OMIT_SELFTEST +int pj_nell_h_selftest (void) {return 0;} +#else + +int pj_nell_h_selftest (void) { + double tolerance_lp = 1e-10; + double tolerance_xy = 1e-7; + + char s_args[] = {"+proj=nell_h +a=6400000 +lat_1=0.5 +lat_2=2"}; + + LP fwd_in[] = { + { 2, 1}, + { 2,-1}, + {-2, 1}, + {-2,-1} + }; + + XY s_fwd_expect[] = { + { 223385.131640952837, 111698.236533561678}, + { 223385.131640952837, -111698.236533561678}, + {-223385.131640952837, 111698.236533561678}, + {-223385.131640952837, -111698.236533561678}, + }; + + XY inv_in[] = { + { 200, 100}, + { 200,-100}, + {-200, 100}, + {-200,-100} + }; + + LP s_inv_expect[] = { + { 0.00179049310989310567, 0.000895246554910125378}, + { 0.00179049310989310567, -0.000895246554910125378}, + {-0.00179049310989310567, 0.000895246554910125378}, + {-0.00179049310989310567, -0.000895246554910125378}, + }; + + return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); +} + + +#endif |