Skip to content

Commit

Permalink
Merge 6dbcf87 into bdb02c7
Browse files Browse the repository at this point in the history
  • Loading branch information
@kbevers
kbevers committed May 24, 2017
2 parents bdb02c7 + 6dbcf87 commit 007b4bf
Show file tree
Hide file tree
Showing 8 changed files with 199 additions and 165 deletions.
204 changes: 103 additions & 101 deletions src/PJ_aea.c
View file
Expand Up @@ -30,8 +30,8 @@
#define PJ_LIB__
#include <projects.h>

# define EPS10 1.e-10
# define TOL7 1.e-7
# 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";
Expand All @@ -42,82 +42,82 @@ PROJ_HEAD(leac, "Lambert Equal Area Conic") "\n\tConic, Sph&Ell\n\tlat_1= south"
# 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 );
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;
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;
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. ? -M_HALFPI : M_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. ? -M_HALFPI : M_HALFPI;
lp.lam = atan2(xy.x, xy.y) / Q->n;
} else {
lp.lam = 0.;
lp.phi = Q->n > 0. ? M_HALFPI : - M_HALFPI;
}
return lp;
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. ? -M_HALFPI : M_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. ? -M_HALFPI : M_HALFPI;
lp.lam = atan2(xy.x, xy.y) / Q->n;
} else {
lp.lam = 0.;
lp.phi = Q->n > 0. ? M_HALFPI : - M_HALFPI;
}
return lp;
}


Expand All @@ -141,47 +141,49 @@ static void freeup (PJ *P) {


static PJ *setup(PJ *P) {
double cosphi, sinphi;
int secant;
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;
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);
if (ml2 == ml1)
return NULL;
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;
}


Expand All @@ -191,8 +193,8 @@ PJ *PROJECTION(aea) {
return freeup_new (P);
P->opaque = Q;

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

Expand All @@ -203,8 +205,8 @@ PJ *PROJECTION(leac) {
return freeup_new (P);
P->opaque = Q;

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

Expand Down
7 changes: 6 additions & 1 deletion src/PJ_eck3.c
View file
Expand Up @@ -24,9 +24,14 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
LP lp = {0.0,0.0};
struct pj_opaque *Q = P->opaque;
double denominator;

lp.phi = xy.y / Q->C_y;
lp.lam = xy.x / (Q->C_x * (Q->A + asqrt(1. - Q->B * lp.phi * lp.phi)));
denominator = (Q->C_x * (Q->A + asqrt(1. - Q->B * lp.phi * lp.phi)));
if ( denominator == 0.0)
lp.lam = HUGE_VAL;
else
lp.lam = xy.x / denominator;
return lp;
}

Expand Down
11 changes: 8 additions & 3 deletions src/PJ_stere.c
View file
Expand Up @@ -34,7 +34,7 @@ static double ssfn_ (double phit, double sinphi, double eccen) {
static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
XY xy = {0.0,0.0};
struct pj_opaque *Q = P->opaque;
double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A, sinphi;
double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A = 0.0, sinphi;

coslam = cos (lp.lam);
sinlam = sin (lp.lam);
Expand All @@ -52,8 +52,13 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
goto xmul; /* but why not just xy.x = A * cosX; break; ? */

case EQUIT:
A = Q->akm1 / (1. + cosX * coslam);
xy.y = A * sinX;
/* avoid zero division */
if (1. + cosX * coslam == 0.0) {
xy.y = HUGE_VAL;
} else {
A = Q->akm1 / (1. + cosX * coslam);
xy.y = A * sinX;
}
xmul:
xy.x = A * cosX;
break;
Expand Down
4 changes: 4 additions & 0 deletions src/pj_ell_set.c
View file
Expand Up @@ -74,6 +74,10 @@ int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
*a = sqrt(*a * b);
*es = 0.;
} else if (pj_param(ctx,pl, "bR_h").i) { /* sphere--harmonic mean */
if ( (*a + b) == 0.0) {
pj_ctx_set_errno(ctx, -20);
goto bomb;
}
*a = 2. * *a * b / (*a + b);
*es = 0.;
} else if ((i = pj_param(ctx,pl, "tR_lat_a").i) || /* sphere--arith. */
Expand Down

0 comments on commit 007b4bf

Please sign in to comment.