Skip to content

Commit

Permalink
Merge 4a563fe into 5f03cc7
Browse files Browse the repository at this point in the history
  • Loading branch information
QuLogic committed May 24, 2017
2 parents 5f03cc7 + 4a563fe commit 2a864e6
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 25 deletions.
35 changes: 18 additions & 17 deletions src/pj_gauss.c
Expand Up @@ -34,7 +34,6 @@ struct GAUSS {
double e;
double ratexp;
};
#define EN ((struct GAUSS *)en)
#define DEL_TOL 1e-14

static double srat(double esinp, double exp) {
Expand All @@ -48,43 +47,45 @@ void *pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)
return (NULL);
es = e * e;
EN->e = e;
en->e = e;
sphi = sin(phi0);
cphi = cos(phi0); cphi *= cphi;
*rc = sqrt(1. - es) / (1. - es * sphi * sphi);
EN->C = sqrt(1. + es * cphi * cphi / (1. - es));
en->C = sqrt(1. + es * cphi * cphi / (1. - es));
if (en->C == 0.0) {
free(en);
return NULL;
}
*chi = asin(sphi / EN->C);
EN->ratexp = 0.5 * EN->C * e;
EN->K = tan(.5 * *chi + M_FORTPI) / (
pow(tan(.5 * phi0 + M_FORTPI), EN->C) *
srat(EN->e * sphi, EN->ratexp) );
*chi = asin(sphi / en->C);
en->ratexp = 0.5 * en->C * e;
en->K = tan(.5 * *chi + M_FORTPI) / (
pow(tan(.5 * phi0 + M_FORTPI), en->C) *
srat(en->e * sphi, en->ratexp) );
return ((void *)en);
}

LP pj_gauss(projCtx ctx, LP elp, const void *en) {
LP pj_gauss(projCtx ctx, LP elp, const void *data) {
struct GAUSS *en = (struct GAUSS *)data;
LP slp;
(void) ctx;

slp.phi = 2. * atan( EN->K *
pow(tan(.5 * elp.phi + M_FORTPI), EN->C) *
srat(EN->e * sin(elp.phi), EN->ratexp) ) - M_HALFPI;
slp.lam = EN->C * (elp.lam);
slp.phi = 2. * atan( en->K *
pow(tan(.5 * elp.phi + M_FORTPI), en->C) *
srat(en->e * sin(elp.phi), en->ratexp) ) - M_HALFPI;
slp.lam = en->C * (elp.lam);
return(slp);
}

LP pj_inv_gauss(projCtx ctx, LP slp, const void *en) {
LP pj_inv_gauss(projCtx ctx, LP slp, const void *data) {
struct GAUSS *en = (struct GAUSS *)data;
LP elp;
double num;
int i;

elp.lam = slp.lam / EN->C;
num = pow(tan(.5 * slp.phi + M_FORTPI)/EN->K, 1./EN->C);
elp.lam = slp.lam / en->C;
num = pow(tan(.5 * slp.phi + M_FORTPI)/en->K, 1./en->C);
for (i = MAX_ITER; i; --i) {
elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))
elp.phi = 2. * atan(num * srat(en->e * sin(slp.phi), -.5 * en->e))
- M_HALFPI;
if (fabs(elp.phi - slp.phi) < DEL_TOL) break;
slp.phi = elp.phi;
Expand Down
17 changes: 9 additions & 8 deletions src/proj_mdist.c
Expand Up @@ -38,7 +38,6 @@ struct MDIST {
double E;
double b[1];
};
#define B ((struct MDIST *)b)
void *
proj_mdist_ini(double es) {
double numf, numfi, twon1, denf, denfi, ens, T, twon;
Expand Down Expand Up @@ -88,28 +87,30 @@ proj_mdist_ini(double es) {
return (b);
}
double
proj_mdist(double phi, double sphi, double cphi, const void *b) {
proj_mdist(double phi, double sphi, double cphi, const void *data) {
struct MDIST *b = (struct MDIST *)data;
double sc, sum, sphi2, D;
int i;

sc = sphi * cphi;
sphi2 = sphi * sphi;
D = phi * B->E - B->es * sc / sqrt(1. - B->es * sphi2);
sum = B->b[i = B->nb];
while (i) sum = B->b[--i] + sphi2 * sum;
D = phi * b->E - b->es * sc / sqrt(1. - b->es * sphi2);
sum = b->b[i = b->nb];
while (i) sum = b->b[--i] + sphi2 * sum;
return(D + sc * sum);
}
double
proj_inv_mdist(projCtx ctx, double dist, const void *b) {
proj_inv_mdist(projCtx ctx, double dist, const void *data) {
struct MDIST *b = (struct MDIST *)data;
double s, t, phi, k;
int i;

k = 1./(1.- B->es);
k = 1./(1.- b->es);
i = MAX_ITER;
phi = dist;
while ( i-- ) {
s = sin(phi);
t = 1. - B->es * s * s;
t = 1. - b->es * s * s;
phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
(t * sqrt(t)) * k;
if (fabs(t) < TOL) /* that is no change */
Expand Down

0 comments on commit 2a864e6

Please sign in to comment.