Permalink
Browse files

Merge f48b423 into 06b2f94

  • Loading branch information...
busstoptaktik committed Nov 14, 2017
2 parents 06b2f94 + f48b423 commit f698b518da866fe0bc3557048d111b2509a72d09
Showing with 76 additions and 140 deletions.
  1. +0 −14 src/PJ_eqdc.c
  2. +0 −16 src/PJ_lcc.c
  3. +69 −78 src/pj_factors.c
  4. +7 −13 src/proj.c
  5. +0 −6 src/proj.h
  6. +0 −3 src/proj_4D_api.c
  7. +0 −10 src/projects.h
View
@@ -54,19 +54,6 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
}
static void special(LP lp, PJ *P, struct FACTORS *fac) {
struct pj_opaque *Q = P->opaque;
double sinphi, cosphi;
sinphi = sin(lp.phi);
cosphi = cos(lp.phi);
fac->code |= IS_ANAL_HK;
fac->h = 1.;
fac->k = Q->n * (Q->c - (Q->ellips ? pj_mlfn(lp.phi, sinphi,
cosphi, Q->en) : lp.phi)) / pj_msfn(sinphi, cosphi, P->es);
}
static void *destructor (PJ *P, int errlev) { /* Destructor */
if (0==P)
return 0;
@@ -124,7 +111,6 @@ PJ *PROJECTION(eqdc) {
P->inv = e_inverse;
P->fwd = e_forward;
P->spc = special;
return P;
}
View
@@ -73,21 +73,6 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
return lp;
}
static void special(LP lp, PJ *P, struct FACTORS *fac) {
struct pj_opaque *Q = P->opaque;
double rho;
if (fabs(fabs(lp.phi) - M_HALFPI) < EPS10) {
if ((lp.phi * Q->n) <= 0.) return;
rho = 0.;
} else
rho = Q->c * (Q->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi),
P->e), Q->n) : pow(tan(M_FORTPI + .5 * lp.phi), -Q->n));
fac->code |= IS_ANAL_HK + IS_ANAL_CONV;
fac->k = fac->h = P->k0 * Q->n * rho /
pj_msfn(sin(lp.phi), cos(lp.phi), P->es);
fac->conv = Q->n * lp.lam;
}
PJ *PROJECTION(lcc) {
double cosphi, sinphi;
@@ -139,7 +124,6 @@ PJ *PROJECTION(lcc) {
P->inv = e_inverse;
P->fwd = e_forward;
P->spc = special;
return P;
}
View
@@ -1,6 +1,8 @@
/* projection scale factors */
#define PJ_LIB__
#include <projects.h>
#include <proj.h>
#include "projects.h"
#include <errno.h>
#ifndef DEFAULT_H
@@ -10,88 +12,77 @@
#define EPS 1.0e-12
int pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
struct DERIVS der;
double cosphi, t, n, r;
int err;
err = proj_errno_reset (P);
/* Check for latitude or longitude overange */
if ((fabs (lp.phi)-M_HALFPI) > EPS || fabs (lp.lam) > 10.) {
proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
return 1;
}
/* Set a reasonable step size for the numerical derivatives */
h = fabs (h);
if (h < EPS)
h = DEFAULT_H;
/* If input latitudes are geocentric, convert to geographic */
if (P->geoc)
lp.phi = atan(P->rone_es * tan(lp.phi));
/* If latitude + one step overshoots the pole, move it slightly inside, */
/* so the numerical derivative still exists */
if (fabs (lp.phi) > (M_HALFPI - h))
lp.phi = lp.phi < 0. ? -(M_HALFPI-h) : (M_HALFPI-h);
der.x_l = 0.0;
der.x_p = 0.0;
der.y_l = 0.0;
der.y_p = 0.0;
/* Longitudinal distance from central meridian */
lp.lam -= P->lam0;
if (!P->over)
lp.lam = adjlon(lp.lam);
/* check for forward and latitude or longitude overange */
if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
pj_ctx_set_errno( P->ctx, -14);
/* Derivatives */
if (pj_deriv (lp, h, P, &(fac->der))) {
proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
return 1;
} else { /* proceed */
errno = pj_errno = 0;
P->ctx->last_errno = 0;
if (h < EPS)
h = DEFAULT_H;
if (fabs(lp.phi) > (M_HALFPI - h))
/* adjust to value around pi/2 where derived still exists*/
lp.phi = lp.phi < 0. ? (-M_HALFPI+h) : (M_HALFPI-h);
else if (P->geoc)
lp.phi = atan(P->rone_es * tan(lp.phi));
lp.lam -= P->lam0; /* compute del lp.lam */
if (!P->over)
lp.lam = adjlon(lp.lam); /* adjust del longitude */
if (P->spc) /* get what projection analytic values */
P->spc(lp, P, fac);
if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
(IS_ANAL_XL_YL+IS_ANAL_XP_YP)) && pj_deriv(lp, h, P, &der))
return 1;
if (!(fac->code & IS_ANAL_XL_YL)) {
fac->der.x_l = der.x_l;
fac->der.y_l = der.y_l;
}
if (!(fac->code & IS_ANAL_XP_YP)) {
fac->der.x_p = der.x_p;
fac->der.y_p = der.y_p;
}
cosphi = cos(lp.phi);
if (!(fac->code & IS_ANAL_HK)) {
fac->h = hypot(fac->der.x_p, fac->der.y_p);
fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
if (P->es != 0.0) {
t = sin(lp.phi);
t = 1. - P->es * t * t;
n = sqrt(t);
fac->h *= t * n / P->one_es;
fac->k *= n;
r = t * t / P->one_es;
} else
r = 1.;
} else if (P->es != 0.0) {
r = sin(lp.phi);
r = 1. - P->es * r * r;
r = r * r / P->one_es;
} else
r = 1.;
/* convergence */
if (!(fac->code & IS_ANAL_CONV)) {
fac->conv = - atan2(fac->der.x_p, fac->der.y_p);
if (fac->code & IS_ANAL_XL_YL)
fac->code |= IS_ANAL_CONV;
}
/* areal scale factor */
fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) * r / cosphi;
/* meridian-parallel angle theta prime */
fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));
/* Tissot ellips axis */
t = fac->k * fac->k + fac->h * fac->h;
fac->a = sqrt(t + 2. * fac->s);
t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
fac->b = 0.5 * (fac->a - t);
fac->a = 0.5 * (fac->a + t);
/* omega */
fac->omega = 2. * aasin(P->ctx,(fac->a - fac->b)/(fac->a + fac->b));
}
/* Scale factors */
cosphi = cos (lp.phi);
fac->h = hypot (fac->der.x_p, fac->der.y_p);
fac->k = hypot (fac->der.x_l, fac->der.y_l) / cosphi;
if (P->es != 0.0) {
t = sin(lp.phi);
t = 1. - P->es * t * t;
n = sqrt(t);
fac->h *= t * n / P->one_es;
fac->k *= n;
r = t * t / P->one_es;
} else
r = 1.;
/* Convergence */
fac->conv = -atan2 (fac->der.x_p, fac->der.y_p);
/* Areal scale factor */
fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) * r / cosphi;
/* Meridian-parallel angle (theta prime) */
fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));
/* Tissot ellipse axis */
t = fac->k * fac->k + fac->h * fac->h;
fac->a = sqrt(t + 2. * fac->s);
t = t - 2. * fac->s;
t = t > 0? sqrt(t): 0;
fac->b = 0.5 * (fac->a - t);
fac->a = 0.5 * (fac->a + t);
/* Angular distortion */
fac->omega = 2. * aasin(P->ctx, (fac->a - fac->b) / (fac->a + fac->b) );
proj_errno_restore (P, err);
return 0;
}
View
@@ -287,21 +287,15 @@ static void vprocess(FILE *fid) {
(void)printf(oform, dat_xy.u); putchar('\n');
(void)fputs("Northing (y): ", stdout);
(void)printf(oform, dat_xy.v); putchar('\n');
(void)printf("Meridian scale (h)%c: %.8f ( %.4g %% error )\n",
facs.code & IS_ANAL_HK ? '*' : ' ', facs.h, (facs.h-1.)*100.);
(void)printf("Parallel scale (k)%c: %.8f ( %.4g %% error )\n",
facs.code & IS_ANAL_HK ? '*' : ' ', facs.k, (facs.k-1.)*100.);
(void)printf("Areal scale (s): %.8f ( %.4g %% error )\n",
facs.s, (facs.s-1.)*100.);
(void)printf("Angular distortion (w): %.3f\n", facs.omega *
RAD_TO_DEG);
(void)printf("Meridian/Parallel angle: %.5f\n",
facs.thetap * RAD_TO_DEG);
(void)printf("Convergence%c: ",facs.code & IS_ANAL_CONV ? '*' : ' ');
(void)printf("Meridian scale (h) : %.8f ( %.4g %% error )\n", facs.h, (facs.h-1.)*100.);
(void)printf("Parallel scale (k) : %.8f ( %.4g %% error )\n", facs.k, (facs.k-1.)*100.);
(void)printf("Areal scale (s): %.8f ( %.4g %% error )\n", facs.s, (facs.s-1.)*100.);
(void)printf("Angular distortion (w): %.3f\n", facs.omega * RAD_TO_DEG);
(void)printf("Meridian/Parallel angle: %.5f\n", facs.thetap * RAD_TO_DEG);
(void)printf("Convergence : ");
(void)fputs(rtodms(pline, facs.conv, 0, 0), stdout);
(void)printf(" [ %.8f ]\n", facs.conv * RAD_TO_DEG);
(void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n",
facs.a, facs.b);
(void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n", facs.a, facs.b);
}
}
View
@@ -289,12 +289,6 @@ union PJ_PAIR {
double v[2]; /* Yes - It's really just a vector! */
};
#define PJ_IS_ANAL_XL_YL 01 /* derivatives of lon analytic */
#define PJ_IS_ANAL_XP_YP 02 /* derivatives of lat analytic */
#define PJ_IS_ANAL_HK 04 /* h and k analytic */
#define PJ_IS_ANAL_CONV 010 /* convergence analytic */
struct PJ_INFO {
char release[64]; /* Release info. Version + date */
char version[64]; /* Full version number */
View
@@ -816,9 +816,6 @@ PJ_FACTORS proj_factors(PJ *P, const LP lp) {
******************************************************************************/
PJ_FACTORS factors;
/* pj_factors rely code being zero */
factors.code = 0;
if (pj_factors(lp, P, 0.0, &factors)) {
/* errno set in pj_factors */
memset(&factors, 0, sizeof(PJ_FACTORS));
View
@@ -257,9 +257,6 @@ struct PJconsts {
PJ_COORD (*fwd4d)(PJ_COORD, PJ *);
PJ_COORD (*inv4d)(PJ_COORD, PJ *);
void (*spc)(LP, PJ *, struct FACTORS *);
void *(*destructor)(PJ *, int);
@@ -460,15 +457,8 @@ struct FACTORS {
double conv; /* convergence */
double s; /* areal scale factor */
double a, b; /* max-min scale error */
int code; /* info as to analytics, see following */
};
#define IS_ANAL_XL_YL 01 /* derivatives of lon analytic */
#define IS_ANAL_XP_YP 02 /* derivatives of lat analytic */
#define IS_ANAL_HK 04 /* h and k analytic */
#define IS_ANAL_CONV 010 /* convergence analytic */
/* datum_type values */
#define PJD_UNKNOWN 0
#define PJD_3PARAM 1

0 comments on commit f698b51

Please sign in to comment.