cpp conversion: move source files in apps/ iso19111/ conversions/ projections/ transformations/ tests/ subdirectories

1 files changed, 163 insertions, 0 deletions

diff --git a/src/projections/PJ_calcofi.cpp b/src/projections/PJ_calcofi.cpp
new file mode 100644
index 00000000..e81e4d2a
--- /dev/null
+++ b/src/projections/PJ_calcofi.cpp
@@ -0,0 +1,163 @@
+#define PJ_LIB__
+
+#include <math.h>
+
+#include "proj.h"
+#include "projects.h"
+#include "proj_api.h"
+
+PROJ_HEAD(calcofi,
+    "Cal Coop Ocean Fish Invest Lines/Stations") "\n\tCyl, Sph&Ell";
+
+
+/* Conversions for the California Cooperative Oceanic Fisheries Investigations
+Line/Station coordinate system following the algorithm of:
+Eber, L.E., and  R.P. Hewitt. 1979. Conversion algorithms for the CalCOFI
+station grid. California Cooperative Oceanic Fisheries Investigations Reports
+20:135-137. (corrected for typographical errors).
+http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf
+They assume 1 unit of CalCOFI Line == 1/5 degree in longitude or
+meridional units at reference point O, and similarly 1 unit of CalCOFI
+Station == 1/15 of a degree at O.
+By convention, CalCOFI Line/Station conversions use Clarke 1866 but we use
+whatever ellipsoid is provided. */
+
+
+#define EPS10 1.e-10
+#define DEG_TO_LINE 5
+#define DEG_TO_STATION 15
+#define LINE_TO_RAD 0.0034906585039886592
+#define STATION_TO_RAD 0.0011635528346628863
+#define PT_O_LINE 80 /* reference point O is at line 80,  */
+#define PT_O_STATION 60 /* station 60,  */
+#define PT_O_LAMBDA -2.1144663887911301 /* lon -121.15 and */
+#define PT_O_PHI 0.59602993955606354 /* lat 34.15 */
+#define ROTATION_ANGLE 0.52359877559829882 /*CalCOFI angle of 30 deg in rad */
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    double oy; /* pt O y value in Mercator */
+    double l1; /* l1 and l2 are distances calculated using trig that sum
+               to the east/west distance between point O and point xy */
+    double l2;
+    double ry; /* r is the point on the same station as o (60) and the same
+               line as xy xy, r, o form a right triangle */
+
+    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return xy;
+    }
+
+    xy.x = lp.lam;
+    xy.y = -log(pj_tsfn(lp.phi, sin(lp.phi), P->e)); /* Mercator transform xy*/
+    oy = -log(pj_tsfn(PT_O_PHI, sin(PT_O_PHI), P->e));
+    l1 = (xy.y - oy) * tan(ROTATION_ANGLE);
+    l2 = -xy.x - l1 + PT_O_LAMBDA;
+    ry = l2 * cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE) + xy.y;
+    ry = pj_phi2(P->ctx, exp(-ry), P->e); /*inverse Mercator*/
+    xy.x = PT_O_LINE - RAD_TO_DEG *
+        (ry - PT_O_PHI) * DEG_TO_LINE / cos(ROTATION_ANGLE);
+    xy.y = PT_O_STATION + RAD_TO_DEG *
+        (ry - lp.phi) * DEG_TO_STATION / sin(ROTATION_ANGLE);
+    /* set a = 1, x0 = 0, and y0 = 0 so that no further unit adjustments
+    are done */
+
+    return xy;
+}
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    double oy;
+    double l1;
+    double l2;
+    double ry;
+    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return xy;
+    }
+    xy.x = lp.lam;
+    xy.y = log(tan(M_FORTPI + .5 * lp.phi));
+    oy = log(tan(M_FORTPI + .5 * PT_O_PHI));
+    l1 = (xy.y - oy) * tan(ROTATION_ANGLE);
+    l2 = -xy.x - l1 + PT_O_LAMBDA;
+    ry = l2 * cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE) + xy.y;
+    ry = M_HALFPI - 2. * atan(exp(-ry));
+    xy.x = PT_O_LINE - RAD_TO_DEG *
+        (ry - PT_O_PHI) * DEG_TO_LINE / cos(ROTATION_ANGLE);
+    xy.y = PT_O_STATION + RAD_TO_DEG *
+        (ry - lp.phi) * DEG_TO_STATION / sin(ROTATION_ANGLE);
+
+    return xy;
+}
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    double ry;     /* y value of point r */
+    double oymctr; /* Mercator-transformed y value of point O */
+    double rymctr; /* Mercator-transformed ry */
+    double xymctr; /* Mercator-transformed xy.y */
+    double l1;
+    double l2;
+
+    ry = PT_O_PHI - LINE_TO_RAD * (xy.x - PT_O_LINE) *
+        cos(ROTATION_ANGLE);
+    lp.phi = ry - STATION_TO_RAD * (xy.y - PT_O_STATION) * sin(ROTATION_ANGLE);
+    oymctr = -log(pj_tsfn(PT_O_PHI, sin(PT_O_PHI), P->e));
+    rymctr = -log(pj_tsfn(ry, sin(ry), P->e));
+    xymctr = -log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
+    l1 = (xymctr - oymctr) * tan(ROTATION_ANGLE);
+    l2 = (rymctr - xymctr) / (cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE));
+    lp.lam = PT_O_LAMBDA - (l1 + l2);
+
+    return lp;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    double ry;
+    double oymctr;
+    double rymctr;
+    double xymctr;
+    double l1;
+    double l2;
+    (void) P;
+
+    ry = PT_O_PHI - LINE_TO_RAD * (xy.x - PT_O_LINE) *
+        cos(ROTATION_ANGLE);
+    lp.phi = ry - STATION_TO_RAD * (xy.y - PT_O_STATION) * sin(ROTATION_ANGLE);
+    oymctr = log(tan(M_FORTPI + .5 * PT_O_PHI));
+    rymctr = log(tan(M_FORTPI + .5 * ry));
+    xymctr = log(tan(M_FORTPI + .5 * lp.phi));
+    l1 = (xymctr - oymctr) * tan(ROTATION_ANGLE);
+    l2 = (rymctr - xymctr) / (cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE));
+    lp.lam = PT_O_LAMBDA - (l1 + l2);
+
+    return lp;
+}
+
+
+PJ *PROJECTION(calcofi) {
+    P->opaque = nullptr;
+
+    /* if the user has specified +lon_0 or +k0 for some reason,
+    we're going to ignore it so that xy is consistent with point O */
+    P->lam0 = 0;
+    P->ra = 1;
+    P->a = 1;
+    P->x0 = 0;
+    P->y0 = 0;
+    P->over = 1;
+
+    if (P->es != 0.0) { /* ellipsoid */
+        P->inv = e_inverse;
+        P->fwd = e_forward;
+    } else { /* sphere */
+        P->inv = s_inverse;
+        P->fwd = s_forward;
+    }
+    return P;
+}