path: root/src/PJ_aea.c

/******************************************************************************
 * Project:  PROJ.4
 * Purpose:  Implementation of the aea (Albers Equal Area) projection.
 *           and the leac (Lambert Equal Area Conic) projection
 * Author:   Gerald Evenden (1995)
 *           Thomas Knudsen (2016) - revise/add regression tests
 *
 ******************************************************************************
 * Copyright (c) 1995, Gerald Evenden
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#define PJ_LIB__
#include <projects.h>

# define EPS10	1.e-10
# define TOL7	1.e-7

PROJ_HEAD(aea, "Albers Equal Area") "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
PROJ_HEAD(leac, "Lambert Equal Area Conic") "\n\tConic, Sph&Ell\n\tlat_1= south";


/* determine latitude angle phi-1 */
# define N_ITER 15
# define EPSILON 1.0e-7
# define TOL 1.0e-10
static double phi1_(double qs, double Te, double Tone_es) {
	int i;
	double Phi, sinpi, cospi, con, com, dphi;

	Phi = asin (.5 * qs);
	if (Te < EPSILON)
		return( Phi );
	i = N_ITER;
	do {
		sinpi = sin (Phi);
		cospi = cos (Phi);
		con = Te * sinpi;
		com = 1. - con * con;
		dphi = .5 * com * com / cospi * (qs / Tone_es -
		   sinpi / com + .5 / Te * log ((1. - con) /
		   (1. + con)));
		Phi += dphi;
	} while (fabs(dphi) > TOL && --i);
	return( i ? Phi : HUGE_VAL );
}


struct pj_opaque {
	double	ec;
	double	n;
	double	c;
	double	dd;
	double	n2;
	double	rho0;
	double	rho;
	double	phi1;
	double	phi2;
	double	*en;
	int		ellips;
};



static XY e_forward (LP lp, PJ *P) {   /* Ellipsoid/spheroid, forward */
    XY xy = {0.0,0.0};
    struct pj_opaque *Q = P->opaque;
	if ((Q->rho = Q->c - (Q->ellips ? Q->n * pj_qsfn(sin(lp.phi),
		P->e, P->one_es) : Q->n2 * sin(lp.phi))) < 0.) F_ERROR
	Q->rho = Q->dd * sqrt(Q->rho);
	xy.x = Q->rho * sin( lp.lam *= Q->n );
	xy.y = Q->rho0 - Q->rho * cos(lp.lam);
	return xy;
}


static LP e_inverse (XY xy, PJ *P) {   /* Ellipsoid/spheroid, inverse */
    LP lp = {0.0,0.0};
    struct pj_opaque *Q = P->opaque;
	if( (Q->rho = hypot(xy.x, xy.y = Q->rho0 - xy.y)) != 0.0 ) {
		if (Q->n < 0.) {
			Q->rho = -Q->rho;
			xy.x = -xy.x;
			xy.y = -xy.y;
		}
		lp.phi =  Q->rho / Q->dd;
		if (Q->ellips) {
			lp.phi = (Q->c - lp.phi * lp.phi) / Q->n;
			if (fabs(Q->ec - fabs(lp.phi)) > TOL7) {
				if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
					I_ERROR
			} else
				lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
		} else if (fabs(lp.phi = (Q->c - lp.phi * lp.phi) / Q->n2) <= 1.)
			lp.phi = asin(lp.phi);
		else
			lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
		lp.lam = atan2(xy.x, xy.y) / Q->n;
	} else {
		lp.lam = 0.;
		lp.phi = Q->n > 0. ? HALFPI : - HALFPI;
	}
	return lp;
}


static void *freeup_new (PJ *P) {                        /* Destructor */
    if (0==P)
        return 0;

    if (0==P->opaque)
        return pj_dealloc (P);

    pj_dealloc (P->opaque->en);
    pj_dealloc (P->opaque);
    return pj_dealloc(P);
}


static void freeup (PJ *P) {
    freeup_new (P);
    return;
}


static PJ *setup(PJ *P) {
	double cosphi, sinphi;
	int secant;
    struct pj_opaque *Q = P->opaque;

    P->inv = e_inverse;
    P->fwd = e_forward;

	if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21);
	Q->n = sinphi = sin(Q->phi1);
	cosphi = cos(Q->phi1);
	secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
	if( (Q->ellips = (P->es > 0.))) {
		double ml1, m1;

		if (!(Q->en = pj_enfn(P->es))) E_ERROR_0;
		m1 = pj_msfn(sinphi, cosphi, P->es);
		ml1 = pj_qsfn(sinphi, P->e, P->one_es);
		if (secant) { /* secant cone */
			double ml2, m2;

			sinphi = sin(Q->phi2);
			cosphi = cos(Q->phi2);
			m2 = pj_msfn(sinphi, cosphi, P->es);
			ml2 = pj_qsfn(sinphi, P->e, P->one_es);
			Q->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
		}
		Q->ec = 1. - .5 * P->one_es * log((1. - P->e) /
			(1. + P->e)) / P->e;
		Q->c = m1 * m1 + Q->n * ml1;
		Q->dd = 1. / Q->n;
		Q->rho0 = Q->dd * sqrt(Q->c - Q->n * pj_qsfn(sin(P->phi0),
			P->e, P->one_es));
	} else {
		if (secant) Q->n = .5 * (Q->n + sin(Q->phi2));
		Q->n2 = Q->n + Q->n;
		Q->c = cosphi * cosphi + Q->n2 * sinphi;
		Q->dd = 1. / Q->n;
		Q->rho0 = Q->dd * sqrt(Q->c - Q->n2 * sin(P->phi0));
	}

	return P;
}


