1 /*
2 The Natural Earth projection was designed by Tom Patterson, US National Park
3 Service, in 2007, using Flex Projector. The shape of the original projection
4 was defined at every 5 degrees and piece-wise cubic spline interpolation was
5 used to compute the complete graticule.
6 The code here uses polynomial functions instead of cubic splines and
7 is therefore much simpler to program. The polynomial approximation was
8 developed by Bojan Savric, in collaboration with Tom Patterson and Bernhard
9 Jenny, Institute of Cartography, ETH Zurich. It slightly deviates from
10 Patterson's original projection by adding additional curvature to meridians
11 where they meet the horizontal pole line. This improvement is by intention
12 and designed in collaboration with Tom Patterson.
13 Port to PROJ.4 by Bernhard Jenny, 6 June 2011
14 */
15 #define PJ_LIB__
16
17 #include <math.h>
18
19 #include "proj.h"
20 #include "proj_internal.h"
21
22 PROJ_HEAD(natearth, "Natural Earth") "\n\tPCyl, Sph";
23
24 #define A0 0.8707
25 #define A1 -0.131979
26 #define A2 -0.013791
27 #define A3 0.003971
28 #define A4 -0.001529
29 #define B0 1.007226
30 #define B1 0.015085
31 #define B2 -0.044475
32 #define B3 0.028874
33 #define B4 -0.005916
34 #define C0 B0
35 #define C1 (3 * B1)
36 #define C2 (7 * B2)
37 #define C3 (9 * B3)
38 #define C4 (11 * B4)
39 #define EPS 1e-11
40 #define MAX_Y (0.8707 * 0.52 * M_PI)
41 /* Not sure at all of the appropriate number for MAX_ITER... */
42 #define MAX_ITER 100
43
44
natearth_s_forward(PJ_LP lp,PJ * P)45 static PJ_XY natearth_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
46 PJ_XY xy = {0.0,0.0};
47 double phi2, phi4;
48 (void) P;
49
50 phi2 = lp.phi * lp.phi;
51 phi4 = phi2 * phi2;
52 xy.x = lp.lam * (A0 + phi2 * (A1 + phi2 * (A2 + phi4 * phi2 * (A3 + phi2 * A4))));
53 xy.y = lp.phi * (B0 + phi2 * (B1 + phi4 * (B2 + B3 * phi2 + B4 * phi4)));
54 return xy;
55 }
56
57
natearth_s_inverse(PJ_XY xy,PJ * P)58 static PJ_LP natearth_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
59 PJ_LP lp = {0.0,0.0};
60 double yc, y2, y4, f, fder;
61 int i;
62 (void) P;
63
64 /* make sure y is inside valid range */
65 if (xy.y > MAX_Y) {
66 xy.y = MAX_Y;
67 } else if (xy.y < -MAX_Y) {
68 xy.y = -MAX_Y;
69 }
70
71 /* latitude */
72 yc = xy.y;
73 for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */
74 y2 = yc * yc;
75 y4 = y2 * y2;
76 f = (yc * (B0 + y2 * (B1 + y4 * (B2 + B3 * y2 + B4 * y4)))) - xy.y;
77 fder = C0 + y2 * (C1 + y4 * (C2 + C3 * y2 + C4 * y4));
78 const double tol = f / fder;
79 yc -= tol;
80 if (fabs(tol) < EPS) {
81 break;
82 }
83 }
84 if( i == 0 )
85 pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
86 lp.phi = yc;
87
88 /* longitude */
89 y2 = yc * yc;
90 lp.lam = xy.x / (A0 + y2 * (A1 + y2 * (A2 + y2 * y2 * y2 * (A3 + y2 * A4))));
91
92 return lp;
93 }
94
95
PROJECTION(natearth)96 PJ *PROJECTION(natearth) {
97 P->es = 0;
98 P->inv = natearth_s_inverse;
99 P->fwd = natearth_s_forward;
100
101 return P;
102 }
103