/src/PROJ/src/projections/stere.cpp

Source (jump to first uncovered line)

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

PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Ell\n\tsouth";

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

namespace { // anonymous namespace
struct pj_stere {
    double phits;
    double sinX1;
    double cosX1;
    double akm1;
    enum Mode mode;
};
} // anonymous namespace

#define sinph0 static_cast<struct pj_stere *>(P->opaque)->sinX1
#define cosph0 static_cast<struct pj_stere *>(P->opaque)->cosX1
#define EPS10 1.e-10
#define TOL 1.e-8
#define NITER 8
#define CONV 1.e-10

static double ssfn_(double phit, double sinphi, double eccen) {
    sinphi *= eccen;
    return (tan(.5 * (M_HALFPI + phit)) *
            pow((1. - sinphi) / (1. + sinphi), .5 * eccen));
}

static PJ_XY stere_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_stere *Q = static_cast<struct pj_stere *>(P->opaque);
    double coslam, sinlam, sinX = 0.0, cosX = 0.0, A = 0.0, sinphi;

    coslam = cos(lp.lam);
    sinlam = sin(lp.lam);
    sinphi = sin(lp.phi);
    if (Q->mode == OBLIQ || Q->mode == EQUIT) {
        const double X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - M_HALFPI;
        sinX = sin(X);
        cosX = cos(X);
    }

    switch (Q->mode) {
    case OBLIQ: {
        const double denom =
            Q->cosX1 * (1. + Q->sinX1 * sinX + Q->cosX1 * cosX * coslam);
        if (denom == 0) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return proj_coord_error().xy;
        }
        A = Q->akm1 / denom;
        xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam);
        xy.x = A * cosX;
        break;
    }

    case EQUIT:
        /* avoid zero division */
        if (1. + cosX * coslam == 0.0) {
            xy.y = HUGE_VAL;
        } else {
            A = Q->akm1 / (1. + cosX * coslam);
            xy.y = A * sinX;
        }
        xy.x = A * cosX;
        break;

    case S_POLE:
        lp.phi = -lp.phi;
        coslam = -coslam;
        sinphi = -sinphi;
        PROJ_FALLTHROUGH;
    case N_POLE:
        if (fabs(lp.phi - M_HALFPI) < 1e-15)
            xy.x = 0;
        else
            xy.x = Q->akm1 * pj_tsfn(lp.phi, sinphi, P->e);
        xy.y = -xy.x * coslam;
        break;
    }

    xy.x = xy.x * sinlam;
    return xy;
}

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

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

    switch (Q->mode) {
    case EQUIT:
        xy.y = 1. + cosphi * coslam;
        goto oblcon;
    case OBLIQ:
        xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
    oblcon:
        if (xy.y <= EPS10) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
        xy.y = Q->akm1 / xy.y;
        xy.x = xy.y * cosphi * sinlam;
        xy.y *= (Q->mode == EQUIT) ? sinphi
                                   : cosph0 * sinphi - sinph0 * cosphi * coslam;
        break;
    case N_POLE:
        coslam = -coslam;
        lp.phi = -lp.phi;
        PROJ_FALLTHROUGH;
    case S_POLE:
        if (fabs(lp.phi - M_HALFPI) < TOL) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
        xy.y = Q->akm1 * tan(M_FORTPI + .5 * lp.phi);
        xy.x = sinlam * xy.y;
        xy.y *= coslam;
        break;
    }
    return xy;
}

