Complete overhaul of Krovak projection. Remove unused and duplicate code. Using standard variables from P and added an opaque object to P which has shared values used both in forward and inverse calculations. Added more background information to comments.

1 files changed, 119 insertions, 134 deletions

diff --git a/src/PJ_krovak.c b/src/PJ_krovak.c
index 9ea0239b..4d1498b5 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -1,6 +1,4 @@
-/******************************************************************************
- * $Id$
- *
+ /*
  * Project:  PROJ.4
  * Purpose:  Implementation of the krovak (Krovak) projection.
  *           Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
@@ -27,177 +25,141 @@
  * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  * SOFTWARE.
+ ******************************************************************************
+ * A description of the (forward) projection is found in:
+ *
+ *      Bohuslav Veverka,
+ *
+ *      KROVAK’S PROJECTION AND ITS USE FOR THE
+ *      CZECH REPUBLIC AND THE SLOVAK REPUBLIC,
+ *
+ *      50 years of the Research Institute of
+ *      and the Slovak Republic Geodesy, Topography and Cartography
+ *
+ * which can be found via the Wayback Machine:
+ *
+ *      https://web.archive.org/web/20150216143806/https://www.vugtk.cz/odis/sborniky/sb2005/Sbornik_50_let_VUGTK/Part_1-Scientific_Contribution/16-Veverka.pdf
+ *
+ * Further info, including the inverse projection, is given by EPSG:
+ *
+ *      Guidance Note 7 part 2
+ *      Coordinate Conversions and Transformations including Formulas
+ *
+ *      http://www.iogp.org/pubs/373-07-2.pdf
+ *
+ * Variable names in this file mostly follows what is used in the
+ * paper by Veverka.
+ *
+ * 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
+ * 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
+ *
  *****************************************************************************/
 
+
 #define PJ_LIB__
 #include <projects.h>
-#include <string.h>
-#include <stdio.h>
 
 PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Ellps.";
 
+#define EPS 1e-15
+#define S45 0.785398163397448  /* 45 deg */
+#define S90 1.570796326794896  /* 90 deg */
+#define UQ  1.04216856380474   /* DU(2, 59, 42, 42.69689) */
+#define S0  1.37008346281555   /* Latitude of pseudo standard parallel 78deg 30'00" N */
 
-/**
-   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
-          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
+struct pj_opaque {
+    double alpha;
+    double k;
+    double n;
+    double rho0;
+    double ad;
+    int czech;
+};
 
-  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
-
- **/
-
-
-static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+static XY e_forward (LP lp, PJ *P) {                /* Ellipsoidal, forward */
+    struct pj_opaque *Q = P->opaque;
     XY xy = {0.0,0.0};
 
-    /* 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;    /* 45deg */
-    s90 = 2 * s45;
-    fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
-
-    /* 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;
-    */
-    a =  1; /* 6377397.155; */
-    /* e2 = P->es;*/       /* 0.006674372230614; */
-    e2 = 0.006674372230614;
-    e = sqrt(e2);
+    double gfi, u, deltav, s, d, eps, rho;
 
-    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
+    gfi = pow ( (1. + P->e * sin(lp.phi)) / (1. - P->e * sin(lp.phi)), Q->alpha * P->e / 2.);
 
-    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.  );
+    u = 2. * (atan(Q->k * pow( tan(lp.phi / 2. + S45), Q->alpha) / gfi)-S45);
+    deltav = -lp.lam * Q->alpha;
 
