/src/PROJ/src/projections/gnom.cpp

Source (jump to first uncovered line)


#include <errno.h>
#include <float.h>
#include <math.h>

#include "geodesic.h"
#include "proj.h"
#include "proj_internal.h"

PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph";

#define EPS10 1.e-10

namespace pj_gnom_ns {
enum Mode { N_POLE = 0, S_POLE = 1, EQUIT = 2, OBLIQ = 3 };
}

namespace { // anonymous namespace
struct pj_gnom_data {
    double sinph0;
    double cosph0;
    enum pj_gnom_ns::Mode mode;
    struct geod_geodesic g;
};
} // anonymous namespace

static PJ_XY gnom_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_gnom_data *Q = static_cast<struct pj_gnom_data *>(P->opaque);
    double coslam, cosphi, sinphi;

    sinphi = sin(lp.phi);
    cosphi = cos(lp.phi);
    coslam = cos(lp.lam);

    switch (Q->mode) {
    case pj_gnom_ns::EQUIT:
        xy.y = cosphi * coslam;
        break;
    case pj_gnom_ns::OBLIQ:
        xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam;
        break;
    case pj_gnom_ns::S_POLE:
        xy.y = -sinphi;
        break;
    case pj_gnom_ns::N_POLE:
        xy.y = sinphi;
        break;
    }

    if (xy.y <= EPS10) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return xy;
    }

    xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam);
    switch (Q->mode) {
    case pj_gnom_ns::EQUIT:
        xy.y *= sinphi;
        break;
    case pj_gnom_ns::OBLIQ:
        xy.y *= Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam;
        break;
    case pj_gnom_ns::N_POLE:
        coslam = -coslam;
        PROJ_FALLTHROUGH;
    case pj_gnom_ns::S_POLE:
        xy.y *= cosphi * coslam;
        break;
    }
    return xy;
}

static PJ_LP gnom_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_gnom_data *Q = static_cast<struct pj_gnom_data *>(P->opaque);
    double rh, cosz, sinz;

    rh = hypot(xy.x, xy.y);
    sinz = sin(lp.phi = atan(rh));
    cosz = sqrt(1. - sinz * sinz);

    if (fabs(rh) <= EPS10) {
        lp.phi = P->phi0;
        lp.lam = 0.;
    } else {
        switch (Q->mode) {
        case pj_gnom_ns::OBLIQ:
            lp.phi = cosz * Q->sinph0 + xy.y * sinz * Q->cosph0 / rh;
            if (fabs(lp.phi) >= 1.)
                lp.phi = lp.phi > 0. ? M_HALFPI : -M_HALFPI;
            else
                lp.phi = asin(lp.phi);
            xy.y = (cosz - Q->sinph0 * sin(lp.phi)) * rh;
            xy.x *= sinz * Q->cosph0;
            break;
        case pj_gnom_ns::EQUIT:
            lp.phi = xy.y * sinz / rh;
            if (fabs(lp.phi) >= 1.)
                lp.phi = lp.phi > 0. ? M_HALFPI : -M_HALFPI;
            else
                lp.phi = asin(lp.phi);
            xy.y = cosz * rh;
            xy.x *= sinz;
            break;
        case pj_gnom_ns::S_POLE:
            lp.phi -= M_HALFPI;
            break;
        case pj_gnom_ns::N_POLE:
            lp.phi = M_HALFPI - lp.phi;
            xy.y = -xy.y;
            break;
        }
        lp.lam = atan2(xy.x, xy.y);
    }
    return lp;
}

static PJ_XY gnom_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_gnom_data *Q = static_cast<struct pj_gnom_data *>(P->opaque);

    double lat0 = P->phi0 / DEG_TO_RAD, lon0 = 0, lat1 = lp.phi / DEG_TO_RAD,
           lon1 = lp.lam / DEG_TO_RAD, azi0, m, M;

    geod_geninverse(&Q->g, lat0, lon0, lat1, lon1, nullptr, &azi0, nullptr, &m,
                    &M, nullptr, nullptr);
    if (M <= 0) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        xy.x = xy.y = HUGE_VAL;
    } else {
        double rho = m / M;
        azi0 *= DEG_TO_RAD;
        xy.x = rho * sin(azi0);
        xy.y = rho * cos(azi0);
    }
    return xy;
}

