Converted ob_tran

2 files changed, 207 insertions, 129 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index afcc2e39..a84a781a 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -361,8 +361,6 @@ int pj_latlong_selftest (void) {return 10000;}
 int pj_lonlat_selftest (void) {return 10000;}
 int pj_longlat_selftest (void) {return 10000;}
 
-int pj_ob_tran_selftest (void) {return 10000;}
-
 int pj_sch_selftest (void) {return 10000;}
 
 int pj_utm_selftest (void) {return 10000;}
diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c
index 4ddba9a4..76222fe9 100644
--- a/src/PJ_ob_tran.c
+++ b/src/PJ_ob_tran.c
@@ -1,145 +1,225 @@
-#define PROJ_PARMS__ \
-	struct PJconsts *link; \
-	double	lamp; \
-	double	cphip, sphip;
 #define PJ_LIB__
 #include <projects.h>
 #include <string.h>
+
+struct pj_opaque {
+    struct PJconsts *link;
+    double  lamp;
+    double  cphip, sphip;
+};
+
 PROJ_HEAD(ob_tran, "General Oblique Transformation") "\n\tMisc Sph"
 "\n\to_proj= plus parameters for projection"
 "\n\to_lat_p= o_lon_p= (new pole) or"
 "\n\to_alpha= o_lon_c= o_lat_c= or"
 "\n\to_lon_1= o_lat_1= o_lon_2= o_lat_2=";
+
 #define TOL 1e-10
-FORWARD(o_forward); /* spheroid */
-	double coslam, sinphi, cosphi;
-
-        (void) xy;
-
-	coslam = cos(lp.lam);
-	sinphi = sin(lp.phi);
-	cosphi = cos(lp.phi);
-	lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), P->sphip * cosphi * coslam +
-		P->cphip * sinphi) + P->lamp);
-	lp.phi = aasin(P->ctx,P->sphip * sinphi - P->cphip * cosphi * coslam);
-	return (P->link->fwd(lp, P->link));
+
+
+static XY o_forward(LP lp, PJ *P) {             /* spheroid */
+    struct pj_opaque *Q = P->opaque;
+    double coslam, sinphi, cosphi;
+
+    coslam = cos(lp.lam);
+    sinphi = sin(lp.phi);
+    cosphi = cos(lp.phi);
+    lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), Q->sphip * cosphi * coslam +
+        Q->cphip * sinphi) + Q->lamp);
+    lp.phi = aasin(P->ctx,Q->sphip * sinphi - Q->cphip * cosphi * coslam);
+
+    return Q->link->fwd(lp, Q->link);
+}
+
+
+static XY t_forward(LP lp, PJ *P) {             /* spheroid */
+    struct pj_opaque *Q = P->opaque;
+    double cosphi, coslam;
+
+    cosphi = cos(lp.phi);
+    coslam = cos(lp.lam);
+    lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), sin(lp.phi)) + Q->lamp);
+    lp.phi = aasin(P->ctx, - cosphi * coslam);
+
+    return Q->link->fwd(lp, Q->link);
 }
-FORWARD(t_forward); /* spheroid */
-	double cosphi, coslam;
 
-        (void) xy;
 
-	cosphi = cos(lp.phi);
-	coslam = cos(lp.lam);
-	lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), sin(lp.phi)) + P->lamp);
-	lp.phi = aasin(P->ctx, - cosphi * coslam);
-	return (P->link->fwd(lp, P->link));
+static LP o_inverse(XY xy, PJ *P) {             /* spheroid */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double coslam, sinphi, cosphi;
+
+    lp = Q->link->inv(xy, Q->link);
+    if (lp.lam != HUGE_VAL) {
+        coslam = cos(lp.lam -= Q->lamp);
+        sinphi = sin(lp.phi);
+        cosphi = cos(lp.phi);
+        lp.phi = aasin(P->ctx,Q->sphip * sinphi + Q->cphip * cosphi * coslam);
+        lp.lam = aatan2(cosphi * sin(lp.lam), Q->sphip * cosphi * coslam -
+            Q->cphip * sinphi);
+    }
+    return lp;
 }
-INVERSE(o_inverse); /* spheroid */
-	double coslam, sinphi, cosphi;
-
-	lp = P->link->inv(xy, P->link);
-	if (lp.lam != HUGE_VAL) {
-		coslam = cos(lp.lam -= P->lamp);
-		sinphi = sin(lp.phi);
-		cosphi = cos(lp.phi);
-		lp.phi = aasin(P->ctx,P->sphip * sinphi + P->cphip * cosphi * coslam);
-		lp.lam = aatan2(cosphi * sin(lp.lam), P->sphip * cosphi * coslam -
-			P->cphip * sinphi);
-	}
-	return (lp);
+
+
+static LP t_inverse(XY xy, PJ *P) {             /* spheroid */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double cosphi, t;
+
+    lp = Q->link->inv(xy, Q->link);
+    if (lp.lam != HUGE_VAL) {
+        cosphi = cos(lp.phi);
+        t = lp.lam - Q->lamp;
+        lp.lam = aatan2(cosphi * sin(t), - sin(lp.phi));
+        lp.phi = aasin(P->ctx,cosphi * cos(t));
+    }
+    return lp;
 }
