/src/proj/src/projections/hatano.cpp
Line | Count | Source |
1 | | |
2 | | |
3 | | #include <math.h> |
4 | | |
5 | | #include "proj.h" |
6 | | #include "proj_internal.h" |
7 | | |
8 | | PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph"; |
9 | | |
10 | 0 | #define NITER 20 |
11 | 0 | #define EPS 1e-7 |
12 | 0 | #define ONETOL 1.000001 |
13 | 0 | #define CN 2.67595 |
14 | 0 | #define CSz 2.43763 |
15 | 0 | #define RCN 0.37369906014686373063 |
16 | 0 | #define RCS 0.41023453108141924738 |
17 | 0 | #define FYCN 1.75859 |
18 | 0 | #define FYCS 1.93052 |
19 | 0 | #define RYCN 0.56863737426006061674 |
20 | 0 | #define RYCS 0.51799515156538134803 |
21 | 0 | #define FXC 0.85 |
22 | 0 | #define RXC 1.17647058823529411764 |
23 | | |
24 | 0 | static PJ_XY hatano_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ |
25 | 0 | PJ_XY xy = {0.0, 0.0}; |
26 | 0 | int i; |
27 | 0 | (void)P; |
28 | |
|
29 | 0 | const double c = sin(lp.phi) * (lp.phi < 0. ? CSz : CN); |
30 | 0 | for (i = NITER; i; --i) { |
31 | 0 | const double th1 = (lp.phi + sin(lp.phi) - c) / (1. + cos(lp.phi)); |
32 | 0 | lp.phi -= th1; |
33 | 0 | if (fabs(th1) < EPS) |
34 | 0 | break; |
35 | 0 | } |
36 | 0 | xy.x = FXC * lp.lam * cos(lp.phi *= .5); |
37 | 0 | xy.y = sin(lp.phi) * (lp.phi < 0. ? FYCS : FYCN); |
38 | |
|
39 | 0 | return xy; |
40 | 0 | } |
41 | | |
42 | 0 | static PJ_LP hatano_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ |
43 | 0 | PJ_LP lp = {0.0, 0.0}; |
44 | 0 | double th; |
45 | |
|
46 | 0 | th = xy.y * (xy.y < 0. ? RYCS : RYCN); |
47 | 0 | if (fabs(th) > 1.) { |
48 | 0 | if (fabs(th) > ONETOL) { |
49 | 0 | proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); |
50 | 0 | return lp; |
51 | 0 | } else { |
52 | 0 | th = th > 0. ? M_HALFPI : -M_HALFPI; |
53 | 0 | } |
54 | 0 | } else { |
55 | 0 | th = asin(th); |
56 | 0 | } |
57 | | |
58 | 0 | lp.lam = RXC * xy.x / cos(th); |
59 | 0 | th += th; |
60 | 0 | lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN); |
61 | 0 | if (fabs(lp.phi) > 1.) { |
62 | 0 | if (fabs(lp.phi) > ONETOL) { |
63 | 0 | proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); |
64 | 0 | return lp; |
65 | 0 | } else { |
66 | 0 | lp.phi = lp.phi > 0. ? M_HALFPI : -M_HALFPI; |
67 | 0 | } |
68 | 0 | } else { |
69 | 0 | lp.phi = asin(lp.phi); |
70 | 0 | } |
71 | | |
72 | 0 | return (lp); |
73 | 0 | } |
74 | | |
75 | 3 | PJ *PJ_PROJECTION(hatano) { |
76 | 3 | P->es = 0.; |
77 | 3 | P->inv = hatano_s_inverse; |
78 | 3 | P->fwd = hatano_s_forward; |
79 | | |
80 | 3 | return P; |
81 | 3 | } |
82 | | |
83 | | #undef NITER |
84 | | #undef EPS |
85 | | #undef ONETOL |
86 | | #undef CN |
87 | | #undef CSz |
88 | | #undef RCN |
89 | | #undef RCS |
90 | | #undef FYCN |
91 | | #undef FYCS |
92 | | #undef RYCN |
93 | | #undef RYCS |
94 | | #undef FXC |
95 | | #undef RXC |