static PJ_LP gnom_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    constexpr int numit_ = 10;
    static const double eps_ = 0.01 * sqrt(DBL_EPSILON);

    PJ_LP lp = {0.0, 0.0};
    struct pj_gnom_data *Q = static_cast<struct pj_gnom_data *>(P->opaque);

    double lat0 = P->phi0 / DEG_TO_RAD, lon0 = 0,
           azi0 = atan2(xy.x, xy.y) / DEG_TO_RAD, // Clockwise from north
        rho = hypot(xy.x, xy.y), s = atan(rho);
    bool little = rho <= 1;
    if (!little)
        rho = 1 / rho;
    struct geod_geodesicline l;
    geod_lineinit(&l, &Q->g, lat0, lon0, azi0,
                  GEOD_LATITUDE | GEOD_LONGITUDE | GEOD_DISTANCE_IN |
                      GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE);

    int count = numit_, trip = 0;
    double lat1 = 0, lon1 = 0;
    while (count--) {
        double m, M;
        geod_genposition(&l, 0, s, &lat1, &lon1, nullptr, &s, &m, &M, nullptr,
                         nullptr);
        if (trip)
            break;
        // If little, solve rho(s) = rho with drho(s)/ds = 1/M^2
        // else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2
        double ds = little ? (m - rho * M) * M : (rho * m - M) * m;
        s -= ds;
        // Reversed test to allow escape with NaNs
        if (!(fabs(ds) >= eps_))
            ++trip;
    }
    if (trip) {
        lp.phi = lat1 * DEG_TO_RAD;
        lp.lam = lon1 * DEG_TO_RAD;
    } else {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        lp.phi = lp.lam = HUGE_VAL;
    }
    return lp;
}

