/src/PROJ/src/projections/imw_p.cpp

Source


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

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

PROJ_HEAD(imw_p, "International Map of the World Polyconic")
"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";

#define TOL 1e-10
#define EPS 1e-10

namespace { // anonymous namespace
enum Mode {
    NONE_IS_ZERO = 0,  /* phi_1 and phi_2 != 0 */
    PHI_1_IS_ZERO = 1, /* phi_1 = 0 */
    PHI_2_IS_ZERO = -1 /* phi_2 = 0 */
};
} // anonymous namespace

namespace { // anonymous namespace
struct pj_imw_p_data {
    double P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
    double phi_1, phi_2, lam_1;
    double *en;
    enum Mode mode;
};
} // anonymous namespace

static int phi12(PJ *P, double *del, double *sig) {
    struct pj_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(P->opaque);
    int err = 0;

    if (!pj_param(P->ctx, P->params, "tlat_1").i) {
        proj_log_error(P, _("Missing parameter: lat_1 should be specified"));
        err = PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    } else if (!pj_param(P->ctx, P->params, "tlat_2").i) {
        proj_log_error(P, _("Missing parameter: lat_2 should be specified"));
        err = PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    } else {
        Q->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
        Q->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
        *del = 0.5 * (Q->phi_2 - Q->phi_1);
        *sig = 0.5 * (Q->phi_2 + Q->phi_1);
        err = (fabs(*del) < EPS || fabs(*sig) < EPS)
                  ? PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE
                  : 0;
        if (err) {
            proj_log_error(
                P, _("Illegal value for lat_1 and lat_2: |lat_1 - lat_2| "
                     "and |lat_1 + lat_2| should be > 0"));
        }
    }
    return err;
}

static PJ_XY loc_for(PJ_LP lp, PJ *P, double *yc) {
    struct pj_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(P->opaque);
    PJ_XY xy;

    if (lp.phi == 0.0) {
        xy.x = lp.lam;
        xy.y = 0.;
    } else {
        double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;

        sp = sin(lp.phi);
        m = pj_mlfn(lp.phi, sp, cos(lp.phi), Q->en);
        xa = Q->Pp + Q->Qp * m;
        ya = Q->P + Q->Q * m;
        R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
        C = sqrt(R * R - xa * xa);
        if (lp.phi < 0.)
            C = -C;
        C += ya - R;
        if (Q->mode == PHI_2_IS_ZERO) {
            xb = lp.lam;
            yb = Q->C2;
        } else {
            t = lp.lam * Q->sphi_2;
            xb = Q->R_2 * sin(t);
            yb = Q->C2 + Q->R_2 * (1. - cos(t));
        }
        if (Q->mode == PHI_1_IS_ZERO) {
            xc = lp.lam;
            *yc = 0.;
        } else {
            t = lp.lam * Q->sphi_1;
            xc = Q->R_1 * sin(t);
            *yc = Q->R_1 * (1. - cos(t));
        }
        D = (xb - xc) / (yb - *yc);
        B = xc + D * (C + R - *yc);
        xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
        if (lp.phi > 0)
            xy.x = -xy.x;
        xy.x = (B + xy.x) / (1. + D * D);
        xy.y = sqrt(R * R - xy.x * xy.x);
        if (lp.phi > 0)
            xy.y = -xy.y;
        xy.y += C + R;
    }
    return xy;
}

static PJ_XY imw_p_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    double yc;
    PJ_XY xy = loc_for(lp, P, &yc);
    return (xy);
}

static PJ_LP imw_p_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(P->opaque);
    PJ_XY t;
    double yc = 0.0;
    int i = 0;
    const int N_MAX_ITER = 1000; /* Arbitrarily chosen number... */

    lp.phi = Q->phi_2;
    lp.lam = xy.x / cos(lp.phi);
    do {
        t = loc_for(lp, P, &yc);
        const double denom = t.y - yc;
        if (denom != 0 || fabs(t.y - xy.y) > TOL) {
            if (denom == 0) {
                proj_errno_set(
                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                return proj_coord_error().lp;
            }
            lp.phi = ((lp.phi - Q->phi_1) * (xy.y - yc) / denom) + Q->phi_1;
        }
        if (t.x != 0 && fabs(t.x - xy.x) > TOL)
            lp.lam = lp.lam * xy.x / t.x;
        i++;
    } while (i < N_MAX_ITER &&
             (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL));

    if (i == N_MAX_ITER) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return proj_coord_error().lp;
    }

    return lp;
}

