/src/proj/src/projections/som.cpp

Source
/******************************************************************************
 * This implements the Space Oblique Mercator (SOM) projection, used by the
 * Multi-angle Imaging SpectroRadiometer (MISR) products, from the NASA EOS
 *Terra platform among others (e.g. ASTER).
 *
 * This code was originally developed for the Landsat SOM projection with the
 * following parameters set for Landsat satellites 1, 2, and 3:
 *
 *   inclination angle = 99.092 degrees
 *   period of revolution = 103.2669323 minutes
 *   ascending longitude = 128.87 degrees - (360 / 251) * path_number
 *
 * or for Landsat satellites greater than 3:
 *
 *   inclination angle = 98.2 degrees
 *   period of revolution = 98.8841202 minutes
 *   ascending longitude = 129.3 degrees - (360 / 233) * path_number
 *
 * For the MISR path based SOM projection, the code is identical to that of
 *Landsat SOM with the following parameter changes:
 *
 *   inclination angle = 98.30382 degrees
 *   period of revolution = 98.88 minutes
 *   ascending longitude = 129.3056 degrees - (360 / 233) * path_number
 *
 * and the following code used for Landsat:
 *
 *   Q->rlm = PI * (1. / 248. + .5161290322580645);
 *
 * changed to:
 *
 *   Q->rlm = 0
 *
 * For the generic SOM projection, the code is identical to the above for MISR
 * except that the following parameters are now taken as input rather than
 *derived from path number:
 *
 *   inclination angle
 *   period of revolution
 *   ascending longitude
 *
 * The change of Q->rlm = 0 is kept.
 *
 *****************************************************************************/
/* based upon Snyder and Linck, USGS-NMD */

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

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

PROJ_HEAD(som, "Space Oblique Mercator")
"\n\tCyl, Sph&Ell\n\tinc_angle= ps_rev= asc_lon= ";
PROJ_HEAD(misrsom, "Space oblique for MISR")
"\n\tCyl, Sph&Ell\n\tpath=";
PROJ_HEAD(lsat, "Space oblique for LANDSAT")
"\n\tCyl, Sph&Ell\n\tlsat= path=";

#define TOL 1e-7

namespace { // anonymous namespace
struct pj_som_data {
    double a2, a4, b, c1, c3;
    double q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
    double alf;
};
} // anonymous namespace

static void seraz0(double lam, double mult, PJ *P) {
    struct pj_som_data *Q = static_cast<struct pj_som_data *>(P->opaque);
    double sdsq, h, s, fc, sd, sq, d_1 = 0;

    lam *= DEG_TO_RAD;
    sd = sin(lam);
    sdsq = sd * sd;
    s = Q->p22 * Q->sa * cos(lam) *
        sqrt((1. + Q->t * sdsq) / ((1. + Q->w * sdsq) * (1. + Q->q * sdsq)));
    d_1 = 1. + Q->q * sdsq;
    h = sqrt((1. + Q->q * sdsq) / (1. + Q->w * sdsq)) *
        ((1. + Q->w * sdsq) / (d_1 * d_1) - Q->p22 * Q->ca);
    sq = sqrt(Q->xj * Q->xj + s * s);
    fc = mult * (h * Q->xj - s * s) / sq;
    Q->b += fc;
    Q->a2 += fc * cos(lam + lam);
    Q->a4 += fc * cos(lam * 4.);
    fc = mult * s * (h + Q->xj) / sq;
    Q->c1 += fc * cos(lam);
    Q->c3 += fc * cos(lam * 3.);
}

