1 #define PJ_LIB__
2 #include <projects.h>
3
4 PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic") "\n\tCyl., Sph.";
5
6 #define NITER 20
7 #define EPS 1e-7
8 #define ONETOL 1.000001
9 #define C 1.70710678118654752440
10 #define RC 0.58578643762690495119
11 #define FYC 1.87475828462269495505
12 #define RYC 0.53340209679417701685
13 #define FXC 0.31245971410378249250
14 #define RXC 3.20041258076506210122
15
16
s_forward(LP lp,PJ * P)17 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
18 XY xy = {0.0,0.0};
19 double th1, c;
20 int i;
21 (void) P;
22
23 c = C * sin(lp.phi);
24 for (i = NITER; i; --i) {
25 lp.phi -= th1 = (sin(.5*lp.phi) + sin(lp.phi) - c) /
26 (.5*cos(.5*lp.phi) + cos(lp.phi));
27 if (fabs(th1) < EPS) break;
28 }
29 xy.x = FXC * lp.lam * (1.0 + 2. * cos(lp.phi)/cos(0.5 * lp.phi));
30 xy.y = FYC * sin(0.5 * lp.phi);
31 return xy;
32 }
33
34
s_inverse(XY xy,PJ * P)35 static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
36 LP lp = {0.0,0.0};
37 double t;
38
39 lp.phi = RYC * xy.y;
40 if (fabs(lp.phi) > 1.) {
41 if (fabs(lp.phi) > ONETOL) I_ERROR
42 else if (lp.phi < 0.) { t = -1.; lp.phi = -M_PI; }
43 else { t = 1.; lp.phi = M_PI; }
44 } else
45 lp.phi = 2. * asin(t = lp.phi);
46 lp.lam = RXC * xy.x / (1. + 2. * cos(lp.phi)/cos(0.5 * lp.phi));
47 lp.phi = RC * (t + sin(lp.phi));
48 if (fabs(lp.phi) > 1.)
49 if (fabs(lp.phi) > ONETOL) I_ERROR
50 else lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
51 else
52 lp.phi = asin(lp.phi);
53 return lp;
54 }
55
56
freeup_new(PJ * P)57 static void *freeup_new (PJ *P) { /* Destructor */
58 if (0==P)
59 return 0;
60
61 return pj_dealloc(P);
62 }
63
freeup(PJ * P)64 static void freeup (PJ *P) {
65 freeup_new (P);
66 return;
67 }
68
69
PROJECTION(mbtfpq)70 PJ *PROJECTION(mbtfpq) {
71
72 P->es = 0.;
73 P->inv = s_inverse;
74 P->fwd = s_forward;
75
76 return P;
77 }
78
79 #ifndef PJ_SELFTEST
pj_mbtfpq_selftest(void)80 int pj_mbtfpq_selftest (void) {return 0;}
81 #else
82
pj_mbtfpq_selftest(void)83 int pj_mbtfpq_selftest (void) {
84 double tolerance_lp = 1e-10;
85 double tolerance_xy = 1e-7;
86
87 char s_args[] = {"+proj=mbtfpq +a=6400000 +lat_1=0.5 +lat_2=2"};
88
89 LP fwd_in[] = {
90 { 2, 1},
91 { 2,-1},
92 {-2, 1},
93 {-2,-1}
94 };
95
96 XY s_fwd_expect[] = {
97 { 209391.854738393013, 119161.040199054827},
98 { 209391.854738393013, -119161.040199054827},
99 {-209391.854738393013, 119161.040199054827},
100 {-209391.854738393013, -119161.040199054827},
101 };
102
103 XY inv_in[] = {
104 { 200, 100},
105 { 200,-100},
106 {-200, 100},
107 {-200,-100}
108 };
109
110 LP s_inv_expect[] = {
111 { 0.00191010555824111571, 0.000839185447792341723},
112 { 0.00191010555824111571, -0.000839185447792341723},
113 {-0.00191010555824111571, 0.000839185447792341723},
114 {-0.00191010555824111571, -0.000839185447792341723},
115 };
116
117 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);
118 }
119
120
121 #endif
122