PJ *PROJECTION(aea) {
    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
    if (0==Q)
        return freeup_new (P);
    P->opaque = Q;

    P->pfree = freeup;
    P->descr = des_aea;

	Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
	Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
    setup(P);
    return P;
}


PJ *PROJECTION(leac) {
    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
    if (0==Q)
        return freeup_new (P);
    P->opaque = Q;

    P->pfree = freeup;
    P->descr = des_leac;
	Q->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
	Q->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? - HALFPI: HALFPI;
    setup (P);
    return P;
}


#ifdef PJ_OMIT_SELFTEST
int pj_aea_selftest (void) {return 10000;}
#else

int pj_aea_selftest (void) {
    double tolerance_lp = 1e-10;
    double tolerance_xy = 1e-7;

    char e_args[] = {"+proj=aea   +ellps=GRS80  +lat_1=0 +lat_2=2"};
    char s_args[] = {"+proj=aea   +a=6400000    +lat_1=0 +lat_2=2"};

    LP fwd_in[] = {
        { 2, 1},
        { 2,-1},
        {-2, 1},
        {-2,-1}
    };

    XY e_fwd_expect[] = {
        {222571.60875710563,  110653.32674302977},
        {222706.30650839131,  -110484.26714439997},
        {-222571.60875710563,  110653.32674302977},
        {-222706.30650839131,  -110484.26714439997},
    };

    XY s_fwd_expect[] = {
        {223334.08517088494,  111780.43188447191},
        {223470.15499168713,  -111610.33943099028},
        {-223334.08517088494,  111780.43188447191},
        {-223470.15499168713,  -111610.33943099028},
    };

    XY inv_in[] = {
        { 200, 100},
        { 200,-100},
        {-200, 100},
        {-200,-100}
    };

    LP e_inv_expect[] = {
        {0.0017966310597749514,  0.00090436885862202158},
        {0.0017966300767030448,  -0.00090437009538581453},
        {-0.0017966310597749514,  0.00090436885862202158},
        {-0.0017966300767030448,  -0.00090437009538581453},
    };

    LP s_inv_expect[] = {
        {0.0017904935979658752,  0.00089524594491375306},
        {0.0017904926216016812,  -0.00089524716502493225},
        {-0.0017904935979658752,  0.00089524594491375306},
        {-0.0017904926216016812,  -0.00089524716502493225},
    };
    return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
}


#endif



#ifdef PJ_OMIT_SELFTEST
int pj_leac_selftest (void) {return 10000;}
#else

int pj_leac_selftest (void) {
    double tolerance_lp = 1e-10;
    double tolerance_xy = 1e-7;

    char e_args[] = {"+proj=leac   +ellps=GRS80  +lat_1=0 +lat_2=2"};
    char s_args[] = {"+proj=leac   +a=6400000    +lat_1=0 +lat_2=2"};

    LP fwd_in[] = {
        { 2, 1},
        { 2,-1},
        {-2, 1},
        {-2,-1}
    };

    XY e_fwd_expect[] = {
        {220685.14054297868,  112983.50088939646},
        {224553.31227982609,  -108128.63674487274},
        {-220685.14054297868,  112983.50088939646},
        {-224553.31227982609,  -108128.63674487274},
    };

    XY s_fwd_expect[] = {
        {221432.86859285168,  114119.45452653214},
        {225331.72412711097,  -109245.82943505641},
        {-221432.86859285168,  114119.45452653214},
        {-225331.72412711097,  -109245.82943505641},
    };

    XY inv_in[] = {
        { 200, 100},
        { 200,-100},
        {-200, 100},
        {-200,-100}
    };

    LP e_inv_expect[] = {
        {0.0017966446840328458,  0.00090435171340223211},
        {0.0017966164523713021,  -0.00090438724081843625},
        {-0.0017966446840328458,  0.00090435171340223211},
        {-0.0017966164523713021,  -0.00090438724081843625},
    };

    LP s_inv_expect[] = {
        {0.0017905070979748127,  0.00089522906964877795},
        {0.001790479121519977,  -0.00089526404022281043},
        {-0.0017905070979748127,  0.00089522906964877795},
        {-0.001790479121519977,  -0.00089526404022281043},
    };

    return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
}











