/src/PROJ/src/projections/aeqd.cpp

Source (jump to first uncovered line)
/******************************************************************************
 * Project:  PROJ.4
 * Purpose:  Implementation of the aeqd (Azimuthal Equidistant) projection.
 * Author:   Gerald Evenden
 *
 ******************************************************************************
 * Copyright (c) 1995, Gerald Evenden
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#include "geodesic.h"
#include "proj.h"
#include "proj_internal.h"
#include <errno.h>
#include <math.h>

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

namespace { // anonymous namespace
struct pj_aeqd_data {
    double sinph0;
    double cosph0;
    double *en;
    double M1;
    double N1;
    double Mp;
    double He;
    double G;
    enum ::pj_aeqd_ns::Mode mode;
    struct geod_geodesic g;
};
} // anonymous namespace

PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam";

#define EPS10 1.e-10
#define TOL 1.e-14

static PJ *destructor(PJ *P, int errlev) { /* Destructor */
    if (nullptr == P)
        return nullptr;

    if (nullptr == P->opaque)
        return pj_default_destructor(P, errlev);

    free(static_cast<struct pj_aeqd_data *>(P->opaque)->en);
    return pj_default_destructor(P, errlev);
}

static PJ_XY e_guam_fwd(PJ_LP lp, PJ *P) { /* Guam elliptical */
    PJ_XY xy = {0.0, 0.0};
    struct pj_aeqd_data *Q = static_cast<struct pj_aeqd_data *>(P->opaque);
    double cosphi, sinphi, t;

    cosphi = cos(lp.phi);
    sinphi = sin(lp.phi);
    t = 1. / sqrt(1. - P->es * sinphi * sinphi);
    xy.x = lp.lam * cosphi * t;
    xy.y = pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->M1 +
           .5 * lp.lam * lp.lam * cosphi * sinphi * t;

    return xy;
}

