1 #define PJ_LIB__
2 #include <projects.h>
3
4 PROJ_HEAD(mbt_fps, "McBryde-Thomas Flat-Pole Sine (No. 2)") "\n\tCyl., Sph.";
5
6 #define MAX_ITER 10
7 #define LOOP_TOL 1e-7
8 #define C1 0.45503
9 #define C2 1.36509
10 #define C3 1.41546
11 #define C_x 0.22248
12 #define C_y 1.44492
13 #define C1_2 0.33333333333333333333333333
14
s_forward(LP lp,PJ * P)15 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
16 XY xy = {0.0,0.0};
17 double k, V, t;
18 int i;
19 (void) P;
20
21 k = C3 * sin(lp.phi);
22 for (i = MAX_ITER; i ; --i) {
23 t = lp.phi / C2;
24 lp.phi -= V = (C1 * sin(t) + sin(lp.phi) - k) /
25 (C1_2 * cos(t) + cos(lp.phi));
26 if (fabs(V) < LOOP_TOL)
27 break;
28 }
29 t = lp.phi / C2;
30 xy.x = C_x * lp.lam * (1. + 3. * cos(lp.phi)/cos(t) );
31 xy.y = C_y * sin(t);
32 return xy;
33 }
34
35
s_inverse(XY xy,PJ * P)36 static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
37 LP lp = {0.0,0.0};
38 double t;
39
40 lp.phi = C2 * (t = aasin(P->ctx,xy.y / C_y));
41 lp.lam = xy.x / (C_x * (1. + 3. * cos(lp.phi)/cos(t)));
42 lp.phi = aasin(P->ctx,(C1 * sin(t) + sin(lp.phi)) / C3);
43 return (lp);
44 }
45
46
freeup_new(PJ * P)47 static void *freeup_new (PJ *P) { /* Destructor */
48 if (0==P)
49 return 0;
50
51 return pj_dealloc(P);
52 }
53
54
freeup(PJ * P)55 static void freeup (PJ *P) {
56 freeup_new (P);
57 return;
58 }
59
60
PROJECTION(mbt_fps)61 PJ *PROJECTION(mbt_fps) {
62
63 P->es = 0;
64 P->inv = s_inverse;
65 P->fwd = s_forward;
66
67 return P;
68 }
69
70 #ifndef PJ_SELFTEST
pj_mbt_fps_selftest(void)71 int pj_mbt_fps_selftest (void) {return 0;}
72 #else
73
pj_mbt_fps_selftest(void)74 int pj_mbt_fps_selftest (void) {
75 double tolerance_lp = 1e-10;
76 double tolerance_xy = 1e-7;
77
78 char s_args[] = {"+proj=mbt_fps +a=6400000 +lat_1=0.5 +lat_2=2"};
79
80 LP fwd_in[] = {
81 { 2, 1},
82 { 2,-1},
83 {-2, 1},
84 {-2,-1}
85 };
86
87 XY s_fwd_expect[] = {
88 { 198798.176129849948, 125512.017254530627},
89 { 198798.176129849948, -125512.017254530627},
90 {-198798.176129849948, 125512.017254530627},
91 {-198798.176129849948, -125512.017254530627},
92 };
93
94 XY inv_in[] = {
95 { 200, 100},
96 { 200,-100},
97 {-200, 100},
98 {-200,-100}
99 };
100
101 LP s_inv_expect[] = {
102 { 0.00201197086238270742, 0.000796711850174446003},
103 { 0.00201197086238270742, -0.000796711850174446003},
104 {-0.00201197086238270742, 0.000796711850174446003},
105 {-0.00201197086238270742, -0.000796711850174446003},
106 };
107
108 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);
109 }
110
111
112 #endif
113