-    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
-
-    k1 = P->k0;
-    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
-    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
-    n = sin(s0);
-    ro0 = k1 * n0 / tan(s0);
-    ad = s90 - uq;
-
-    /* Transformation */
-    gfi = pow ( ((1. + e * sin(lp.phi)) /
-                (1. - e * sin(lp.phi))) , (alfa * e / 2.));
-
-    u = 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
-
-    deltav = - lp.lam * alfa;
-
-    s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
+    s = asin(cos(Q->ad) * sin(u) + sin(Q->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;
+    eps = Q->n * d;
+    rho = Q->rho0 * pow(tan(S0 / 2. + S45) , Q->n) / pow(tan(s / 2. + S45) , Q->n);
+
+    xy.y = rho * cos(eps);
+    xy.x = rho * sin(eps);
 
-    if( !pj_param(P->ctx, P->params, "tczech").i ) {
-        xy.y *= -1.0;
-        xy.x *= -1.0;
-    }
+    xy.y *= Q->czech;
+    xy.x *= Q->czech;
 
     return xy;
 }
 
 
-static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+static LP e_inverse (XY xy, PJ *P) {                /* Ellipsoidal, inverse */
+    struct pj_opaque *Q = P->opaque;
     LP lp = {0.0,0.0};
 
-    /* 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;
+    double u, deltav, s, d, eps, rho, fi1, xy0;
     int ok;
 
-    s45 = 0.785398163397448;    /* 45deg */
-    s90 = 2 * s45;
-    fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
-
-
-   /* 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;
-   */
-    a = 1; /* 6377397.155; */
-    /* e2 = P->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 = P->k0;
-    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
-    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 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(P->ctx, P->params, "tczech").i ) {
-        xy.x *= -1.0;
-        xy.y *= -1.0;
-    }
+    xy.x *= Q->czech;
+    xy.y *= Q->czech;
 
-    ro = sqrt(xy.x * xy.x + xy.y * xy.y);
+    rho = 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));
+    d = eps / sin(S0);
+    s = 2. * (atan(  pow(Q->rho0 / rho, 1. / Q->n) * tan(S0 / 2. + S45)) - S45);
+
+    u = asin(cos(Q->ad) * sin(s) - sin(Q->ad) * cos(s) * cos(d));
     deltav = asin(cos(s) * sin(d) / cos(u));
 
-    lp.lam = P->lam0 - deltav / alfa;
+    lp.lam = P->lam0 - deltav / Q->alpha;
 
     /* ITERATION FOR lp.phi */
     fi1 = u;
 
     ok = 0;
     do {
-        lp.phi = 2. * ( atan( pow( k, -1. / alfa)  *
-                              pow( tan(u / 2. + s45) , 1. / alfa)  *
-                              pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
-                            )  - s45);
+        lp.phi = 2. * ( atan( pow( Q->k, -1. / Q->alpha)  *
+                              pow( tan(u / 2. + S45) , 1. / Q->alpha)  *
+                              pow( (1. + P->e * sin(fi1)) / (1. - P->e * sin(fi1)) , P->e / 2.)
+                            )  - S45);
 
-        if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1;
+        if (fabs(fi1 - lp.phi) < EPS) ok=1;
         fi1 = lp.phi;
    } while (ok==0);
 
@@ -207,10 +169,13 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
 }
 
 
-static void *freeup_new (PJ *P) {                       /* Destructor */
+static void *freeup_new (PJ *P) {                   /* Destructor */
     if (0==P)
         return 0;
+    if (0==P->opaque)
+        return pj_dealloc(P);
 
+    pj_dealloc(P->opaque);
     return pj_dealloc(P);
 }
 
@@ -221,6 +186,12 @@ static void freeup (PJ *P) {
 
 
 PJ *PROJECTION(krovak) {
+    double u0, n0, g;
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
     /* we want Bessel as fixed ellipsoid */
     P->a = 6377397.155;
     P->e = sqrt(P->es = 0.006674372230614);
@@ -239,7 +210,20 @@ PJ *PROJECTION(krovak) {
     if (!pj_param(P->ctx, P->params, "tk").i)
             P->k0 = 0.9999;
 
-    /* always the same */
+    Q->czech = 1;
+    if( !pj_param(P->ctx, P->params, "tczech").i )
+        Q->czech = -1;
+
+    /* Set up shared parameters between forward and inverse */
+    Q->alpha = sqrt(1. + (P->es * pow(cos(P->phi0), 4)) / (1. - P->es));
+    u0 = asin(sin(P->phi0) / Q->alpha);
+    g = pow( (1. + P->e * sin(P->phi0)) / (1. - P->e * sin(P->phi0)) , Q->alpha * P->e / 2. );
+    Q->k = tan( u0 / 2. + S45) / pow  (tan(P->phi0 / 2. + S45) , Q->alpha) * g;
+    n0 = sqrt(1. - P->es) / (1. - P->es * pow(sin(P->phi0), 2));
+    Q->n = sin(S0);
+    Q->rho0 = P->k0 * n0 / tan(S0);
+    Q->ad = S90 - UQ;
+
     P->inv = e_inverse;
     P->fwd = e_forward;
 
@@ -255,7 +239,8 @@ int pj_krovak_selftest (void) {
     double tolerance_lp = 1e-10;
     double tolerance_xy = 1e-7;
 
-    char e_args[] = {"+proj=krovak   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    /* No need to specify an ellipsoid as the projection is hard-coded to Bessel */
+    char e_args[] = {"+proj=krovak"};
 
     LP fwd_in[] = {
         { 2, 1},