static PJ_XY aeqd_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_aeqd_data *Q = static_cast<struct pj_aeqd_data *>(P->opaque);
    double coslam, cosphi, sinphi, rho;
    double azi1, azi2, s12;
    double lat1, lon1, lat2, lon2;

    coslam = cos(lp.lam);
    switch (Q->mode) {
    case pj_aeqd_ns::N_POLE:
        coslam = -coslam;
        PROJ_FALLTHROUGH;
    case pj_aeqd_ns::S_POLE:
        cosphi = cos(lp.phi);
        sinphi = sin(lp.phi);
        rho = fabs(Q->Mp - pj_mlfn(lp.phi, sinphi, cosphi, Q->en));
        xy.x = rho * sin(lp.lam);
        xy.y = rho * coslam;
        break;
    case pj_aeqd_ns::EQUIT:
    case pj_aeqd_ns::OBLIQ:
        if (fabs(lp.lam) < EPS10 && fabs(lp.phi - P->phi0) < EPS10) {
            xy.x = xy.y = 0.;
            break;
        }

        lat1 = P->phi0 / DEG_TO_RAD;
        lon1 = 0;
        lat2 = lp.phi / DEG_TO_RAD;
        lon2 = lp.lam / DEG_TO_RAD;

        geod_inverse(&Q->g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
        azi1 *= DEG_TO_RAD;
        xy.x = s12 * sin(azi1);
        xy.y = s12 * cos(azi1);
        break;
    }
    return xy;
}

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

    if (Q->mode == pj_aeqd_ns::EQUIT) {
        const double cosphi = cos(lp.phi);
        const double sinphi = sin(lp.phi);
        const double coslam = cos(lp.lam);
        const double sinlam = sin(lp.lam);

        xy.y = cosphi * coslam;

        if (fabs(fabs(xy.y) - 1.) < TOL) {
            if (xy.y < 0.) {
                proj_errno_set(
                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                return xy;
            } else
                return aeqd_e_forward(lp, P);
        } else {
            xy.y = acos(xy.y);
            xy.y /= sin(xy.y);
            xy.x = xy.y * cosphi * sinlam;
            xy.y *= sinphi;
        }
    } else if (Q->mode == pj_aeqd_ns::OBLIQ) {
        const double cosphi = cos(lp.phi);
        const double sinphi = sin(lp.phi);
        const double coslam = cos(lp.lam);
        const double sinlam = sin(lp.lam);
        const double cosphi_x_coslam = cosphi * coslam;

        xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi_x_coslam;

        if (fabs(fabs(xy.y) - 1.) < TOL) {
            if (xy.y < 0.) {
                proj_errno_set(
                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                return xy;
            } else
                return aeqd_e_forward(lp, P);
        } else {
            xy.y = acos(xy.y);
            xy.y /= sin(xy.y);
            xy.x = xy.y * cosphi * sinlam;
            xy.y *= Q->cosph0 * sinphi - Q->sinph0 * cosphi_x_coslam;
        }
    } else {
        double coslam = cos(lp.lam);
        double sinlam = sin(lp.lam);
        if (Q->mode == pj_aeqd_ns::N_POLE) {
            lp.phi = -lp.phi;
            coslam = -coslam;
        }
        if (fabs(lp.phi - M_HALFPI) < EPS10) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
        xy.y = (M_HALFPI + lp.phi);
        xy.x = xy.y * sinlam;
        xy.y *= coslam;
    }
    return xy;
}

static PJ_LP e_guam_inv(PJ_XY xy, PJ *P) { /* Guam elliptical */
    PJ_LP lp = {0.0, 0.0};
    struct pj_aeqd_data *Q = static_cast<struct pj_aeqd_data *>(P->opaque);
    double x2, t = 0.0;
    int i;

    x2 = 0.5 * xy.x * xy.x;
    lp.phi = P->phi0;
    for (i = 0; i < 3; ++i) {
        t = P->e * sin(lp.phi);
        t = sqrt(1. - t * t);
        lp.phi = pj_inv_mlfn(Q->M1 + xy.y - x2 * tan(lp.phi) * t, Q->en);
    }
    lp.lam = xy.x * t / cos(lp.phi);
    return lp;
}

static PJ_LP aeqd_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_aeqd_data *Q = static_cast<struct pj_aeqd_data *>(P->opaque);
    double azi1, azi2, s12, lat1, lon1, lat2, lon2;

    if ((s12 = hypot(xy.x, xy.y)) < EPS10) {
        lp.phi = P->phi0;
        lp.lam = 0.;
        return (lp);
    }
    if (Q->mode == pj_aeqd_ns::OBLIQ || Q->mode == pj_aeqd_ns::EQUIT) {
        lat1 = P->phi0 / DEG_TO_RAD;
        lon1 = 0;
        azi1 = atan2(xy.x, xy.y) / DEG_TO_RAD; // Clockwise from north
        geod_direct(&Q->g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2);
        lp.phi = lat2 * DEG_TO_RAD;
        lp.lam = lon2 * DEG_TO_RAD;
    } else { /* Polar */
        lp.phi = pj_inv_mlfn(
            Q->mode == pj_aeqd_ns::N_POLE ? Q->Mp - s12 : Q->Mp + s12, Q->en);
        lp.lam = atan2(xy.x, Q->mode == pj_aeqd_ns::N_POLE ? -xy.y : xy.y);
    }
    return lp;
}