static PJ_LP stere_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_stere *Q = static_cast<struct pj_stere *>(P->opaque);
    double cosphi, sinphi, tp = 0.0, phi_l = 0.0, rho, halfe = 0.0,
                           halfpi = 0.0;

    rho = hypot(xy.x, xy.y);

    switch (Q->mode) {
    case OBLIQ:
    case EQUIT:
        tp = 2. * atan2(rho * Q->cosX1, Q->akm1);
        cosphi = cos(tp);
        sinphi = sin(tp);
        if (rho == 0.0)
            phi_l = asin(cosphi * Q->sinX1);
        else
            phi_l = asin(cosphi * Q->sinX1 + (xy.y * sinphi * Q->cosX1 / rho));

        tp = tan(.5 * (M_HALFPI + phi_l));
        xy.x *= sinphi;
        xy.y = rho * Q->cosX1 * cosphi - xy.y * Q->sinX1 * sinphi;
        halfpi = M_HALFPI;
        halfe = .5 * P->e;
        break;
    case N_POLE:
        xy.y = -xy.y;
        PROJ_FALLTHROUGH;
    case S_POLE:
        tp = -rho / Q->akm1;
        phi_l = M_HALFPI - 2. * atan(tp);
        halfpi = -M_HALFPI;
        halfe = -.5 * P->e;
        break;
    }

    for (int i = NITER; i > 0; --i) {
        sinphi = P->e * sin(phi_l);
        lp.phi =
            2. * atan(tp * pow((1. + sinphi) / (1. - sinphi), halfe)) - halfpi;
        if (fabs(phi_l - lp.phi) < CONV) {
            if (Q->mode == S_POLE)
                lp.phi = -lp.phi;
            lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);
            return lp;
        }
        phi_l = lp.phi;
    }

    proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
    return lp;
}

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

    const double rh = hypot(xy.x, xy.y);
    c = 2. * atan(rh / Q->akm1);
    sinc = sin(c);
    cosc = cos(c);
    lp.lam = 0.;

    switch (Q->mode) {
    case EQUIT:
        if (fabs(rh) <= EPS10)
            lp.phi = 0.;
        else
            lp.phi = asin(xy.y * sinc / rh);
        if (cosc != 0. || xy.x != 0.)
            lp.lam = atan2(xy.x * sinc, cosc * rh);
        break;
    case OBLIQ:
        if (fabs(rh) <= EPS10)
            lp.phi = P->phi0;
        else
            lp.phi = asin(cosc * sinph0 + xy.y * sinc * cosph0 / rh);
        c = cosc - sinph0 * sin(lp.phi);
        if (c != 0. || xy.x != 0.)
            lp.lam = atan2(xy.x * sinc * cosph0, c * rh);
        break;
    case N_POLE:
        xy.y = -xy.y;
        PROJ_FALLTHROUGH;
    case S_POLE:
        if (fabs(rh) <= EPS10)
            lp.phi = P->phi0;
        else
            lp.phi = asin(Q->mode == S_POLE ? -cosc : cosc);
        lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);
        break;
    }
    return lp;
}

