Converted lcca

2 files changed, 146 insertions, 61 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index fb26487c..78021a0e 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -366,7 +366,6 @@ int pj_latlon_selftest (void) {return 10000;}
 int pj_latlong_selftest (void) {return 10000;}
 int pj_lonlat_selftest (void) {return 10000;}
 int pj_longlat_selftest (void) {return 10000;}
-int pj_lcca_selftest (void) {return 10000;}
 
 int pj_lee_os_selftest (void) {return 10000;}
 int pj_loxim_selftest (void) {return 10000;}
diff --git a/src/PJ_lcca.c b/src/PJ_lcca.c
index 320d52db..eef46110 100644
--- a/src/PJ_lcca.c
+++ b/src/PJ_lcca.c
@@ -1,70 +1,156 @@
-/* PROJ.4 Cartographic Projection System 
+/* PROJ.4 Cartographic Projection System
 */
-#define MAX_ITER 10
-#define DEL_TOL 1e-12
-#define PROJ_PARMS__ \
-	double	*en; \
-	double	r0, l, M0; \
-	double	C;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
 
 PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative")
-	"\n\tConic, Sph&Ell\n\tlat_0=";
+    "\n\tConic, Sph&Ell\n\tlat_0=";
+
+#define MAX_ITER 10
+#define DEL_TOL 1e-12
+
+struct pj_opaque {
+    double *en;
+    double r0, l, M0;
+    double C;
+};
+
+
+static double fS(double S, double C) {        /* func to compute dr */
+
+    return S * ( 1. + S * S * C);
+}
+
+
+static double fSp(double S, double C) {       /* deriv of fs */
+
+    return 1. + 3.* S * S * C;
+}
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double S, r, dr;
+
+    S = pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), Q->en) - Q->M0;
+    dr = fS(S, Q->C);
+    r = Q->r0 - dr;
+    xy.x = P->k0 * (r * sin( lp.lam *= Q->l ) );
+    xy.y = P->k0 * (Q->r0 - r * cos(lp.lam) );
+    return xy;
+}
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double theta, dr, S, dif;
+    int i;
+
+    xy.x /= P->k0;
+    xy.y /= P->k0;
+    theta = atan2(xy.x , Q->r0 - xy.y);
+    dr = xy.y - xy.x * tan(0.5 * theta);
+    lp.lam = theta / Q->l;
+    S = dr;
+    for (i = MAX_ITER; i ; --i) {
+        S -= (dif = (fS(S, Q->C) - dr) / fSp(S, Q->C));
+        if (fabs(dif) < DEL_TOL) break;
+    }
+    if (!i) I_ERROR
+    lp.phi = pj_inv_mlfn(P->ctx, S + Q->M0, P->es, Q->en);
 
-	static double /* func to compute dr */
-fS(double S, double C) {
-		return(S * ( 1. + S * S * C));
+    return lp;
 }
-	static double /* deriv of fs */
-fSp(double S, double C) {
-	return(1. + 3.* S * S * C);
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    pj_dealloc (P->opaque->en);
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
 }
-FORWARD(e_forward); /* ellipsoid */
-	double S, r, dr;
-	
-	S = pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), P->en) - P->M0;
-	dr = fS(S, P->C);
-	r = P->r0 - dr;
-	xy.x = P->k0 * (r * sin( lp.lam *= P->l ) );
-	xy.y = P->k0 * (P->r0 - r * cos(lp.lam) );
-	return (xy);
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-INVERSE(e_inverse); /* ellipsoid & spheroid */
-	double theta, dr, S, dif;
-	int i;
-
-	xy.x /= P->k0;
-	xy.y /= P->k0;
-	theta = atan2(xy.x , P->r0 - xy.y);
-	dr = xy.y - xy.x * tan(0.5 * theta);
-	lp.lam = theta / P->l;
-	S = dr;
-	for (i = MAX_ITER; i ; --i) {
-		S -= (dif = (fS(S, P->C) - dr) / fSp(S, P->C));
-		if (fabs(dif) < DEL_TOL) break;
-	}
-	if (!i) I_ERROR
-	lp.phi = pj_inv_mlfn(P->ctx, S + P->M0, P->es, P->en);
-	return (lp);
+
+
+PJ *PROJECTION(lcca) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+    double s2p0, N0, R0, tan0;
+
+    (Q->en = pj_enfn(P->es));
+    if (!Q->en) E_ERROR_0;
+    if (!pj_param(P->ctx, P->params, "tlat_0").i) E_ERROR(50);
+    if (P->phi0 == 0.) E_ERROR(51);
+    Q->l = sin(P->phi0);
+    Q->M0 = pj_mlfn(P->phi0, Q->l, cos(P->phi0), Q->en);
+    s2p0 = Q->l * Q->l;
+    R0 = 1. / (1. - P->es * s2p0);
+    N0 = sqrt(R0);
+    R0 *= P->one_es * N0;
+    tan0 = tan(P->phi0);
+    Q->r0 = N0 / tan0;
+    Q->C = 1. / (6. * R0 * N0);
+
+    P->inv = e_inverse;
+    P->fwd = e_forward;
+
+    return P;
 }
-FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
-ENTRY0(lcca)
-	double s2p0, N0, R0, tan0;
-
-	if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-	if (!pj_param(P->ctx, P->params, "tlat_0").i) E_ERROR(50);
-	if (P->phi0 == 0.) E_ERROR(51);
-	P->l = sin(P->phi0);
-	P->M0 = pj_mlfn(P->phi0, P->l, cos(P->phi0), P->en);
-	s2p0 = P->l * P->l;
-	R0 = 1. / (1. - P->es * s2p0);
-	N0 = sqrt(R0);
-	R0 *= P->one_es * N0;
-	tan0 = tan(P->phi0);
-	P->r0 = N0 / tan0;
-	P->C = 1. / (6. * R0 * N0);
-	P->inv = e_inverse;
-	P->fwd = e_forward;
-ENDENTRY(P)
 
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_lcca_selftest (void) {return 0;}
+#else
+
+int pj_lcca_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=lcca   +ellps=GRS80  +lat_0=1 +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222605.285770237417,  67.8060072715846616},
+        { 222740.037637936533, -221125.539829601563},
+        {-222605.285770237417,  67.8060072715846616},
+        {-222740.037637936533, -221125.539829601563},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.00179690290525662526, 1.00090436621350798},
+        { 0.00179690192174008037, 0.999095632791497268},
+        {-0.00179690290525662526, 1.00090436621350798},
+        {-0.00179690192174008037, 0.999095632791497268},
+    };
+
+    return pj_generic_selftest (e_args, 0, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, 0, inv_in, e_inv_expect, 0);
+}
+
+
+#endif