/src/PROJ/src/projections/mod_ster.cpp

Source (jump to first uncovered line)
/* based upon Snyder and Linck, USGS-NMD */

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

PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)";
PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)";
PROJ_HEAD(gs48, "Modified Stereographic of 48 U.S.") "\n\tAzi(mod)";
PROJ_HEAD(alsk, "Modified Stereographic of Alaska") "\n\tAzi(mod)";
PROJ_HEAD(gs50, "Modified Stereographic of 50 U.S.") "\n\tAzi(mod)";

#define EPSLN 1e-12

namespace { // anonymous namespace
struct pj_mod_ster_data {
    const COMPLEX *zcoeff;
    double cchio, schio;
    int n;
};
} // anonymous namespace

static PJ_XY mod_ster_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_mod_ster_data *Q =
        static_cast<struct pj_mod_ster_data *>(P->opaque);
    double sinlon, coslon, esphi, chi, schi, cchi, s;
    COMPLEX p;

    sinlon = sin(lp.lam);
    coslon = cos(lp.lam);
    esphi = P->e * sin(lp.phi);
    chi = 2. * atan(tan((M_HALFPI + lp.phi) * .5) *
                    pow((1. - esphi) / (1. + esphi), P->e * .5)) -
          M_HALFPI;
    schi = sin(chi);
    cchi = cos(chi);
    const double denom = 1. + Q->schio * schi + Q->cchio * cchi * coslon;
    if (denom == 0) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return xy;
    }
    s = 2. / denom;
    p.r = s * cchi * sinlon;
    p.i = s * (Q->cchio * schi - Q->schio * cchi * coslon);
    p = pj_zpoly1(p, Q->zcoeff, Q->n);
    xy.x = p.r;
    xy.y = p.i;

    return xy;
}

static PJ_LP mod_ster_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_mod_ster_data *Q =
        static_cast<struct pj_mod_ster_data *>(P->opaque);
    int nn;
    COMPLEX p, fxy, fpxy, dp;
    double den, rh = 0.0, z, sinz = 0.0, cosz = 0.0, chi, phi = 0.0, esphi;

    p.r = xy.x;
    p.i = xy.y;
    for (nn = 20; nn; --nn) {
        fxy = pj_zpolyd1(p, Q->zcoeff, Q->n, &fpxy);
        fxy.r -= xy.x;
        fxy.i -= xy.y;
        den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
        dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
        dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
        p.r += dp.r;
        p.i += dp.i;
        if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
            break;
    }
    if (nn) {
        rh = hypot(p.r, p.i);
        z = 2. * atan(.5 * rh);
        sinz = sin(z);
        cosz = cos(z);
        if (fabs(rh) <= EPSLN) {
            /* if we end up here input coordinates were (0,0).
             * pj_inv() adds P->lam0 to lp.lam, this way we are
             * sure to get the correct offset */
            lp.lam = 0.0;
            lp.phi = P->phi0;
            return lp;
        }
        chi = aasin(P->ctx, cosz * Q->schio + p.i * sinz * Q->cchio / rh);
        phi = chi;
        for (nn = 20; nn; --nn) {
            double dphi;
            esphi = P->e * sin(phi);
            dphi = 2. * atan(tan((M_HALFPI + chi) * .5) *
                             pow((1. + esphi) / (1. - esphi), P->e * .5)) -
                   M_HALFPI - phi;
            phi += dphi;
            if (fabs(dphi) <= EPSLN)
                break;
        }
    }
    if (nn) {
        lp.phi = phi;
        lp.lam =
            atan2(p.r * sinz, rh * Q->cchio * cosz - p.i * Q->schio * sinz);
    } else
        lp.lam = lp.phi = HUGE_VAL;
    return lp;
}

