1 /*
2 The Natural Earth II projection was designed by Tom Patterson, US National
3 Park Service, in 2012, using Flex Projector. The polynomial equation was
4 developed by Bojan Savric and Bernhard Jenny, College of Earth, Ocean,
5 and Atmospheric Sciences, Oregon State University.
6 Port to PROJ.4 by Bojan Savric, 4 April 2016
7 */
8 #define PJ_LIB__
9 #include <projects.h>
10
11 PROJ_HEAD(natearth2, "Natural Earth 2") "\n\tPCyl., Sph.";
12
13 #define A0 0.84719
14 #define A1 -0.13063
15 #define A2 -0.04515
16 #define A3 0.05494
17 #define A4 -0.02326
18 #define A5 0.00331
19 #define B0 1.01183
20 #define B1 -0.02625
21 #define B2 0.01926
22 #define B3 -0.00396
23 #define C0 B0
24 #define C1 (9 * B1)
25 #define C2 (11 * B2)
26 #define C3 (13 * B3)
27 #define EPS 1e-11
28 #define MAX_Y (0.84719 * 0.535117535153096 * M_PI)
29
30
s_forward(LP lp,PJ * P)31 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
32 XY xy = {0.0,0.0};
33 double phi2, phi4, phi6;
34 (void) P;
35
36 phi2 = lp.phi * lp.phi;
37 phi4 = phi2 * phi2;
38 phi6 = phi2 * phi4;
39
40 xy.x = lp.lam * (A0 + A1 * phi2 + phi6 * phi6 * (A2 + A3 * phi2 + A4 * phi4 + A5 * phi6));
41 xy.y = lp.phi * (B0 + phi4 * phi4 * (B1 + B2 * phi2 + B3 * phi4));
42 return xy;
43 }
44
45
s_inverse(XY xy,PJ * P)46 static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
47 LP lp = {0.0,0.0};
48 double yc, tol, y2, y4, y6, f, fder;
49 (void) P;
50
51 /* make sure y is inside valid range */
52 if (xy.y > MAX_Y) {
53 xy.y = MAX_Y;
54 } else if (xy.y < -MAX_Y) {
55 xy.y = -MAX_Y;
56 }
57
58 /* latitude */
59 yc = xy.y;
60 for (;;) { /* Newton-Raphson */
61 y2 = yc * yc;
62 y4 = y2 * y2;
63 f = (yc * (B0 + y4 * y4 * (B1 + B2 * y2 + B3 * y4))) - xy.y;
64 fder = C0 + y4 * y4 * (C1 + C2 * y2 + C3 * y4);
65 yc -= tol = f / fder;
66 if (fabs(tol) < EPS) {
67 break;
68 }
69 }
70 lp.phi = yc;
71
72 /* longitude */
73 y2 = yc * yc;
74 y4 = y2 * y2;
75 y6 = y2 * y4;
76
77 lp.lam = xy.x / (A0 + A1 * y2 + y6 * y6 * (A2 + A3 * y2 + A4 * y4 + A5 * y6));
78
79 return lp;
80 }
81
82
freeup_new(PJ * P)83 static void *freeup_new (PJ *P) { /* Destructor */
84 if (0==P)
85 return 0;
86
87 return pj_dealloc(P);
88 }
89
90
freeup(PJ * P)91 static void freeup (PJ *P) {
92 freeup_new (P);
93 return;
94 }
95
96
PROJECTION(natearth2)97 PJ *PROJECTION(natearth2) {
98 P->es = 0;
99 P->inv = s_inverse;
100 P->fwd = s_forward;
101
102 return P;
103 }
104
105 #ifndef PJ_SELFTEST
pj_natearth2_selftest(void)106 int pj_natearth2_selftest (void) {return 0;}
107 #else
108
pj_natearth2_selftest(void)109 int pj_natearth2_selftest (void) {
110 double tolerance_lp = 1e-10;
111 double tolerance_xy = 1e-7;
112
113 char s_args[] = {"+proj=natearth2 +a=6400000 +lat_1=0.5 +lat_2=2"};
114
115 LP fwd_in[] = {
116 { 2, 1},
117 { 2,-1},
118 {-2, 1},
119 {-2,-1}
120 };
121
122 XY s_fwd_expect[] = {
123 { 189255.172934730799, 113022.495810907014},
124 { 189255.172934730799, -113022.495810907014},
125 {-189255.172934730799, 113022.495810907014},
126 {-189255.172934730799, -113022.495810907014},
127 };
128
129 XY inv_in[] = {
130 { 200, 100},
131 { 200,-100},
132 {-200, 100},
133 {-200,-100}
134 };
135
136 LP s_inv_expect[] = {
137 { 0.00211344929691056112, 0.000884779612080993237},
138 { 0.00211344929691056112, -0.000884779612080993237},
139 {-0.00211344929691056112, 0.000884779612080993237},
140 {-0.00211344929691056112, -0.000884779612080993237},
141 };
142
143 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);
144 }
145
146
147 #endif
148