Coverage Report

Created: 2025-06-13 06:29

/src/proj/src/projections/sterea.cpp
Line
Count
Source (jump to first uncovered line)
1
/*
2
** libproj -- library of cartographic projections
3
**
4
** Copyright (c) 2003   Gerald I. Evenden
5
*/
6
/*
7
** Permission is hereby granted, free of charge, to any person obtaining
8
** a copy of this software and associated documentation files (the
9
** "Software"), to deal in the Software without restriction, including
10
** without limitation the rights to use, copy, modify, merge, publish,
11
** distribute, sublicense, and/or sell copies of the Software, and to
12
** permit persons to whom the Software is furnished to do so, subject to
13
** the following conditions:
14
**
15
** The above copyright notice and this permission notice shall be
16
** included in all copies or substantial portions of the Software.
17
**
18
** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
19
** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20
** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
21
** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
22
** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23
** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24
** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25
*/
26
27
#include "proj.h"
28
#include "proj_internal.h"
29
#include <errno.h>
30
#include <math.h>
31
32
namespace { // anonymous namespace
33
struct pj_opaque {
34
    double phic0;
35
    double cosc0, sinc0;
36
    double R2;
37
    void *en;
38
};
39
} // anonymous namespace
40
41
PROJ_HEAD(sterea, "Oblique Stereographic Alternative") "\n\tAzimuthal, Sph&Ell";
42
43
0
static PJ_XY sterea_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
44
0
    PJ_XY xy = {0.0, 0.0};
45
0
    struct pj_opaque *Q = static_cast<struct pj_opaque *>(P->opaque);
46
0
    double cosc, sinc, cosl, k;
47
48
0
    lp = pj_gauss(P->ctx, lp, Q->en);
49
0
    sinc = sin(lp.phi);
50
0
    cosc = cos(lp.phi);
51
0
    cosl = cos(lp.lam);
52
0
    const double denom = 1. + Q->sinc0 * sinc + Q->cosc0 * cosc * cosl;
53
0
    if (denom == 0.0) {
54
0
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
55
0
        return proj_coord_error().xy;
56
0
    }
57
0
    k = P->k0 * Q->R2 / denom;
58
0
    xy.x = k * cosc * sin(lp.lam);
59
0
    xy.y = k * (Q->cosc0 * sinc - Q->sinc0 * cosc * cosl);
60
0
    return xy;
61
0
}
62
63
0
static PJ_LP sterea_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
64
0
    PJ_LP lp = {0.0, 0.0};
65
0
    struct pj_opaque *Q = static_cast<struct pj_opaque *>(P->opaque);
66
0
    double rho, c, sinc, cosc;
67
68
0
    xy.x /= P->k0;
69
0
    xy.y /= P->k0;
70
0
    if ((rho = hypot(xy.x, xy.y)) != 0.0) {
71
0
        c = 2. * atan2(rho, Q->R2);
72
0
        sinc = sin(c);
73
0
        cosc = cos(c);
74
0
        lp.phi = asin(cosc * Q->sinc0 + xy.y * sinc * Q->cosc0 / rho);
75
0
        lp.lam =
76
0
            atan2(xy.x * sinc, rho * Q->cosc0 * cosc - xy.y * Q->sinc0 * sinc);
77
0
    } else {
78
0
        lp.phi = Q->phic0;
79
0
        lp.lam = 0.;
80
0
    }
81
0
    return pj_inv_gauss(P->ctx, lp, Q->en);
82
0
}
83
84
53
static PJ *destructor(PJ *P, int errlev) {
85
53
    if (nullptr == P)
86
0
        return nullptr;
87
88
53
    if (nullptr == P->opaque)
89
0
        return pj_default_destructor(P, errlev);
90
91
53
    free(static_cast<struct pj_opaque *>(P->opaque)->en);
92
53
    return pj_default_destructor(P, errlev);
93
53
}
94
95
53
PJ *PJ_PROJECTION(sterea) {
96
53
    double R;
97
53
    struct pj_opaque *Q =
98
53
        static_cast<struct pj_opaque *>(calloc(1, sizeof(struct pj_opaque)));
99
100
53
    if (nullptr == Q)
101
0
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
102
53
    P->opaque = Q;
103
104
53
    Q->en = pj_gauss_ini(P->e, P->phi0, &(Q->phic0), &R);
105
53
    if (nullptr == Q->en)
106
0
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
107
108
53
    Q->sinc0 = sin(Q->phic0);
109
53
    Q->cosc0 = cos(Q->phic0);
110
53
    Q->R2 = 2. * R;
111
112
53
    P->inv = sterea_e_inverse;
113
53
    P->fwd = sterea_e_forward;
114
53
    P->destructor = destructor;
115
116
53
    return P;
117
53
}