Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
  • Loading branch information
Fil committed Sep 21, 2018
1 parent 6b0567b commit fce77eb
Showing 1 changed file with 24 additions and 21 deletions.
45 changes: 24 additions & 21 deletions src/PJ_bertin1953.c
Original file line number Diff line number Diff line change
Expand Up @@ -31,44 +31,45 @@ static XY s_forward (LP lp, PJ *P) {
XY xy = {0.0,0.0};
struct pj_opaque *Q = P->opaque;

// Projection constants
/* Projection constants */
double fu = 1.4, k = 12., w = 1.68;

// Aliases
/* Aliases */
double lambda = lp.lam, phi = lp.phi;

// Variable
/* Variable */
double d;

// Apply rotation
/* Apply rotation */
lambda += Q->deltaLambda;

double cosphi = cos(phi),
x = cos(lambda) * cosphi,
y = sin(lambda) * cosphi,
z = sin(phi),
z0 = z * Q->cosDeltaPhi + x * Q->sinDeltaPhi;
double cosphi, x, y, z, z0;
cosphi = cos(phi);
x = cos(lambda) * cosphi;
y = sin(lambda) * cosphi;
z = sin(phi);
z0 = z * Q->cosDeltaPhi + x * Q->sinDeltaPhi;
lambda = atan2(y * Q->cosDeltaGamma - z0 * Q->sinDeltaGamma,
x * Q->cosDeltaPhi - z * Q->sinDeltaPhi);
z0 = z0 * Q->cosDeltaGamma + y * Q->sinDeltaGamma;
phi = asin(z0);

lambda = adjlon(lambda);

// Adjust pre-projection
/* Adjust pre-projection */
if (lambda + phi < -fu) {
double u = (lambda - phi + 1.6) * (lambda + phi + fu) / 8.;
lambda += u;
phi -= 0.8 * u * sin(phi + M_PI / 2.);
d = (lambda - phi + 1.6) * (lambda + phi + fu) / 8.;
lambda += d;
phi -= 0.8 * d * sin(phi + M_PI / 2.);
}

// Project with Hammer (1.68,2)
/* Project with Hammer (1.68,2) */
cosphi = cos(phi);
d = sqrt(2./(1. + cosphi * cos(lambda / 2.)));
xy.x = w * d * cosphi * sin(lambda / 2.);
xy.y = d * sin(phi);

// Adjust post-projection
/* Adjust post-projection */
d = (1. - cos(lambda * phi)) / k;
if (xy.y < 0.) {
xy.x *= 1. + d;
Expand All @@ -87,12 +88,14 @@ PJ *PROJECTION(bertin1953) {
return pj_default_destructor (P, ENOMEM);
P->opaque = Q;

// force +lon_0 = -16.5
double deltaLambda = P->lam0 -16.5 / 180. * M_PI;
double deltaLambda, deltaPhi, deltaGamma;

/* force +lon_0 = -16.5 */
deltaLambda = P->lam0 -16.5 / 180. * M_PI;

// force +lat_0=-42
double deltaPhi = -42. / 180. * M_PI,
deltaGamma = 0. / 180. * M_PI;
/* force +lat_0=-42 */
deltaPhi = -42. / 180. * M_PI;
deltaGamma = 0. / 180. * M_PI;

Q->deltaLambda = deltaLambda;
Q->cosDeltaPhi = cos(deltaPhi);
Expand All @@ -102,7 +105,7 @@ PJ *PROJECTION(bertin1953) {

P->es = 0.;
P->fwd = s_forward;
// P->inv = s_inverse;
/* P->inv = s_inverse; */

return P;
}

0 comments on commit fce77eb

Please sign in to comment.