static PJ *mod_ster_setup(PJ *P) { /* general initialization */
    struct pj_mod_ster_data *Q =
        static_cast<struct pj_mod_ster_data *>(P->opaque);
    double esphi, chio;

    if (P->es != 0.0) {
        esphi = P->e * sin(P->phi0);
        chio = 2. * atan(tan((M_HALFPI + P->phi0) * .5) *
                         pow((1. - esphi) / (1. + esphi), P->e * .5)) -
               M_HALFPI;
    } else
        chio = P->phi0;
    Q->schio = sin(chio);
    Q->cchio = cos(chio);
    P->inv = mod_ster_e_inverse;
    P->fwd = mod_ster_e_forward;

    return P;
}

/* Miller Oblated Stereographic */
PJ *PJ_PROJECTION(mil_os) {
    static const COMPLEX AB[] = {{0.924500, 0.}, {0., 0.}, {0.019430, 0.}};

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

    Q->n = 2;
    P->lam0 = DEG_TO_RAD * 20.;
    P->phi0 = DEG_TO_RAD * 18.;
    Q->zcoeff = AB;
    P->es = 0.;

    return mod_ster_setup(P);
}

/* Lee Oblated Stereographic */
PJ *PJ_PROJECTION(lee_os) {
    static const COMPLEX AB[] = {
        {0.721316, 0.}, {0., 0.}, {-0.0088162, -0.00617325}};

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

    Q->n = 2;
    P->lam0 = DEG_TO_RAD * -165.;
    P->phi0 = DEG_TO_RAD * -10.;
    Q->zcoeff = AB;
    P->es = 0.;

    return mod_ster_setup(P);
}

PJ *PJ_PROJECTION(gs48) {
    static const COMPLEX /* 48 United States */
        AB[] = {
            {0.98879, 0.}, {0., 0.}, {-0.050909, 0.}, {0., 0.}, {0.075528, 0.}};

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

    Q->n = 4;
    P->lam0 = DEG_TO_RAD * -96.;
    P->phi0 = DEG_TO_RAD * 39.;
    Q->zcoeff = AB;
    P->es = 0.;
    P->a = 6370997.;

    return mod_ster_setup(P);
}

PJ *PJ_PROJECTION(alsk) {
    static const COMPLEX ABe[] = {
        /* Alaska ellipsoid */
        {.9945303, 0.},         {.0052083, -.0027404}, {.0072721, .0048181},
        {-.0151089, -.1932526}, {.0642675, -.1381226}, {.3582802, -.2884586},
    };

    static const COMPLEX ABs[] = {/* Alaska sphere */
                                  {.9972523, 0.},        {.0052513, -.0041175},
                                  {.0074606, .0048125},  {-.0153783, -.1968253},
                                  {.0636871, -.1408027}, {.3660976, -.2937382}};

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

    Q->n = 5;
    P->lam0 = DEG_TO_RAD * -152.;
    P->phi0 = DEG_TO_RAD * 64.;
    if (P->es != 0.0) { /* fixed ellipsoid/sphere */
        Q->zcoeff = ABe;
        P->a = 6378206.4;
        P->e = sqrt(P->es = 0.00676866);
    } else {
        Q->zcoeff = ABs;
        P->a = 6370997.;
    }

    return mod_ster_setup(P);
}

PJ *PJ_PROJECTION(gs50) {
    static const COMPLEX ABe[] = {
        /* GS50 ellipsoid */
        {.9827497, 0.},         {.0210669, .0053804},  {-.1031415, -.0571664},
        {-.0323337, -.0322847}, {.0502303, .1211983},  {.0251805, .0895678},
        {-.0012315, -.1416121}, {.0072202, -.1317091}, {-.0194029, .0759677},
        {-.0210072, .0834037}};

    static const COMPLEX ABs[] = {
        /* GS50 sphere */
        {.9842990, 0.},         {.0211642, .0037608},  {-.1036018, -.0575102},
        {-.0329095, -.0320119}, {.0499471, .1223335},  {.0260460, .0899805},
        {.0007388, -.1435792},  {.0075848, -.1334108}, {-.0216473, .0776645},
        {-.0225161, .0853673}};

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

    Q->n = 9;
    P->lam0 = DEG_TO_RAD * -120.;
    P->phi0 = DEG_TO_RAD * 45.;
    if (P->es != 0.0) { /* fixed ellipsoid/sphere */
        Q->zcoeff = ABe;
        P->a = 6378206.4;
        P->e = sqrt(P->es = 0.00676866);
    } else {
        Q->zcoeff = ABs;
        P->a = 6370997.;
    }

    return mod_ster_setup(P);
}

