/src/PROJ/src/projections/eqdc.cpp

Source (jump to first uncovered line)


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

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

namespace { // anonymous namespace
struct pj_eqdc_data {
    double phi1;
    double phi2;
    double n;
    double rho;
    double rho0;
    double c;
    double *en;
    int ellips;
};
} // anonymous namespace

PROJ_HEAD(eqdc, "Equidistant Conic")
"\n\tConic, Sph&Ell\n\tlat_1= lat_2=";
#define EPS10 1.e-10

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

    Q->rho =
        Q->c -
        (Q->ellips ? pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), Q->en) : lp.phi);
    const double lam_mul_n = lp.lam * Q->n;
    xy.x = Q->rho * sin(lam_mul_n);
    xy.y = Q->rho0 - Q->rho * cos(lam_mul_n);

    return xy;
}

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

    if ((Q->rho = hypot(xy.x, xy.y = Q->rho0 - xy.y)) != 0.0) {
        if (Q->n < 0.) {
            Q->rho = -Q->rho;
            xy.x = -xy.x;
            xy.y = -xy.y;
        }
        lp.phi = Q->c - Q->rho;
        if (Q->ellips)
            lp.phi = pj_inv_mlfn(lp.phi, Q->en);
        lp.lam = atan2(xy.x, xy.y) / Q->n;
    } else {
        lp.lam = 0.;
        lp.phi = Q->n > 0. ? M_HALFPI : -M_HALFPI;
    }
    return lp;
}

static PJ *pj_eqdc_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_eqdc_data *>(P->opaque)->en);
    return pj_default_destructor(P, errlev);
}