static PJ_LP aeqd_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_aeqd_data *Q = static_cast<struct pj_aeqd_data *>(P->opaque);
    double cosc, c_rh, sinc;

    c_rh = hypot(xy.x, xy.y);
    if (c_rh > M_PI) {
        if (c_rh - EPS10 > M_PI) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return lp;
        }
        c_rh = M_PI;
    } else if (c_rh < EPS10) {
        lp.phi = P->phi0;
        lp.lam = 0.;
        return (lp);
    }
    if (Q->mode == pj_aeqd_ns::OBLIQ || Q->mode == pj_aeqd_ns::EQUIT) {
        sinc = sin(c_rh);
        cosc = cos(c_rh);
        if (Q->mode == pj_aeqd_ns::EQUIT) {
            lp.phi = aasin(P->ctx, xy.y * sinc / c_rh);
            xy.x *= sinc;
            xy.y = cosc * c_rh;
        } else {
            lp.phi = aasin(P->ctx,
                           cosc * Q->sinph0 + xy.y * sinc * Q->cosph0 / c_rh);
            xy.y = (cosc - Q->sinph0 * sin(lp.phi)) * c_rh;
            xy.x *= sinc * Q->cosph0;
        }
        lp.lam = xy.y == 0. ? 0. : atan2(xy.x, xy.y);
    } else if (Q->mode == pj_aeqd_ns::N_POLE) {
        lp.phi = M_HALFPI - c_rh;
        lp.lam = atan2(xy.x, -xy.y);
    } else {
        lp.phi = c_rh - M_HALFPI;
        lp.lam = atan2(xy.x, xy.y);
    }
    return lp;
}

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

    geod_init(&Q->g, 1, P->f);

    if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) {
        Q->mode = P->phi0 < 0. ? pj_aeqd_ns::S_POLE : pj_aeqd_ns::N_POLE;
        Q->sinph0 = P->phi0 < 0. ? -1. : 1.;
        Q->cosph0 = 0.;
    } else if (fabs(P->phi0) < EPS10) {
        Q->mode = pj_aeqd_ns::EQUIT;
        Q->sinph0 = 0.;
        Q->cosph0 = 1.;
    } else {
        Q->mode = pj_aeqd_ns::OBLIQ;
        Q->sinph0 = sin(P->phi0);
        Q->cosph0 = cos(P->phi0);
    }
    if (P->es == 0.0) {
        P->inv = aeqd_s_inverse;
        P->fwd = aeqd_s_forward;
    } else {
        if (!(Q->en = pj_enfn(P->n)))
            return pj_default_destructor(P, 0);
        if (pj_param(P->ctx, P->params, "bguam").i) {
            Q->M1 = pj_mlfn(P->phi0, Q->sinph0, Q->cosph0, Q->en);
            P->inv = e_guam_inv;
            P->fwd = e_guam_fwd;
        } else {
            switch (Q->mode) {
            case pj_aeqd_ns::N_POLE:
                Q->Mp = pj_mlfn(M_HALFPI, 1., 0., Q->en);
                break;
            case pj_aeqd_ns::S_POLE:
                Q->Mp = pj_mlfn(-M_HALFPI, -1., 0., Q->en);
                break;
            case pj_aeqd_ns::EQUIT:
            case pj_aeqd_ns::OBLIQ:
                Q->N1 = 1. / sqrt(1. - P->es * Q->sinph0 * Q->sinph0);
                Q->He = P->e / sqrt(P->one_es);
                Q->G = Q->sinph0 * Q->He;
                Q->He *= Q->cosph0;
                break;
            }
            P->inv = aeqd_e_inverse;
            P->fwd = aeqd_e_forward;
        }
    }

    return P;
}

#undef EPS10
#undef TOL

Coverage Report

Created: 2025-07-12 06:32