-INVERSE(t_inverse); /* spheroid */
-	double cosphi, t;
-
-	lp = P->link->inv(xy, P->link);
-	if (lp.lam != HUGE_VAL) {
-		cosphi = cos(lp.phi);
-		t = lp.lam - P->lamp;
-		lp.lam = aatan2(cosphi * sin(t), - sin(lp.phi));
-		lp.phi = aasin(P->ctx,cosphi * cos(t));
-	}
-	return (lp);
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    if (P->opaque->link)
+        return pj_dealloc (P->opaque->link);
+
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
 }
-FREEUP;
-	if (P) {
-		if (P->link)
-			(*(P->link->pfree))(P->link);
-		pj_dalloc(P);
-	}
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-ENTRY1(ob_tran, link)
-	int i;
-	double phip;
-	char *name, *s;
-
-	/* get name of projection to be translated */
-	if (!(name = pj_param(P->ctx, P->params, "so_proj").s)) E_ERROR(-26);
-	for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ;
-	if (!s || !(P->link = (*pj_list[i].proj)(0))) E_ERROR(-37);
-	/* copy existing header into new */
-	P->es = 0.; /* force to spherical */
-	P->link->params = P->params;
-        P->link->ctx = P->ctx;
-	P->link->over = P->over;
-	P->link->geoc = P->geoc;
-	P->link->a = P->a;
-	P->link->es = P->es;
-	P->link->ra = P->ra;
-	P->link->lam0 = P->lam0;
-	P->link->phi0 = P->phi0;
-	P->link->x0 = P->x0;
-	P->link->y0 = P->y0;
-	P->link->k0 = P->k0;
-	/* force spherical earth */
-	P->link->one_es = P->link->rone_es = 1.;
-	P->link->es = P->link->e = 0.;
-	if (!(P->link = pj_list[i].proj(P->link))) {
-		freeup(P);
-		return 0;
-	}
-	if (pj_param(P->ctx, P->params, "to_alpha").i) {
-		double lamc, phic, alpha;
-
-		lamc	= pj_param(P->ctx, P->params, "ro_lon_c").f;
-		phic	= pj_param(P->ctx, P->params, "ro_lat_c").f;
-		alpha	= pj_param(P->ctx, P->params, "ro_alpha").f;
+
+
+PJ *PROJECTION(ob_tran) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    int i;
+    double phip;
+    char *name, *s;
+
+    /* get name of projection to be translated */
+    if (!(name = pj_param(P->ctx, P->params, "so_proj").s)) E_ERROR(-26);
+    for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ;
+    if (!s || !(Q->link = (*pj_list[i].proj)(0))) E_ERROR(-37);
+    /* copy existing header into new */
+    P->es = 0.; /* force to spherical */
+    Q->link->params = P->params;
+    Q->link->ctx = P->ctx;
+    Q->link->over = P->over;
+    Q->link->geoc = P->geoc;
+    Q->link->a = P->a;
+    Q->link->es = P->es;
+    Q->link->ra = P->ra;
+    Q->link->lam0 = P->lam0;
+    Q->link->phi0 = P->phi0;
+    Q->link->x0 = P->x0;
+    Q->link->y0 = P->y0;
+    Q->link->k0 = P->k0;
+    /* force spherical earth */
+    Q->link->one_es = Q->link->rone_es = 1.;
+    Q->link->es = Q->link->e = 0.;
+    if (!(Q->link = pj_list[i].proj(Q->link))) {
+        return freeup_new(P);
+    }
+    if (pj_param(P->ctx, P->params, "to_alpha").i) {
+        double lamc, phic, alpha;
+
+        lamc    = pj_param(P->ctx, P->params, "ro_lon_c").f;
+        phic    = pj_param(P->ctx, P->params, "ro_lat_c").f;
+        alpha   = pj_param(P->ctx, P->params, "ro_alpha").f;
 /*
-		if (fabs(phic) <= TOL ||
-			fabs(fabs(phic) - HALFPI) <= TOL ||
-			fabs(fabs(alpha) - HALFPI) <= TOL)
+        if (fabs(phic) <= TOL ||
+            fabs(fabs(phic) - HALFPI) <= TOL ||
+            fabs(fabs(alpha) - HALFPI) <= TOL)
 */
-		if (fabs(fabs(phic) - HALFPI) <= TOL)
-			E_ERROR(-32);
-		P->lamp = lamc + aatan2(-cos(alpha), -sin(alpha) * sin(phic));
-		phip = aasin(P->ctx,cos(phic) * sin(alpha));
-	} else if (pj_param(P->ctx, P->params, "to_lat_p").i) { /* specified new pole */
-		P->lamp = pj_param(P->ctx, P->params, "ro_lon_p").f;
-		phip = pj_param(P->ctx, P->params, "ro_lat_p").f;
-	} else { /* specified new "equator" points */
-		double lam1, lam2, phi1, phi2, con;
-
-		lam1 = pj_param(P->ctx, P->params, "ro_lon_1").f;
-		phi1 = pj_param(P->ctx, P->params, "ro_lat_1").f;
-		lam2 = pj_param(P->ctx, P->params, "ro_lon_2").f;
-		phi2 = pj_param(P->ctx, P->params, "ro_lat_2").f;
-		if (fabs(phi1 - phi2) <= TOL ||
-			(con = fabs(phi1)) <= TOL ||
-			fabs(con - HALFPI) <= TOL ||
-			fabs(fabs(phi2) - HALFPI) <= TOL) E_ERROR(-33);
-		P->lamp = atan2(cos(phi1) * sin(phi2) * cos(lam1) -
-			sin(phi1) * cos(phi2) * cos(lam2),
-			sin(phi1) * cos(phi2) * sin(lam2) -
-			cos(phi1) * sin(phi2) * sin(lam1));
-		phip = atan(-cos(P->lamp - lam1) / tan(phi1));
-	}
-	if (fabs(phip) > TOL) { /* oblique */
-		P->cphip = cos(phip);
-		P->sphip = sin(phip);
-		P->fwd = o_forward;
-		P->inv = P->link->inv ? o_inverse : 0;
-	} else { /* transverse */
-		P->fwd = t_forward;
-		P->inv = P->link->inv ? t_inverse : 0;
-	}
-ENDENTRY(P)
+        if (fabs(fabs(phic) - HALFPI) <= TOL)
+            E_ERROR(-32);
+        Q->lamp = lamc + aatan2(-cos(alpha), -sin(alpha) * sin(phic));
+        phip = aasin(P->ctx,cos(phic) * sin(alpha));
+    } else if (pj_param(P->ctx, P->params, "to_lat_p").i) { /* specified new pole */
+        Q->lamp = pj_param(P->ctx, P->params, "ro_lon_p").f;
+        phip = pj_param(P->ctx, P->params, "ro_lat_p").f;
+    } else { /* specified new "equator" points */
+        double lam1, lam2, phi1, phi2, con;
+
+        lam1 = pj_param(P->ctx, P->params, "ro_lon_1").f;
+        phi1 = pj_param(P->ctx, P->params, "ro_lat_1").f;
+        lam2 = pj_param(P->ctx, P->params, "ro_lon_2").f;
+        phi2 = pj_param(P->ctx, P->params, "ro_lat_2").f;
+        if (fabs(phi1 - phi2) <= TOL ||
+            (con = fabs(phi1)) <= TOL ||
+            fabs(con - HALFPI) <= TOL ||
+            fabs(fabs(phi2) - HALFPI) <= TOL) E_ERROR(-33);
+        Q->lamp = atan2(cos(phi1) * sin(phi2) * cos(lam1) -
+            sin(phi1) * cos(phi2) * cos(lam2),
+            sin(phi1) * cos(phi2) * sin(lam2) -
+            cos(phi1) * sin(phi2) * sin(lam1));
+        phip = atan(-cos(Q->lamp - lam1) / tan(phi1));
+    }
+    if (fabs(phip) > TOL) { /* oblique */
+        Q->cphip = cos(phip);
+        Q->sphip = sin(phip);
+        P->fwd = o_forward;
+        P->inv = Q->link->inv ? o_inverse : 0;
+    } else { /* transverse */
+        P->fwd = t_forward;
+        P->inv = Q->link->inv ? t_inverse : 0;
+    }
+
+    return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_ob_tran_selftest (void) {return 0;}
+#else
+
+int pj_ob_tran_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=ob_tran +o_proj=latlon +o_lon_p=20 +o_lat_p=20 +lon_0=180"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        {-2.6856872138416592, 1.2374302350496296},
+        {-2.6954069748943286, 1.2026833954513816},
+        {-2.8993663925401947, 1.2374302350496296},
+        {-2.8896466314875244, 1.2026833954513816},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 121.5518748407577,  -2.5361001573966084},
+        { 63.261184340201858,  17.585319578673531},
+        {-141.10073322351622,  26.091712304855108},
+        {-65.862385598848391,  51.830295078417215},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+#endif