static void xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
    double F;

    *sp = sin(phi);
    *R = 1. / (tan(phi) * sqrt(1. - P->es * *sp * *sp));
    F = static_cast<struct pj_imw_p_data *>(P->opaque)->lam_1 * *sp;
    *y = *R * (1 - cos(F));
    *x = *R * sin(F);
}

static PJ *pj_imw_p_destructor(PJ *P, int errlev) {
    if (nullptr == P)
        return nullptr;

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

    if (static_cast<struct pj_imw_p_data *>(P->opaque)->en)
        free(static_cast<struct pj_imw_p_data *>(P->opaque)->en);

    return pj_default_destructor(P, errlev);
}

PJ *PJ_PROJECTION(imw_p) {
    double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
    int err;
    struct pj_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(
        calloc(1, sizeof(struct pj_imw_p_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    if (!(Q->en = pj_enfn(P->n)))
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    if ((err = phi12(P, &del, &sig)) != 0) {
        return pj_imw_p_destructor(P, err);
    }
    if (Q->phi_2 < Q->phi_1) { /* make sure P->phi_1 most southerly */
        del = Q->phi_1;
        Q->phi_1 = Q->phi_2;
        Q->phi_2 = del;
    }
    if (pj_param(P->ctx, P->params, "tlon_1").i)
        Q->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
    else { /* use predefined based upon latitude */
        sig = fabs(sig * RAD_TO_DEG);
        if (sig <= 60)
            sig = 2.;
        else if (sig <= 76)
            sig = 4.;
        else
            sig = 8.;
        Q->lam_1 = sig * DEG_TO_RAD;
    }
    Q->mode = NONE_IS_ZERO;
    if (Q->phi_1 != 0.0)
        xy(P, Q->phi_1, &x1, &y1, &Q->sphi_1, &Q->R_1);
    else {
        Q->mode = PHI_1_IS_ZERO;
        y1 = 0.;
        x1 = Q->lam_1;
    }
    if (Q->phi_2 != 0.0)
        xy(P, Q->phi_2, &x2, &T2, &Q->sphi_2, &Q->R_2);
    else {
        Q->mode = PHI_2_IS_ZERO;
        T2 = 0.;
        x2 = Q->lam_1;
    }
    m1 = pj_mlfn(Q->phi_1, Q->sphi_1, cos(Q->phi_1), Q->en);
    m2 = pj_mlfn(Q->phi_2, Q->sphi_2, cos(Q->phi_2), Q->en);
    t = m2 - m1;
    s = x2 - x1;
    y2 = sqrt(t * t - s * s) + y1;
    Q->C2 = y2 - T2;
    t = 1. / t;
    Q->P = (m2 * y1 - m1 * y2) * t;
    Q->Q = (y2 - y1) * t;
    Q->Pp = (m2 * x1 - m1 * x2) * t;
    Q->Qp = (x2 - x1) * t;

    P->fwd = imw_p_e_forward;
    P->inv = imw_p_e_inverse;
    P->destructor = pj_imw_p_destructor;

    return P;
}

#undef TOL
#undef EPS

Coverage Report

Created: 2026-06-09 07:13

Line	Count	Source
1
2
3		#include <errno.h>
4		#include <math.h>
5
6		#include "proj.h"
7		#include "proj_internal.h"
8
9		PROJ_HEAD(imw_p, "International Map of the World Polyconic")
10		"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
11
12	0	#define TOL 1e-10
13	599	#define EPS 1e-10
14
15		namespace { // anonymous namespace
16		enum Mode {
17		NONE_IS_ZERO = 0, /* phi_1 and phi_2 != 0 */
18		PHI_1_IS_ZERO = 1, /* phi_1 = 0 */
19		PHI_2_IS_ZERO = -1 /* phi_2 = 0 */
20		};
21		} // anonymous namespace
22
23		namespace { // anonymous namespace
24		struct pj_imw_p_data {
25		double P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
26		double phi_1, phi_2, lam_1;
27		double *en;
28		enum Mode mode;
29		};
30		} // anonymous namespace
31
32	216	static int phi12(PJ P, double del, double *sig) {
33	216	struct pj_imw_p_data Q = static_cast<struct pj_imw_p_data >(P->opaque);
34	216	int err = 0;
35
36	216	if (!pj_param(P->ctx, P->params, "tlat_1").i) {
37	1	proj_log_error(P, _("Missing parameter: lat_1 should be specified"));
38	1	err = PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
39	215	} else if (!pj_param(P->ctx, P->params, "tlat_2").i) {
40	3	proj_log_error(P, _("Missing parameter: lat_2 should be specified"));
41	3	err = PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
42	212	} else {
43	212	Q->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
44	212	Q->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
45	212	del = 0.5 (Q->phi_2 - Q->phi_1);
46	212	sig = 0.5 (Q->phi_2 + Q->phi_1);
47	212	err = (fabs(del) < EPS \|\| fabs(sig) < EPS)
48	212	? PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE
49	212	: 0;
50	212	if (err) {
51	40	proj_log_error(
52	40	P, _("Illegal value for lat_1 and lat_2: \|lat_1 - lat_2\| "
53	40	"and \|lat_1 + lat_2\| should be > 0"));
54	40	}
55	212	}
56	216	return err;
57	216	}
58
59	0	static PJ_XY loc_for(PJ_LP lp, PJ P, double yc) {
60	0	struct pj_imw_p_data Q = static_cast<struct pj_imw_p_data >(P->opaque);
61	0	PJ_XY xy;
62
63	0	if (lp.phi == 0.0) {
64	0	xy.x = lp.lam;
65	0	xy.y = 0.;
66	0	} else {
67	0	double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
68
69	0	sp = sin(lp.phi);
70	0	m = pj_mlfn(lp.phi, sp, cos(lp.phi), Q->en);
71	0	xa = Q->Pp + Q->Qp * m;
72	0	ya = Q->P + Q->Q * m;
73	0	R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
74	0	C = sqrt(R * R - xa * xa);
75	0	if (lp.phi < 0.)
76	0	C = -C;
77	0	C += ya - R;
78	0	if (Q->mode == PHI_2_IS_ZERO) {
79	0	xb = lp.lam;
80	0	yb = Q->C2;
81	0	} else {
82	0	t = lp.lam * Q->sphi_2;
83	0	xb = Q->R_2 * sin(t);
84	0	yb = Q->C2 + Q->R_2 * (1. - cos(t));
85	0	}
86	0	if (Q->mode == PHI_1_IS_ZERO) {
87	0	xc = lp.lam;
88	0	*yc = 0.;
89	0	} else {
90	0	t = lp.lam * Q->sphi_1;
91	0	xc = Q->R_1 * sin(t);
92	0	yc = Q->R_1 (1. - cos(t));
93	0	}
94	0	D = (xb - xc) / (yb - *yc);
95	0	B = xc + D * (C + R - *yc);
96	0	xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
97	0	if (lp.phi > 0)
98	0	xy.x = -xy.x;
99	0	xy.x = (B + xy.x) / (1. + D * D);
100	0	xy.y = sqrt(R * R - xy.x * xy.x);
101	0	if (lp.phi > 0)
102	0	xy.y = -xy.y;
103	0	xy.y += C + R;
104	0	}
105	0	return xy;
106	0	}
107
108	0	static PJ_XY imw_p_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
109	0	double yc;
110	0	PJ_XY xy = loc_for(lp, P, &yc);
111	0	return (xy);
112	0	}
113
114	0	static PJ_LP imw_p_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
115	0	PJ_LP lp = {0.0, 0.0};
116	0	struct pj_imw_p_data Q = static_cast<struct pj_imw_p_data >(P->opaque);
117	0	PJ_XY t;
118	0	double yc = 0.0;
119	0	int i = 0;
120	0	const int N_MAX_ITER = 1000; /* Arbitrarily chosen number... */
121
122	0	lp.phi = Q->phi_2;
123	0	lp.lam = xy.x / cos(lp.phi);
124	0	do {
125	0	t = loc_for(lp, P, &yc);
126	0	const double denom = t.y - yc;
127	0	if (denom != 0 \|\| fabs(t.y - xy.y) > TOL) {
128	0	if (denom == 0) {
129	0	proj_errno_set(
130	0	P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
131	0	return proj_coord_error().lp;
132	0	}
133	0	lp.phi = ((lp.phi - Q->phi_1) * (xy.y - yc) / denom) + Q->phi_1;
134	0	}
135	0	if (t.x != 0 && fabs(t.x - xy.x) > TOL)
136	0	lp.lam = lp.lam * xy.x / t.x;
137	0	i++;
138	0	} while (i < N_MAX_ITER &&
139	0	(fabs(t.x - xy.x) > TOL \|\| fabs(t.y - xy.y) > TOL));
140
141	0	if (i == N_MAX_ITER) {
142	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
143	0	return proj_coord_error().lp;
144	0	}
145
146	0	return lp;
147	0	}
148
149	172	static void xy(PJ P, double phi, double x, double y, double sp, double *R) {
150	172	double F;
151
152	172	*sp = sin(phi);
153	172	R = 1. / (tan(phi) sqrt(1. - P->es * sp *sp));
154	172	F = static_cast<struct pj_imw_p_data >(P->opaque)->lam_1 *sp;
155	172	y = R * (1 - cos(F));
156	172	x = R * sin(F);
157	172	}
158
159	216	static PJ pj_imw_p_destructor(PJ P, int errlev) {
160	216	if (nullptr == P)
161	0	return nullptr;
162
163	216	if (nullptr == P->opaque)
164	0	return pj_default_destructor(P, errlev);
165
166	216	if (static_cast<struct pj_imw_p_data *>(P->opaque)->en)
167	216	free(static_cast<struct pj_imw_p_data *>(P->opaque)->en);
168
169	216	return pj_default_destructor(P, errlev);
170	216	}
171
172	216	PJ *PJ_PROJECTION(imw_p) {
173	216	double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
174	216	int err;
175	216	struct pj_imw_p_data Q = static_cast<struct pj_imw_p_data >(
176	216	calloc(1, sizeof(struct pj_imw_p_data)));
177	216	if (nullptr == Q)
178	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
179	216	P->opaque = Q;
180
181	216	if (!(Q->en = pj_enfn(P->n)))
182	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
183	216	if ((err = phi12(P, &del, &sig)) != 0) {
184	44	return pj_imw_p_destructor(P, err);
185	44	}
186	172	if (Q->phi_2 < Q->phi_1) { /* make sure P->phi_1 most southerly */
187	52	del = Q->phi_1;
188	52	Q->phi_1 = Q->phi_2;
189	52	Q->phi_2 = del;
190	52	}
191	172	if (pj_param(P->ctx, P->params, "tlon_1").i)
192	3	Q->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
193	169	else { /* use predefined based upon latitude */
194	169	sig = fabs(sig * RAD_TO_DEG);
195	169	if (sig <= 60)
196	100	sig = 2.;
197	69	else if (sig <= 76)
198	43	sig = 4.;
199	26	else
200	26	sig = 8.;
201	169	Q->lam_1 = sig * DEG_TO_RAD;
202	169	}
203	172	Q->mode = NONE_IS_ZERO;
204	172	if (Q->phi_1 != 0.0)
205	120	xy(P, Q->phi_1, &x1, &y1, &Q->sphi_1, &Q->R_1);
206	52	else {
207	52	Q->mode = PHI_1_IS_ZERO;
208	52	y1 = 0.;
209	52	x1 = Q->lam_1;
210	52	}
211	172	if (Q->phi_2 != 0.0)
212	52	xy(P, Q->phi_2, &x2, &T2, &Q->sphi_2, &Q->R_2);
213	120	else {
214	120	Q->mode = PHI_2_IS_ZERO;
215	120	T2 = 0.;
216	120	x2 = Q->lam_1;
217	120	}
218	172	m1 = pj_mlfn(Q->phi_1, Q->sphi_1, cos(Q->phi_1), Q->en);
219	172	m2 = pj_mlfn(Q->phi_2, Q->sphi_2, cos(Q->phi_2), Q->en);
220	172	t = m2 - m1;
221	172	s = x2 - x1;
222	172	y2 = sqrt(t * t - s * s) + y1;
223	172	Q->C2 = y2 - T2;
224	172	t = 1. / t;
225	172	Q->P = (m2 * y1 - m1 * y2) * t;
226	172	Q->Q = (y2 - y1) * t;
227	172	Q->Pp = (m2 * x1 - m1 * x2) * t;
228	172	Q->Qp = (x2 - x1) * t;
229
230	172	P->fwd = imw_p_e_forward;
231	172	P->inv = imw_p_e_inverse;
232	172	P->destructor = pj_imw_p_destructor;
233
234	172	return P;
235	216	}
236
237		#undef TOL
238		#undef EPS