Converted poly

2 files changed, 146 insertions, 29 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 63e2dfa9..831890fa 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -365,7 +365,6 @@ int pj_lonlat_selftest (void) {return 10000;}
 int pj_longlat_selftest (void) {return 10000;}
 
 int pj_ob_tran_selftest (void) {return 10000;}
-int pj_poly_selftest (void) {return 10000;}
 int pj_putp2_selftest (void) {return 10000;}
 int pj_putp3_selftest (void) {return 10000;}
 int pj_putp3p_selftest (void) {return 10000;}
diff --git a/src/PJ_poly.c b/src/PJ_poly.c
index ef5b8905..50c0ee5f 100644
--- a/src/PJ_poly.c
+++ b/src/PJ_poly.c
@@ -1,42 +1,67 @@
-#define PROJ_PARMS__ \
-	double ml0; \
-	double *en;
 #define PJ_LIB__
 #include <projects.h>
+
 PROJ_HEAD(poly, "Polyconic (American)")
 	"\n\tConic, Sph&Ell";
+
+struct pj_opaque {
+	double ml0; \
+	double *en;
+};
+
 #define TOL	1e-10
 #define CONV	1e-10
 #define N_ITER	10
 #define I_ITER 20
 #define ITOL 1.e-12
-FORWARD(e_forward); /* ellipsoid */
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double  ms, sp, cp;
 
-	if (fabs(lp.phi) <= TOL) { xy.x = lp.lam; xy.y = -P->ml0; }
-	else {
+	if (fabs(lp.phi) <= TOL) {
+        xy.x = lp.lam;
+        xy.y = -Q->ml0;
+    } else {
 		sp = sin(lp.phi);
 		ms = fabs(cp = cos(lp.phi)) > TOL ? pj_msfn(sp, cp, P->es) / sp : 0.;
 		xy.x = ms * sin(lp.lam *= sp);
-		xy.y = (pj_mlfn(lp.phi, sp, cp, P->en) - P->ml0) + ms * (1. - cos(lp.lam));
+		xy.y = (pj_mlfn(lp.phi, sp, cp, Q->en) - Q->ml0) + ms * (1. - cos(lp.lam));
 	}
-	return (xy);
+
+	return xy;
 }
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double  cot, E;
 
-	if (fabs(lp.phi) <= TOL) { xy.x = lp.lam; xy.y = P->ml0; }
-	else {
+	if (fabs(lp.phi) <= TOL) {
+        xy.x = lp.lam;
+        xy.y = Q->ml0;
+    } else {
 		cot = 1. / tan(lp.phi);
 		xy.x = sin(E = lp.lam * sin(lp.phi)) * cot;
 		xy.y = lp.phi - P->phi0 + cot * (1. - cos(E));
 	}
-	return (xy);
+
+	return xy;
 }
-INVERSE(e_inverse); /* ellipsoid */
-	xy.y += P->ml0;
-	if (fabs(xy.y) <= TOL) { lp.lam = xy.x; lp.phi = 0.; }
-	else {
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+	xy.y += Q->ml0;
+	if (fabs(xy.y) <= TOL) {
+        lp.lam = xy.x;
+        lp.phi = 0.;
+    } else {
 		double r, c, sp, cp, s2ph, ml, mlb, mlp, dPhi;
 		int i;
 
@@ -47,7 +72,7 @@ INVERSE(e_inverse); /* ellipsoid */
 			if (fabs(cp) < ITOL)
 				I_ERROR;
 			c = sp * (mlp = sqrt(1. - P->es * sp * sp)) / cp;
-			ml = pj_mlfn(lp.phi, sp, cp, P->en);
+			ml = pj_mlfn(lp.phi, sp, cp, Q->en);
 			mlb = ml * ml + r;
 			mlp = P->one_es / (mlp * mlp * mlp);
 			lp.phi += ( dPhi =
@@ -62,14 +87,20 @@ INVERSE(e_inverse); /* ellipsoid */
 		c = sin(lp.phi);
 		lp.lam = asin(xy.x * tan(lp.phi) * sqrt(1. - P->es * c * c)) / sin(lp.phi);
 	}
-	return (lp);
+
+	return lp;
 }
-INVERSE(s_inverse); /* spheroid */
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
 	double B, dphi, tp;
 	int i;
 
-	if (fabs(xy.y = P->phi0 + xy.y) <= TOL) { lp.lam = xy.x; lp.phi = 0.; }
-	else {
+	if (fabs(xy.y = P->phi0 + xy.y) <= TOL) {
+        lp.lam = xy.x;
+        lp.phi = 0.;
+    } else {
 		lp.phi = xy.y;
 		B = xy.x * xy.x + xy.y * xy.y;
 		i = N_ITER;
@@ -82,18 +113,105 @@ INVERSE(s_inverse); /* spheroid */
 		if (! i) I_ERROR;
 		lp.lam = asin(xy.x * tan(lp.phi)) / sin(lp.phi);
 	}
-	return (lp);
+
+	return lp;
 }
-FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
-ENTRY1(poly, en)
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+    if (P->opaque->en)
+        pj_dealloc (P->opaque->en);
+    pj_dealloc (P->opaque);
+
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(poly) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
 	if (P->es) {
-		if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-		P->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en);
+		if (!(Q->en = pj_enfn(P->es))) E_ERROR_0;
+		Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
 		P->inv = e_inverse;
 		P->fwd = e_forward;
 	} else {
-		P->ml0 = -P->phi0;
+		Q->ml0 = -P->phi0;
 		P->inv = s_inverse;
 		P->fwd = s_forward;
 	}
-ENDENTRY(P)
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_poly_selftest (void) {return 0;}
+#else
+
+int pj_poly_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=poly   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=poly   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222605.285770237475,  110642.194561440483},
+        { 222605.285770237475, -110642.194561440483},
+        {-222605.285770237475,  110642.194561440483},
+        {-222605.285770237475, -110642.194561440483},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223368.105210218986,  111769.110491224754},
+        { 223368.105210218986, -111769.110491224754},
+        {-223368.105210218986,  111769.110491224754},
+        {-223368.105210218986, -111769.110491224754},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.00179663056846135222,  0.000904369476631838518},
+        { 0.00179663056846135222, -0.000904369476631838518},
+        {-0.00179663056846135222,  0.000904369476631838518},
+        {-0.00179663056846135222, -0.000904369476631838518},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931100023887,  0.000895246554454779222},
+        { 0.0017904931100023887, -0.000895246554454779222},
+        {-0.0017904931100023887,  0.000895246554454779222},
+        {-0.0017904931100023887, -0.000895246554454779222},
+    };
+
+    return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
+}
+
+
+#endif