PJ *PJ_PROJECTION(gnom) {
    struct pj_gnom_data *Q = static_cast<struct pj_gnom_data *>(
        calloc(1, sizeof(struct pj_gnom_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    if (P->es == 0) {
        if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) {
            Q->mode = P->phi0 < 0. ? pj_gnom_ns::S_POLE : pj_gnom_ns::N_POLE;
        } else if (fabs(P->phi0) < EPS10) {
            Q->mode = pj_gnom_ns::EQUIT;
        } else {
            Q->mode = pj_gnom_ns::OBLIQ;
            Q->sinph0 = sin(P->phi0);
            Q->cosph0 = cos(P->phi0);
        }

        P->inv = gnom_s_inverse;
        P->fwd = gnom_s_forward;
    } else {
        geod_init(&Q->g, 1, P->f);

        P->inv = gnom_e_inverse;
        P->fwd = gnom_e_forward;
    }
    P->es = 0.;

    return P;
}

Coverage Report

Created: 2025-07-23 06:58

Line	Count	Source (jump to first uncovered line)
1
2
3		#include <errno.h>
4		#include <float.h>
5		#include <math.h>
6
7		#include "geodesic.h"
8		#include "proj.h"
9		#include "proj_internal.h"
10
11		PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph";
12
13	33	#define EPS10 1.e-10
14
15		namespace pj_gnom_ns {
16		enum Mode { N_POLE = 0, S_POLE = 1, EQUIT = 2, OBLIQ = 3 };
17		}
18
19		namespace { // anonymous namespace
20		struct pj_gnom_data {
21		double sinph0;
22		double cosph0;
23		enum pj_gnom_ns::Mode mode;
24		struct geod_geodesic g;
25		};
26		} // anonymous namespace
27
28	0	static PJ_XY gnom_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
29	0	PJ_XY xy = {0.0, 0.0};
30	0	struct pj_gnom_data Q = static_cast<struct pj_gnom_data >(P->opaque);
31	0	double coslam, cosphi, sinphi;
32
33	0	sinphi = sin(lp.phi);
34	0	cosphi = cos(lp.phi);
35	0	coslam = cos(lp.lam);
36
37	0	switch (Q->mode) {
38	0	case pj_gnom_ns::EQUIT:
39	0	xy.y = cosphi * coslam;
40	0	break;
41	0	case pj_gnom_ns::OBLIQ:
42	0	xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam;
43	0	break;
44	0	case pj_gnom_ns::S_POLE:
45	0	xy.y = -sinphi;
46	0	break;
47	0	case pj_gnom_ns::N_POLE:
48	0	xy.y = sinphi;
49	0	break;
50	0	}
51
52	0	if (xy.y <= EPS10) {
53	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
54	0	return xy;
55	0	}
56
57	0	xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam);
58	0	switch (Q->mode) {
59	0	case pj_gnom_ns::EQUIT:
60	0	xy.y *= sinphi;
61	0	break;
62	0	case pj_gnom_ns::OBLIQ:
63	0	xy.y = Q->cosph0 sinphi - Q->sinph0 * cosphi * coslam;
64	0	break;
65	0	case pj_gnom_ns::N_POLE:
66	0	coslam = -coslam;
67	0	PROJ_FALLTHROUGH;
68	0	case pj_gnom_ns::S_POLE:
69	0	xy.y = cosphi coslam;
70	0	break;
71	0	}
72	0	return xy;
73	0	}
74
75	0	static PJ_LP gnom_s_inverse(PJ_XY xy, PJ P) { / Spheroidal, inverse */
76	0	PJ_LP lp = {0.0, 0.0};
77	0	struct pj_gnom_data Q = static_cast<struct pj_gnom_data >(P->opaque);
78	0	double rh, cosz, sinz;
79
80	0	rh = hypot(xy.x, xy.y);
81	0	sinz = sin(lp.phi = atan(rh));
82	0	cosz = sqrt(1. - sinz * sinz);
83
84	0	if (fabs(rh) <= EPS10) {
85	0	lp.phi = P->phi0;
86	0	lp.lam = 0.;
87	0	} else {
88	0	switch (Q->mode) {
89	0	case pj_gnom_ns::OBLIQ:
90	0	lp.phi = cosz * Q->sinph0 + xy.y * sinz * Q->cosph0 / rh;
91	0	if (fabs(lp.phi) >= 1.)
92	0	lp.phi = lp.phi > 0. ? M_HALFPI : -M_HALFPI;
93	0	else
94	0	lp.phi = asin(lp.phi);
95	0	xy.y = (cosz - Q->sinph0 * sin(lp.phi)) * rh;
96	0	xy.x = sinz Q->cosph0;
97	0	break;
98	0	case pj_gnom_ns::EQUIT:
99	0	lp.phi = xy.y * sinz / rh;
100	0	if (fabs(lp.phi) >= 1.)
101	0	lp.phi = lp.phi > 0. ? M_HALFPI : -M_HALFPI;
102	0	else
103	0	lp.phi = asin(lp.phi);
104	0	xy.y = cosz * rh;
105	0	xy.x *= sinz;
106	0	break;
107	0	case pj_gnom_ns::S_POLE:
108	0	lp.phi -= M_HALFPI;
109	0	break;
110	0	case pj_gnom_ns::N_POLE:
111	0	lp.phi = M_HALFPI - lp.phi;
112	0	xy.y = -xy.y;
113	0	break;
114	0	}
115	0	lp.lam = atan2(xy.x, xy.y);
116	0	}
117	0	return lp;
118	0	}
119
120	15.3k	static PJ_XY gnom_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
121	15.3k	PJ_XY xy = {0.0, 0.0};
122	15.3k	struct pj_gnom_data Q = static_cast<struct pj_gnom_data >(P->opaque);
123
124	15.3k	double lat0 = P->phi0 / DEG_TO_RAD, lon0 = 0, lat1 = lp.phi / DEG_TO_RAD,
125	15.3k	lon1 = lp.lam / DEG_TO_RAD, azi0, m, M;
126
127	15.3k	geod_geninverse(&Q->g, lat0, lon0, lat1, lon1, nullptr, &azi0, nullptr, &m,
128	15.3k	&M, nullptr, nullptr);
129	15.3k	if (M <= 0) {
130	6.78k	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
131	6.78k	xy.x = xy.y = HUGE_VAL;
132	8.58k	} else {
133	8.58k	double rho = m / M;
134	8.58k	azi0 *= DEG_TO_RAD;
135	8.58k	xy.x = rho * sin(azi0);
136	8.58k	xy.y = rho * cos(azi0);
137	8.58k	}
138	15.3k	return xy;
139	15.3k	}
140
141	0	static PJ_LP gnom_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
142	0	constexpr int numit_ = 10;
143	0	static const double eps_ = 0.01 * sqrt(DBL_EPSILON);
144
145	0	PJ_LP lp = {0.0, 0.0};
146	0	struct pj_gnom_data Q = static_cast<struct pj_gnom_data >(P->opaque);
147
148	0	double lat0 = P->phi0 / DEG_TO_RAD, lon0 = 0,
149	0	azi0 = atan2(xy.x, xy.y) / DEG_TO_RAD, // Clockwise from north
150	0	rho = hypot(xy.x, xy.y), s = atan(rho);
151	0	bool little = rho <= 1;
152	0	if (!little)
153	0	rho = 1 / rho;
154	0	struct geod_geodesicline l;
155	0	geod_lineinit(&l, &Q->g, lat0, lon0, azi0,
156	0	GEOD_LATITUDE \| GEOD_LONGITUDE \| GEOD_DISTANCE_IN \|
157	0	GEOD_REDUCEDLENGTH \| GEOD_GEODESICSCALE);
158
159	0	int count = numit_, trip = 0;
160	0	double lat1 = 0, lon1 = 0;
161	0	while (count--) {
162	0	double m, M;
163	0	geod_genposition(&l, 0, s, &lat1, &lon1, nullptr, &s, &m, &M, nullptr,
164	0	nullptr);
165	0	if (trip)
166	0	break;
167		// If little, solve rho(s) = rho with drho(s)/ds = 1/M^2
168		// else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2
169	0	double ds = little ? (m - rho * M) * M : (rho * m - M) * m;
170	0	s -= ds;
171		// Reversed test to allow escape with NaNs
172	0	if (!(fabs(ds) >= eps_))
173	0	++trip;
174	0	}
175	0	if (trip) {
176	0	lp.phi = lat1 * DEG_TO_RAD;
177	0	lp.lam = lon1 * DEG_TO_RAD;
178	0	} else {
179	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
180	0	lp.phi = lp.lam = HUGE_VAL;
181	0	}
182	0	return lp;
183	0	}
184
185	265	PJ *PJ_PROJECTION(gnom) {
186	265	struct pj_gnom_data Q = static_cast<struct pj_gnom_data >(
187	265	calloc(1, sizeof(struct pj_gnom_data)));
188	265	if (nullptr == Q)
189	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
190	265	P->opaque = Q;
191
192	265	if (P->es == 0) {
193	18	if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) {
194	3	Q->mode = P->phi0 < 0. ? pj_gnom_ns::S_POLE : pj_gnom_ns::N_POLE;
195	15	} else if (fabs(P->phi0) < EPS10) {
196	9	Q->mode = pj_gnom_ns::EQUIT;
197	9	} else {
198	6	Q->mode = pj_gnom_ns::OBLIQ;
199	6	Q->sinph0 = sin(P->phi0);
200	6	Q->cosph0 = cos(P->phi0);
201	6	}
202
203	18	P->inv = gnom_s_inverse;
204	18	P->fwd = gnom_s_forward;
205	247	} else {
206	247	geod_init(&Q->g, 1, P->f);
207
208	247	P->inv = gnom_e_inverse;
209	247	P->fwd = gnom_e_forward;
210	247	}
211	265	P->es = 0.;
212
213	265	return P;
214	265	}