Line	Count	Source (jump to first uncovered line)
1		/******************************************************************************
2		* Project: PROJ.4
3		* Purpose: Implementation of the aeqd (Azimuthal Equidistant) projection.
4		* Author: Gerald Evenden
5		*
6		******************************************************************************
7		* Copyright (c) 1995, Gerald Evenden
8		*
9		* Permission is hereby granted, free of charge, to any person obtaining a
10		* copy of this software and associated documentation files (the "Software"),
11		* to deal in the Software without restriction, including without limitation
12		* the rights to use, copy, modify, merge, publish, distribute, sublicense,
13		* and/or sell copies of the Software, and to permit persons to whom the
14		* Software is furnished to do so, subject to the following conditions:
15		*
16		* The above copyright notice and this permission notice shall be included
17		* in all copies or substantial portions of the Software.
18		*
19		* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20		* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21		* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22		* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23		* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24		* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25		* DEALINGS IN THE SOFTWARE.
26		*****************************************************************************/
27
28		#include "geodesic.h"
29		#include "proj.h"
30		#include "proj_internal.h"
31		#include <errno.h>
32		#include <math.h>
33
34		namespace pj_aeqd_ns {
35		enum Mode { N_POLE = 0, S_POLE = 1, EQUIT = 2, OBLIQ = 3 };
36		}
37
38		namespace { // anonymous namespace
39		struct pj_aeqd_data {
40		double sinph0;
41		double cosph0;
42		double *en;
43		double M1;
44		double N1;
45		double Mp;
46		double He;
47		double G;
48		enum ::pj_aeqd_ns::Mode mode;
49		struct geod_geodesic g;
50		};
51		} // anonymous namespace
52
53		PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam";
54
55	10.8k	#define EPS10 1.e-10
56	0	#define TOL 1.e-14
57
58	310	static PJ destructor(PJ P, int errlev) { /* Destructor */
59	310	if (nullptr == P)
60	0	return nullptr;
61
62	310	if (nullptr == P->opaque)
63	0	return pj_default_destructor(P, errlev);
64
65	310	free(static_cast<struct pj_aeqd_data *>(P->opaque)->en);
66	310	return pj_default_destructor(P, errlev);
67	310	}
68
69	5.20k	static PJ_XY e_guam_fwd(PJ_LP lp, PJ P) { / Guam elliptical */
70	5.20k	PJ_XY xy = {0.0, 0.0};
71	5.20k	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(P->opaque);
72	5.20k	double cosphi, sinphi, t;
73
74	5.20k	cosphi = cos(lp.phi);
75	5.20k	sinphi = sin(lp.phi);
76	5.20k	t = 1. / sqrt(1. - P->es * sinphi * sinphi);
77	5.20k	xy.x = lp.lam * cosphi * t;
78	5.20k	xy.y = pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->M1 +
79	5.20k	.5 * lp.lam * lp.lam * cosphi * sinphi * t;
80
81	5.20k	return xy;
82	5.20k	}
83
84	5.12k	static PJ_XY aeqd_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
85	5.12k	PJ_XY xy = {0.0, 0.0};
86	5.12k	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(P->opaque);
87	5.12k	double coslam, cosphi, sinphi, rho;
88	5.12k	double azi1, azi2, s12;
89	5.12k	double lat1, lon1, lat2, lon2;
90
91	5.12k	coslam = cos(lp.lam);
92	5.12k	switch (Q->mode) {
93	0	case pj_aeqd_ns::N_POLE:
94	0	coslam = -coslam;
95	0	PROJ_FALLTHROUGH;
96	0	case pj_aeqd_ns::S_POLE:
97	0	cosphi = cos(lp.phi);
98	0	sinphi = sin(lp.phi);
99	0	rho = fabs(Q->Mp - pj_mlfn(lp.phi, sinphi, cosphi, Q->en));
100	0	xy.x = rho * sin(lp.lam);
101	0	xy.y = rho * coslam;
102	0	break;
103	5.12k	case pj_aeqd_ns::EQUIT:
104	5.12k	case pj_aeqd_ns::OBLIQ:
105	5.12k	if (fabs(lp.lam) < EPS10 && fabs(lp.phi - P->phi0) < EPS10) {
106	0	xy.x = xy.y = 0.;
107	0	break;
108	0	}
109
110	5.12k	lat1 = P->phi0 / DEG_TO_RAD;
111	5.12k	lon1 = 0;
112	5.12k	lat2 = lp.phi / DEG_TO_RAD;
113	5.12k	lon2 = lp.lam / DEG_TO_RAD;
114
115	5.12k	geod_inverse(&Q->g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
116	5.12k	azi1 *= DEG_TO_RAD;
117	5.12k	xy.x = s12 * sin(azi1);
118	5.12k	xy.y = s12 * cos(azi1);
119	5.12k	break;
120	5.12k	}
121	5.12k	return xy;
122	5.12k	}
123
124	0	static PJ_XY aeqd_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
125	0	PJ_XY xy = {0.0, 0.0};
126	0	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(P->opaque);
127
128	0	if (Q->mode == pj_aeqd_ns::EQUIT) {
129	0	const double cosphi = cos(lp.phi);
130	0	const double sinphi = sin(lp.phi);
131	0	const double coslam = cos(lp.lam);
132	0	const double sinlam = sin(lp.lam);
133
134	0	xy.y = cosphi * coslam;
135
136	0	if (fabs(fabs(xy.y) - 1.) < TOL) {
137	0	if (xy.y < 0.) {
138	0	proj_errno_set(
139	0	P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
140	0	return xy;
141	0	} else
142	0	return aeqd_e_forward(lp, P);
143	0	} else {
144	0	xy.y = acos(xy.y);
145	0	xy.y /= sin(xy.y);
146	0	xy.x = xy.y * cosphi * sinlam;
147	0	xy.y *= sinphi;
148	0	}
149	0	} else if (Q->mode == pj_aeqd_ns::OBLIQ) {
150	0	const double cosphi = cos(lp.phi);
151	0	const double sinphi = sin(lp.phi);
152	0	const double coslam = cos(lp.lam);
153	0	const double sinlam = sin(lp.lam);
154	0	const double cosphi_x_coslam = cosphi * coslam;
155
156	0	xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi_x_coslam;
157
158	0	if (fabs(fabs(xy.y) - 1.) < TOL) {
159	0	if (xy.y < 0.) {
160	0	proj_errno_set(
161	0	P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
162	0	return xy;
163	0	} else
164	0	return aeqd_e_forward(lp, P);
165	0	} else {
166	0	xy.y = acos(xy.y);
167	0	xy.y /= sin(xy.y);
168	0	xy.x = xy.y * cosphi * sinlam;
169	0	xy.y = Q->cosph0 sinphi - Q->sinph0 * cosphi_x_coslam;
170	0	}
171	0	} else {
172	0	double coslam = cos(lp.lam);
173	0	double sinlam = sin(lp.lam);
174	0	if (Q->mode == pj_aeqd_ns::N_POLE) {
175	0	lp.phi = -lp.phi;
176	0	coslam = -coslam;
177	0	}
178	0	if (fabs(lp.phi - M_HALFPI) < EPS10) {
179	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
180	0	return xy;
181	0	}
182	0	xy.y = (M_HALFPI + lp.phi);
183	0	xy.x = xy.y * sinlam;
184	0	xy.y *= coslam;
185	0	}
186	0	return xy;
187	0	}
188
189	0	static PJ_LP e_guam_inv(PJ_XY xy, PJ P) { / Guam elliptical */
190	0	PJ_LP lp = {0.0, 0.0};
191	0	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(P->opaque);
192	0	double x2, t = 0.0;
193	0	int i;
194
195	0	x2 = 0.5 * xy.x * xy.x;
196	0	lp.phi = P->phi0;
197	0	for (i = 0; i < 3; ++i) {
198	0	t = P->e * sin(lp.phi);
199	0	t = sqrt(1. - t * t);
200	0	lp.phi = pj_inv_mlfn(Q->M1 + xy.y - x2 * tan(lp.phi) * t, Q->en);
201	0	}
202	0	lp.lam = xy.x * t / cos(lp.phi);
203	0	return lp;
204	0	}
205
206	0	static PJ_LP aeqd_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
207	0	PJ_LP lp = {0.0, 0.0};
208	0	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(P->opaque);
209	0	double azi1, azi2, s12, lat1, lon1, lat2, lon2;
210
211	0	if ((s12 = hypot(xy.x, xy.y)) < EPS10) {
212	0	lp.phi = P->phi0;
213	0	lp.lam = 0.;
214	0	return (lp);
215	0	}
216	0	if (Q->mode == pj_aeqd_ns::OBLIQ \|\| Q->mode == pj_aeqd_ns::EQUIT) {
217	0	lat1 = P->phi0 / DEG_TO_RAD;
218	0	lon1 = 0;
219	0	azi1 = atan2(xy.x, xy.y) / DEG_TO_RAD; // Clockwise from north
220	0	geod_direct(&Q->g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2);
221	0	lp.phi = lat2 * DEG_TO_RAD;
222	0	lp.lam = lon2 * DEG_TO_RAD;
223	0	} else { /* Polar */
224	0	lp.phi = pj_inv_mlfn(
225	0	Q->mode == pj_aeqd_ns::N_POLE ? Q->Mp - s12 : Q->Mp + s12, Q->en);
226	0	lp.lam = atan2(xy.x, Q->mode == pj_aeqd_ns::N_POLE ? -xy.y : xy.y);
227	0	}
228	0	return lp;
229	0	}
230
231	0	static PJ_LP aeqd_s_inverse(PJ_XY xy, PJ P) { / Spheroidal, inverse */
232	0	PJ_LP lp = {0.0, 0.0};
233	0	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(P->opaque);
234	0	double cosc, c_rh, sinc;
235
236	0	c_rh = hypot(xy.x, xy.y);
237	0	if (c_rh > M_PI) {
238	0	if (c_rh - EPS10 > M_PI) {
239	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
240	0	return lp;
241	0	}
242	0	c_rh = M_PI;
243	0	} else if (c_rh < EPS10) {
244	0	lp.phi = P->phi0;
245	0	lp.lam = 0.;
246	0	return (lp);
247	0	}
248	0	if (Q->mode == pj_aeqd_ns::OBLIQ \|\| Q->mode == pj_aeqd_ns::EQUIT) {
249	0	sinc = sin(c_rh);
250	0	cosc = cos(c_rh);
251	0	if (Q->mode == pj_aeqd_ns::EQUIT) {
252	0	lp.phi = aasin(P->ctx, xy.y * sinc / c_rh);
253	0	xy.x *= sinc;
254	0	xy.y = cosc * c_rh;
255	0	} else {
256	0	lp.phi = aasin(P->ctx,
257	0	cosc * Q->sinph0 + xy.y * sinc * Q->cosph0 / c_rh);
258	0	xy.y = (cosc - Q->sinph0 * sin(lp.phi)) * c_rh;
259	0	xy.x = sinc Q->cosph0;
260	0	}
261	0	lp.lam = xy.y == 0. ? 0. : atan2(xy.x, xy.y);
262	0	} else if (Q->mode == pj_aeqd_ns::N_POLE) {
263	0	lp.phi = M_HALFPI - c_rh;
264	0	lp.lam = atan2(xy.x, -xy.y);
265	0	} else {
266	0	lp.phi = c_rh - M_HALFPI;
267	0	lp.lam = atan2(xy.x, xy.y);
268	0	}
269	0	return lp;
270	0	}
271
272	310	PJ *PJ_PROJECTION(aeqd) {
273	310	struct pj_aeqd_data Q = static_cast<struct pj_aeqd_data >(
274	310	calloc(1, sizeof(struct pj_aeqd_data)));
275	310	if (nullptr == Q)
276	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
277	310	P->opaque = Q;
278	310	P->destructor = destructor;
279
280	310	geod_init(&Q->g, 1, P->f);
281
282	310	if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) {
283	13	Q->mode = P->phi0 < 0. ? pj_aeqd_ns::S_POLE : pj_aeqd_ns::N_POLE;
284	13	Q->sinph0 = P->phi0 < 0. ? -1. : 1.;
285	13	Q->cosph0 = 0.;
286	297	} else if (fabs(P->phi0) < EPS10) {
287	250	Q->mode = pj_aeqd_ns::EQUIT;
288	250	Q->sinph0 = 0.;
289	250	Q->cosph0 = 1.;
290	250	} else {
291	47	Q->mode = pj_aeqd_ns::OBLIQ;
292	47	Q->sinph0 = sin(P->phi0);
293	47	Q->cosph0 = cos(P->phi0);
294	47	}
295	310	if (P->es == 0.0) {
296	7	P->inv = aeqd_s_inverse;
297	7	P->fwd = aeqd_s_forward;
298	303	} else {
299	303	if (!(Q->en = pj_enfn(P->n)))
300	0	return pj_default_destructor(P, 0);
301	303	if (pj_param(P->ctx, P->params, "bguam").i) {
302	56	Q->M1 = pj_mlfn(P->phi0, Q->sinph0, Q->cosph0, Q->en);
303	56	P->inv = e_guam_inv;
304	56	P->fwd = e_guam_fwd;
305	247	} else {
306	247	switch (Q->mode) {
307	2	case pj_aeqd_ns::N_POLE:
308	2	Q->Mp = pj_mlfn(M_HALFPI, 1., 0., Q->en);
309	2	break;
310	11	case pj_aeqd_ns::S_POLE:
311	11	Q->Mp = pj_mlfn(-M_HALFPI, -1., 0., Q->en);
312	11	break;
313	194	case pj_aeqd_ns::EQUIT:
314	234	case pj_aeqd_ns::OBLIQ:
315	234	Q->N1 = 1. / sqrt(1. - P->es * Q->sinph0 * Q->sinph0);
316	234	Q->He = P->e / sqrt(P->one_es);
317	234	Q->G = Q->sinph0 * Q->He;
318	234	Q->He *= Q->cosph0;
319	234	break;
320	247	}
321	247	P->inv = aeqd_e_inverse;
322	247	P->fwd = aeqd_e_forward;
323	247	}
324	303	}
325
326	310	return P;
327	310	}
328
329		#undef EPS10
330		#undef TOL