1 #define PJ_LIB__
2
3 #include <errno.h>
4 #include <math.h>
5
6 #include "proj.h"
7 #include "proj_internal.h"
8
9 PROJ_HEAD(moll, "Mollweide") "\n\tPCyl, Sph";
10 PROJ_HEAD(wag4, "Wagner IV") "\n\tPCyl, Sph";
11 PROJ_HEAD(wag5, "Wagner V") "\n\tPCyl, Sph";
12
13 #define MAX_ITER 10
14 #define LOOP_TOL 1e-7
15
16 namespace { // anonymous namespace
17 struct pj_opaque {
18 double C_x, C_y, C_p;
19 };
20 } // anonymous namespace
21
22
moll_s_forward(PJ_LP lp,PJ * P)23 static PJ_XY moll_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
24 PJ_XY xy = {0.0,0.0};
25 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
26 int i;
27
28 const double k = Q->C_p * sin(lp.phi);
29 for (i = MAX_ITER; i ; --i) {
30 const double V = (lp.phi + sin(lp.phi) - k) /
31 (1. + cos(lp.phi));
32 lp.phi -= V;
33 if (fabs(V) < LOOP_TOL)
34 break;
35 }
36 if (!i)
37 lp.phi = (lp.phi < 0.) ? -M_HALFPI : M_HALFPI;
38 else
39 lp.phi *= 0.5;
40 xy.x = Q->C_x * lp.lam * cos(lp.phi);
41 xy.y = Q->C_y * sin(lp.phi);
42 return xy;
43 }
44
45
moll_s_inverse(PJ_XY xy,PJ * P)46 static PJ_LP moll_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
47 PJ_LP lp = {0.0,0.0};
48 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
49 lp.phi = aasin(P->ctx, xy.y / Q->C_y);
50 lp.lam = xy.x / (Q->C_x * cos(lp.phi));
51 if (fabs(lp.lam) < M_PI) {
52 lp.phi += lp.phi;
53 lp.phi = aasin(P->ctx, (lp.phi + sin(lp.phi)) / Q->C_p);
54 } else {
55 lp.lam = lp.phi = HUGE_VAL;
56 }
57 return lp;
58 }
59
60
setup(PJ * P,double p)61 static PJ * setup(PJ *P, double p) {
62 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
63 double r, sp, p2 = p + p;
64
65 P->es = 0;
66 sp = sin(p);
67 r = sqrt(M_TWOPI * sp / (p2 + sin(p2)));
68
69 Q->C_x = 2. * r / M_PI;
70 Q->C_y = r / sp;
71 Q->C_p = p2 + sin(p2);
72
73 P->inv = moll_s_inverse;
74 P->fwd = moll_s_forward;
75 return P;
76 }
77
78
PROJECTION(moll)79 PJ *PROJECTION(moll) {
80 struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
81 if (nullptr==Q)
82 return pj_default_destructor (P, ENOMEM);
83 P->opaque = Q;
84
85 return setup(P, M_HALFPI);
86 }
87
88
PROJECTION(wag4)89 PJ *PROJECTION(wag4) {
90 struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
91 if (nullptr==Q)
92 return pj_default_destructor (P, ENOMEM);
93 P->opaque = Q;
94
95 return setup(P, M_PI/3.);
96 }
97
PROJECTION(wag5)98 PJ *PROJECTION(wag5) {
99 struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
100 if (nullptr==Q)
101 return pj_default_destructor (P, ENOMEM);
102 P->opaque = Q;
103
104 P->es = 0;
105 Q->C_x = 0.90977;
106 Q->C_y = 1.65014;
107 Q->C_p = 3.00896;
108
109 P->inv = moll_s_inverse;
110 P->fwd = moll_s_forward;
111
112 return P;
113 }
114