/src/PROJ/src/projections/mbtfpq.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | |
2 | | |
3 | | #include <math.h> |
4 | | |
5 | | #include "proj.h" |
6 | | #include "proj_internal.h" |
7 | | |
8 | | PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic") "\n\tCyl, Sph"; |
9 | | |
10 | 0 | #define NITER 20 |
11 | 0 | #define EPS 1e-7 |
12 | 0 | #define ONETOL 1.000001 |
13 | 0 | #define C 1.70710678118654752440 |
14 | 0 | #define RC 0.58578643762690495119 |
15 | 0 | #define FYC 1.87475828462269495505 |
16 | 0 | #define RYC 0.53340209679417701685 |
17 | 0 | #define FXC 0.31245971410378249250 |
18 | 0 | #define RXC 3.20041258076506210122 |
19 | | |
20 | 0 | static PJ_XY mbtfpq_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ |
21 | 0 | PJ_XY xy = {0.0, 0.0}; |
22 | 0 | (void)P; |
23 | |
|
24 | 0 | const double c = C * sin(lp.phi); |
25 | 0 | for (int i = NITER; i; --i) { |
26 | 0 | const double th1 = (sin(.5 * lp.phi) + sin(lp.phi) - c) / |
27 | 0 | (.5 * cos(.5 * lp.phi) + cos(lp.phi)); |
28 | 0 | lp.phi -= th1; |
29 | 0 | if (fabs(th1) < EPS) |
30 | 0 | break; |
31 | 0 | } |
32 | 0 | xy.x = FXC * lp.lam * (1.0 + 2. * cos(lp.phi) / cos(0.5 * lp.phi)); |
33 | 0 | xy.y = FYC * sin(0.5 * lp.phi); |
34 | 0 | return xy; |
35 | 0 | } |
36 | | |
37 | 0 | static PJ_LP mbtfpq_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ |
38 | 0 | PJ_LP lp = {0.0, 0.0}; |
39 | 0 | double t; |
40 | |
|
41 | 0 | lp.phi = RYC * xy.y; |
42 | 0 | if (fabs(lp.phi) > 1.) { |
43 | 0 | if (fabs(lp.phi) > ONETOL) { |
44 | 0 | proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); |
45 | 0 | return lp; |
46 | 0 | } else if (lp.phi < 0.) { |
47 | 0 | t = -1.; |
48 | 0 | lp.phi = -M_PI; |
49 | 0 | } else { |
50 | 0 | t = 1.; |
51 | 0 | lp.phi = M_PI; |
52 | 0 | } |
53 | 0 | } else { |
54 | 0 | t = lp.phi; |
55 | 0 | lp.phi = 2. * asin(lp.phi); |
56 | 0 | } |
57 | 0 | lp.lam = RXC * xy.x / (1. + 2. * cos(lp.phi) / cos(0.5 * lp.phi)); |
58 | 0 | lp.phi = RC * (t + sin(lp.phi)); |
59 | 0 | if (fabs(lp.phi) > 1.) |
60 | 0 | if (fabs(lp.phi) > ONETOL) { |
61 | 0 | proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); |
62 | 0 | return lp; |
63 | 0 | } else |
64 | 0 | lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI; |
65 | 0 | else |
66 | 0 | lp.phi = asin(lp.phi); |
67 | 0 | return lp; |
68 | 0 | } |
69 | | |
70 | 0 | PJ *PJ_PROJECTION(mbtfpq) { |
71 | |
|
72 | 0 | P->es = 0.; |
73 | 0 | P->inv = mbtfpq_s_inverse; |
74 | 0 | P->fwd = mbtfpq_s_forward; |
75 | |
|
76 | 0 | return P; |
77 | 0 | } |
78 | | |
79 | | #undef NITER |
80 | | #undef EPS |
81 | | #undef ONETOL |
82 | | #undef C |
83 | | #undef RC |
84 | | #undef FYC |
85 | | #undef RYC |
86 | | #undef FXC |
87 | | #undef RXC |