/src/proj/src/ell_set.cpp

Source
/* set ellipsoid parameters a and es */

#include <math.h>
#include <stddef.h>
#include <string.h>

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

/* Prototypes of the pj_ellipsoid helper functions */
static int ellps_ellps(PJ *P);
static int ellps_size(PJ *P);
static int ellps_shape(PJ *P);
static int ellps_spherification(PJ *P);

static paralist *pj_get_param(paralist *list, const char *key);
static char *pj_param_value(paralist *list);
static const PJ_ELLPS *pj_find_ellps(const char *name);

/***************************************************************************************/
int pj_ellipsoid(PJ *P) {
    /****************************************************************************************
        This is a replacement for the classic PROJ pj_ell_set function. The main
    difference is that pj_ellipsoid augments the PJ object with a copy of the
    exact tags used to define its related ellipsoid.

        This makes it possible to let a new PJ object inherit the geometrical
    properties of an existing one.

        A complete ellipsoid definition comprises a size (primary) and a shape
    (secondary) parameter.

        Size parameters supported are:
            R, defining the radius of a spherical planet
            a, defining the semimajor axis of an ellipsoidal planet

        Shape parameters supported are:
            rf, the reverse flattening of the ellipsoid
            f,  the flattening of the ellipsoid
            es, the eccentricity squared
            e,  the eccentricity
            b,  the semiminor axis

        The ellps=xxx parameter provides both size and shape for a number of
    built in ellipsoid definitions.

        The ellipsoid definition may be augmented with a spherification flag,
    turning the ellipsoid into a sphere with features defined by the ellipsoid.

        Spherification parameters supported are:
            R_A, which gives a sphere with the same surface area as the
            ellipsoid R_V, which gives a sphere with the same volume as the
    ellipsoid

            R_a, which gives a sphere with R = (a + b)/2   (arithmetic mean)
            R_g, which gives a sphere with R = sqrt(a*b)   (geometric mean)
            R_h, which gives a sphere with R = 2*a*b/(a+b) (harmonic mean)

            R_lat_a=phi, which gives a sphere with R being the arithmetic mean
            of of the corresponding ellipsoid at latitude phi.
            R_lat_g=phi, which gives
            a sphere with R being the geometric mean of of the corresponding
    ellipsoid at latitude phi.

            R_C, which gives a sphere with the radius of the conformal sphere
            at phi0.

        If R is given as size parameter, any shape and spherification parameters
        given are ignored.

        If size and shape are given as ellps=xxx, later shape and size
    parameters are are taken into account as modifiers for the built in
    ellipsoid definition.

        While this may seem strange, it is in accordance with historical PROJ
        behavior. It can e.g. be used to define coordinates on the ellipsoid
        scaled to unit semimajor axis by specifying "+ellps=xxx +a=1"

    ****************************************************************************************/
    int err = proj_errno_reset(P);
    const char *empty = {""};

    free(P->def_size);
    P->def_size = nullptr;
    free(P->def_shape);
    P->def_shape = nullptr;
    free(P->def_spherification);
    P->def_spherification = nullptr;
    free(P->def_ellps);
    P->def_ellps = nullptr;

    /* Specifying R overrules everything */
    if (pj_get_param(P->params, "R")) {
        if (0 != ellps_size(P))
            return 1;
        pj_calc_ellipsoid_params(P, P->a, 0);
        if (proj_errno(P))
            return 1;
        return proj_errno_restore(P, err);
    }

    /* If an ellps argument is specified, start by using that */
    if (0 != ellps_ellps(P))
        return 1;

    /* We may overwrite the size */
    if (0 != ellps_size(P))
        return 2;

    /* We may also overwrite the shape */
    if (0 != ellps_shape(P))
        return 3;

    /* When we're done with it, we compute all related ellipsoid parameters */
    pj_calc_ellipsoid_params(P, P->a, P->es);

    /* And finally, we may turn it into a sphere */
    if (0 != ellps_spherification(P))
        return 4;

    proj_log_trace(P, "pj_ellipsoid - final: a=%.3f f=1/%7.3f, errno=%d", P->a,
                   P->f != 0 ? 1 / P->f : 0, proj_errno(P));
    proj_log_trace(P, "pj_ellipsoid - final: %s %s %s %s",
                   P->def_size ? P->def_size : empty,
                   P->def_shape ? P->def_shape : empty,
                   P->def_spherification ? P->def_spherification : empty,
                   P->def_ellps ? P->def_ellps : empty);

    if (proj_errno(P))
        return 5;

    /* success */
    return proj_errno_restore(P, err);
}

