/src/proj/src/projections/putp6.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | |
2 | | |
3 | | #include <errno.h> |
4 | | #include <math.h> |
5 | | |
6 | | #include "proj.h" |
7 | | #include "proj_internal.h" |
8 | | |
9 | | namespace { // anonymous namespace |
10 | | struct pj_putp6 { |
11 | | double C_x, C_y, A, B, D; |
12 | | }; |
13 | | } // anonymous namespace |
14 | | |
15 | | PROJ_HEAD(putp6, "Putnins P6") "\n\tPCyl, Sph"; |
16 | | PROJ_HEAD(putp6p, "Putnins P6'") "\n\tPCyl, Sph"; |
17 | | |
18 | 0 | #define EPS 1e-10 |
19 | 0 | #define NITER 10 |
20 | 0 | #define CON_POLE 1.732050807568877 /* sqrt(3) */ |
21 | | |
22 | 0 | static PJ_XY putp6_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ |
23 | 0 | PJ_XY xy = {0.0, 0.0}; |
24 | 0 | struct pj_putp6 *Q = static_cast<struct pj_putp6 *>(P->opaque); |
25 | 0 | int i; |
26 | |
|
27 | 0 | const double p = Q->B * sin(lp.phi); |
28 | 0 | lp.phi *= 1.10265779; |
29 | 0 | for (i = NITER; i; --i) { |
30 | 0 | const double r = sqrt(1. + lp.phi * lp.phi); |
31 | 0 | const double V = |
32 | 0 | ((Q->A - r) * lp.phi - log(lp.phi + r) - p) / (Q->A - 2. * r); |
33 | 0 | lp.phi -= V; |
34 | 0 | if (fabs(V) < EPS) |
35 | 0 | break; |
36 | 0 | } |
37 | 0 | double sqrt_1_plus_phi2; |
38 | 0 | if (!i) { |
39 | | // Note: it seems this case is rarely reached as from experimenting, |
40 | | // i seems to be >= 6 |
41 | 0 | lp.phi = p < 0. ? -CON_POLE : CON_POLE; |
42 | | // the formula of the else case would also work, but this makes |
43 | | // some cppcheck versions happier. |
44 | 0 | sqrt_1_plus_phi2 = 2; |
45 | 0 | } else { |
46 | 0 | sqrt_1_plus_phi2 = sqrt(1. + lp.phi * lp.phi); |
47 | 0 | } |
48 | 0 | xy.x = Q->C_x * lp.lam * (Q->D - sqrt_1_plus_phi2); |
49 | 0 | xy.y = Q->C_y * lp.phi; |
50 | |
|
51 | 0 | return xy; |
52 | 0 | } |
53 | | |
54 | 0 | static PJ_LP putp6_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ |
55 | 0 | PJ_LP lp = {0.0, 0.0}; |
56 | 0 | struct pj_putp6 *Q = static_cast<struct pj_putp6 *>(P->opaque); |
57 | 0 | double r; |
58 | |
|
59 | 0 | lp.phi = xy.y / Q->C_y; |
60 | 0 | r = sqrt(1. + lp.phi * lp.phi); |
61 | 0 | lp.lam = xy.x / (Q->C_x * (Q->D - r)); |
62 | 0 | lp.phi = aasin(P->ctx, ((Q->A - r) * lp.phi - log(lp.phi + r)) / Q->B); |
63 | |
|
64 | 0 | return lp; |
65 | 0 | } |
66 | | |
67 | 94 | PJ *PJ_PROJECTION(putp6) { |
68 | 94 | struct pj_putp6 *Q = |
69 | 94 | static_cast<struct pj_putp6 *>(calloc(1, sizeof(struct pj_putp6))); |
70 | 94 | if (nullptr == Q) |
71 | 0 | return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
72 | 94 | P->opaque = Q; |
73 | | |
74 | 94 | Q->C_x = 1.01346; |
75 | 94 | Q->C_y = 0.91910; |
76 | 94 | Q->A = 4.; |
77 | 94 | Q->B = 2.1471437182129378784; |
78 | 94 | Q->D = 2.; |
79 | | |
80 | 94 | P->es = 0.; |
81 | 94 | P->inv = putp6_s_inverse; |
82 | 94 | P->fwd = putp6_s_forward; |
83 | | |
84 | 94 | return P; |
85 | 94 | } |
86 | | |
87 | 12 | PJ *PJ_PROJECTION(putp6p) { |
88 | 12 | struct pj_putp6 *Q = |
89 | 12 | static_cast<struct pj_putp6 *>(calloc(1, sizeof(struct pj_putp6))); |
90 | 12 | if (nullptr == Q) |
91 | 0 | return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
92 | 12 | P->opaque = Q; |
93 | | |
94 | 12 | Q->C_x = 0.44329; |
95 | 12 | Q->C_y = 0.80404; |
96 | 12 | Q->A = 6.; |
97 | 12 | Q->B = 5.61125; |
98 | 12 | Q->D = 3.; |
99 | | |
100 | 12 | P->es = 0.; |
101 | 12 | P->inv = putp6_s_inverse; |
102 | 12 | P->fwd = putp6_s_forward; |
103 | | |
104 | 12 | return P; |
105 | 12 | } |
106 | | |
107 | | #undef EPS |
108 | | #undef NITER |
109 | | #undef CON_POLE |