mirror of
https://github.com/boostorg/geometry.git
synced 2026-02-02 21:02:13 +00:00
[projections] ONLY: fixed withspace
This commit is contained in:
@@ -72,17 +72,17 @@ namespace boost { namespace geometry { namespace projections
|
||||
double en[EN_SIZE];
|
||||
int ellips;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/* determine latitude angle phi-1 */
|
||||
inline 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 );
|
||||
|
||||
@@ -68,11 +68,11 @@ namespace boost { namespace geometry { namespace projections
|
||||
double G;
|
||||
int mode;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -92,7 +92,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, cosphi, sinphi, rho, s, H, H2, c, Az, t, ct, st, cA, sA;
|
||||
|
||||
|
||||
coslam = cos(lp_lon);
|
||||
cosphi = cos(lp_lat);
|
||||
sinphi = sin(lp_lat);
|
||||
@@ -133,7 +133,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double c, Az, cosAz, A, B, D, E, F, psi, t;
|
||||
|
||||
|
||||
if ((c = boost::math::hypot(xy_x, xy_y)) < EPS10) {
|
||||
lp_lat = this->m_par.phi0;
|
||||
lp_lon = 0.;
|
||||
@@ -183,7 +183,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double cosphi, sinphi, t;
|
||||
|
||||
|
||||
cosphi = cos(lp_lat);
|
||||
sinphi = sin(lp_lat);
|
||||
t = 1. / sqrt(1. - this->m_par.es * sinphi * sinphi);
|
||||
@@ -196,7 +196,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double x2, t;
|
||||
int i;
|
||||
|
||||
|
||||
x2 = 0.5 * xy_x * xy_x;
|
||||
lp_lat = this->m_par.phi0;
|
||||
for (i = 0; i < 3; ++i) {
|
||||
@@ -226,7 +226,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, cosphi, sinphi;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
@@ -239,7 +239,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
oblcon:
|
||||
if (fabs(fabs(xy_y) - 1.) < TOL)
|
||||
if (xy_y < 0.)
|
||||
throw proj_exception();
|
||||
throw proj_exception();
|
||||
else
|
||||
xy_x = xy_y = 0.;
|
||||
else {
|
||||
@@ -264,7 +264,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double cosc, c_rh, sinc;
|
||||
|
||||
|
||||
if ((c_rh = boost::math::hypot(xy_x, xy_y)) > PI) {
|
||||
if (c_rh - EPS10 > PI) throw proj_exception();;
|
||||
c_rh = PI;
|
||||
|
||||
@@ -61,10 +61,10 @@ namespace boost { namespace geometry { namespace projections
|
||||
int mode;
|
||||
int no_cut; /* do not cut at hemisphere limit */
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -84,7 +84,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double sinlam, coslam, cosphi, sinphi, t, s, Krho, cosz;
|
||||
|
||||
|
||||
sinlam = sin(lp_lon);
|
||||
coslam = cos(lp_lon);
|
||||
switch (this->m_proj_parm.mode) {
|
||||
|
||||
@@ -53,10 +53,10 @@ namespace boost { namespace geometry { namespace projections
|
||||
double cosphi1;
|
||||
int mode;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -76,7 +76,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double c, d;
|
||||
|
||||
|
||||
if((d = acos(cos(lp_lat) * cos(c = 0.5 * lp_lon)))) {/* basic Aitoff */
|
||||
xy_x = 2. * d * cos(lp_lat) * sin(c) * (xy_y = 1. / sin(d));
|
||||
xy_y *= d * sin(lp_lat);
|
||||
|
||||
@@ -66,7 +66,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double t, c1, c, x1, x12, y1, y12;
|
||||
|
||||
|
||||
t = tan(.5 * lp_lat);
|
||||
c1 = sqrt(1. - t * t);
|
||||
c = 1. + c1 * cos(lp_lon *= .5);
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double ax, f;
|
||||
|
||||
|
||||
xy_y = this->m_proj_parm.bacn ? HALFPI * sin(lp_lat) : lp_lat;
|
||||
if ((ax = fabs(lp_lon)) >= EPS) {
|
||||
if (this->m_proj_parm.ortl && ax >= HALFPI)
|
||||
|
||||
@@ -90,7 +90,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
|
||||
int tag;
|
||||
|
||||
|
||||
cphi = cos(lp_lat);
|
||||
sphi = sin(lp_lat);
|
||||
cdlam = cos(sdlam = lamB - lp_lon);
|
||||
@@ -140,8 +140,8 @@ namespace boost { namespace geometry { namespace projections
|
||||
xy_y += (tag ? -r : r) * cos(t);
|
||||
if (this->m_proj_parm.noskew) {
|
||||
t = xy_x;
|
||||
xy_x = -xy_x * cAzc - xy_y * sAzc;
|
||||
xy_y = -xy_y * cAzc + t * sAzc;
|
||||
xy_x = -xy_x * cAzc - xy_y * sAzc;
|
||||
xy_y = -xy_y * cAzc + t * sAzc;
|
||||
}
|
||||
}
|
||||
|
||||
@@ -149,11 +149,11 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
|
||||
int neg, i;
|
||||
|
||||
|
||||
if (this->m_proj_parm.noskew) {
|
||||
t = xy_x;
|
||||
xy_x = -xy_x * cAzc + xy_y * sAzc;
|
||||
xy_y = -xy_y * cAzc - t * sAzc;
|
||||
xy_x = -xy_x * cAzc + xy_y * sAzc;
|
||||
xy_y = -xy_y * cAzc - t * sAzc;
|
||||
}
|
||||
if( (neg = (xy_x < 0.)) ) {
|
||||
xy_y = rhoc - xy_y;
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double theta, th1, c;
|
||||
int i;
|
||||
|
||||
|
||||
theta = lp_lat;
|
||||
if (fabs(fabs(lp_lat) - HALFPI) < EPS)
|
||||
xy_x = 0.;
|
||||
|
||||
@@ -76,7 +76,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double rh, E, c;
|
||||
|
||||
|
||||
rh = this->m_proj_parm.am1 + this->m_proj_parm.m1 - pj_mlfn(lp_lat, E = sin(lp_lat), c = cos(lp_lat), this->m_proj_parm.en);
|
||||
E = c * lp_lon / (rh * sqrt(1. - this->m_par.es * E * E));
|
||||
xy_x = rh * sin(E);
|
||||
@@ -86,7 +86,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double s, rh;
|
||||
|
||||
|
||||
rh = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.am1 - xy_y);
|
||||
lp_lat = pj_inv_mlfn(this->m_proj_parm.am1 + this->m_proj_parm.m1 - rh, this->m_par.es, this->m_proj_parm.en);
|
||||
if ((s = fabs(lp_lat)) < HALFPI) {
|
||||
@@ -117,7 +117,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double E, rh;
|
||||
|
||||
|
||||
rh = this->m_proj_parm.cphi1 + this->m_proj_parm.phi1 - lp_lat;
|
||||
if (fabs(rh) > EPS10) {
|
||||
xy_x = rh * sin(E = lp_lon * cos(lp_lat) / rh);
|
||||
@@ -129,7 +129,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double rh;
|
||||
|
||||
|
||||
rh = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.cphi1 - xy_y);
|
||||
lp_lat = this->m_proj_parm.cphi1 + this->m_proj_parm.phi1 - rh;
|
||||
if (fabs(lp_lat) > HALFPI) throw proj_exception();;
|
||||
|
||||
@@ -103,7 +103,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double ph1;
|
||||
|
||||
|
||||
ph1 = pj_inv_mlfn(this->m_proj_parm.m0 + xy_y, this->m_par.es, this->m_proj_parm.en);
|
||||
this->m_proj_parm.tn = tan(ph1); this->m_proj_parm.t = this->m_proj_parm.tn * this->m_proj_parm.tn;
|
||||
this->m_proj_parm.n = sin(ph1);
|
||||
|
||||
@@ -108,7 +108,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
if ((t = fabs(xy_y *= this->m_par.k0)) - EPS <= 1.) {
|
||||
if (t >= 1.)
|
||||
lp_lat = xy_y < 0. ? -HALFPI : HALFPI;
|
||||
|
||||
@@ -71,7 +71,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
vect(double dphi, double c1, double s1, double c2, double s2, double dlam) {
|
||||
VECT v;
|
||||
double cdl, dp, dl;
|
||||
|
||||
|
||||
cdl = cos(dlam);
|
||||
if (fabs(dphi) > 1. || fabs(dlam) > 1.)
|
||||
v.r = aacos(s1 * s2 + c1 * c2 * cdl);
|
||||
@@ -111,7 +111,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
double sinphi, cosphi, a;
|
||||
VECT v[3];
|
||||
int i, j;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
for (i = 0; i < 3; ++i) { /* dist/azimiths from control */
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double p, V, s, c;
|
||||
int i;
|
||||
|
||||
|
||||
p = C_p * sin(lp_lat);
|
||||
V = lp_lat * lp_lat;
|
||||
lp_lat *= 0.895168 + V * ( 0.0218849 + V * 0.00826809 );
|
||||
@@ -97,7 +97,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double c;
|
||||
|
||||
|
||||
lp_lat = aasin(xy_y / C_y);
|
||||
lp_lon = xy_x / (C_x * (1. + (c = cos(lp_lat))));
|
||||
lp_lat = aasin((lp_lat + sin(lp_lat) * (c + 2.)) / C_p);
|
||||
|
||||
@@ -107,7 +107,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fac(Geographic lp, Factors &fac) const
|
||||
{
|
||||
double sinphi, cosphi;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
this->m_fac.code |= IS_ANAL_HK;
|
||||
@@ -128,7 +128,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
proj_parm.phi2 = pj_param(par.params, "rlat_2").f;
|
||||
if (fabs(proj_parm.phi1 + proj_parm.phi2) < EPS10) throw proj_exception(-21);
|
||||
pj_enfn(par.es, proj_parm.en);
|
||||
|
||||
|
||||
proj_parm.n = sinphi = sin(proj_parm.phi1);
|
||||
cosphi = cos(proj_parm.phi1);
|
||||
secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= EPS10;
|
||||
|
||||
@@ -72,7 +72,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
t = cos(lp_lat);
|
||||
xy_x = lp_lon * t / (this->m_proj_parm.n + this->m_proj_parm.n1 * t);
|
||||
xy_y = this->m_proj_parm.n * lp_lat + this->m_proj_parm.n1 * sin(lp_lat);
|
||||
@@ -82,7 +82,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double V;
|
||||
int i;
|
||||
|
||||
|
||||
if (this->m_proj_parm.n) {
|
||||
lp_lat = xy_y;
|
||||
for (i = MAX_ITER; i ; --i) {
|
||||
|
||||
@@ -47,10 +47,10 @@ namespace boost { namespace geometry { namespace projections
|
||||
#ifndef DOXYGEN_NO_DETAIL
|
||||
namespace detail { namespace geocent{
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -84,7 +84,6 @@ namespace boost { namespace geometry { namespace projections
|
||||
void setup_geocent(Parameters& par)
|
||||
{
|
||||
par.is_geocent = 1;
|
||||
|
||||
par.x0 = 0.0;
|
||||
par.y0 = 0.0;
|
||||
// par.inv = inverse;
|
||||
|
||||
@@ -59,8 +59,8 @@ namespace boost { namespace geometry { namespace projections
|
||||
std::string sweep_axis;
|
||||
int flip_axis;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -80,7 +80,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double r, Vx, Vy, Vz, tmp;
|
||||
|
||||
|
||||
/* Calculation of geocentric latitude. */
|
||||
lp_lat = atan (this->m_proj_parm.radius_p2 * tan (lp_lat));
|
||||
/* Calculation of the three components of the vector from satellite to
|
||||
@@ -109,7 +109,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double Vx, Vy, Vz, a, b, det, k;
|
||||
|
||||
|
||||
/* Setting three components of vector from satellite to position.*/
|
||||
Vx = -1.0;
|
||||
if(this->m_proj_parm.flip_axis)
|
||||
@@ -157,7 +157,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double Vx, Vy, Vz, tmp;
|
||||
|
||||
|
||||
/* Calculation of the three components of the vector from satellite to
|
||||
** position on earth surface (lon,lat).*/
|
||||
tmp = cos(lp_lat);
|
||||
@@ -183,7 +183,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double Vx, Vy, Vz, a, b, det, k;
|
||||
|
||||
|
||||
/* Setting three components of vector from satellite to position.*/
|
||||
Vx = -1.0;
|
||||
if(this->m_proj_parm.flip_axis)
|
||||
|
||||
@@ -68,7 +68,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double t = lp_lat * lp_lat;
|
||||
|
||||
|
||||
xy_y = lp_lat * (1. + t * C12);
|
||||
xy_x = lp_lon * (1. - Cp * t);
|
||||
t = lp_lon * lp_lon;
|
||||
|
||||
@@ -77,7 +77,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double s, c;
|
||||
|
||||
|
||||
xy_y = pj_mlfn(lp_lat, s = sin(lp_lat), c = cos(lp_lat), this->m_proj_parm.en);
|
||||
xy_x = lp_lon * c / sqrt(1. - this->m_par.es * s * s);
|
||||
}
|
||||
@@ -85,7 +85,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double s;
|
||||
|
||||
|
||||
if ((s = fabs(lp_lat = pj_inv_mlfn(xy_y, this->m_par.es, this->m_proj_parm.en))) < HALFPI) {
|
||||
s = sin(lp_lat);
|
||||
lp_lon = xy_x * sqrt(1. - this->m_par.es * s * s) / cos(lp_lat);
|
||||
@@ -119,7 +119,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
else {
|
||||
double k, V;
|
||||
int i;
|
||||
|
||||
|
||||
k = this->m_proj_parm.n * sin(lp_lat);
|
||||
for (i = MAX_ITER; i ; --i) {
|
||||
lp_lat -= V = (this->m_proj_parm.m * lp_lat + sin(lp_lat) - k) /
|
||||
@@ -172,7 +172,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
void setup_sinu(Parameters& par, par_gn_sinu& proj_parm)
|
||||
{
|
||||
pj_enfn(par.es, proj_parm.en);
|
||||
|
||||
|
||||
if (par.es) {
|
||||
// par.inv = e_inverse;
|
||||
// par.fwd = e_forward;
|
||||
|
||||
@@ -77,7 +77,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, cosphi, sinphi;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
@@ -115,7 +115,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double rh, cosz, sinz;
|
||||
|
||||
|
||||
rh = boost::math::hypot(xy_x, xy_y);
|
||||
sinz = sin(lp_lat = atan(rh));
|
||||
cosz = sqrt(1. - sinz * sinz);
|
||||
|
||||
@@ -59,7 +59,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
double XS;
|
||||
double YS;
|
||||
};
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
|
||||
@@ -71,7 +71,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double cosphi, d;
|
||||
|
||||
|
||||
d = sqrt(2./(1. + (cosphi = cos(lp_lat)) * cos(lp_lon *= this->m_proj_parm.w)));
|
||||
xy_x = this->m_proj_parm.m * d * cosphi * sin(lp_lon);
|
||||
xy_y = this->m_proj_parm.rm * d * sin(lp_lat);
|
||||
|
||||
@@ -79,7 +79,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double th1, c;
|
||||
int i;
|
||||
|
||||
|
||||
c = sin(lp_lat) * (lp_lat < 0. ? CS_ : CN_);
|
||||
for (i = NITER; i; --i) {
|
||||
lp_lat -= th1 = (lp_lat + sin(lp_lat) - c) / (1. + cos(lp_lat));
|
||||
@@ -92,7 +92,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double th;
|
||||
|
||||
|
||||
th = xy_y * ( xy_y < 0. ? RYCS : RYCN);
|
||||
if (fabs(th) > 1.)
|
||||
if (fabs(th) > ONETOL) throw proj_exception();
|
||||
|
||||
@@ -63,7 +63,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline int
|
||||
phi12(Parameters& par, par_imw_p& proj_parm, double *del, double *sig) {
|
||||
int err = 0;
|
||||
|
||||
|
||||
if (!pj_param(par.params, "tlat_1").i ||
|
||||
!pj_param(par.params, "tlat_2").i) {
|
||||
err = -41;
|
||||
@@ -80,13 +80,13 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline PXY
|
||||
loc_for(double const& lp_lam, double const& lp_phi, const Parameters& par, par_imw_p const& proj_parm, double *yc) {
|
||||
PXY xy;
|
||||
|
||||
|
||||
if (! lp_phi) {
|
||||
xy.x = lp_lam;
|
||||
xy.y = 0.;
|
||||
} else {
|
||||
double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
|
||||
|
||||
|
||||
sp = sin(lp_phi);
|
||||
m = pj_mlfn(lp_phi, sp, cos(lp_phi), proj_parm.en);
|
||||
xa = proj_parm.Pp + proj_parm.Qp * m;
|
||||
@@ -124,19 +124,19 @@ namespace boost { namespace geometry { namespace projections
|
||||
}
|
||||
return (xy);
|
||||
}
|
||||
|
||||
|
||||
template <typename Parameters>
|
||||
inline void
|
||||
xy(Parameters& par, par_imw_p& proj_parm, double phi, double *x, double *y, double *sp, double *R) {
|
||||
double F;
|
||||
|
||||
|
||||
*sp = sin(phi);
|
||||
*R = 1./(tan(phi) * sqrt(1. - par.es * *sp * *sp ));
|
||||
F = proj_parm.lam_1 * *sp;
|
||||
*y = *R * (1 - cos(F));
|
||||
*x = *R * sin(F);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -164,7 +164,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
PXY t;
|
||||
double yc = 0;
|
||||
|
||||
|
||||
lp_lat = this->m_proj_parm.phi_2;
|
||||
lp_lon = xy_x / cos(lp_lat);
|
||||
do {
|
||||
|
||||
@@ -51,36 +51,36 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double C_x;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
NOTES: According to EPSG the full Krovak projection method should have
|
||||
the following parameters. Within PROJ.4 the azimuth, and pseudo
|
||||
standard parallel are hardcoded in the algorithm and can't be
|
||||
standard parallel are hardcoded in the algorithm and can't be
|
||||
altered from outside. The others all have defaults to match the
|
||||
common usage with Krovak projection.
|
||||
|
||||
|
||||
lat_0 = latitude of centre of the projection
|
||||
|
||||
|
||||
lon_0 = longitude of centre of the projection
|
||||
|
||||
|
||||
** = azimuth (true) of the centre line passing through the centre of the projection
|
||||
|
||||
|
||||
** = latitude of pseudo standard parallel
|
||||
|
||||
|
||||
k = scale factor on the pseudo standard parallel
|
||||
|
||||
|
||||
x_0 = False Easting of the centre of the projection at the apex of the cone
|
||||
|
||||
|
||||
y_0 = False Northing of the centre of the projection at the apex of the cone
|
||||
|
||||
|
||||
**/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -100,17 +100,17 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
/* calculate xy from lat/lon */
|
||||
|
||||
|
||||
/* Constants, identical to inverse transform function */
|
||||
double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
|
||||
double gfi, u, fi0, deltav, s, d, eps, ro;
|
||||
|
||||
|
||||
|
||||
|
||||
s45 = 0.785398163397448; /* 45 DEG */
|
||||
s90 = 2 * s45;
|
||||
fi0 = this->m_par.phi0; /* Latitude of projection centre 49 DEG 30' */
|
||||
|
||||
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
|
||||
|
||||
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
|
||||
be set to 1 here.
|
||||
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
|
||||
e2=0.006674372230614;
|
||||
@@ -119,67 +119,67 @@ namespace boost { namespace geometry { namespace projections
|
||||
/* e2 = this->m_par.es;*/ /* 0.006674372230614; */
|
||||
e2 = 0.006674372230614;
|
||||
e = sqrt(e2);
|
||||
|
||||
|
||||
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
|
||||
|
||||
|
||||
uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */
|
||||
u0 = asin(sin(fi0) / alfa);
|
||||
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. );
|
||||
|
||||
|
||||
k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g;
|
||||
|
||||
|
||||
k1 = this->m_par.k0;
|
||||
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
|
||||
s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78 DEG 30'00" N */
|
||||
n = sin(s0);
|
||||
ro0 = k1 * n0 / tan(s0);
|
||||
ad = s90 - uq;
|
||||
|
||||
|
||||
/* Transformation */
|
||||
|
||||
|
||||
gfi =pow ( ((1. + e * sin(lp_lat)) /
|
||||
(1. - e * sin(lp_lat))) , (alfa * e / 2.));
|
||||
|
||||
|
||||
u= 2. * (atan(k * pow( tan(lp_lat / 2. + s45), alfa) / gfi)-s45);
|
||||
|
||||
|
||||
deltav = - lp_lon * alfa;
|
||||
|
||||
|
||||
s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
|
||||
d = asin(cos(u) * sin(deltav) / cos(s));
|
||||
eps = n * d;
|
||||
ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n) ;
|
||||
|
||||
|
||||
/* x and y are reverted! */
|
||||
xy_y = ro * cos(eps) / a;
|
||||
xy_x = ro * sin(eps) / a;
|
||||
|
||||
|
||||
if( !pj_param(this->m_par.params, "tczech").i )
|
||||
{
|
||||
xy_y *= -1.0;
|
||||
xy_x *= -1.0;
|
||||
}
|
||||
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
/* calculate lat/lon from xy */
|
||||
|
||||
|
||||
/* Constants, identisch wie in der Umkehrfunktion */
|
||||
double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
|
||||
double u, deltav, s, d, eps, ro, fi1, xy0;
|
||||
int ok;
|
||||
|
||||
|
||||
s45 = 0.785398163397448; /* 45 DEG */
|
||||
s90 = 2 * s45;
|
||||
fi0 = this->m_par.phi0; /* Latitude of projection centre 49 DEG 30' */
|
||||
|
||||
|
||||
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
|
||||
|
||||
|
||||
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
|
||||
be set to 1 here.
|
||||
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
|
||||
e2=0.006674372230614;
|
||||
@@ -188,47 +188,47 @@ namespace boost { namespace geometry { namespace projections
|
||||
/* e2 = this->m_par.es; */ /* 0.006674372230614; */
|
||||
e2 = 0.006674372230614;
|
||||
e = sqrt(e2);
|
||||
|
||||
|
||||
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
|
||||
uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */
|
||||
u0 = asin(sin(fi0) / alfa);
|
||||
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. );
|
||||
|
||||
|
||||
k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g;
|
||||
|
||||
|
||||
k1 = this->m_par.k0;
|
||||
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
|
||||
s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78 DEG 30'00" N */
|
||||
n = sin(s0);
|
||||
ro0 = k1 * n0 / tan(s0);
|
||||
ad = s90 - uq;
|
||||
|
||||
|
||||
|
||||
|
||||
/* Transformation */
|
||||
/* revert y, x*/
|
||||
xy0=xy_x;
|
||||
xy_x=xy_y;
|
||||
xy_y=xy0;
|
||||
|
||||
|
||||
if( !pj_param(this->m_par.params, "tczech").i )
|
||||
{
|
||||
xy_x *= -1.0;
|
||||
xy_y *= -1.0;
|
||||
}
|
||||
|
||||
|
||||
ro = sqrt(xy_x * xy_x + xy_y * xy_y);
|
||||
eps = atan2(xy_y, xy_x);
|
||||
d = eps / sin(s0);
|
||||
s = 2. * (atan( pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
|
||||
|
||||
|
||||
u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
|
||||
deltav = asin(cos(s) * sin(d) / cos(u));
|
||||
|
||||
|
||||
lp_lon = this->m_par.lam0 - deltav / alfa;
|
||||
|
||||
|
||||
/* ITERATION FOR lp_lat */
|
||||
fi1 = u;
|
||||
|
||||
|
||||
ok = 0;
|
||||
do
|
||||
{
|
||||
@@ -236,18 +236,18 @@ namespace boost { namespace geometry { namespace projections
|
||||
pow( tan(u / 2. + s45) , 1. / alfa) *
|
||||
pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
|
||||
) - s45);
|
||||
|
||||
|
||||
if (fabs(fi1 - lp_lat) < 0.000000000000001) ok=1;
|
||||
fi1 = lp_lat;
|
||||
|
||||
|
||||
}
|
||||
while (ok==0);
|
||||
|
||||
|
||||
lp_lon -= this->m_par.lam0;
|
||||
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
// Krovak
|
||||
@@ -259,14 +259,12 @@ namespace boost { namespace geometry { namespace projections
|
||||
* here Latitude Truescale */
|
||||
ts = pj_param(par.params, "rlat_ts").f;
|
||||
proj_parm.C_x = ts;
|
||||
|
||||
/* we want Bessel as fixed ellipsoid */
|
||||
par.a = 6377397.155;
|
||||
par.e = sqrt(par.es = 0.006674372230614);
|
||||
/* if latitude of projection center is not set, use 49d30'N */
|
||||
if (!pj_param(par.params, "tlat_0").i)
|
||||
par.phi0 = 0.863937979737193;
|
||||
|
||||
/* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
|
||||
/* that will correspond to using longitudes relative to greenwich */
|
||||
/* as input and output, instead of lat/long relative to Ferro */
|
||||
@@ -277,7 +275,6 @@ namespace boost { namespace geometry { namespace projections
|
||||
par.k0 = 0.9999;
|
||||
/* always the same */
|
||||
// par.inv = e_inverse;
|
||||
|
||||
// par.fwd = e_forward;
|
||||
}
|
||||
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double V1, V2, ps, sinps, cosps, sinps2, cosps2, I1, I2, I3, I4, I5, I6,
|
||||
x2, y2, t;
|
||||
|
||||
|
||||
V1 = this->m_proj_parm.A * log( tan(FORTPI + .5 * lp_lat) );
|
||||
t = this->m_par.e * sin(lp_lat);
|
||||
V2 = .5 * this->m_par.e * this->m_proj_parm.A * log ((1. + t)/(1. - t));
|
||||
@@ -104,7 +104,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
double x2, y2, V1, V2, V3, V4, t, t2, ps, pe, tpe, s,
|
||||
I7, I8, I9, I10, I11, d, Re;
|
||||
int i;
|
||||
|
||||
|
||||
x2 = xy_x * xy_x;
|
||||
y2 = xy_y * xy_y;
|
||||
V1 = 3. * xy_x * y2 - xy_x * x2;
|
||||
|
||||
@@ -88,7 +88,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
|
||||
|
||||
|
||||
coslam = cos(lp_lon);
|
||||
sinlam = sin(lp_lon);
|
||||
sinphi = sin(lp_lat);
|
||||
@@ -121,7 +121,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
goto eqcon;
|
||||
break;
|
||||
case EQUIT:
|
||||
xy_y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * this->m_proj_parm.ymf;
|
||||
xy_y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * this->m_proj_parm.ymf;
|
||||
eqcon:
|
||||
xy_x = this->m_proj_parm.xmf * b * cosb * sinlam;
|
||||
break;
|
||||
@@ -139,7 +139,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double cCe, sCe, q, rho, ab=0.0;
|
||||
|
||||
|
||||
switch (this->m_proj_parm.mode) {
|
||||
case EQUIT:
|
||||
case OBLIQ:
|
||||
@@ -197,7 +197,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, cosphi, sinphi;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
@@ -228,7 +228,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double cosz=0.0, rh, sinz=0.0;
|
||||
|
||||
|
||||
rh = boost::math::hypot(xy_x, xy_y);
|
||||
if ((lp_lat = rh * .5 ) > 1.) throw proj_exception();;
|
||||
lp_lat = 2. * asin(lp_lat);
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double v, c;
|
||||
|
||||
|
||||
if (fabs(fabs(lp_lat) - HALFPI) < TOL) {
|
||||
xy_x = 0;
|
||||
xy_y = lp_lat < 0 ? -2. : 2.;
|
||||
|
||||
@@ -75,7 +75,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double l2, p2;
|
||||
|
||||
|
||||
l2 = lp_lon * lp_lon;
|
||||
p2 = lp_lat * lp_lat;
|
||||
xy_x = lp_lon * (a10 + p2 * (a12 + l2 * a32 + p2 * a14));
|
||||
|
||||
@@ -49,9 +49,9 @@ namespace boost { namespace geometry { namespace projections
|
||||
#ifndef DOXYGEN_NO_DETAIL
|
||||
namespace detail { namespace latlong{
|
||||
|
||||
|
||||
|
||||
/* very loosely based upon DMA code by Bradford W. Drew */
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -69,14 +69,14 @@ namespace boost { namespace geometry { namespace projections
|
||||
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
|
||||
|
||||
xy_x = lp_lon / this->m_par.a;
|
||||
xy_y = lp_lat / this->m_par.a;
|
||||
}
|
||||
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
|
||||
|
||||
lp_lat = xy_y * this->m_par.a;
|
||||
lp_lon = xy_x * this->m_par.a;
|
||||
}
|
||||
|
||||
@@ -56,8 +56,8 @@ namespace boost { namespace geometry { namespace projections
|
||||
double r0, l, M0;
|
||||
double C;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
inline double /* func to compute dr */
|
||||
fS(double S, double C) {
|
||||
return(S * ( 1. + S * S * C));
|
||||
@@ -85,7 +85,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double S, r, dr;
|
||||
|
||||
|
||||
S = pj_mlfn(lp_lat, sin(lp_lat), cos(lp_lat), this->m_proj_parm.en) - this->m_proj_parm.M0;
|
||||
dr = fS(S, this->m_proj_parm.C);
|
||||
r = this->m_proj_parm.r0 - dr;
|
||||
@@ -97,7 +97,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double theta, dr, S, dif;
|
||||
int i;
|
||||
|
||||
|
||||
xy_x /= this->m_par.k0;
|
||||
xy_y /= this->m_par.k0;
|
||||
theta = atan2(xy_x , this->m_proj_parm.r0 - xy_y);
|
||||
|
||||
@@ -60,14 +60,14 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void
|
||||
seraz0(double lam, double mult, Parameters& par, par_lsat& proj_parm) {
|
||||
double sdsq, h, s, fc, sd, sq, d__1;
|
||||
|
||||
|
||||
lam *= DEG_TO_RAD;
|
||||
sd = sin(lam);
|
||||
sdsq = sd * sd;
|
||||
s = proj_parm.p22 * proj_parm.sa * cos(lam) * sqrt((1. + proj_parm.t * sdsq) / ((
|
||||
1. + proj_parm.w * sdsq) * (1. + proj_parm.q * sdsq)));
|
||||
d__1 = 1. + proj_parm.q * sdsq;
|
||||
h = sqrt((1. + proj_parm.q * sdsq) / (1. + proj_parm.w * sdsq)) * ((1. +
|
||||
h = sqrt((1. + proj_parm.q * sdsq) / (1. + proj_parm.w * sdsq)) * ((1. +
|
||||
proj_parm.w * sdsq) / (d__1 * d__1) - proj_parm.p22 * proj_parm.ca);
|
||||
sq = sqrt(proj_parm.xj * proj_parm.xj + s * s);
|
||||
proj_parm.b += fc = mult * (h * proj_parm.xj - s * s) / sq;
|
||||
@@ -98,7 +98,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
int l, nn;
|
||||
double lamt, xlam, sdsq, c, d, s, lamdp, phidp, lampp, tanph,
|
||||
lamtp, cl, sd, sp, fac, sav, tanphi;
|
||||
|
||||
|
||||
if (lp_lat > HALFPI)
|
||||
lp_lat = HALFPI;
|
||||
else if (lp_lat < -HALFPI)
|
||||
@@ -131,7 +131,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
}
|
||||
if (l) {
|
||||
sp = sin(lp_lat);
|
||||
phidp = aasin((this->m_par.one_es * this->m_proj_parm.ca * sp - this->m_proj_parm.sa * cos(lp_lat) *
|
||||
phidp = aasin((this->m_par.one_es * this->m_proj_parm.ca * sp - this->m_proj_parm.sa * cos(lp_lat) *
|
||||
sin(lamt)) / sqrt(1. - this->m_par.es * sp * sp));
|
||||
tanph = log(tan(FORTPI + .5 * phidp));
|
||||
sd = sin(lamdp);
|
||||
@@ -150,7 +150,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
int nn;
|
||||
double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
|
||||
|
||||
|
||||
lamdp = xy_x / this->m_proj_parm.b;
|
||||
nn = 50;
|
||||
do {
|
||||
@@ -165,7 +165,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
lamdp /= this->m_proj_parm.b;
|
||||
} while (fabs(lamdp - sav) >= TOL && --nn);
|
||||
sl = sin(lamdp);
|
||||
fac = exp(sqrt(1. + s * s / this->m_proj_parm.xj / this->m_proj_parm.xj) * (xy_y -
|
||||
fac = exp(sqrt(1. + s * s / this->m_proj_parm.xj / this->m_proj_parm.xj) * (xy_y -
|
||||
this->m_proj_parm.c1 * sl - this->m_proj_parm.c3 * sin(lamdp * 3.)));
|
||||
phidp = 2. * (atan(fac) - FORTPI);
|
||||
dd = sl * sl;
|
||||
@@ -173,9 +173,9 @@ namespace boost { namespace geometry { namespace projections
|
||||
lamdp -= TOL;
|
||||
spp = sin(phidp);
|
||||
sppsq = spp * spp;
|
||||
lamt = atan(((1. - sppsq * this->m_par.rone_es) * tan(lamdp) *
|
||||
lamt = atan(((1. - sppsq * this->m_par.rone_es) * tan(lamdp) *
|
||||
this->m_proj_parm.ca - spp * this->m_proj_parm.sa * sqrt((1. + this->m_proj_parm.q * dd) * (
|
||||
1. - sppsq) - sppsq * this->m_proj_parm.u) / cos(lamdp)) / (1. - sppsq
|
||||
1. - sppsq) - sppsq * this->m_proj_parm.u) / cos(lamdp)) / (1. - sppsq
|
||||
* (1. + this->m_proj_parm.u)));
|
||||
sl = lamt >= 0. ? 1. : -1.;
|
||||
scl = cos(lamdp) >= 0. ? 1. : -1;
|
||||
|
||||
@@ -74,7 +74,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double k, V, t;
|
||||
int i;
|
||||
|
||||
|
||||
k = C3 * sin(lp_lat);
|
||||
for (i = MAX_ITER; i ; --i) {
|
||||
t = lp_lat / C2;
|
||||
@@ -91,7 +91,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
lp_lat = C2 * (t = aasin(xy_y / C_y));
|
||||
lp_lon = xy_x / (C_x * (1. + 3. * cos(lp_lat)/cos(t)));
|
||||
lp_lat = aasin((C1 * sin(t) + sin(lp_lat)) / C3);
|
||||
|
||||
@@ -75,7 +75,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double th1, c;
|
||||
int i;
|
||||
|
||||
|
||||
c = C * sin(lp_lat);
|
||||
for (i = NITER; i; --i) {
|
||||
lp_lat -= th1 = (sin(.5*lp_lat) + sin(lp_lat) - c) /
|
||||
@@ -89,7 +89,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
lp_lat = RYC * xy_y;
|
||||
if (fabs(lp_lat) > 1.) {
|
||||
if (fabs(lp_lat) > ONETOL) throw proj_exception();
|
||||
|
||||
@@ -57,7 +57,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
int n;
|
||||
};
|
||||
/* based upon Snyder and Linck, USGS-NMD */
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -78,7 +78,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double sinlon, coslon, esphi, chi, schi, cchi, s;
|
||||
COMPLEX p;
|
||||
|
||||
|
||||
sinlon = sin(lp_lon);
|
||||
coslon = cos(lp_lon);
|
||||
esphi = this->m_par.e * sin(lp_lat);
|
||||
@@ -99,7 +99,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
int nn;
|
||||
COMPLEX p, fxy, fpxy, dp;
|
||||
double den, rh = 0, z, sinz = 0, cosz = 0, chi, phi = 0, dphi, esphi;
|
||||
|
||||
|
||||
p.r = xy_x;
|
||||
p.i = xy_y;
|
||||
for (nn = 20; nn ;--nn) {
|
||||
@@ -137,7 +137,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
}
|
||||
if (nn) {
|
||||
lp_lat = phi;
|
||||
lp_lon = atan2(p.r * sinz, rh * this->m_proj_parm.cchio * cosz - p.i *
|
||||
lp_lon = atan2(p.r * sinz, rh * this->m_proj_parm.cchio * cosz - p.i *
|
||||
this->m_proj_parm.schio * sinz);
|
||||
} else
|
||||
lp_lon = lp_lat = HUGE_VAL;
|
||||
|
||||
@@ -75,7 +75,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double k, V;
|
||||
int i;
|
||||
|
||||
|
||||
k = this->m_proj_parm.C_p * sin(lp_lat);
|
||||
for (i = MAX_ITER; i ; --i) {
|
||||
lp_lat -= V = (lp_lat + sin(lp_lat) - k) /
|
||||
|
||||
@@ -68,7 +68,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double k, V;
|
||||
int i;
|
||||
|
||||
|
||||
k = 2. * sin(lp_lat);
|
||||
V = lp_lat * lp_lat;
|
||||
lp_lat *= 1.00371 + V * (-0.0935382 + V * -0.011412);
|
||||
|
||||
@@ -74,7 +74,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double V, c, p;
|
||||
int i;
|
||||
|
||||
|
||||
p = 0.5 * xy_y;
|
||||
for (i = NITER; i ; --i) {
|
||||
c = cos(0.5 * lp_lat);
|
||||
|
||||
@@ -79,7 +79,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
xy_y = lp_lat;
|
||||
} else {
|
||||
double tb, c, d, m, n, r2, sp;
|
||||
|
||||
|
||||
tb = HALFPI / lp_lon - lp_lon / HALFPI;
|
||||
c = lp_lat / HALFPI;
|
||||
d = (1 - c * c)/((sp = sin(lp_lat)) - c);
|
||||
|
||||
@@ -89,7 +89,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, cosphi, sinphi;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
@@ -126,7 +126,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
}
|
||||
if (this->m_proj_parm.tilt) {
|
||||
double yt, ba;
|
||||
|
||||
|
||||
yt = xy_y * this->m_proj_parm.cg + xy_x * this->m_proj_parm.sg;
|
||||
ba = 1. / (yt * this->m_proj_parm.sw * this->m_proj_parm.h + this->m_proj_parm.cw);
|
||||
xy_x = (xy_x * this->m_proj_parm.cg - xy_y * this->m_proj_parm.sg) * this->m_proj_parm.cw * ba;
|
||||
@@ -137,10 +137,10 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double rh, cosz, sinz;
|
||||
|
||||
|
||||
if (this->m_proj_parm.tilt) {
|
||||
double bm, bq, yt;
|
||||
|
||||
|
||||
yt = 1./(this->m_proj_parm.pn1 - xy_y * this->m_proj_parm.sw);
|
||||
bm = this->m_proj_parm.pn1 * xy_x * yt;
|
||||
bq = this->m_proj_parm.pn1 * xy_y * this->m_proj_parm.cw * yt;
|
||||
|
||||
@@ -53,10 +53,10 @@ namespace boost { namespace geometry { namespace projections
|
||||
static const int Ntpsi = 9;
|
||||
static const int Ntphi = 8;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
static COMPLEX
|
||||
bf[] = {
|
||||
{.7557853228, 0.0},
|
||||
@@ -90,7 +90,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
COMPLEX p;
|
||||
double *C;
|
||||
int i;
|
||||
|
||||
|
||||
lp_lat = (lp_lat - this->m_par.phi0) * RAD_TO_SEC5;
|
||||
for (p.r = *(C = tpsi + (i = Ntpsi)); i ; --i)
|
||||
p.r = *--C + lp_lat * p.r;
|
||||
@@ -106,7 +106,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
int nn, i;
|
||||
COMPLEX p, f, fp, dp;
|
||||
double den, *C;
|
||||
|
||||
|
||||
p.r = xy_y;
|
||||
p.i = xy_x;
|
||||
for (nn = 20; nn ;--nn) {
|
||||
|
||||
@@ -78,9 +78,9 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, sinphi, cosphi;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
coslam = cos(lp_lon);
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
@@ -93,7 +93,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double coslam, sinphi, cosphi;
|
||||
|
||||
|
||||
m_proj_parm.link->inv(xy_x, xy_y, lp_lon, lp_lat);
|
||||
if (lp_lon != HUGE_VAL) {
|
||||
coslam = cos(lp_lon -= this->m_proj_parm.lamp);
|
||||
@@ -124,9 +124,9 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double cosphi, coslam;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
lp_lon = adjlon(aatan2(cosphi * sin(lp_lon), sin(lp_lat)) + this->m_proj_parm.lamp);
|
||||
@@ -137,7 +137,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double cosphi, t;
|
||||
|
||||
|
||||
m_proj_parm.link->inv(xy_x, xy_y, lp_lon, lp_lat);
|
||||
if (lp_lon != HUGE_VAL) {
|
||||
cosphi = cos(lp_lat);
|
||||
@@ -154,8 +154,8 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
int i;
|
||||
double phip;
|
||||
|
||||
|
||||
|
||||
|
||||
Parameters pj;
|
||||
/* copy existing header into new */
|
||||
par.es = 0.;
|
||||
@@ -175,7 +175,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
pj.one_es = pj.rone_es = 1.;
|
||||
pj.es = pj.e = 0.;
|
||||
pj.name = pj_param(par.params, "so_proj").s;
|
||||
|
||||
|
||||
factory<Geographic, Cartesian, Parameters> fac;
|
||||
if (create)
|
||||
{
|
||||
|
||||
@@ -75,7 +75,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
xy_y = sin(lp_lon);
|
||||
/*
|
||||
xy_x = atan2((tan(lp_lat) * this->m_proj_parm.cosphi + this->m_proj_parm.sinphi * xy_y) , cos(lp_lon));
|
||||
@@ -91,7 +91,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double t, s;
|
||||
|
||||
|
||||
xy_y /= this->m_proj_parm.rok;
|
||||
xy_x /= this->m_proj_parm.rtk;
|
||||
t = sqrt(1. - xy_y * xy_y);
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double Az, M, N, cp, sp, cl, shz;
|
||||
|
||||
|
||||
cp = cos(lp_lat);
|
||||
sp = sin(lp_lat);
|
||||
cl = cos(lp_lon);
|
||||
@@ -88,7 +88,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double N, M, xp, yp, z, Az, cz, sz, cAz;
|
||||
|
||||
|
||||
N = this->m_proj_parm.hn * aasin(xy_y * this->m_proj_parm.rn);
|
||||
M = this->m_proj_parm.hm * aasin(xy_x * this->m_proj_parm.rm * cos(N * this->m_proj_parm.two_r_n) / cos(N));
|
||||
xp = 2. * sin(M);
|
||||
|
||||
@@ -61,8 +61,8 @@ namespace boost { namespace geometry { namespace projections
|
||||
double v_pole_n, v_pole_s, u_0;
|
||||
int no_rot;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -82,7 +82,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double Q, S, T, U, V, temp, u, v;
|
||||
|
||||
|
||||
if (fabs(fabs(lp_lat) - HALFPI) > EPS) {
|
||||
Q = this->m_proj_parm.E / pow(pj_tsfn(lp_lat, sin(lp_lat), this->m_par.e), this->m_proj_parm.B);
|
||||
temp = 1. / Q;
|
||||
@@ -95,7 +95,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
v = 0.5 * this->m_proj_parm.ArB * log((1. - U)/(1. + U));
|
||||
temp = cos(this->m_proj_parm.B * lp_lon);
|
||||
u = (fabs(temp) < TOL) ? this->m_proj_parm.AB * lp_lon :
|
||||
this->m_proj_parm.ArB * atan2((S * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam) , temp);
|
||||
this->m_proj_parm.ArB * atan2((S * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam) , temp);
|
||||
} else {
|
||||
v = lp_lat > 0 ? this->m_proj_parm.v_pole_n : this->m_proj_parm.v_pole_s;
|
||||
u = this->m_proj_parm.ArB * lp_lat;
|
||||
@@ -113,7 +113,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double u, v, Qp, Sp, Tp, Vp, Up;
|
||||
|
||||
|
||||
if (this->m_proj_parm.no_rot) {
|
||||
v = xy_y;
|
||||
u = xy_x;
|
||||
@@ -153,7 +153,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
gamma = pj_param(par.params, "rgamma").f;
|
||||
if (alp || gam) {
|
||||
lamc = pj_param(par.params, "rlonc").f;
|
||||
no_off =
|
||||
no_off =
|
||||
/* For libproj4 compatability */
|
||||
pj_param(par.params, "tno_off").i
|
||||
/* for backward compatibility */
|
||||
|
||||
@@ -77,7 +77,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, cosphi, sinphi;
|
||||
|
||||
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
switch (this->m_proj_parm.mode) {
|
||||
@@ -100,12 +100,12 @@ namespace boost { namespace geometry { namespace projections
|
||||
xy_x = cosphi * sin(lp_lon);
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double rh, cosc, sinc;
|
||||
|
||||
|
||||
if ((sinc = (rh = boost::math::hypot(xy_x, xy_y))) > 1.) {
|
||||
if ((sinc - 1.) > EPS10) throw proj_exception();;
|
||||
sinc = 1.;
|
||||
@@ -145,7 +145,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
// Orthographic
|
||||
|
||||
@@ -78,7 +78,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double ms, sp, cp;
|
||||
|
||||
|
||||
if (fabs(lp_lat) <= TOL) { xy_x = lp_lon; xy_y = -this->m_proj_parm.ml0; }
|
||||
else {
|
||||
sp = sin(lp_lat);
|
||||
@@ -95,7 +95,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
else {
|
||||
double r, c, sp, cp, s2ph, ml, mlb, mlp, dPhi;
|
||||
int i;
|
||||
|
||||
|
||||
r = xy_y * xy_y + xy_x * xy_x;
|
||||
for (lp_lat = xy_y, i = I_ITER; i ; --i) {
|
||||
sp = sin(lp_lat);
|
||||
@@ -139,7 +139,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double cot, E;
|
||||
|
||||
|
||||
if (fabs(lp_lat) <= TOL) { xy_x = lp_lon; xy_y = this->m_proj_parm.ml0; }
|
||||
else {
|
||||
cot = 1. / tan(lp_lat);
|
||||
@@ -152,7 +152,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double B, dphi, tp;
|
||||
int i;
|
||||
|
||||
|
||||
if (fabs(xy_y = this->m_par.phi0 + xy_y) <= TOL) { lp_lon = xy_x; lp_lat = 0.; }
|
||||
else {
|
||||
lp_lat = xy_y;
|
||||
|
||||
@@ -72,7 +72,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double p, c, s, V;
|
||||
int i;
|
||||
|
||||
|
||||
p = C_p * sin(lp_lat);
|
||||
s = lp_lat * lp_lat;
|
||||
lp_lat *= 0.615709 + s * ( 0.00909953 + s * 0.0046292 );
|
||||
@@ -93,7 +93,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double c;
|
||||
|
||||
|
||||
lp_lat = aasin(xy_y / C_y);
|
||||
lp_lon = xy_x / (C_x * ((c = cos(lp_lat)) - 0.5));
|
||||
lp_lat = aasin((lp_lat + sin(lp_lat) * (c - 1.)) / C_p);
|
||||
|
||||
@@ -75,7 +75,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double p, r, V;
|
||||
int i;
|
||||
|
||||
|
||||
p = this->m_proj_parm.B * sin(lp_lat);
|
||||
lp_lat *= 1.10265779;
|
||||
for (i = NITER; i ; --i) {
|
||||
@@ -94,7 +94,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double r;
|
||||
|
||||
|
||||
lp_lat = xy_y / this->m_proj_parm.C_y;
|
||||
r = sqrt(1. + lp_lat * lp_lat);
|
||||
lp_lon = xy_x / (this->m_proj_parm.C_x * (this->m_proj_parm.D - r));
|
||||
|
||||
@@ -55,12 +55,12 @@ namespace boost { namespace geometry { namespace projections
|
||||
static const double ONEEPS = 1.000001;
|
||||
static const double EPS = 1e-8;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
struct COEFS {
|
||||
double c0, c1, c2, c3;
|
||||
};
|
||||
|
||||
|
||||
static const struct COEFS X[] = {
|
||||
{1, 2.2199e-17, -7.15515e-05, 3.1103e-06},
|
||||
{0.9986, -0.000482243, -2.4897e-05, -1.3309e-06},
|
||||
@@ -127,7 +127,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
int i;
|
||||
double dphi;
|
||||
|
||||
|
||||
i = int_floor((dphi = fabs(lp_lat)) * C1);
|
||||
if (i >= NODES) i = NODES - 1;
|
||||
dphi = RAD_TO_DEG * (dphi - RC1 * i);
|
||||
@@ -141,7 +141,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
int i;
|
||||
double t, t1;
|
||||
struct COEFS T;
|
||||
|
||||
|
||||
lp_lon = xy_x / FXC;
|
||||
lp_lat = fabs(xy_y / FYC);
|
||||
if (lp_lat >= 1.) { /* simple pathologic cases */
|
||||
|
||||
@@ -76,7 +76,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double s, al, cp, sp, al2, s2;
|
||||
|
||||
|
||||
cp = cos(lp_lat);
|
||||
sp = sin(lp_lat);
|
||||
s = proj_mdist(lp_lat, sp, cp, this->m_proj_parm.en) - this->m_proj_parm.s0;
|
||||
@@ -93,7 +93,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double s, al, x = xy_x / this->m_par.k0, y = xy_y / this->m_par.k0, x2, y2;;
|
||||
|
||||
|
||||
x2 = x * x;
|
||||
y2 = y * y;
|
||||
al = x*(1.-this->m_proj_parm.C1*y2+x2*(this->m_proj_parm.C2+this->m_proj_parm.C3*y-this->m_proj_parm.C4*x2+this->m_proj_parm.C5*y2-this->m_proj_parm.C7*x2*y)
|
||||
@@ -113,7 +113,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double N0, es2, t, t2, R_R0_2, R_R0_4;
|
||||
proj_mdist_ini(par.es, proj_parm.en);
|
||||
|
||||
|
||||
es2 = sin(par.phi0);
|
||||
proj_parm.s0 = proj_mdist(par.phi0, es2, cos(par.phi0), proj_parm.en);
|
||||
t = 1. - (es2 = par.es * es2 * es2);
|
||||
|
||||
@@ -74,7 +74,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double fa;
|
||||
|
||||
|
||||
if (this->m_proj_parm.mode)
|
||||
fa = tan(lp_lon * this->m_proj_parm.fxb) * this->m_proj_parm.fxa;
|
||||
else
|
||||
|
||||
@@ -72,7 +72,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
phi12(Parameters& par, par_sconics& proj_parm, double *del) {
|
||||
double p1, p2;
|
||||
int err = 0;
|
||||
|
||||
|
||||
if (!pj_param(par.params, "tlat_1").i ||
|
||||
!pj_param(par.params, "tlat_2").i) {
|
||||
err = -41;
|
||||
@@ -105,7 +105,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double rho;
|
||||
|
||||
|
||||
switch (this->m_proj_parm.type) {
|
||||
case MURD2:
|
||||
rho = this->m_proj_parm.rho_c + tan(this->m_proj_parm.sig - lp_lat);
|
||||
@@ -124,7 +124,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double rho;
|
||||
|
||||
|
||||
rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho_0 - xy_y);
|
||||
if (this->m_proj_parm.n < 0.) {
|
||||
rho = - rho;
|
||||
@@ -180,7 +180,6 @@ namespace boost { namespace geometry { namespace projections
|
||||
proj_parm.n = sin(proj_parm.sig) * sin(del) / del;
|
||||
del *= 0.5;
|
||||
proj_parm.rho_c = del / (tan(del) * tan(proj_parm.sig)) + proj_parm.sig;
|
||||
|
||||
proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
|
||||
break;
|
||||
case PCONIC:
|
||||
|
||||
@@ -72,7 +72,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double phip, lamp, phipp, lampp, sp, cp;
|
||||
|
||||
|
||||
sp = this->m_par.e * sin(lp_lat);
|
||||
phip = 2.* atan( exp( this->m_proj_parm.c * (
|
||||
log(tan(FORTPI + 0.5 * lp_lat)) - this->m_proj_parm.hlf_e * log((1. + sp)/(1. - sp)))
|
||||
@@ -89,7 +89,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double phip, lamp, phipp, lampp, cp, esp, con, delp;
|
||||
int i;
|
||||
|
||||
|
||||
phipp = 2. * (atan(exp(xy_y / this->m_proj_parm.kR)) - FORTPI);
|
||||
lampp = xy_x / this->m_proj_parm.kR;
|
||||
cp = cos(phipp);
|
||||
|
||||
@@ -90,7 +90,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double coslam, sinlam, sinX=0.0, cosX=0.0, X, A, sinphi;
|
||||
|
||||
|
||||
coslam = cos(lp_lon);
|
||||
sinlam = sin(lp_lon);
|
||||
sinphi = sin(lp_lat);
|
||||
@@ -126,7 +126,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
|
||||
int i;
|
||||
|
||||
|
||||
rho = boost::math::hypot(xy_x, xy_y);
|
||||
switch (this->m_proj_parm.mode) {
|
||||
case OBLIQ:
|
||||
@@ -137,7 +137,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
phi_l = asin(cosphi * this->m_proj_parm.sinX1);
|
||||
else
|
||||
phi_l = asin(cosphi * this->m_proj_parm.sinX1 + (xy_y * sinphi * this->m_proj_parm.cosX1 / rho));
|
||||
|
||||
|
||||
tp = tan(.5 * (HALFPI + phi_l));
|
||||
xy_x *= sinphi;
|
||||
xy_y = rho * this->m_proj_parm.cosX1 * cosphi - xy_y * this->m_proj_parm.sinX1* sinphi;
|
||||
@@ -185,7 +185,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double sinphi, cosphi, coslam, sinlam;
|
||||
|
||||
|
||||
sinphi = sin(lp_lat);
|
||||
cosphi = cos(lp_lat);
|
||||
coslam = cos(lp_lon);
|
||||
@@ -216,7 +216,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double c, rh, sinc, cosc;
|
||||
|
||||
|
||||
sinc = sin(c = 2. * atan((rh = boost::math::hypot(xy_x, xy_y)) / this->m_proj_parm.akm1));
|
||||
cosc = cos(c);
|
||||
lp_lon = 0.;
|
||||
|
||||
@@ -59,9 +59,9 @@ namespace boost { namespace geometry { namespace projections
|
||||
double R2;
|
||||
gauss::GAUSS en;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// template class, using CRTP to implement forward/inverse
|
||||
template <typename Geographic, typename Cartesian, typename Parameters>
|
||||
@@ -81,7 +81,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double cosc, sinc, cosl_, k;
|
||||
|
||||
|
||||
detail::gauss::gauss(m_proj_parm.en, lp_lon, lp_lat);
|
||||
sinc = sin(lp_lat);
|
||||
cosc = cos(lp_lat);
|
||||
@@ -94,7 +94,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double rho, c, sinc, cosc;
|
||||
|
||||
|
||||
xy_x /= this->m_par.k0;
|
||||
xy_y /= this->m_par.k0;
|
||||
if((rho = boost::math::hypot(xy_x, xy_y))) {
|
||||
@@ -116,7 +116,6 @@ namespace boost { namespace geometry { namespace projections
|
||||
template <typename Parameters>
|
||||
void setup_sterea(Parameters& par, par_sterea& proj_parm)
|
||||
{
|
||||
|
||||
double R;
|
||||
proj_parm.en = detail::gauss::gauss_ini(par.e, par.phi0, proj_parm.phic0, R);
|
||||
proj_parm.sinc0 = sin(proj_parm.phic0);
|
||||
|
||||
@@ -72,7 +72,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double c;
|
||||
|
||||
|
||||
xy_x = this->m_proj_parm.C_x * lp_lon * cos(lp_lat);
|
||||
xy_y = this->m_proj_parm.C_y;
|
||||
lp_lat *= this->m_proj_parm.C_p;
|
||||
@@ -89,7 +89,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double c;
|
||||
|
||||
|
||||
xy_y /= this->m_proj_parm.C_y;
|
||||
c = cos(lp_lat = this->m_proj_parm.tan_mode ? atan(xy_y) : aasin(xy_y));
|
||||
lp_lat /= this->m_proj_parm.C_p;
|
||||
|
||||
@@ -71,7 +71,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double b, bt;
|
||||
|
||||
|
||||
b = cos(lp_lat) * sin(lp_lon);
|
||||
if ((bt = 1. - b * b) < EPS10) throw proj_exception();;
|
||||
xy_x = b / sqrt(bt);
|
||||
|
||||
@@ -76,7 +76,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
xy_y = xy_y * this->m_proj_parm.rk0 + this->m_par.phi0;
|
||||
xy_x *= this->m_par.k0;
|
||||
t = sqrt(1. - xy_x * xy_x);
|
||||
|
||||
@@ -86,12 +86,12 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double al, als, n, cosphi, sinphi, t;
|
||||
|
||||
|
||||
/*
|
||||
* Fail if our longitude is more than 90 degrees from the
|
||||
* central meridian since the results are essentially garbage.
|
||||
* Fail if our longitude is more than 90 degrees from the
|
||||
* central meridian since the results are essentially garbage.
|
||||
* Is error -20 really an appropriate return value?
|
||||
*
|
||||
*
|
||||
* http://trac.osgeo.org/proj/ticket/5
|
||||
*/
|
||||
if( lp_lon < -HALFPI || lp_lon > HALFPI )
|
||||
@@ -101,7 +101,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
throw proj_exception(-14 );
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
sinphi = sin(lp_lat); cosphi = cos(lp_lat);
|
||||
t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.;
|
||||
t *= t;
|
||||
@@ -125,7 +125,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double n, con, cosphi, d, ds, sinphi, t;
|
||||
|
||||
|
||||
lp_lat = pj_inv_mlfn(this->m_proj_parm.ml0 + xy_y / this->m_par.k0, this->m_par.es, this->m_proj_parm.en);
|
||||
if (fabs(lp_lat) >= HALFPI) {
|
||||
lp_lat = xy_y < 0. ? -HALFPI : HALFPI;
|
||||
@@ -172,12 +172,12 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double b, cosphi;
|
||||
|
||||
|
||||
/*
|
||||
* Fail if our longitude is more than 90 degrees from the
|
||||
* central meridian since the results are essentially garbage.
|
||||
* Fail if our longitude is more than 90 degrees from the
|
||||
* central meridian since the results are essentially garbage.
|
||||
* Is error -20 really an appropriate return value?
|
||||
*
|
||||
*
|
||||
* http://trac.osgeo.org/proj/ticket/5
|
||||
*/
|
||||
if( lp_lon < -HALFPI || lp_lon > HALFPI )
|
||||
@@ -187,7 +187,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
throw proj_exception(-14 );
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
b = (cosphi = cos(lp_lat)) * sin(lp_lon);
|
||||
if (fabs(fabs(b) - 1.) <= EPS10) throw proj_exception();;
|
||||
xy_x = this->m_proj_parm.ml0 * log((1. + b) / (1. - b));
|
||||
@@ -203,7 +203,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double h, g;
|
||||
|
||||
|
||||
h = exp(xy_x / this->m_proj_parm.esp);
|
||||
g = .5 * (h - 1. / h);
|
||||
h = cos(this->m_par.phi0 + xy_y / this->m_proj_parm.esp);
|
||||
|
||||
@@ -71,7 +71,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double t, z1, z2, dl1, dl2, sp, cp;
|
||||
|
||||
|
||||
sp = sin(lp_lat);
|
||||
cp = cos(lp_lat);
|
||||
z1 = aacos(this->m_proj_parm.sp1 * sp + this->m_proj_parm.cp1 * cp * cos(dl1 = lp_lon + this->m_proj_parm.dlam2));
|
||||
@@ -88,7 +88,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double cz1, cz2, s, d, cp, sp;
|
||||
|
||||
|
||||
cz1 = cos(boost::math::hypot(xy_y, xy_x + this->m_proj_parm.hz0));
|
||||
cz2 = cos(boost::math::hypot(xy_y, xy_x - this->m_proj_parm.hz0));
|
||||
s = cz1 + cz2;
|
||||
|
||||
@@ -70,7 +70,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double t;
|
||||
|
||||
|
||||
t = lp_lat = aasin(this->m_proj_parm.n * sin(lp_lat));
|
||||
xy_x = this->m_proj_parm.m * lp_lon * cos(lp_lat);
|
||||
t *= t;
|
||||
|
||||
@@ -73,7 +73,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double al, al2, g, g2, p2;
|
||||
|
||||
|
||||
p2 = fabs(lp_lat / HALFPI);
|
||||
if ((p2 - TOL) > 1.) throw proj_exception();;
|
||||
if (p2 > 1.)
|
||||
@@ -107,7 +107,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
|
||||
{
|
||||
double t, c0, c1, c2, c3, al, r2, r, m, d, ay, x2, y2;
|
||||
|
||||
|
||||
x2 = xy_x * xy_x;
|
||||
if ((ay = fabs(xy_y)) < TOL) {
|
||||
lp_lat = 0.;
|
||||
|
||||
@@ -72,7 +72,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double x1, at, bt, ct;
|
||||
|
||||
|
||||
bt = fabs(TWORPI * lp_lat);
|
||||
if ((ct = 1. - bt * bt) < 0.)
|
||||
ct = 0.;
|
||||
|
||||
@@ -67,7 +67,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double x1, t, bt, ct, ft, bt2, ct2, dt, dt2;
|
||||
|
||||
|
||||
if (fabs(lp_lat) < TOL) {
|
||||
xy_x = lp_lon;
|
||||
xy_y = 0.;
|
||||
|
||||
@@ -65,7 +65,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
|
||||
{
|
||||
double theta, ct, D;
|
||||
|
||||
|
||||
theta = asin(xy_y = 0.90630778703664996 * sin(lp_lat));
|
||||
xy_x = 2.66723 * (ct = cos(theta)) * sin(lp_lon /= 3.);
|
||||
xy_y *= 1.24104 * (D = 1/(sqrt(0.5 * (1 + ct * cos(lp_lon)))));
|
||||
|
||||
@@ -74,7 +74,7 @@ namespace boost { namespace geometry { namespace projections
|
||||
{
|
||||
double k, V;
|
||||
int i;
|
||||
|
||||
|
||||
xy_y = lp_lat * TWO_D_PI;
|
||||
k = PI * sin(lp_lat);
|
||||
lp_lat *= 1.8;
|
||||
|
||||
Reference in New Issue
Block a user