PJ *PJ_PROJECTION(eqdc) {
    double cosphi, sinphi;
    int secant;

    struct pj_eqdc_data *Q = static_cast<struct pj_eqdc_data *>(
        calloc(1, sizeof(struct pj_eqdc_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_eqdc_destructor;

    Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
    Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;

    if (fabs(Q->phi1) > M_HALFPI) {
        proj_log_error(P,
                       _("Invalid value for lat_1: |lat_1| should be <= 90°"));
        return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    if (fabs(Q->phi2) > M_HALFPI) {
        proj_log_error(P,
                       _("Invalid value for lat_2: |lat_2| should be <= 90°"));
        return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    if (fabs(Q->phi1 + Q->phi2) < EPS10) {
        proj_log_error(P, _("Invalid value for lat_1 and lat_2: |lat_1 + "
                            "lat_2| should be > 0"));
        return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    if (!(Q->en = pj_enfn(P->n)))
        return pj_eqdc_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);

    sinphi = sin(Q->phi1);
    Q->n = sinphi;
    cosphi = cos(Q->phi1);
    secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
    Q->ellips = (P->es > 0.);
    if (Q->ellips) {
        double ml1, m1;

        m1 = pj_msfn(sinphi, cosphi, P->es);
        ml1 = pj_mlfn(Q->phi1, sinphi, cosphi, Q->en);
        if (secant) { /* secant cone */
            sinphi = sin(Q->phi2);
            cosphi = cos(Q->phi2);
            const double ml2 = pj_mlfn(Q->phi2, sinphi, cosphi, Q->en);
            if (ml1 == ml2) {
                proj_log_error(P, _("Eccentricity too close to 1"));
                return pj_eqdc_destructor(
                    P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
            }
            Q->n = (m1 - pj_msfn(sinphi, cosphi, P->es)) / (ml2 - ml1);
            if (Q->n == 0) {
                // Not quite, but es is very close to 1...
                proj_log_error(P, _("Invalid value for eccentricity"));
                return pj_eqdc_destructor(
                    P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
            }
        }
        Q->c = ml1 + m1 / Q->n;
        Q->rho0 = Q->c - pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
    } else {
        if (secant)
            Q->n = (cosphi - cos(Q->phi2)) / (Q->phi2 - Q->phi1);
        if (Q->n == 0) {
            proj_log_error(P, _("Invalid value for lat_1 and lat_2: lat_1 + "
                                "lat_2 should be > 0"));
            return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        Q->c = Q->phi1 + cos(Q->phi1) / Q->n;
        Q->rho0 = Q->c - P->phi0;
    }

    P->inv = eqdc_e_inverse;
    P->fwd = eqdc_e_forward;

    return P;
}

#undef EPS10

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		#include <math.h>
9
10		namespace { // anonymous namespace
11		struct pj_eqdc_data {
12		double phi1;
13		double phi2;
14		double n;
15		double rho;
16		double rho0;
17		double c;
18		double *en;
19		int ellips;
20		};
21		} // anonymous namespace
22
23		PROJ_HEAD(eqdc, "Equidistant Conic")
24		"\n\tConic, Sph&Ell\n\tlat_1= lat_2=";
25	890	#define EPS10 1.e-10
26
27	0	static PJ_XY eqdc_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
28	0	PJ_XY xy = {0.0, 0.0};
29	0	struct pj_eqdc_data Q = static_cast<struct pj_eqdc_data >(P->opaque);
30
31	0	Q->rho =
32	0	Q->c -
33	0	(Q->ellips ? pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), Q->en) : lp.phi);
34	0	const double lam_mul_n = lp.lam * Q->n;
35	0	xy.x = Q->rho * sin(lam_mul_n);
36	0	xy.y = Q->rho0 - Q->rho * cos(lam_mul_n);
37
38	0	return xy;
39	0	}
40
41	0	static PJ_LP eqdc_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
42	0	PJ_LP lp = {0.0, 0.0};
43	0	struct pj_eqdc_data Q = static_cast<struct pj_eqdc_data >(P->opaque);
44
45	0	if ((Q->rho = hypot(xy.x, xy.y = Q->rho0 - xy.y)) != 0.0) {
46	0	if (Q->n < 0.) {
47	0	Q->rho = -Q->rho;
48	0	xy.x = -xy.x;
49	0	xy.y = -xy.y;
50	0	}
51	0	lp.phi = Q->c - Q->rho;
52	0	if (Q->ellips)
53	0	lp.phi = pj_inv_mlfn(lp.phi, Q->en);
54	0	lp.lam = atan2(xy.x, xy.y) / Q->n;
55	0	} else {
56	0	lp.lam = 0.;
57	0	lp.phi = Q->n > 0. ? M_HALFPI : -M_HALFPI;
58	0	}
59	0	return lp;
60	0	}
61
62	456	static PJ pj_eqdc_destructor(PJ P, int errlev) { /* Destructor */
63	456	if (nullptr == P)
64	0	return nullptr;
65
66	456	if (nullptr == P->opaque)
67	0	return pj_default_destructor(P, errlev);
68
69	456	free(static_cast<struct pj_eqdc_data *>(P->opaque)->en);
70	456	return pj_default_destructor(P, errlev);
71	456	}
72
73	456	PJ *PJ_PROJECTION(eqdc) {
74	456	double cosphi, sinphi;
75	456	int secant;
76
77	456	struct pj_eqdc_data Q = static_cast<struct pj_eqdc_data >(
78	456	calloc(1, sizeof(struct pj_eqdc_data)));
79	456	if (nullptr == Q)
80	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
81	456	P->opaque = Q;
82	456	P->destructor = pj_eqdc_destructor;
83
84	456	Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
85	456	Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
86
87	456	if (fabs(Q->phi1) > M_HALFPI) {
88	3	proj_log_error(P,
89	3	_("Invalid value for lat_1: \|lat_1\| should be <= 90°"));
90	3	return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
91	3	}
92
93	453	if (fabs(Q->phi2) > M_HALFPI) {
94	3	proj_log_error(P,
95	3	_("Invalid value for lat_2: \|lat_2\| should be <= 90°"));
96	3	return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
97	3	}
98	450	if (fabs(Q->phi1 + Q->phi2) < EPS10) {
99	10	proj_log_error(P, _("Invalid value for lat_1 and lat_2: \|lat_1 + "
100	10	"lat_2\| should be > 0"));
101	10	return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
102	10	}
103
104	440	if (!(Q->en = pj_enfn(P->n)))
105	0	return pj_eqdc_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
106
107	440	sinphi = sin(Q->phi1);
108	440	Q->n = sinphi;
109	440	cosphi = cos(Q->phi1);
110	440	secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
111	440	Q->ellips = (P->es > 0.);
112	440	if (Q->ellips) {
113	215	double ml1, m1;
114
115	215	m1 = pj_msfn(sinphi, cosphi, P->es);
116	215	ml1 = pj_mlfn(Q->phi1, sinphi, cosphi, Q->en);
117	215	if (secant) { /* secant cone */
118	139	sinphi = sin(Q->phi2);
119	139	cosphi = cos(Q->phi2);
120	139	const double ml2 = pj_mlfn(Q->phi2, sinphi, cosphi, Q->en);
121	139	if (ml1 == ml2) {
122	0	proj_log_error(P, _("Eccentricity too close to 1"));
123	0	return pj_eqdc_destructor(
124	0	P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
125	0	}
126	139	Q->n = (m1 - pj_msfn(sinphi, cosphi, P->es)) / (ml2 - ml1);
127	139	if (Q->n == 0) {
128		// Not quite, but es is very close to 1...
129	3	proj_log_error(P, _("Invalid value for eccentricity"));
130	3	return pj_eqdc_destructor(
131	3	P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
132	3	}
133	139	}
134	212	Q->c = ml1 + m1 / Q->n;
135	212	Q->rho0 = Q->c - pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
136	225	} else {
137	225	if (secant)
138	149	Q->n = (cosphi - cos(Q->phi2)) / (Q->phi2 - Q->phi1);
139	225	if (Q->n == 0) {
140	0	proj_log_error(P, _("Invalid value for lat_1 and lat_2: lat_1 + "
141	0	"lat_2 should be > 0"));
142	0	return pj_eqdc_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
143	0	}
144	225	Q->c = Q->phi1 + cos(Q->phi1) / Q->n;
145	225	Q->rho0 = Q->c - P->phi0;
146	225	}
147
148	437	P->inv = eqdc_e_inverse;
149	437	P->fwd = eqdc_e_forward;
150
151	437	return P;
152	440	}
153
154		#undef EPS10

Coverage Report

Created: 2024-02-25 06:14