src/map/krovak.cpp

   1 #include "ellipsoid.h"
   2 #include "krovak.h"
   3
   4 Krovak::Krovak(const Ellipsoid *ellipsoid, double standardParallel,
   5   double azimuth, double scale, double centerLatitude, double longitudeOrigin,
   6   double falseEasting, double falseNorthing)
   7 {
   8         double phiC = deg2rad(centerLatitude);
   9         double sinPhiC = sin(phiC);
  10         double sinPhiC2 = sinPhiC * sinPhiC;
  11         double cosPhiC = cos(phiC);
  12         double cosPhiC2 = cosPhiC * cosPhiC;
  13         double cosPhiC4 = cosPhiC2 * cosPhiC2;
  14         double alphaC = deg2rad(azimuth);
  15
  16         _phiP = deg2rad(standardParallel);
  17         _e = sqrt(ellipsoid->es());
  18         _A = ellipsoid->radius() * sqrt(1.0 - ellipsoid->es())
  19           / (1.0 - ellipsoid->es() * sinPhiC2);
  20         _B = sqrt(1.0 + (ellipsoid->es() * cosPhiC4 / (1.0 - ellipsoid->es())));
  21         double gamma0 = asin(sinPhiC / _B);
  22         _t0 = tan(M_PI_4 + gamma0 / 2.0) * pow((1.0 + _e * sinPhiC) /
  23           (1.0 - _e * sinPhiC), _e*_B / 2.0) / pow(tan(M_PI_4 + phiC/2.0), _B);
  24         _n = sin(_phiP);
  25         _r0 = scale * _A / tan(_phiP);
  26         _FE = falseEasting;
  27         _FN = falseNorthing;
  28         _cosAlphaC = cos(alphaC);
  29         _sinAlphaC = sin(alphaC);
  30         _lambda0 = deg2rad(longitudeOrigin);
  31 }
  32
  33 PointD Krovak::ll2xy(const Coordinates &c) const
  34 {
  35         double phi = deg2rad(c.lat());
  36         double lambda = deg2rad(c.lon());
  37         double eSinPhi = _e * sin(phi);
  38         double U = 2.0 * (atan(_t0 * pow(tan(phi/2.0 + M_PI_4), _B)
  39           / pow((1.0 + eSinPhi) / (1.0 - eSinPhi), _e * _B/2.0)) - M_PI_4);
  40         double cosU = cos(U);
  41         double V = _B * (_lambda0 - lambda);
  42         double T = asin(_cosAlphaC * sin(U) + _sinAlphaC * cosU * cos(V));
  43         double D = asin(cosU * sin(V) / cos(T));
  44         double theta = _n * D;
  45         double r = _r0 * pow(tan(M_PI_4 + _phiP/2.0), _n)
  46           / pow(tan(T/2.0 + M_PI_4), _n);
  47
  48         return PointD(r * sin(theta) + _FE, r * cos(theta) + _FN);
  49 }
  50
  51 Coordinates Krovak::xy2ll(const PointD &p) const
  52 {
  53         double Xp = p.y() - _FN;
  54         double Yp = p.x() - _FE;
  55         double Xp2 = Xp * Xp;
  56         double Yp2 = Yp * Yp;
  57         double r = sqrt(Xp2 + Yp2);
  58         double theta = atan(Yp / Xp);
  59         double D = theta / sin(_phiP);
  60         double T = 2.0 * (atan(pow(_r0 / r, 1.0/_n) * tan(M_PI_4 + _phiP/2.0))
  61           - M_PI_4);
  62         double U = asin(_cosAlphaC * sin(T) - _sinAlphaC * cos(T) * cos(D));
  63         double V = asin(cos(T) * sin(D) / cos(U));
  64         double phi = U;
  65         for (int i = 0; i < 3; i++)
  66                 phi = 2.0 * (atan(pow(_t0, -1.0/_B) * pow(tan(U/2.0 + M_PI_4), 1.0/_B)
  67                   * pow((1.0 + _e * sin(phi))/(1.0 - _e * sin(phi)), _e/2.0)) - M_PI_4);
  68
  69         return Coordinates(rad2deg(_lambda0 - V/_B), rad2deg(phi));
  70 }