/***************************************************************************************/
static int ellps_ellps(PJ *P) {
    /***************************************************************************************/
    const PJ_ELLPS *ellps;
    paralist *par = nullptr;
    char *name;
    int err;

    /* Sail home if ellps=xxx is not specified */
    par = pj_get_param(P->params, "ellps");
    if (nullptr == par)
        return 0;

    /* Then look up the right size and shape parameters from the builtin list */
    if (strlen(par->param) < 7) {
        proj_log_error(P, _("Invalid value for +ellps"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    name = par->param + 6;
    ellps = pj_find_ellps(name);
    if (nullptr == ellps) {
        proj_log_error(P, _("Unrecognized value for +ellps"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    /* Now, get things ready for ellps_size/ellps_shape, make them do their
     * thing
     */
    err = proj_errno_reset(P);

    paralist *new_params = pj_mkparam(ellps->major);
    if (nullptr == new_params)
        return proj_errno_set(P, PROJ_ERR_OTHER /*ENOMEM*/);
    new_params->next = pj_mkparam(ellps->ell);
    if (nullptr == new_params->next) {
        free(new_params);
        return proj_errno_set(P, PROJ_ERR_OTHER /*ENOMEM*/);
    }
    paralist *old_params = P->params;
    P->params = new_params;

    {
        PJ empty_PJ;
        pj_inherit_ellipsoid_def(&empty_PJ, P);
    }
    const bool errorOnSizeOrShape = (ellps_size(P) || ellps_shape(P));
    P->params = old_params;
    free(new_params->next);
    free(new_params);
    if (errorOnSizeOrShape)
        return proj_errno_set(P, PROJ_ERR_OTHER /*ENOMEM*/);
    if (proj_errno(P))
        return proj_errno(P);

    /* Finally update P and sail home */
    P->def_ellps = pj_strdup(par->param);
    par->used = 1;

    return proj_errno_restore(P, err);
}

/***************************************************************************************/
static int ellps_size(PJ *P) {
    /***************************************************************************************/
    paralist *par = nullptr;
    int a_was_set = 0;

    free(P->def_size);
    P->def_size = nullptr;

    /* A size parameter *must* be given, but may have been given as ellps prior
     */
    if (P->a != 0)
        a_was_set = 1;

    /* Check which size key is specified */
    par = pj_get_param(P->params, "R");
    if (nullptr == par)
        par = pj_get_param(P->params, "a");
    if (nullptr == par) {
        if (a_was_set)
            return 0;
        if (P->need_ellps)
            proj_log_error(P, _("Major axis not given"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    P->def_size = pj_strdup(par->param);
    par->used = 1;
    P->a = pj_atof(pj_param_value(par));
    if (P->a <= 0) {
        proj_log_error(P, _("Invalid value for major axis"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    if (HUGE_VAL == P->a) {
        proj_log_error(P, _("Invalid value for major axis"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    if ('R' == par->param[0]) {
        P->es = P->f = P->e = P->rf = 0;
        P->b = P->a;
    }
    return 0;
}

/***************************************************************************************/
static int ellps_shape(PJ *P) {
    /***************************************************************************************/
    const char *keys[] = {"rf", "f", "es", "e", "b"};
    paralist *par = nullptr;
    size_t i, len;

    par = nullptr;
    len = sizeof(keys) / sizeof(char *);

    free(P->def_shape);
    P->def_shape = nullptr;

    /* Check which shape key is specified */
    for (i = 0; i < len; i++) {
        par = pj_get_param(P->params, keys[i]);
        if (par)
            break;
    }

    /* Not giving a shape parameter means selecting a sphere, unless shape */
    /* has been selected previously via ellps=xxx */
    if (nullptr == par && P->es != 0)
        return 0;
    if (nullptr == par && P->es == 0) {
        P->es = P->f = 0;
        P->b = P->a;
        return 0;
    }

    P->def_shape = pj_strdup(par->param);
    par->used = 1;
    P->es = P->f = P->b = P->e = P->rf = 0;

    switch (i) {

    /* reverse flattening, rf */
    case 0:
        P->rf = pj_atof(pj_param_value(par));
        if (HUGE_VAL == P->rf || P->rf <= 0) {
            proj_log_error(P, _("Invalid value for rf. Should be > 0"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        P->f = 1 / P->rf;
        P->es = 2 * P->f - P->f * P->f;
        break;

    /* flattening, f */
    case 1:
        P->f = pj_atof(pj_param_value(par));
        if (HUGE_VAL == P->f || P->f < 0) {
            proj_log_error(P, _("Invalid value for f. Should be >= 0"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }

        P->rf = P->f != 0.0 ? 1.0 / P->f : HUGE_VAL;
        P->es = 2 * P->f - P->f * P->f;
        break;

    /* eccentricity squared, es */
    case 2:
        P->es = pj_atof(pj_param_value(par));
        if (HUGE_VAL == P->es || P->es < 0 || P->es >= 1) {
            proj_log_error(P,
                           _("Invalid value for es. Should be in [0,1[ range"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        break;

    /* eccentricity, e */
    case 3:
        P->e = pj_atof(pj_param_value(par));
        if (HUGE_VAL == P->e || P->e < 0 || P->e >= 1) {
            proj_log_error(P,
                           _("Invalid value for e. Should be in [0,1[ range"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        P->es = P->e * P->e;
        break;

    /* semiminor axis, b */
    case 4:
        P->b = pj_atof(pj_param_value(par));
        if (HUGE_VAL == P->b || P->b <= 0) {
            proj_log_error(P, _("Invalid value for b. Should be > 0"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        if (P->b == P->a)
            break;
        // coverity[division_by_zero]
        P->f = (P->a - P->b) / P->a;
        P->es = 2 * P->f - P->f * P->f;
        break;
    default:
        // shouldn't happen
        return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    }

    // Written that way to catch NaN
    if (!(P->es >= 0)) {
        proj_log_error(P, _("Invalid eccentricity"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    return 0;
}

/* series coefficients for calculating ellipsoid-equivalent spheres */
static const double SIXTH = 1 / 6.;
static const double RA4 = 17 / 360.;
static const double RA6 = 67 / 3024.;
static const double RV4 = 5 / 72.;
static const double RV6 = 55 / 1296.;

/***************************************************************************************/
static int ellps_spherification(PJ *P) {
    /***************************************************************************************/
    const char *keys[] = {"R_A", "R_V",     "R_a",     "R_g",
                          "R_h", "R_lat_a", "R_lat_g", "R_C"};
    size_t len, i;
    paralist *par = nullptr;

    double t;
    char *v, *endp;

    len = sizeof(keys) / sizeof(char *);

    /* Check which spherification key is specified */
    for (i = 0; i < len; i++) {
        par = pj_get_param(P->params, keys[i]);
        if (par)
            break;
    }

    /* No spherification specified? Then we're done */
    if (i == len)
        return 0;

    /* Store definition */
    P->def_spherification = pj_strdup(par->param);
    par->used = 1;

    switch (i) {

    /* R_A - a sphere with same area as ellipsoid */
    case 0:
        P->a *= 1. - P->es * (SIXTH + P->es * (RA4 + P->es * RA6));
        break;

    /* R_V - a sphere with same volume as ellipsoid */
    case 1:
        P->a *= 1. - P->es * (SIXTH + P->es * (RV4 + P->es * RV6));
        break;

    /* R_a - a sphere with R = the arithmetic mean of the ellipsoid */
    case 2:
        P->a = (P->a + P->b) / 2;
        break;

    /* R_g - a sphere with R = the geometric mean of the ellipsoid */
    case 3:
        P->a = sqrt(P->a * P->b);
        break;

    /* R_h - a sphere with R = the harmonic mean of the ellipsoid */
    case 4:
        if (P->a + P->b == 0)
            return proj_errno_set(
                P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        P->a = (2 * P->a * P->b) / (P->a + P->b);
        break;

    /* R_lat_a - a sphere with R = the arithmetic mean of the ellipsoid at given
     * latitude */
    case 5:
    /* R_lat_g - a sphere with R = the geometric  mean of the ellipsoid at given
     * latitude */
    case 6:
        v = pj_param_value(par);
        t = proj_dmstor(v, &endp);
        if (fabs(t) > M_HALFPI) {
            proj_log_error(
                P, _("Invalid value for lat_g. |lat_g| should be <= 90°"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        t = sin(t);
        t = 1 - P->es * t * t;
        if (t == 0.) {
            proj_log_error(P, _("Invalid eccentricity"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        if (i == 5) /* arithmetic */
            P->a *= (1. - P->es + t) / (2 * t * sqrt(t));
        else /* geometric */
            P->a *= sqrt(1 - P->es) / t;
        break;

    /* R_C - a sphere with R = radius of conformal sphere, taken at a
     * latitude that is phi0 (note: at least for mercator. for other
     * projection methods, this could be phi1)
     * Formula from IOGP Publication 373-7-2 – Geomatics Guidance Note number 7,
     * part 2 1.1 Ellipsoid parameters
     */
    case 7:
        t = sin(P->phi0);
        t = 1 - P->es * t * t;
        if (t == 0.) {
            proj_log_error(P, _("Invalid eccentricity"));
            return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
        P->a *= sqrt(1 - P->es) / t;
        break;
    }

    if (P->a <= 0.) {
        proj_log_error(P, _("Invalid or missing major axis"));
        return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    /* Clean up the ellipsoidal parameters to reflect the sphere */
    P->es = P->e = P->f = 0;
    P->rf = HUGE_VAL;
    P->b = P->a;
    pj_calc_ellipsoid_params(P, P->a, 0);

    return 0;
}

/* locate parameter in list */
static paralist *pj_get_param(paralist *list, const char *key) {
    size_t l = strlen(key);
    while (list && !(0 == strncmp(list->param, key, l) &&
                     (0 == list->param[l] || list->param[l] == '=')))
        list = list->next;
    return list;
}

static char *pj_param_value(paralist *list) {
    char *key, *value;
    if (nullptr == list)
        return nullptr;

    key = list->param;
    value = strchr(key, '=');

    /* a flag (i.e. a key without value) has its own name (key) as value */
    return value ? value + 1 : key;
}

static const PJ_ELLPS *pj_find_ellps(const char *name) {
    int i;
    const char *s;
    const PJ_ELLPS *ellps;

    if (nullptr == name)
        return nullptr;

    ellps = proj_list_ellps();

    /* Search through internal ellipsoid list for name */
    for (i = 0; (s = ellps[i].id) && strcmp(name, s); ++i)
        ;
    if (nullptr == s)
        return nullptr;
    return ellps + i;
}

/**************************************************************************************/
void pj_inherit_ellipsoid_def(const PJ *src, PJ *dst) {
    /***************************************************************************************
        Brute force copy the ellipsoidal parameters from src to dst.  This code
    was written before the actual ellipsoid setup parameters were kept available
    in the PJ->def_xxx elements.
    ***************************************************************************************/

    /* The linear parameters */
    dst->a = src->a;
    dst->b = src->b;
    dst->ra = src->ra;
    dst->rb = src->rb;

    /* The eccentricities */
    dst->alpha = src->alpha;
    dst->e = src->e;
    dst->es = src->es;
    dst->e2 = src->e2;
    dst->e2s = src->e2s;
    dst->e3 = src->e3;
    dst->e3s = src->e3s;
    dst->one_es = src->one_es;
    dst->rone_es = src->rone_es;

    /* The flattenings */
    dst->f = src->f;
    dst->f2 = src->f2;
    dst->n = src->n;
    dst->rf = src->rf;
    dst->rf2 = src->rf2;
    dst->rn = src->rn;

    /* This one's for GRS80 */
    dst->J = src->J;

    /* es and a before any +proj related adjustment */
    dst->es_orig = src->es_orig;
    dst->a_orig = src->a_orig;
}

/***************************************************************************************/
int pj_calc_ellipsoid_params(PJ *P, double a, double es) {
    /****************************************************************************************
        Calculate a large number of ancillary ellipsoidal parameters, in
    addition to the two traditional PROJ defining parameters: Semimajor axis, a,
    and the eccentricity squared, es.

        Most of these parameters are fairly cheap to compute in comparison to
    the overall effort involved in initializing a PJ object. They may, however,
    take a substantial part of the time taken in computing an individual point
    transformation.

        So by providing them up front, we can amortize the (already modest) cost
    over all transformations carried out over the entire lifetime of a PJ
    object, rather than incur that cost for every single transformation.

        Most of the parameter calculations here are based on the "angular
    eccentricity", i.e. the angle, measured from the semiminor axis, of a line
    going from the north pole to one of the foci of the ellipsoid - or in other
    words: The arc sine of the eccentricity.

        The formulae used are mostly taken from:

        Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
        Columbus, Ohio:  Dept. of Geodetic Science
        and Surveying, Ohio State University.

    ****************************************************************************************/

    P->a = a;
    P->es = es;

    /* Compute some ancillary ellipsoidal parameters */
    if (P->e == 0)
        P->e = sqrt(P->es); /* eccentricity */
    P->alpha = asin(P->e);  /* angular eccentricity */

    /* second eccentricity */
    P->e2 = tan(P->alpha);
    P->e2s = P->e2 * P->e2;

    /* third eccentricity */
    P->e3 = (0 != P->alpha)
                ? sin(P->alpha) / sqrt(2 - sin(P->alpha) * sin(P->alpha))
                : 0;
    P->e3s = P->e3 * P->e3;

    /* flattening */
    if (0 == P->f)
        P->f = 1 - cos(P->alpha); /* = 1 - sqrt (1 - PIN->es); */
    if (!(P->f >= 0.0 && P->f < 1.0)) {
        proj_log_error(P, _("Invalid eccentricity"));
        proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    }
    P->rf = P->f != 0.0 ? 1.0 / P->f : HUGE_VAL;

    /* second flattening */
    P->f2 = (cos(P->alpha) != 0) ? 1 / cos(P->alpha) - 1 : 0;
    P->rf2 = P->f2 != 0.0 ? 1 / P->f2 : HUGE_VAL;

    /* third flattening */
    P->n = pow(tan(P->alpha / 2), 2);
    P->rn = P->n != 0.0 ? 1 / P->n : HUGE_VAL;

    /* ...and a few more */
    if (0 == P->b)
        P->b = (1 - P->f) * P->a;
    P->rb = 1. / P->b;
    P->ra = 1. / P->a;

    P->one_es = 1. - P->es;
    if (P->one_es == 0.) {
        proj_log_error(P, _("Invalid eccentricity"));
        proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    }

    P->rone_es = 1. / P->one_es;

    return 0;
}

/**************************************************************************************/
int pj_ell_set(PJ_CONTEXT *ctx, paralist *pl, double *a, double *es) {
    /***************************************************************************************
        Initialize ellipsoidal parameters by emulating the original ellipsoid
    setup function by Gerald Evenden, through a call to pj_ellipsoid
    ***************************************************************************************/
    PJ B;
    int ret;

    B.ctx = ctx;
    B.params = pl;

    ret = pj_ellipsoid(&B);

    free(B.def_size);
    free(B.def_shape);
    free(B.def_spherification);
    free(B.def_ellps);

    if (ret)
        return ret;

    *a = B.a;
    *es = B.es;
    return 0;
}

Coverage Report

Created: 2026-07-14 06:21

Line	Count	Source
1		/* set ellipsoid parameters a and es */
2
3		#include <math.h>
4		#include <stddef.h>
5		#include <string.h>
6
7		#include "proj.h"
8		#include "proj_internal.h"
9
10		/* Prototypes of the pj_ellipsoid helper functions */
11		static int ellps_ellps(PJ *P);
12		static int ellps_size(PJ *P);
13		static int ellps_shape(PJ *P);
14		static int ellps_spherification(PJ *P);
15
16		static paralist pj_get_param(paralist list, const char *key);
17		static char pj_param_value(paralist list);
18		static const PJ_ELLPS pj_find_ellps(const char name);
19
20		/***************************************************************************************/
21	0	int pj_ellipsoid(PJ *P) {
22		/****************************************************************************************
23		This is a replacement for the classic PROJ pj_ell_set function. The main
24		difference is that pj_ellipsoid augments the PJ object with a copy of the
25		exact tags used to define its related ellipsoid.
26
27		This makes it possible to let a new PJ object inherit the geometrical
28		properties of an existing one.
29
30		A complete ellipsoid definition comprises a size (primary) and a shape
31		(secondary) parameter.
32
33		Size parameters supported are:
34		R, defining the radius of a spherical planet
35		a, defining the semimajor axis of an ellipsoidal planet
36
37		Shape parameters supported are:
38		rf, the reverse flattening of the ellipsoid
39		f, the flattening of the ellipsoid
40		es, the eccentricity squared
41		e, the eccentricity
42		b, the semiminor axis
43
44		The ellps=xxx parameter provides both size and shape for a number of
45		built in ellipsoid definitions.
46
47		The ellipsoid definition may be augmented with a spherification flag,
48		turning the ellipsoid into a sphere with features defined by the ellipsoid.
49
50		Spherification parameters supported are:
51		R_A, which gives a sphere with the same surface area as the
52		ellipsoid R_V, which gives a sphere with the same volume as the
53		ellipsoid
54
55		R_a, which gives a sphere with R = (a + b)/2 (arithmetic mean)
56		R_g, which gives a sphere with R = sqrt(a*b) (geometric mean)
57		R_h, which gives a sphere with R = 2ab/(a+b) (harmonic mean)
58
59		R_lat_a=phi, which gives a sphere with R being the arithmetic mean
60		of of the corresponding ellipsoid at latitude phi.
61		R_lat_g=phi, which gives
62		a sphere with R being the geometric mean of of the corresponding
63		ellipsoid at latitude phi.
64
65		R_C, which gives a sphere with the radius of the conformal sphere
66		at phi0.
67
68		If R is given as size parameter, any shape and spherification parameters
69		given are ignored.
70
71		If size and shape are given as ellps=xxx, later shape and size
72		parameters are are taken into account as modifiers for the built in
73		ellipsoid definition.
74
75		While this may seem strange, it is in accordance with historical PROJ
76		behavior. It can e.g. be used to define coordinates on the ellipsoid
77		scaled to unit semimajor axis by specifying "+ellps=xxx +a=1"
78
79		****************************************************************************************/
80	0	int err = proj_errno_reset(P);
81	0	const char *empty = {""};
82
83	0	free(P->def_size);
84	0	P->def_size = nullptr;
85	0	free(P->def_shape);
86	0	P->def_shape = nullptr;
87	0	free(P->def_spherification);
88	0	P->def_spherification = nullptr;
89	0	free(P->def_ellps);
90	0	P->def_ellps = nullptr;
91
92		/* Specifying R overrules everything */
93	0	if (pj_get_param(P->params, "R")) {
94	0	if (0 != ellps_size(P))
95	0	return 1;
96	0	pj_calc_ellipsoid_params(P, P->a, 0);
97	0	if (proj_errno(P))
98	0	return 1;
99	0	return proj_errno_restore(P, err);
100	0	}
101
102		/* If an ellps argument is specified, start by using that */
103	0	if (0 != ellps_ellps(P))
104	0	return 1;
105
106		/* We may overwrite the size */
107	0	if (0 != ellps_size(P))
108	0	return 2;
109
110		/* We may also overwrite the shape */
111	0	if (0 != ellps_shape(P))
112	0	return 3;
113
114		/* When we're done with it, we compute all related ellipsoid parameters */
115	0	pj_calc_ellipsoid_params(P, P->a, P->es);
116
117		/* And finally, we may turn it into a sphere */
118	0	if (0 != ellps_spherification(P))
119	0	return 4;
120
121	0	proj_log_trace(P, "pj_ellipsoid - final: a=%.3f f=1/%7.3f, errno=%d", P->a,
122	0	P->f != 0 ? 1 / P->f : 0, proj_errno(P));
123	0	proj_log_trace(P, "pj_ellipsoid - final: %s %s %s %s",
124	0	P->def_size ? P->def_size : empty,
125	0	P->def_shape ? P->def_shape : empty,
126	0	P->def_spherification ? P->def_spherification : empty,
127	0	P->def_ellps ? P->def_ellps : empty);
128
129	0	if (proj_errno(P))
130	0	return 5;
131
132		/* success */
133	0	return proj_errno_restore(P, err);
134	0	}
135
136		/***************************************************************************************/
137	0	static int ellps_ellps(PJ *P) {
138		/***************************************************************************************/
139	0	const PJ_ELLPS *ellps;
140	0	paralist *par = nullptr;
141	0	char *name;
142	0	int err;
143
144		/* Sail home if ellps=xxx is not specified */
145	0	par = pj_get_param(P->params, "ellps");
146	0	if (nullptr == par)
147	0	return 0;
148
149		/* Then look up the right size and shape parameters from the builtin list */
150	0	if (strlen(par->param) < 7) {
151	0	proj_log_error(P, _("Invalid value for +ellps"));
152	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
153	0	}
154	0	name = par->param + 6;
155	0	ellps = pj_find_ellps(name);
156	0	if (nullptr == ellps) {
157	0	proj_log_error(P, _("Unrecognized value for +ellps"));
158	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
159	0	}
160
161		/* Now, get things ready for ellps_size/ellps_shape, make them do their
162		* thing
163		*/
164	0	err = proj_errno_reset(P);
165
166	0	paralist *new_params = pj_mkparam(ellps->major);
167	0	if (nullptr == new_params)
168	0	return proj_errno_set(P, PROJ_ERR_OTHER /ENOMEM/);
169	0	new_params->next = pj_mkparam(ellps->ell);
170	0	if (nullptr == new_params->next) {
171	0	free(new_params);
172	0	return proj_errno_set(P, PROJ_ERR_OTHER /ENOMEM/);
173	0	}
174	0	paralist *old_params = P->params;
175	0	P->params = new_params;
176
177	0	{
178	0	PJ empty_PJ;
179	0	pj_inherit_ellipsoid_def(&empty_PJ, P);
180	0	}
181	0	const bool errorOnSizeOrShape = (ellps_size(P) \|\| ellps_shape(P));
182	0	P->params = old_params;
183	0	free(new_params->next);
184	0	free(new_params);
185	0	if (errorOnSizeOrShape)
186	0	return proj_errno_set(P, PROJ_ERR_OTHER /ENOMEM/);
187	0	if (proj_errno(P))
188	0	return proj_errno(P);
189
190		/* Finally update P and sail home */
191	0	P->def_ellps = pj_strdup(par->param);
192	0	par->used = 1;
193
194	0	return proj_errno_restore(P, err);
195	0	}
196
197		/***************************************************************************************/
198	0	static int ellps_size(PJ *P) {
199		/***************************************************************************************/
200	0	paralist *par = nullptr;
201	0	int a_was_set = 0;
202
203	0	free(P->def_size);
204	0	P->def_size = nullptr;
205
206		/* A size parameter must be given, but may have been given as ellps prior
207		*/
208	0	if (P->a != 0)
209	0	a_was_set = 1;
210
211		/* Check which size key is specified */
212	0	par = pj_get_param(P->params, "R");
213	0	if (nullptr == par)
214	0	par = pj_get_param(P->params, "a");
215	0	if (nullptr == par) {
216	0	if (a_was_set)
217	0	return 0;
218	0	if (P->need_ellps)
219	0	proj_log_error(P, _("Major axis not given"));
220	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
221	0	}
222
223	0	P->def_size = pj_strdup(par->param);
224	0	par->used = 1;
225	0	P->a = pj_atof(pj_param_value(par));
226	0	if (P->a <= 0) {
227	0	proj_log_error(P, _("Invalid value for major axis"));
228	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
229	0	}
230	0	if (HUGE_VAL == P->a) {
231	0	proj_log_error(P, _("Invalid value for major axis"));
232	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
233	0	}
234
235	0	if ('R' == par->param[0]) {
236	0	P->es = P->f = P->e = P->rf = 0;
237	0	P->b = P->a;
238	0	}
239	0	return 0;
240	0	}
241
242		/***************************************************************************************/
243	0	static int ellps_shape(PJ *P) {
244		/***************************************************************************************/
245	0	const char *keys[] = {"rf", "f", "es", "e", "b"};
246	0	paralist *par = nullptr;
247	0	size_t i, len;
248
249	0	par = nullptr;
250	0	len = sizeof(keys) / sizeof(char *);
251
252	0	free(P->def_shape);
253	0	P->def_shape = nullptr;
254
255		/* Check which shape key is specified */
256	0	for (i = 0; i < len; i++) {
257	0	par = pj_get_param(P->params, keys[i]);
258	0	if (par)
259	0	break;
260	0	}
261
262		/* Not giving a shape parameter means selecting a sphere, unless shape */
263		/* has been selected previously via ellps=xxx */
264	0	if (nullptr == par && P->es != 0)
265	0	return 0;
266	0	if (nullptr == par && P->es == 0) {
267	0	P->es = P->f = 0;
268	0	P->b = P->a;
269	0	return 0;
270	0	}
271
272	0	P->def_shape = pj_strdup(par->param);
273	0	par->used = 1;
274	0	P->es = P->f = P->b = P->e = P->rf = 0;
275
276	0	switch (i) {
277
278		/* reverse flattening, rf */
279	0	case 0:
280	0	P->rf = pj_atof(pj_param_value(par));
281	0	if (HUGE_VAL == P->rf \|\| P->rf <= 0) {
282	0	proj_log_error(P, _("Invalid value for rf. Should be > 0"));
283	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
284	0	}
285	0	P->f = 1 / P->rf;
286	0	P->es = 2 * P->f - P->f * P->f;
287	0	break;
288
289		/* flattening, f */
290	0	case 1:
291	0	P->f = pj_atof(pj_param_value(par));
292	0	if (HUGE_VAL == P->f \|\| P->f < 0) {
293	0	proj_log_error(P, _("Invalid value for f. Should be >= 0"));
294	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
295	0	}
296
297	0	P->rf = P->f != 0.0 ? 1.0 / P->f : HUGE_VAL;
298	0	P->es = 2 * P->f - P->f * P->f;
299	0	break;
300
301		/* eccentricity squared, es */
302	0	case 2:
303	0	P->es = pj_atof(pj_param_value(par));
304	0	if (HUGE_VAL == P->es \|\| P->es < 0 \|\| P->es >= 1) {
305	0	proj_log_error(P,
306	0	_("Invalid value for es. Should be in [0,1[ range"));
307	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
308	0	}
309	0	break;
310
311		/* eccentricity, e */
312	0	case 3:
313	0	P->e = pj_atof(pj_param_value(par));
314	0	if (HUGE_VAL == P->e \|\| P->e < 0 \|\| P->e >= 1) {
315	0	proj_log_error(P,
316	0	_("Invalid value for e. Should be in [0,1[ range"));
317	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
318	0	}
319	0	P->es = P->e * P->e;
320	0	break;
321
322		/* semiminor axis, b */
323	0	case 4:
324	0	P->b = pj_atof(pj_param_value(par));
325	0	if (HUGE_VAL == P->b \|\| P->b <= 0) {
326	0	proj_log_error(P, _("Invalid value for b. Should be > 0"));
327	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
328	0	}
329	0	if (P->b == P->a)
330	0	break;
331		// coverity[division_by_zero]
332	0	P->f = (P->a - P->b) / P->a;
333	0	P->es = 2 * P->f - P->f * P->f;
334	0	break;
335	0	default:
336		// shouldn't happen
337	0	return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
338	0	}
339
340		// Written that way to catch NaN
341	0	if (!(P->es >= 0)) {
342	0	proj_log_error(P, _("Invalid eccentricity"));
343	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
344	0	}
345	0	return 0;
346	0	}
347
348		/* series coefficients for calculating ellipsoid-equivalent spheres */
349		static const double SIXTH = 1 / 6.;
350		static const double RA4 = 17 / 360.;
351		static const double RA6 = 67 / 3024.;
352		static const double RV4 = 5 / 72.;
353		static const double RV6 = 55 / 1296.;
354
355		/***************************************************************************************/
356	0	static int ellps_spherification(PJ *P) {
357		/***************************************************************************************/
358	0	const char *keys[] = {"R_A", "R_V", "R_a", "R_g",
359	0	"R_h", "R_lat_a", "R_lat_g", "R_C"};
360	0	size_t len, i;
361	0	paralist *par = nullptr;
362
363	0	double t;
364	0	char v, endp;
365
366	0	len = sizeof(keys) / sizeof(char *);
367
368		/* Check which spherification key is specified */
369	0	for (i = 0; i < len; i++) {
370	0	par = pj_get_param(P->params, keys[i]);
371	0	if (par)
372	0	break;
373	0	}
374
375		/* No spherification specified? Then we're done */
376	0	if (i == len)
377	0	return 0;
378
379		/* Store definition */
380	0	P->def_spherification = pj_strdup(par->param);
381	0	par->used = 1;
382
383	0	switch (i) {
384
385		/* R_A - a sphere with same area as ellipsoid */
386	0	case 0:
387	0	P->a = 1. - P->es (SIXTH + P->es * (RA4 + P->es * RA6));
388	0	break;
389
390		/* R_V - a sphere with same volume as ellipsoid */
391	0	case 1:
392	0	P->a = 1. - P->es (SIXTH + P->es * (RV4 + P->es * RV6));
393	0	break;
394
395		/* R_a - a sphere with R = the arithmetic mean of the ellipsoid */
396	0	case 2:
397	0	P->a = (P->a + P->b) / 2;
398	0	break;
399
400		/* R_g - a sphere with R = the geometric mean of the ellipsoid */
401	0	case 3:
402	0	P->a = sqrt(P->a * P->b);
403	0	break;
404
405		/* R_h - a sphere with R = the harmonic mean of the ellipsoid */
406	0	case 4:
407	0	if (P->a + P->b == 0)
408	0	return proj_errno_set(
409	0	P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
410	0	P->a = (2 * P->a * P->b) / (P->a + P->b);
411	0	break;
412
413		/* R_lat_a - a sphere with R = the arithmetic mean of the ellipsoid at given
414		* latitude */
415	0	case 5:
416		/* R_lat_g - a sphere with R = the geometric mean of the ellipsoid at given
417		* latitude */
418	0	case 6:
419	0	v = pj_param_value(par);
420	0	t = proj_dmstor(v, &endp);
421	0	if (fabs(t) > M_HALFPI) {
422	0	proj_log_error(
423	0	P, _("Invalid value for lat_g. \|lat_g\| should be <= 90°"));
424	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
425	0	}
426	0	t = sin(t);
427	0	t = 1 - P->es * t * t;
428	0	if (t == 0.) {
429	0	proj_log_error(P, _("Invalid eccentricity"));
430	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
431	0	}
432	0	if (i == 5) /* arithmetic */
433	0	P->a = (1. - P->es + t) / (2 t * sqrt(t));
434	0	else /* geometric */
435	0	P->a *= sqrt(1 - P->es) / t;
436	0	break;
437
438		/* R_C - a sphere with R = radius of conformal sphere, taken at a
439		* latitude that is phi0 (note: at least for mercator. for other
440		* projection methods, this could be phi1)
441		* Formula from IOGP Publication 373-7-2 – Geomatics Guidance Note number 7,
442		* part 2 1.1 Ellipsoid parameters
443		*/
444	0	case 7:
445	0	t = sin(P->phi0);
446	0	t = 1 - P->es * t * t;
447	0	if (t == 0.) {
448	0	proj_log_error(P, _("Invalid eccentricity"));
449	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
450	0	}
451	0	P->a *= sqrt(1 - P->es) / t;
452	0	break;
453	0	}
454
455	0	if (P->a <= 0.) {
456	0	proj_log_error(P, _("Invalid or missing major axis"));
457	0	return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
458	0	}
459
460		/* Clean up the ellipsoidal parameters to reflect the sphere */
461	0	P->es = P->e = P->f = 0;
462	0	P->rf = HUGE_VAL;
463	0	P->b = P->a;
464	0	pj_calc_ellipsoid_params(P, P->a, 0);
465
466	0	return 0;
467	0	}
468
469		/* locate parameter in list */
470	0	static paralist pj_get_param(paralist list, const char *key) {
471	0	size_t l = strlen(key);
472	0	while (list && !(0 == strncmp(list->param, key, l) &&
473	0	(0 == list->param[l] \|\| list->param[l] == '=')))
474	0	list = list->next;
475	0	return list;
476	0	}
477
478	0	static char pj_param_value(paralist list) {
479	0	char key, value;
480	0	if (nullptr == list)
481	0	return nullptr;
482
483	0	key = list->param;
484	0	value = strchr(key, '=');
485
486		/* a flag (i.e. a key without value) has its own name (key) as value */
487	0	return value ? value + 1 : key;
488	0	}
489
490	0	static const PJ_ELLPS pj_find_ellps(const char name) {
491	0	int i;
492	0	const char *s;
493	0	const PJ_ELLPS *ellps;
494
495	0	if (nullptr == name)
496	0	return nullptr;
497
498	0	ellps = proj_list_ellps();
499
500		/* Search through internal ellipsoid list for name */
501	0	for (i = 0; (s = ellps[i].id) && strcmp(name, s); ++i)
502	0	;
503	0	if (nullptr == s)
504	0	return nullptr;
505	0	return ellps + i;
506	0	}
507
508		/**************************************************************************************/
509	0	void pj_inherit_ellipsoid_def(const PJ src, PJ dst) {
510		/***************************************************************************************
511		Brute force copy the ellipsoidal parameters from src to dst. This code
512		was written before the actual ellipsoid setup parameters were kept available
513		in the PJ->def_xxx elements.
514		***************************************************************************************/
515
516		/* The linear parameters */
517	0	dst->a = src->a;
518	0	dst->b = src->b;
519	0	dst->ra = src->ra;
520	0	dst->rb = src->rb;
521
522		/* The eccentricities */
523	0	dst->alpha = src->alpha;
524	0	dst->e = src->e;
525	0	dst->es = src->es;
526	0	dst->e2 = src->e2;
527	0	dst->e2s = src->e2s;
528	0	dst->e3 = src->e3;
529	0	dst->e3s = src->e3s;
530	0	dst->one_es = src->one_es;
531	0	dst->rone_es = src->rone_es;
532
533		/* The flattenings */
534	0	dst->f = src->f;
535	0	dst->f2 = src->f2;
536	0	dst->n = src->n;
537	0	dst->rf = src->rf;
538	0	dst->rf2 = src->rf2;
539	0	dst->rn = src->rn;
540
541		/* This one's for GRS80 */
542	0	dst->J = src->J;
543
544		/* es and a before any +proj related adjustment */
545	0	dst->es_orig = src->es_orig;
546	0	dst->a_orig = src->a_orig;
547	0	}
548
549		/***************************************************************************************/
550	0	int pj_calc_ellipsoid_params(PJ *P, double a, double es) {
551		/****************************************************************************************
552		Calculate a large number of ancillary ellipsoidal parameters, in
553		addition to the two traditional PROJ defining parameters: Semimajor axis, a,
554		and the eccentricity squared, es.
555
556		Most of these parameters are fairly cheap to compute in comparison to
557		the overall effort involved in initializing a PJ object. They may, however,
558		take a substantial part of the time taken in computing an individual point
559		transformation.
560
561		So by providing them up front, we can amortize the (already modest) cost
562		over all transformations carried out over the entire lifetime of a PJ
563		object, rather than incur that cost for every single transformation.
564
565		Most of the parameter calculations here are based on the "angular
566		eccentricity", i.e. the angle, measured from the semiminor axis, of a line
567		going from the north pole to one of the foci of the ellipsoid - or in other
568		words: The arc sine of the eccentricity.
569
570		The formulae used are mostly taken from:
571
572		Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
573		Columbus, Ohio: Dept. of Geodetic Science
574		and Surveying, Ohio State University.
575
576		****************************************************************************************/
577
578	0	P->a = a;
579	0	P->es = es;
580
581		/* Compute some ancillary ellipsoidal parameters */
582	0	if (P->e == 0)
583	0	P->e = sqrt(P->es); /* eccentricity */
584	0	P->alpha = asin(P->e); /* angular eccentricity */
585
586		/* second eccentricity */
587	0	P->e2 = tan(P->alpha);
588	0	P->e2s = P->e2 * P->e2;
589
590		/* third eccentricity */
591	0	P->e3 = (0 != P->alpha)
592	0	? sin(P->alpha) / sqrt(2 - sin(P->alpha) * sin(P->alpha))
593	0	: 0;
594	0	P->e3s = P->e3 * P->e3;
595
596		/* flattening */
597	0	if (0 == P->f)
598	0	P->f = 1 - cos(P->alpha); /* = 1 - sqrt (1 - PIN->es); */
599	0	if (!(P->f >= 0.0 && P->f < 1.0)) {
600	0	proj_log_error(P, _("Invalid eccentricity"));
601	0	proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
602	0	return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
603	0	}
604	0	P->rf = P->f != 0.0 ? 1.0 / P->f : HUGE_VAL;
605
606		/* second flattening */
607	0	P->f2 = (cos(P->alpha) != 0) ? 1 / cos(P->alpha) - 1 : 0;
608	0	P->rf2 = P->f2 != 0.0 ? 1 / P->f2 : HUGE_VAL;
609
610		/* third flattening */
611	0	P->n = pow(tan(P->alpha / 2), 2);
612	0	P->rn = P->n != 0.0 ? 1 / P->n : HUGE_VAL;
613
614		/* ...and a few more */
615	0	if (0 == P->b)
616	0	P->b = (1 - P->f) * P->a;
617	0	P->rb = 1. / P->b;
618	0	P->ra = 1. / P->a;
619
620	0	P->one_es = 1. - P->es;
621	0	if (P->one_es == 0.) {
622	0	proj_log_error(P, _("Invalid eccentricity"));
623	0	proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
624	0	return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
625	0	}
626
627	0	P->rone_es = 1. / P->one_es;
628
629	0	return 0;
630	0	}
631
632		/**************************************************************************************/
633	0	int pj_ell_set(PJ_CONTEXT ctx, paralist pl, double a, double es) {
634		/***************************************************************************************
635		Initialize ellipsoidal parameters by emulating the original ellipsoid
636		setup function by Gerald Evenden, through a call to pj_ellipsoid
637		***************************************************************************************/
638	0	PJ B;
639	0	int ret;
640
641	0	B.ctx = ctx;
642	0	B.params = pl;
643
644	0	ret = pj_ellipsoid(&B);
645
646	0	free(B.def_size);
647	0	free(B.def_shape);
648	0	free(B.def_spherification);
649	0	free(B.def_ellps);
650
651	0	if (ret)
652	0	return ret;
653
654	0	*a = B.a;
655	0	*es = B.es;
656	0	return 0;
657	0	}