#undef EPSLN

Coverage Report

Created: 2025-07-11 06:33

Line	Count	Source (jump to first uncovered line)
1		/* based upon Snyder and Linck, USGS-NMD */
2
3		#include "proj.h"
4		#include "proj_internal.h"
5		#include <errno.h>
6		#include <math.h>
7
8		PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)";
9		PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)";
10		PROJ_HEAD(gs48, "Modified Stereographic of 48 U.S.") "\n\tAzi(mod)";
11		PROJ_HEAD(alsk, "Modified Stereographic of Alaska") "\n\tAzi(mod)";
12		PROJ_HEAD(gs50, "Modified Stereographic of 50 U.S.") "\n\tAzi(mod)";
13
14	0	#define EPSLN 1e-12
15
16		namespace { // anonymous namespace
17		struct pj_mod_ster_data {
18		const COMPLEX *zcoeff;
19		double cchio, schio;
20		int n;
21		};
22		} // anonymous namespace
23
24	3.27k	static PJ_XY mod_ster_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
25	3.27k	PJ_XY xy = {0.0, 0.0};
26	3.27k	struct pj_mod_ster_data *Q =
27	3.27k	static_cast<struct pj_mod_ster_data *>(P->opaque);
28	3.27k	double sinlon, coslon, esphi, chi, schi, cchi, s;
29	3.27k	COMPLEX p;
30
31	3.27k	sinlon = sin(lp.lam);
32	3.27k	coslon = cos(lp.lam);
33	3.27k	esphi = P->e * sin(lp.phi);
34	3.27k	chi = 2. * atan(tan((M_HALFPI + lp.phi) * .5) *
35	3.27k	pow((1. - esphi) / (1. + esphi), P->e * .5)) -
36	3.27k	M_HALFPI;
37	3.27k	schi = sin(chi);
38	3.27k	cchi = cos(chi);
39	3.27k	const double denom = 1. + Q->schio * schi + Q->cchio * cchi * coslon;
40	3.27k	if (denom == 0) {
41	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
42	0	return xy;
43	0	}
44	3.27k	s = 2. / denom;
45	3.27k	p.r = s * cchi * sinlon;
46	3.27k	p.i = s * (Q->cchio * schi - Q->schio * cchi * coslon);
47	3.27k	p = pj_zpoly1(p, Q->zcoeff, Q->n);
48	3.27k	xy.x = p.r;
49	3.27k	xy.y = p.i;
50
51	3.27k	return xy;
52	3.27k	}
53
54	0	static PJ_LP mod_ster_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
55	0	PJ_LP lp = {0.0, 0.0};
56	0	struct pj_mod_ster_data *Q =
57	0	static_cast<struct pj_mod_ster_data *>(P->opaque);
58	0	int nn;
59	0	COMPLEX p, fxy, fpxy, dp;
60	0	double den, rh = 0.0, z, sinz = 0.0, cosz = 0.0, chi, phi = 0.0, esphi;
61
62	0	p.r = xy.x;
63	0	p.i = xy.y;
64	0	for (nn = 20; nn; --nn) {
65	0	fxy = pj_zpolyd1(p, Q->zcoeff, Q->n, &fpxy);
66	0	fxy.r -= xy.x;
67	0	fxy.i -= xy.y;
68	0	den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
69	0	dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
70	0	dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
71	0	p.r += dp.r;
72	0	p.i += dp.i;
73	0	if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
74	0	break;
75	0	}
76	0	if (nn) {
77	0	rh = hypot(p.r, p.i);
78	0	z = 2. * atan(.5 * rh);
79	0	sinz = sin(z);
80	0	cosz = cos(z);
81	0	if (fabs(rh) <= EPSLN) {
82		/* if we end up here input coordinates were (0,0).
83		* pj_inv() adds P->lam0 to lp.lam, this way we are
84		* sure to get the correct offset */
85	0	lp.lam = 0.0;
86	0	lp.phi = P->phi0;
87	0	return lp;
88	0	}
89	0	chi = aasin(P->ctx, cosz * Q->schio + p.i * sinz * Q->cchio / rh);
90	0	phi = chi;
91	0	for (nn = 20; nn; --nn) {
92	0	double dphi;
93	0	esphi = P->e * sin(phi);
94	0	dphi = 2. * atan(tan((M_HALFPI + chi) * .5) *
95	0	pow((1. + esphi) / (1. - esphi), P->e * .5)) -
96	0	M_HALFPI - phi;
97	0	phi += dphi;
98	0	if (fabs(dphi) <= EPSLN)
99	0	break;
100	0	}
101	0	}
102	0	if (nn) {
103	0	lp.phi = phi;
104	0	lp.lam =
105	0	atan2(p.r * sinz, rh * Q->cchio * cosz - p.i * Q->schio * sinz);
106	0	} else
107	0	lp.lam = lp.phi = HUGE_VAL;
108	0	return lp;
109	0	}
110
111	195	static PJ mod_ster_setup(PJ P) { /* general initialization */
112	195	struct pj_mod_ster_data *Q =
113	195	static_cast<struct pj_mod_ster_data *>(P->opaque);
114	195	double esphi, chio;
115
116	195	if (P->es != 0.0) {
117	138	esphi = P->e * sin(P->phi0);
118	138	chio = 2. * atan(tan((M_HALFPI + P->phi0) * .5) *
119	138	pow((1. - esphi) / (1. + esphi), P->e * .5)) -
120	138	M_HALFPI;
121	138	} else
122	57	chio = P->phi0;
123	195	Q->schio = sin(chio);
124	195	Q->cchio = cos(chio);
125	195	P->inv = mod_ster_e_inverse;
126	195	P->fwd = mod_ster_e_forward;
127
128	195	return P;
129	195	}
130
131		/* Miller Oblated Stereographic */
132	2	PJ *PJ_PROJECTION(mil_os) {
133	2	static const COMPLEX AB[] = {{0.924500, 0.}, {0., 0.}, {0.019430, 0.}};
134
135	2	struct pj_mod_ster_data Q = static_cast<struct pj_mod_ster_data >(
136	2	calloc(1, sizeof(struct pj_mod_ster_data)));
137	2	if (nullptr == Q)
138	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
139	2	P->opaque = Q;
140
141	2	Q->n = 2;
142	2	P->lam0 = DEG_TO_RAD * 20.;
143	2	P->phi0 = DEG_TO_RAD * 18.;
144	2	Q->zcoeff = AB;
145	2	P->es = 0.;
146
147	2	return mod_ster_setup(P);
148	2	}
149
150		/* Lee Oblated Stereographic */
151	19	PJ *PJ_PROJECTION(lee_os) {
152	19	static const COMPLEX AB[] = {
153	19	{0.721316, 0.}, {0., 0.}, {-0.0088162, -0.00617325}};
154
155	19	struct pj_mod_ster_data Q = static_cast<struct pj_mod_ster_data >(
156	19	calloc(1, sizeof(struct pj_mod_ster_data)));
157	19	if (nullptr == Q)
158	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
159	19	P->opaque = Q;
160
161	19	Q->n = 2;
162	19	P->lam0 = DEG_TO_RAD * -165.;
163	19	P->phi0 = DEG_TO_RAD * -10.;
164	19	Q->zcoeff = AB;
165	19	P->es = 0.;
166
167	19	return mod_ster_setup(P);
168	19	}
169
170	20	PJ *PJ_PROJECTION(gs48) {
171	20	static const COMPLEX /* 48 United States */
172	20	AB[] = {
173	20	{0.98879, 0.}, {0., 0.}, {-0.050909, 0.}, {0., 0.}, {0.075528, 0.}};
174
175	20	struct pj_mod_ster_data Q = static_cast<struct pj_mod_ster_data >(
176	20	calloc(1, sizeof(struct pj_mod_ster_data)));
177	20	if (nullptr == Q)
178	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
179	20	P->opaque = Q;
180
181	20	Q->n = 4;
182	20	P->lam0 = DEG_TO_RAD * -96.;
183	20	P->phi0 = DEG_TO_RAD * 39.;
184	20	Q->zcoeff = AB;
185	20	P->es = 0.;
186	20	P->a = 6370997.;
187
188	20	return mod_ster_setup(P);
189	20	}
190
191	136	PJ *PJ_PROJECTION(alsk) {
192	136	static const COMPLEX ABe[] = {
193		/* Alaska ellipsoid */
194	136	{.9945303, 0.}, {.0052083, -.0027404}, {.0072721, .0048181},
195	136	{-.0151089, -.1932526}, {.0642675, -.1381226}, {.3582802, -.2884586},
196	136	};
197
198	136	static const COMPLEX ABs[] = {/* Alaska sphere */
199	136	{.9972523, 0.}, {.0052513, -.0041175},
200	136	{.0074606, .0048125}, {-.0153783, -.1968253},
201	136	{.0636871, -.1408027}, {.3660976, -.2937382}};
202
203	136	struct pj_mod_ster_data Q = static_cast<struct pj_mod_ster_data >(
204	136	calloc(1, sizeof(struct pj_mod_ster_data)));
205	136	if (nullptr == Q)
206	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
207	136	P->opaque = Q;
208
209	136	Q->n = 5;
210	136	P->lam0 = DEG_TO_RAD * -152.;
211	136	P->phi0 = DEG_TO_RAD * 64.;
212	136	if (P->es != 0.0) { /* fixed ellipsoid/sphere */
213	123	Q->zcoeff = ABe;
214	123	P->a = 6378206.4;
215	123	P->e = sqrt(P->es = 0.00676866);
216	123	} else {
217	13	Q->zcoeff = ABs;
218	13	P->a = 6370997.;
219	13	}
220
221	136	return mod_ster_setup(P);
222	136	}
223
224	18	PJ *PJ_PROJECTION(gs50) {
225	18	static const COMPLEX ABe[] = {
226		/* GS50 ellipsoid */
227	18	{.9827497, 0.}, {.0210669, .0053804}, {-.1031415, -.0571664},
228	18	{-.0323337, -.0322847}, {.0502303, .1211983}, {.0251805, .0895678},
229	18	{-.0012315, -.1416121}, {.0072202, -.1317091}, {-.0194029, .0759677},
230	18	{-.0210072, .0834037}};
231
232	18	static const COMPLEX ABs[] = {
233		/* GS50 sphere */
234	18	{.9842990, 0.}, {.0211642, .0037608}, {-.1036018, -.0575102},
235	18	{-.0329095, -.0320119}, {.0499471, .1223335}, {.0260460, .0899805},
236	18	{.0007388, -.1435792}, {.0075848, -.1334108}, {-.0216473, .0776645},
237	18	{-.0225161, .0853673}};
238
239	18	struct pj_mod_ster_data Q = static_cast<struct pj_mod_ster_data >(
240	18	calloc(1, sizeof(struct pj_mod_ster_data)));
241	18	if (nullptr == Q)
242	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
243	18	P->opaque = Q;
244
245	18	Q->n = 9;
246	18	P->lam0 = DEG_TO_RAD * -120.;
247	18	P->phi0 = DEG_TO_RAD * 45.;
248	18	if (P->es != 0.0) { /* fixed ellipsoid/sphere */
249	15	Q->zcoeff = ABe;
250	15	P->a = 6378206.4;
251	15	P->e = sqrt(P->es = 0.00676866);
252	15	} else {
253	3	Q->zcoeff = ABs;
254	3	P->a = 6370997.;
255	3	}
256
257	18	return mod_ster_setup(P);
258	18	}
259
260		#undef EPSLN