static PJ *stere_setup(PJ *P) { /* general initialization */
    double t;
    struct pj_stere *Q = static_cast<struct pj_stere *>(P->opaque);

    if (fabs((t = fabs(P->phi0)) - M_HALFPI) < EPS10)
        Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
    else
        Q->mode = t > EPS10 ? OBLIQ : EQUIT;
    Q->phits = fabs(Q->phits);

    if (P->es != 0.0) {
        double X;

        switch (Q->mode) {
        case N_POLE:
        case S_POLE:
            if (fabs(Q->phits - M_HALFPI) < EPS10)
                Q->akm1 =
                    2. * P->k0 /
                    sqrt(pow(1 + P->e, 1 + P->e) * pow(1 - P->e, 1 - P->e));
            else {
                t = sin(Q->phits);
                Q->akm1 = cos(Q->phits) / pj_tsfn(Q->phits, t, P->e);
                t *= P->e;
                Q->akm1 /= sqrt(1. - t * t);
            }
            break;
        case EQUIT:
        case OBLIQ:
            t = sin(P->phi0);
            X = 2. * atan(ssfn_(P->phi0, t, P->e)) - M_HALFPI;
            t *= P->e;
            Q->akm1 = 2. * P->k0 * cos(P->phi0) / sqrt(1. - t * t);
            Q->sinX1 = sin(X);
            Q->cosX1 = cos(X);
            break;
        }
        P->inv = stere_e_inverse;
        P->fwd = stere_e_forward;
    } else {
        switch (Q->mode) {
        case OBLIQ:
            sinph0 = sin(P->phi0);
            cosph0 = cos(P->phi0);
            PROJ_FALLTHROUGH;
        case EQUIT:
            Q->akm1 = 2. * P->k0;
            break;
        case S_POLE:
        case N_POLE:
            Q->akm1 = fabs(Q->phits - M_HALFPI) >= EPS10
                          ? cos(Q->phits) / tan(M_FORTPI - .5 * Q->phits)
                          : 2. * P->k0;
            break;
        }

        P->inv = stere_s_inverse;
        P->fwd = stere_s_forward;
    }
    return P;
}

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

    Q->phits = pj_param(P->ctx, P->params, "tlat_ts").i
                   ? pj_param(P->ctx, P->params, "rlat_ts").f
                   : M_HALFPI;

    return stere_setup(P);
}

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

    /* International Ellipsoid */
    P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? -M_HALFPI : M_HALFPI;
    if (P->es == 0.0) {
        proj_log_error(
            P,
            _("Invalid value for es: only ellipsoidal formulation supported"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    P->k0 = .994;
    P->x0 = 2000000.;
    P->y0 = 2000000.;
    Q->phits = M_HALFPI;
    P->lam0 = 0.;

    return stere_setup(P);
}

#undef sinph0
#undef cosph0
#undef EPS10
#undef TOL
#undef NITER
#undef CONV

Coverage Report

Created: 2025-07-23 06:58

Line	Count	Source (jump to first uncovered line)
1
2		#include "proj.h"
3		#include "proj_internal.h"
4		#include <errno.h>
5		#include <math.h>
6
7		PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
8		PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Ell\n\tsouth";
9
10		namespace { // anonymous namespace
11		enum Mode { S_POLE = 0, N_POLE = 1, OBLIQ = 2, EQUIT = 3 };
12		} // anonymous namespace
13
14		namespace { // anonymous namespace
15		struct pj_stere {
16		double phits;
17		double sinX1;
18		double cosX1;
19		double akm1;
20		enum Mode mode;
21		};
22		} // anonymous namespace
23
24	4	#define sinph0 static_cast<struct pj_stere *>(P->opaque)->sinX1
25	4	#define cosph0 static_cast<struct pj_stere *>(P->opaque)->cosX1
26	2.24k	#define EPS10 1.e-10
27	0	#define TOL 1.e-8
28	0	#define NITER 8
29	0	#define CONV 1.e-10
30
31	17.1k	static double ssfn_(double phit, double sinphi, double eccen) {
32	17.1k	sinphi *= eccen;
33	17.1k	return (tan(.5 * (M_HALFPI + phit)) *
34	17.1k	pow((1. - sinphi) / (1. + sinphi), .5 * eccen));
35	17.1k	}
36
37	29.7k	static PJ_XY stere_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
38	29.7k	PJ_XY xy = {0.0, 0.0};
39	29.7k	struct pj_stere Q = static_cast<struct pj_stere >(P->opaque);
40	29.7k	double coslam, sinlam, sinX = 0.0, cosX = 0.0, A = 0.0, sinphi;
41
42	29.7k	coslam = cos(lp.lam);
43	29.7k	sinlam = sin(lp.lam);
44	29.7k	sinphi = sin(lp.phi);
45	29.7k	if (Q->mode == OBLIQ \|\| Q->mode == EQUIT) {
46	16.5k	const double X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - M_HALFPI;
47	16.5k	sinX = sin(X);
48	16.5k	cosX = cos(X);
49	16.5k	}
50
51	29.7k	switch (Q->mode) {
52	7.56k	case OBLIQ: {
53	7.56k	const double denom =
54	7.56k	Q->cosX1 * (1. + Q->sinX1 * sinX + Q->cosX1 * cosX * coslam);
55	7.56k	if (denom == 0) {
56	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
57	0	return proj_coord_error().xy;
58	0	}
59	7.56k	A = Q->akm1 / denom;
60	7.56k	xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam);
61	7.56k	xy.x = A * cosX;
62	7.56k	break;
63	7.56k	}
64
65	8.98k	case EQUIT:
66		/* avoid zero division */
67	8.98k	if (1. + cosX * coslam == 0.0) {
68	0	xy.y = HUGE_VAL;
69	8.98k	} else {
70	8.98k	A = Q->akm1 / (1. + cosX * coslam);
71	8.98k	xy.y = A * sinX;
72	8.98k	}
73	8.98k	xy.x = A * cosX;
74	8.98k	break;
75
76	2.60k	case S_POLE:
77	2.60k	lp.phi = -lp.phi;
78	2.60k	coslam = -coslam;
79	2.60k	sinphi = -sinphi;
80	2.60k	PROJ_FALLTHROUGH;
81	13.1k	case N_POLE:
82	13.1k	if (fabs(lp.phi - M_HALFPI) < 1e-15)
83	0	xy.x = 0;
84	13.1k	else
85	13.1k	xy.x = Q->akm1 * pj_tsfn(lp.phi, sinphi, P->e);
86	13.1k	xy.y = -xy.x * coslam;
87	13.1k	break;
88	29.7k	}
89
90	29.7k	xy.x = xy.x * sinlam;
91	29.7k	return xy;
92	29.7k	}
93
94	0	static PJ_XY stere_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
95	0	PJ_XY xy = {0.0, 0.0};
96	0	struct pj_stere Q = static_cast<struct pj_stere >(P->opaque);
97	0	double sinphi, cosphi, coslam, sinlam;
98
99	0	sinphi = sin(lp.phi);
100	0	cosphi = cos(lp.phi);
101	0	coslam = cos(lp.lam);
102	0	sinlam = sin(lp.lam);
103
104	0	switch (Q->mode) {
105	0	case EQUIT:
106	0	xy.y = 1. + cosphi * coslam;
107	0	goto oblcon;
108	0	case OBLIQ:
109	0	xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
110	0	oblcon:
111	0	if (xy.y <= EPS10) {
112	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
113	0	return xy;
114	0	}
115	0	xy.y = Q->akm1 / xy.y;
116	0	xy.x = xy.y * cosphi * sinlam;
117	0	xy.y *= (Q->mode == EQUIT) ? sinphi
118	0	: cosph0 * sinphi - sinph0 * cosphi * coslam;
119	0	break;
120	0	case N_POLE:
121	0	coslam = -coslam;
122	0	lp.phi = -lp.phi;
123	0	PROJ_FALLTHROUGH;
124	0	case S_POLE:
125	0	if (fabs(lp.phi - M_HALFPI) < TOL) {
126	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
127	0	return xy;
128	0	}
129	0	xy.y = Q->akm1 * tan(M_FORTPI + .5 * lp.phi);
130	0	xy.x = sinlam * xy.y;
131	0	xy.y *= coslam;
132	0	break;
133	0	}
134	0	return xy;
135	0	}
136
137	0	static PJ_LP stere_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
138	0	PJ_LP lp = {0.0, 0.0};
139	0	struct pj_stere Q = static_cast<struct pj_stere >(P->opaque);
140	0	double cosphi, sinphi, tp = 0.0, phi_l = 0.0, rho, halfe = 0.0,
141	0	halfpi = 0.0;
142
143	0	rho = hypot(xy.x, xy.y);
144
145	0	switch (Q->mode) {
146	0	case OBLIQ:
147	0	case EQUIT:
148	0	tp = 2. * atan2(rho * Q->cosX1, Q->akm1);
149	0	cosphi = cos(tp);
150	0	sinphi = sin(tp);
151	0	if (rho == 0.0)
152	0	phi_l = asin(cosphi * Q->sinX1);
153	0	else
154	0	phi_l = asin(cosphi * Q->sinX1 + (xy.y * sinphi * Q->cosX1 / rho));
155
156	0	tp = tan(.5 * (M_HALFPI + phi_l));
157	0	xy.x *= sinphi;
158	0	xy.y = rho * Q->cosX1 * cosphi - xy.y * Q->sinX1 * sinphi;
159	0	halfpi = M_HALFPI;
160	0	halfe = .5 * P->e;
161	0	break;
162	0	case N_POLE:
163	0	xy.y = -xy.y;
164	0	PROJ_FALLTHROUGH;
165	0	case S_POLE:
166	0	tp = -rho / Q->akm1;
167	0	phi_l = M_HALFPI - 2. * atan(tp);
168	0	halfpi = -M_HALFPI;
169	0	halfe = -.5 * P->e;
170	0	break;
171	0	}
172
173	0	for (int i = NITER; i > 0; --i) {
174	0	sinphi = P->e * sin(phi_l);
175	0	lp.phi =
176	0	2. * atan(tp * pow((1. + sinphi) / (1. - sinphi), halfe)) - halfpi;
177	0	if (fabs(phi_l - lp.phi) < CONV) {
178	0	if (Q->mode == S_POLE)
179	0	lp.phi = -lp.phi;
180	0	lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);
181	0	return lp;
182	0	}
183	0	phi_l = lp.phi;
184	0	}
185
186	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
187	0	return lp;
188	0	}
189
190	0	static PJ_LP stere_s_inverse(PJ_XY xy, PJ P) { / Spheroidal, inverse */
191	0	PJ_LP lp = {0.0, 0.0};
192	0	struct pj_stere Q = static_cast<struct pj_stere >(P->opaque);
193	0	double c, sinc, cosc;
194
195	0	const double rh = hypot(xy.x, xy.y);
196	0	c = 2. * atan(rh / Q->akm1);
197	0	sinc = sin(c);
198	0	cosc = cos(c);
199	0	lp.lam = 0.;
200
201	0	switch (Q->mode) {
202	0	case EQUIT:
203	0	if (fabs(rh) <= EPS10)
204	0	lp.phi = 0.;
205	0	else
206	0	lp.phi = asin(xy.y * sinc / rh);
207	0	if (cosc != 0. \|\| xy.x != 0.)
208	0	lp.lam = atan2(xy.x * sinc, cosc * rh);
209	0	break;
210	0	case OBLIQ:
211	0	if (fabs(rh) <= EPS10)
212	0	lp.phi = P->phi0;
213	0	else
214	0	lp.phi = asin(cosc * sinph0 + xy.y * sinc * cosph0 / rh);
215	0	c = cosc - sinph0 * sin(lp.phi);
216	0	if (c != 0. \|\| xy.x != 0.)
217	0	lp.lam = atan2(xy.x * sinc * cosph0, c * rh);
218	0	break;
219	0	case N_POLE:
220	0	xy.y = -xy.y;
221	0	PROJ_FALLTHROUGH;
222	0	case S_POLE:
223	0	if (fabs(rh) <= EPS10)
224	0	lp.phi = P->phi0;
225	0	else
226	0	lp.phi = asin(Q->mode == S_POLE ? -cosc : cosc);
227	0	lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);
228	0	break;
229	0	}
230	0	return lp;
231	0	}
232
233	1.12k	static PJ stere_setup(PJ P) { /* general initialization */
234	1.12k	double t;
235	1.12k	struct pj_stere Q = static_cast<struct pj_stere >(P->opaque);
236
237	1.12k	if (fabs((t = fabs(P->phi0)) - M_HALFPI) < EPS10)
238	447	Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
239	673	else
240	673	Q->mode = t > EPS10 ? OBLIQ : EQUIT;
241	1.12k	Q->phits = fabs(Q->phits);
242
243	1.12k	if (P->es != 0.0) {
244	1.06k	double X;
245
246	1.06k	switch (Q->mode) {
247	380	case N_POLE:
248	433	case S_POLE:
249	433	if (fabs(Q->phits - M_HALFPI) < EPS10)
250	417	Q->akm1 =
251	417	2. * P->k0 /
252	417	sqrt(pow(1 + P->e, 1 + P->e) * pow(1 - P->e, 1 - P->e));
253	16	else {
254	16	t = sin(Q->phits);
255	16	Q->akm1 = cos(Q->phits) / pj_tsfn(Q->phits, t, P->e);
256	16	t *= P->e;
257	16	Q->akm1 /= sqrt(1. - t * t);
258	16	}
259	433	break;
260	518	case EQUIT:
261	632	case OBLIQ:
262	632	t = sin(P->phi0);
263	632	X = 2. * atan(ssfn_(P->phi0, t, P->e)) - M_HALFPI;
264	632	t *= P->e;
265	632	Q->akm1 = 2. * P->k0 * cos(P->phi0) / sqrt(1. - t * t);
266	632	Q->sinX1 = sin(X);
267	632	Q->cosX1 = cos(X);
268	632	break;
269	1.06k	}
270	1.06k	P->inv = stere_e_inverse;
271	1.06k	P->fwd = stere_e_forward;
272	1.06k	} else {
273	55	switch (Q->mode) {
274	4	case OBLIQ:
275	4	sinph0 = sin(P->phi0);
276	4	cosph0 = cos(P->phi0);
277	4	PROJ_FALLTHROUGH;
278	41	case EQUIT:
279	41	Q->akm1 = 2. * P->k0;
280	41	break;
281	0	case S_POLE:
282	14	case N_POLE:
283	14	Q->akm1 = fabs(Q->phits - M_HALFPI) >= EPS10
284	14	? cos(Q->phits) / tan(M_FORTPI - .5 * Q->phits)
285	14	: 2. * P->k0;
286	14	break;
287	55	}
288
289	55	P->inv = stere_s_inverse;
290	55	P->fwd = stere_s_forward;
291	55	}
292	1.12k	return P;
293	1.12k	}
294
295	825	PJ *PJ_PROJECTION(stere) {
296	825	struct pj_stere *Q =
297	825	static_cast<struct pj_stere *>(calloc(1, sizeof(struct pj_stere)));
298	825	if (nullptr == Q)
299	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
300	825	P->opaque = Q;
301
302	825	Q->phits = pj_param(P->ctx, P->params, "tlat_ts").i
303	825	? pj_param(P->ctx, P->params, "rlat_ts").f
304	825	: M_HALFPI;
305
306	825	return stere_setup(P);
307	825	}
308
309	298	PJ *PJ_PROJECTION(ups) {
310	298	struct pj_stere *Q =
311	298	static_cast<struct pj_stere *>(calloc(1, sizeof(struct pj_stere)));
312	298	if (nullptr == Q)
313	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
314	298	P->opaque = Q;
315
316		/* International Ellipsoid */
317	298	P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? -M_HALFPI : M_HALFPI;
318	298	if (P->es == 0.0) {
319	3	proj_log_error(
320	3	P,
321	3	_("Invalid value for es: only ellipsoidal formulation supported"));
322	3	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
323	3	}
324	295	P->k0 = .994;
325	295	P->x0 = 2000000.;
326	295	P->y0 = 2000000.;
327	295	Q->phits = M_HALFPI;
328	295	P->lam0 = 0.;
329
330	295	return stere_setup(P);
331	298	}
332
333		#undef sinph0
334		#undef cosph0
335		#undef EPS10
336		#undef TOL
337		#undef NITER
338		#undef CONV