/***********************************************************************
                          SELFTEST STUBS
************************************************************************

Selftest stubs temporarily placed here.

To be removed as real selftest functions are added to the projection
source files

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

int pj_aeqd_selftest     (void)     {return 10000;}
int pj_alsk_selftest     (void)     {return 10000;}
int pj_eck3_selftest     (void)     {return 10000;}
int pj_eck4_selftest     (void)     {return 10000;}
int pj_eck5_selftest     (void)     {return 10000;}
int pj_eck6_selftest     (void)     {return 10000;}
int pj_eqc_selftest (void) {return 10000;}
int pj_eqdc_selftest (void) {return 10000;}
int pj_etmerc_selftest (void) {return 10000;}
int pj_fahey_selftest (void) {return 10000;}

int pj_fouc_s_selftest (void) {return 10000;}
int pj_gall_selftest (void) {return 10000;}
int pj_geos_selftest (void) {return 10000;}
int pj_geocent_selftest (void) {return 10000;}
int pj_gins8_selftest (void) {return 10000;}
int pj_gn_sinu_selftest (void) {return 10000;}
int pj_gnom_selftest (void) {return 10000;}
int pj_goode_selftest (void) {return 10000;}
int pj_gs48_selftest (void) {return 10000;}
int pj_gs50_selftest (void) {return 10000;}
int pj_hammer_selftest (void) {return 10000;}
int pj_hatano_selftest (void) {return 10000;}
int pj_healpix_selftest (void) {return 10000;}
int pj_rhealpix_selftest (void) {return 10000;}
int pj_igh_selftest (void) {return 10000;}
int pj_imw_p_selftest (void) {return 10000;}
int pj_isea_selftest (void) {return 10000;}

int pj_kav7_selftest (void) {return 10000;}
int pj_krovak_selftest (void) {return 10000;}
int pj_labrd_selftest (void) {return 10000;}
int pj_laea_selftest (void) {return 10000;}
int pj_lagrng_selftest (void) {return 10000;}
int pj_larr_selftest (void) {return 10000;}
int pj_lask_selftest (void) {return 10000;}
int pj_latlon_selftest (void) {return 10000;}
int pj_latlong_selftest (void) {return 10000;}
int pj_lonlat_selftest (void) {return 10000;}
int pj_longlat_selftest (void) {return 10000;}
int pj_lcc_selftest (void) {return 10000;}
int pj_lcca_selftest (void) {return 10000;}

int pj_lee_os_selftest (void) {return 10000;}
int pj_loxim_selftest (void) {return 10000;}
int pj_lsat_selftest (void) {return 10000;}

int pj_mbt_fps_selftest (void) {return 10000;}
int pj_mbtfpp_selftest (void) {return 10000;}
int pj_mbtfpq_selftest (void) {return 10000;}
int pj_mbtfps_selftest (void) {return 10000;}
int pj_merc_selftest (void) {return 10000;}
int pj_mil_os_selftest (void) {return 10000;}
int pj_mill_selftest (void) {return 10000;}
int pj_misrsom_selftest (void) {return 10000;}
int pj_moll_selftest (void) {return 10000;}
int pj_natearth_selftest (void) {return 10000;}
int pj_natearth2_selftest (void) {return 10000;}
int pj_nell_selftest (void) {return 10000;}
int pj_nell_h_selftest (void) {return 10000;}
int pj_nicol_selftest (void) {return 10000;}
int pj_nsper_selftest (void) {return 10000;}
int pj_nzmg_selftest (void) {return 10000;}
int pj_ob_tran_selftest (void) {return 10000;}
int pj_ocea_selftest (void) {return 10000;}
int pj_oea_selftest (void) {return 10000;}
int pj_omerc_selftest (void) {return 10000;}
int pj_ortho_selftest (void) {return 10000;}
int pj_patterson_selftest (void) {return 10000;}
int pj_poly_selftest (void) {return 10000;}
int pj_putp1_selftest (void) {return 10000;}
int pj_putp2_selftest (void) {return 10000;}
int pj_putp3_selftest (void) {return 10000;}
int pj_putp3p_selftest (void) {return 10000;}
int pj_putp4p_selftest (void) {return 10000;}
int pj_putp5_selftest (void) {return 10000;}
int pj_putp5p_selftest (void) {return 10000;}
int pj_putp6_selftest (void) {return 10000;}
int pj_putp6p_selftest (void) {return 10000;}

int pj_qsc_selftest (void) {return 10000;}
int pj_robin_selftest (void) {return 10000;}
int pj_rouss_selftest (void) {return 10000;}
int pj_rpoly_selftest (void) {return 10000;}
int pj_sch_selftest (void) {return 10000;}
int pj_sinu_selftest (void) {return 10000;}

int pj_gstmerc_selftest (void) {return 10000;}
int pj_tpers_selftest (void) {return 10000;}

int pj_utm_selftest (void) {return 10000;}
int pj_vandg_selftest (void) {return 10000;}
int pj_vandg2_selftest (void) {return 10000;}
int pj_vandg3_selftest (void) {return 10000;}
int pj_vandg4_selftest (void) {return 10000;}
int pj_wag4_selftest (void) {return 10000;}
int pj_wag5_selftest (void) {return 10000;}
int pj_wag6_selftest (void) {return 10000;}
int pj_weren_selftest (void) {return 10000;}
#endif