static PJ_XY som_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_som_data *Q = static_cast<struct pj_som_data *>(P->opaque);
    int l, nn;
    double lamt = 0.0, xlam, sdsq, c, d, s, lamdp = 0.0, phidp, lampp, tanph;
    double lamtp, cl, sd, sp, sav, tanphi;

    if (lp.phi > M_HALFPI)
        lp.phi = M_HALFPI;
    else if (lp.phi < -M_HALFPI)
        lp.phi = -M_HALFPI;
    if (lp.phi >= 0.)
        lampp = M_HALFPI;
    else
        lampp = M_PI_HALFPI;
    tanphi = tan(lp.phi);
    for (nn = 0;;) {
        double fac;
        sav = lampp;
        lamtp = lp.lam + Q->p22 * lampp;
        cl = cos(lamtp);
        if (cl < 0)
            fac = lampp + sin(lampp) * M_HALFPI;
        else
            fac = lampp - sin(lampp) * M_HALFPI;
        for (l = 50; l >= 0; --l) {
            lamt = lp.lam + Q->p22 * sav;
            c = cos(lamt);
            if (fabs(c) < TOL)
                lamt -= TOL;
            xlam = (P->one_es * tanphi * Q->sa + sin(lamt) * Q->ca) / c;
            lamdp = atan(xlam) + fac;
            if (fabs(fabs(sav) - fabs(lamdp)) < TOL)
                break;
            sav = lamdp;
        }
        if (!l || ++nn >= 3 || (lamdp > Q->rlm && lamdp < Q->rlm2))
            break;
        if (lamdp <= Q->rlm)
            lampp = M_TWOPI_HALFPI;
        else if (lamdp >= Q->rlm2)
            lampp = M_HALFPI;
    }
    if (l) {
        sp = sin(lp.phi);
        phidp = aasin(
            P->ctx, (P->one_es * Q->ca * sp - Q->sa * cos(lp.phi) * sin(lamt)) /
                        sqrt(1. - P->es * sp * sp));
        tanph = log(tan(M_FORTPI + .5 * phidp));
        sd = sin(lamdp);
        sdsq = sd * sd;
        s = Q->p22 * Q->sa * cos(lamdp) *
            sqrt((1. + Q->t * sdsq) /
                 ((1. + Q->w * sdsq) * (1. + Q->q * sdsq)));
        d = sqrt(Q->xj * Q->xj + s * s);
        xy.x = Q->b * lamdp + Q->a2 * sin(2. * lamdp) +
               Q->a4 * sin(lamdp * 4.) - tanph * s / d;
        xy.y = Q->c1 * sd + Q->c3 * sin(lamdp * 3.) + tanph * Q->xj / d;
    } else
        xy.x = xy.y = HUGE_VAL;
    return xy;
}

static PJ_LP som_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_som_data *Q = static_cast<struct pj_som_data *>(P->opaque);
    int nn;
    double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;

    lamdp = xy.x / Q->b;
    nn = 50;
    do {
        sav = lamdp;
        sd = sin(lamdp);
        sdsq = sd * sd;
        s = Q->p22 * Q->sa * cos(lamdp) *
            sqrt((1. + Q->t * sdsq) /
                 ((1. + Q->w * sdsq) * (1. + Q->q * sdsq)));
        lamdp = xy.x + xy.y * s / Q->xj - Q->a2 * sin(2. * lamdp) -
                Q->a4 * sin(lamdp * 4.) -
                s / Q->xj * (Q->c1 * sin(lamdp) + Q->c3 * sin(lamdp * 3.));
        lamdp /= Q->b;
    } while (fabs(lamdp - sav) >= TOL && --nn);
    sl = sin(lamdp);
    fac = exp(sqrt(1. + s * s / Q->xj / Q->xj) *
              (xy.y - Q->c1 * sl - Q->c3 * sin(lamdp * 3.)));
    phidp = 2. * (atan(fac) - M_FORTPI);
    dd = sl * sl;
    if (fabs(cos(lamdp)) < TOL)
        lamdp -= TOL;
    spp = sin(phidp);
    sppsq = spp * spp;
    const double denom = 1. - sppsq * (1. + Q->u);
    if (denom == 0.0) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return proj_coord_error().lp;
    }
    lamt = atan(
        ((1. - sppsq * P->rone_es) * tan(lamdp) * Q->ca -
         spp * Q->sa * sqrt((1. + Q->q * dd) * (1. - sppsq) - sppsq * Q->u) /
             cos(lamdp)) /
        denom);
    sl = lamt >= 0. ? 1. : -1.;
    scl = cos(lamdp) >= 0. ? 1. : -1;
    lamt -= M_HALFPI * (1. - scl) * sl;
    lp.lam = lamt - Q->p22 * lamdp;
    if (fabs(Q->sa) < TOL)
        lp.phi =
            aasin(P->ctx, spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
    else
        lp.phi = atan((tan(lamdp) * cos(lamt) - Q->ca * sin(lamt)) /
                      (P->one_es * Q->sa));
    return lp;
}

static PJ *som_setup(PJ *P) {
    double esc, ess, lam;
    struct pj_som_data *Q = static_cast<struct pj_som_data *>(P->opaque);
    Q->sa = sin(Q->alf);
    Q->ca = cos(Q->alf);
    if (fabs(Q->ca) < 1e-9)
        Q->ca = 1e-9;
    esc = P->es * Q->ca * Q->ca;
    ess = P->es * Q->sa * Q->sa;
    Q->w = (1. - esc) * P->rone_es;
    Q->w = Q->w * Q->w - 1.;
    Q->q = ess * P->rone_es;
    Q->t = ess * (2. - P->es) * P->rone_es * P->rone_es;
    Q->u = esc * P->rone_es;
    Q->xj = P->one_es * P->one_es * P->one_es;
    Q->rlm2 = Q->rlm + M_TWOPI;
    Q->a2 = Q->a4 = Q->b = Q->c1 = Q->c3 = 0.;
    seraz0(0., 1., P);
    for (lam = 9.; lam <= 81.0001; lam += 18.)
        seraz0(lam, 4., P);
    for (lam = 18; lam <= 72.0001; lam += 18.)
        seraz0(lam, 2., P);
    seraz0(90., 1., P);
    Q->a2 /= 30.;
    Q->a4 /= 60.;
    Q->b /= 30.;
    Q->c1 /= 15.;
    Q->c3 /= 45.;

    P->inv = som_e_inverse;
    P->fwd = som_e_forward;

    return P;
}

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

    // ascending longitude (radians)
    P->lam0 = pj_param(P->ctx, P->params, "rasc_lon").f;
    if (P->lam0 < -M_TWOPI || P->lam0 > M_TWOPI) {
        proj_log_error(P,
                       _("Invalid value for ascending longitude: should be in "
                         "[-2pi, 2pi] range"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    // inclination angle (radians)
    Q->alf = pj_param(P->ctx, P->params, "rinc_angle").f;
    if (Q->alf < 0 || Q->alf > M_PI) {
        proj_log_error(P, _("Invalid value for inclination angle: should be in "
                            "[0, pi] range"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    // period of revolution (day / rev)
    Q->p22 = pj_param(P->ctx, P->params, "dps_rev").f;
    if (Q->p22 < 0) {
        proj_log_error(P, _("Number of days per rotation should be positive"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    Q->rlm = 0;

    return som_setup(P);
}

PJ *PJ_PROJECTION(misrsom) {
    int path;

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

    path = pj_param(P->ctx, P->params, "ipath").i;
    if (path <= 0 || path > 233) {
        proj_log_error(
            P, _("Invalid value for path: path should be in [1, 233] range"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    P->lam0 = DEG_TO_RAD * 129.3056 - M_TWOPI / 233. * path;
    Q->alf = 98.30382 * DEG_TO_RAD;
    Q->p22 = 98.88 / 1440.0;

    Q->rlm = 0;

    return som_setup(P);
}

PJ *PJ_PROJECTION(lsat) {
    int land, path;
    struct pj_som_data *Q = static_cast<struct pj_som_data *>(
        calloc(1, sizeof(struct pj_som_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    land = pj_param(P->ctx, P->params, "ilsat").i;
    if (land <= 0 || land > 5) {
        proj_log_error(
            P, _("Invalid value for lsat: lsat should be in [1, 5] range"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    path = pj_param(P->ctx, P->params, "ipath").i;
    const int maxPathVal = (land <= 3 ? 251 : 233);
    if (path <= 0 || path > maxPathVal) {
        proj_log_error(
            P, _("Invalid value for path: path should be in [1, %d] range"),
            maxPathVal);
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    if (land <= 3) {
        P->lam0 = DEG_TO_RAD * 128.87 - M_TWOPI / 251. * path;
        Q->p22 = 103.2669323;
        Q->alf = DEG_TO_RAD * 99.092;
    } else {
        P->lam0 = DEG_TO_RAD * 129.3 - M_TWOPI / 233. * path;
        Q->p22 = 98.8841202;
        Q->alf = DEG_TO_RAD * 98.2;
    }
    Q->p22 /= 1440.;

    Q->rlm = M_PI * (1. / 248. + .5161290322580645);

    return som_setup(P);
}

#undef TOL

Coverage Report

Created: 2025-12-31 06:48

Line	Count	Source
1		/******************************************************************************
2		* This implements the Space Oblique Mercator (SOM) projection, used by the
3		* Multi-angle Imaging SpectroRadiometer (MISR) products, from the NASA EOS
4		*Terra platform among others (e.g. ASTER).
5		*
6		* This code was originally developed for the Landsat SOM projection with the
7		* following parameters set for Landsat satellites 1, 2, and 3:
8		*
9		* inclination angle = 99.092 degrees
10		* period of revolution = 103.2669323 minutes
11		* ascending longitude = 128.87 degrees - (360 / 251) * path_number
12		*
13		* or for Landsat satellites greater than 3:
14		*
15		* inclination angle = 98.2 degrees
16		* period of revolution = 98.8841202 minutes
17		* ascending longitude = 129.3 degrees - (360 / 233) * path_number
18		*
19		* For the MISR path based SOM projection, the code is identical to that of
20		*Landsat SOM with the following parameter changes:
21		*
22		* inclination angle = 98.30382 degrees
23		* period of revolution = 98.88 minutes
24		* ascending longitude = 129.3056 degrees - (360 / 233) * path_number
25		*
26		* and the following code used for Landsat:
27		*
28		* Q->rlm = PI * (1. / 248. + .5161290322580645);
29		*
30		* changed to:
31		*
32		* Q->rlm = 0
33		*
34		* For the generic SOM projection, the code is identical to the above for MISR
35		* except that the following parameters are now taken as input rather than
36		*derived from path number:
37		*
38		* inclination angle
39		* period of revolution
40		* ascending longitude
41		*
42		* The change of Q->rlm = 0 is kept.
43		*
44		*****************************************************************************/
45		/* based upon Snyder and Linck, USGS-NMD */
46
47		#include <errno.h>
48		#include <math.h>
49
50		#include "proj.h"
51		#include "proj_internal.h"
52
53		PROJ_HEAD(som, "Space Oblique Mercator")
54		"\n\tCyl, Sph&Ell\n\tinc_angle= ps_rev= asc_lon= ";
55		PROJ_HEAD(misrsom, "Space oblique for MISR")
56		"\n\tCyl, Sph&Ell\n\tpath=";
57		PROJ_HEAD(lsat, "Space oblique for LANDSAT")
58		"\n\tCyl, Sph&Ell\n\tlsat= path=";
59
60	0	#define TOL 1e-7
61
62		namespace { // anonymous namespace
63		struct pj_som_data {
64		double a2, a4, b, c1, c3;
65		double q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
66		double alf;
67		};
68		} // anonymous namespace
69
70	132	static void seraz0(double lam, double mult, PJ *P) {
71	132	struct pj_som_data Q = static_cast<struct pj_som_data >(P->opaque);
72	132	double sdsq, h, s, fc, sd, sq, d_1 = 0;
73
74	132	lam *= DEG_TO_RAD;
75	132	sd = sin(lam);
76	132	sdsq = sd * sd;
77	132	s = Q->p22 * Q->sa * cos(lam) *
78	132	sqrt((1. + Q->t * sdsq) / ((1. + Q->w * sdsq) * (1. + Q->q * sdsq)));
79	132	d_1 = 1. + Q->q * sdsq;
80	132	h = sqrt((1. + Q->q * sdsq) / (1. + Q->w * sdsq)) *
81	132	((1. + Q->w * sdsq) / (d_1 * d_1) - Q->p22 * Q->ca);
82	132	sq = sqrt(Q->xj * Q->xj + s * s);
83	132	fc = mult * (h * Q->xj - s * s) / sq;
84	132	Q->b += fc;
85	132	Q->a2 += fc * cos(lam + lam);
86	132	Q->a4 += fc * cos(lam * 4.);
87	132	fc = mult * s * (h + Q->xj) / sq;
88	132	Q->c1 += fc * cos(lam);
89	132	Q->c3 += fc * cos(lam * 3.);
90	132	}
91
92	0	static PJ_XY som_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
93	0	PJ_XY xy = {0.0, 0.0};
94	0	struct pj_som_data Q = static_cast<struct pj_som_data >(P->opaque);
95	0	int l, nn;
96	0	double lamt = 0.0, xlam, sdsq, c, d, s, lamdp = 0.0, phidp, lampp, tanph;
97	0	double lamtp, cl, sd, sp, sav, tanphi;
98
99	0	if (lp.phi > M_HALFPI)
100	0	lp.phi = M_HALFPI;
101	0	else if (lp.phi < -M_HALFPI)
102	0	lp.phi = -M_HALFPI;
103	0	if (lp.phi >= 0.)
104	0	lampp = M_HALFPI;
105	0	else
106	0	lampp = M_PI_HALFPI;
107	0	tanphi = tan(lp.phi);
108	0	for (nn = 0;;) {
109	0	double fac;
110	0	sav = lampp;
111	0	lamtp = lp.lam + Q->p22 * lampp;
112	0	cl = cos(lamtp);
113	0	if (cl < 0)
114	0	fac = lampp + sin(lampp) * M_HALFPI;
115	0	else
116	0	fac = lampp - sin(lampp) * M_HALFPI;
117	0	for (l = 50; l >= 0; --l) {
118	0	lamt = lp.lam + Q->p22 * sav;
119	0	c = cos(lamt);
120	0	if (fabs(c) < TOL)
121	0	lamt -= TOL;
122	0	xlam = (P->one_es * tanphi * Q->sa + sin(lamt) * Q->ca) / c;
123	0	lamdp = atan(xlam) + fac;
124	0	if (fabs(fabs(sav) - fabs(lamdp)) < TOL)
125	0	break;
126	0	sav = lamdp;
127	0	}
128	0	if (!l \|\| ++nn >= 3 \|\| (lamdp > Q->rlm && lamdp < Q->rlm2))
129	0	break;
130	0	if (lamdp <= Q->rlm)
131	0	lampp = M_TWOPI_HALFPI;
132	0	else if (lamdp >= Q->rlm2)
133	0	lampp = M_HALFPI;
134	0	}
135	0	if (l) {
136	0	sp = sin(lp.phi);
137	0	phidp = aasin(
138	0	P->ctx, (P->one_es * Q->ca * sp - Q->sa * cos(lp.phi) * sin(lamt)) /
139	0	sqrt(1. - P->es * sp * sp));
140	0	tanph = log(tan(M_FORTPI + .5 * phidp));
141	0	sd = sin(lamdp);
142	0	sdsq = sd * sd;
143	0	s = Q->p22 * Q->sa * cos(lamdp) *
144	0	sqrt((1. + Q->t * sdsq) /
145	0	((1. + Q->w * sdsq) * (1. + Q->q * sdsq)));
146	0	d = sqrt(Q->xj * Q->xj + s * s);
147	0	xy.x = Q->b * lamdp + Q->a2 * sin(2. * lamdp) +
148	0	Q->a4 * sin(lamdp * 4.) - tanph * s / d;
149	0	xy.y = Q->c1 * sd + Q->c3 * sin(lamdp * 3.) + tanph * Q->xj / d;
150	0	} else
151	0	xy.x = xy.y = HUGE_VAL;
152	0	return xy;
153	0	}
154
155	0	static PJ_LP som_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
156	0	PJ_LP lp = {0.0, 0.0};
157	0	struct pj_som_data Q = static_cast<struct pj_som_data >(P->opaque);
158	0	int nn;
159	0	double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
160
161	0	lamdp = xy.x / Q->b;
162	0	nn = 50;
163	0	do {
164	0	sav = lamdp;
165	0	sd = sin(lamdp);
166	0	sdsq = sd * sd;
167	0	s = Q->p22 * Q->sa * cos(lamdp) *
168	0	sqrt((1. + Q->t * sdsq) /
169	0	((1. + Q->w * sdsq) * (1. + Q->q * sdsq)));
170	0	lamdp = xy.x + xy.y * s / Q->xj - Q->a2 * sin(2. * lamdp) -
171	0	Q->a4 * sin(lamdp * 4.) -
172	0	s / Q->xj * (Q->c1 * sin(lamdp) + Q->c3 * sin(lamdp * 3.));
173	0	lamdp /= Q->b;
174	0	} while (fabs(lamdp - sav) >= TOL && --nn);
175	0	sl = sin(lamdp);
176	0	fac = exp(sqrt(1. + s * s / Q->xj / Q->xj) *
177	0	(xy.y - Q->c1 * sl - Q->c3 * sin(lamdp * 3.)));
178	0	phidp = 2. * (atan(fac) - M_FORTPI);
179	0	dd = sl * sl;
180	0	if (fabs(cos(lamdp)) < TOL)
181	0	lamdp -= TOL;
182	0	spp = sin(phidp);
183	0	sppsq = spp * spp;
184	0	const double denom = 1. - sppsq * (1. + Q->u);
185	0	if (denom == 0.0) {
186	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
187	0	return proj_coord_error().lp;
188	0	}
189	0	lamt = atan(
190	0	((1. - sppsq * P->rone_es) * tan(lamdp) * Q->ca -
191	0	spp * Q->sa * sqrt((1. + Q->q * dd) * (1. - sppsq) - sppsq * Q->u) /
192	0	cos(lamdp)) /
193	0	denom);
194	0	sl = lamt >= 0. ? 1. : -1.;
195	0	scl = cos(lamdp) >= 0. ? 1. : -1;
196	0	lamt -= M_HALFPI * (1. - scl) * sl;
197	0	lp.lam = lamt - Q->p22 * lamdp;
198	0	if (fabs(Q->sa) < TOL)
199	0	lp.phi =
200	0	aasin(P->ctx, spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
201	0	else
202	0	lp.phi = atan((tan(lamdp) * cos(lamt) - Q->ca * sin(lamt)) /
203	0	(P->one_es * Q->sa));
204	0	return lp;
205	0	}
206
207	12	static PJ som_setup(PJ P) {
208	12	double esc, ess, lam;
209	12	struct pj_som_data Q = static_cast<struct pj_som_data >(P->opaque);
210	12	Q->sa = sin(Q->alf);
211	12	Q->ca = cos(Q->alf);
212	12	if (fabs(Q->ca) < 1e-9)
213	0	Q->ca = 1e-9;
214	12	esc = P->es * Q->ca * Q->ca;
215	12	ess = P->es * Q->sa * Q->sa;
216	12	Q->w = (1. - esc) * P->rone_es;
217	12	Q->w = Q->w * Q->w - 1.;
218	12	Q->q = ess * P->rone_es;
219	12	Q->t = ess * (2. - P->es) * P->rone_es * P->rone_es;
220	12	Q->u = esc * P->rone_es;
221	12	Q->xj = P->one_es * P->one_es * P->one_es;
222	12	Q->rlm2 = Q->rlm + M_TWOPI;
223	12	Q->a2 = Q->a4 = Q->b = Q->c1 = Q->c3 = 0.;
224	12	seraz0(0., 1., P);
225	72	for (lam = 9.; lam <= 81.0001; lam += 18.)
226	60	seraz0(lam, 4., P);
227	60	for (lam = 18; lam <= 72.0001; lam += 18.)
228	48	seraz0(lam, 2., P);
229	12	seraz0(90., 1., P);
230	12	Q->a2 /= 30.;
231	12	Q->a4 /= 60.;
232	12	Q->b /= 30.;
233	12	Q->c1 /= 15.;
234	12	Q->c3 /= 45.;
235
236	12	P->inv = som_e_inverse;
237	12	P->fwd = som_e_forward;
238
239	12	return P;
240	12	}
241
242	12	PJ *PJ_PROJECTION(som) {
243	12	struct pj_som_data Q = static_cast<struct pj_som_data >(
244	12	calloc(1, sizeof(struct pj_som_data)));
245	12	if (nullptr == Q)
246	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
247	12	P->opaque = Q;
248
249		// ascending longitude (radians)
250	12	P->lam0 = pj_param(P->ctx, P->params, "rasc_lon").f;
251	12	if (P->lam0 < -M_TWOPI \|\| P->lam0 > M_TWOPI) {
252	0	proj_log_error(P,
253	0	_("Invalid value for ascending longitude: should be in "
254	0	"[-2pi, 2pi] range"));
255	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
256	0	}
257
258		// inclination angle (radians)
259	12	Q->alf = pj_param(P->ctx, P->params, "rinc_angle").f;
260	12	if (Q->alf < 0 \|\| Q->alf > M_PI) {
261	0	proj_log_error(P, _("Invalid value for inclination angle: should be in "
262	0	"[0, pi] range"));
263	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
264	0	}
265
266		// period of revolution (day / rev)
267	12	Q->p22 = pj_param(P->ctx, P->params, "dps_rev").f;
268	12	if (Q->p22 < 0) {
269	0	proj_log_error(P, _("Number of days per rotation should be positive"));
270	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
271	0	}
272
273	12	Q->rlm = 0;
274
275	12	return som_setup(P);
276	12	}
277
278	0	PJ *PJ_PROJECTION(misrsom) {
279	0	int path;
280
281	0	struct pj_som_data Q = static_cast<struct pj_som_data >(
282	0	calloc(1, sizeof(struct pj_som_data)));
283	0	if (nullptr == Q)
284	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
285	0	P->opaque = Q;
286
287	0	path = pj_param(P->ctx, P->params, "ipath").i;
288	0	if (path <= 0 \|\| path > 233) {
289	0	proj_log_error(
290	0	P, _("Invalid value for path: path should be in [1, 233] range"));
291	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
292	0	}
293
294	0	P->lam0 = DEG_TO_RAD * 129.3056 - M_TWOPI / 233. * path;
295	0	Q->alf = 98.30382 * DEG_TO_RAD;
296	0	Q->p22 = 98.88 / 1440.0;
297
298	0	Q->rlm = 0;
299
300	0	return som_setup(P);
301	0	}
302
303	4	PJ *PJ_PROJECTION(lsat) {
304	4	int land, path;
305	4	struct pj_som_data Q = static_cast<struct pj_som_data >(
306	4	calloc(1, sizeof(struct pj_som_data)));
307	4	if (nullptr == Q)
308	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
309	4	P->opaque = Q;
310
311	4	land = pj_param(P->ctx, P->params, "ilsat").i;
312	4	if (land <= 0 \|\| land > 5) {
313	4	proj_log_error(
314	4	P, _("Invalid value for lsat: lsat should be in [1, 5] range"));
315	4	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
316	4	}
317
318	0	path = pj_param(P->ctx, P->params, "ipath").i;
319	0	const int maxPathVal = (land <= 3 ? 251 : 233);
320	0	if (path <= 0 \|\| path > maxPathVal) {
321	0	proj_log_error(
322	0	P, _("Invalid value for path: path should be in [1, %d] range"),
323	0	maxPathVal);
324	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
325	0	}
326
327	0	if (land <= 3) {
328	0	P->lam0 = DEG_TO_RAD * 128.87 - M_TWOPI / 251. * path;
329	0	Q->p22 = 103.2669323;
330	0	Q->alf = DEG_TO_RAD * 99.092;
331	0	} else {
332	0	P->lam0 = DEG_TO_RAD * 129.3 - M_TWOPI / 233. * path;
333	0	Q->p22 = 98.8841202;
334	0	Q->alf = DEG_TO_RAD * 98.2;
335	0	}
336	0	Q->p22 /= 1440.;
337
338	0	Q->rlm = M_PI * (1. / 248. + .5161290322580645);
339
340	0	return som_setup(P);
341	0	}
342
343		#undef TOL