Convert somerc, stere, ups, sterea, kav5, fouc, mbt_s,, qua_aut

i.e. the files
PJ_somerc.c, PJ_stere.c, PJ_sterea.c, PJ_sts.c

47 projections converted so far - 95 to go...

5 files changed, 1038 insertions, 378 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 0822efd6..8cca2386 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -366,7 +366,7 @@ int pj_eqc_selftest (void) {return 10000;}
 int pj_eqdc_selftest (void) {return 10000;}
 int pj_etmerc_selftest (void) {return 10000;}
 int pj_fahey_selftest (void) {return 10000;}
-int pj_fouc_selftest (void) {return 10000;}
+
 int pj_fouc_s_selftest (void) {return 10000;}
 int pj_gall_selftest (void) {return 10000;}
 int pj_geos_selftest (void) {return 10000;}
@@ -384,7 +384,7 @@ int pj_rhealpix_selftest (void) {return 10000;}
 int pj_igh_selftest (void) {return 10000;}
 int pj_imw_p_selftest (void) {return 10000;}
 int pj_isea_selftest (void) {return 10000;}
-int pj_kav5_selftest (void) {return 10000;}
+
 int pj_kav7_selftest (void) {return 10000;}
 int pj_krovak_selftest (void) {return 10000;}
 int pj_labrd_selftest (void) {return 10000;}
@@ -402,7 +402,7 @@ int pj_lcca_selftest (void) {return 10000;}
 int pj_lee_os_selftest (void) {return 10000;}
 int pj_loxim_selftest (void) {return 10000;}
 int pj_lsat_selftest (void) {return 10000;}
-int pj_mbt_s_selftest (void) {return 10000;}
+
 int pj_mbt_fps_selftest (void) {return 10000;}
 int pj_mbtfpp_selftest (void) {return 10000;}
 int pj_mbtfpq_selftest (void) {return 10000;}
@@ -435,19 +435,17 @@ int pj_putp5_selftest (void) {return 10000;}
 int pj_putp5p_selftest (void) {return 10000;}
 int pj_putp6_selftest (void) {return 10000;}
 int pj_putp6p_selftest (void) {return 10000;}
-int pj_qua_aut_selftest (void) {return 10000;}
+
 int pj_qsc_selftest (void) {return 10000;}
 int pj_robin_selftest (void) {return 10000;}
 int pj_rouss_selftest (void) {return 10000;}
 int pj_rpoly_selftest (void) {return 10000;}
 int pj_sch_selftest (void) {return 10000;}
 int pj_sinu_selftest (void) {return 10000;}
-int pj_somerc_selftest (void) {return 10000;}
-int pj_stere_selftest (void) {return 10000;}
-int pj_sterea_selftest (void) {return 10000;}
+
 int pj_gstmerc_selftest (void) {return 10000;}
 int pj_tpers_selftest (void) {return 10000;}
-int pj_ups_selftest (void) {return 10000;}
+
 int pj_utm_selftest (void) {return 10000;}
 int pj_vandg_selftest (void) {return 10000;}
 int pj_vandg2_selftest (void) {return 10000;}
diff --git a/src/PJ_somerc.c b/src/PJ_somerc.c
index f7d1f119..e7659876 100644
--- a/src/PJ_somerc.c
+++ b/src/PJ_somerc.c
@@ -1,66 +1,165 @@
-#define PROJ_PARMS__ \
-	double	K, c, hlf_e, kR, cosp0, sinp0;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
+
 PROJ_HEAD(somerc, "Swiss. Obl. Mercator") "\n\tCyl, Ell\n\tFor CH1903";
-#define EPS	1.e-10
+
+struct pj_opaque {
+    double K, c, hlf_e, kR, cosp0, sinp0;
+};
+
+#define EPS 1.e-10
 #define NITER 6
-FORWARD(e_forward);
-	double phip, lamp, phipp, lampp, sp, cp;
-
-	sp = P->e * sin(lp.phi);
-	phip = 2.* atan( exp( P->c * (
-		log(tan(FORTPI + 0.5 * lp.phi)) - P->hlf_e * log((1. + sp)/(1. - sp)))
-		+ P->K)) - HALFPI;
-	lamp = P->c * lp.lam;
-	cp = cos(phip);
-	phipp = aasin(P->ctx,P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp));
-	lampp = aasin(P->ctx,cp * sin(lamp) / cos(phipp));
-	xy.x = P->kR * lampp;
-	xy.y = P->kR * log(tan(FORTPI + 0.5 * phipp));
-	return (xy);
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0, 0.0};
+    double phip, lamp, phipp, lampp, sp, cp;
+    struct pj_opaque *Q = P->opaque;
+
+    sp = P->e * sin (lp.phi);
+    phip = 2.* atan ( exp ( Q->c * (
+        log (tan (FORTPI + 0.5 * lp.phi)) - Q->hlf_e * log ((1. + sp)/(1. - sp)))
+        + Q->K)) - HALFPI;
+    lamp = Q->c * lp.lam;
+    cp = cos(phip);
+    phipp = aasin (P->ctx, Q->cosp0 * sin (phip) - Q->sinp0 * cp * cos (lamp));
+    lampp = aasin (P->ctx, cp * sin (lamp) / cos (phipp));
+    xy.x = Q->kR * lampp;
+    xy.y = Q->kR * log (tan (FORTPI + 0.5 * phipp));
+    return xy;
 }
-INVERSE(e_inverse); /* ellipsoid & spheroid */
-	double phip, lamp, phipp, lampp, cp, esp, con, delp;
-	int i;
-
-	phipp = 2. * (atan(exp(xy.y / P->kR)) - FORTPI);
-	lampp = xy.x / P->kR;
-	cp = cos(phipp);
-	phip = aasin(P->ctx,P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp));
-	lamp = aasin(P->ctx,cp * sin(lampp) / cos(phip));
-	con = (P->K - log(tan(FORTPI + 0.5 * phip)))/P->c;
-	for (i = NITER; i ; --i) {
-		esp = P->e * sin(phip);
-		delp = (con + log(tan(FORTPI + 0.5 * phip)) - P->hlf_e *
-			log((1. + esp)/(1. - esp)) ) *
-			(1. - esp * esp) * cos(phip) * P->rone_es;
-		phip -= delp;
-		if (fabs(delp) < EPS)
-			break;
-	}
-	if (i) {
-		lp.phi = phip;
-		lp.lam = lamp / P->c;
-	} else
-		I_ERROR
-	return (lp);
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double phip, lamp, phipp, lampp, cp, esp, con, delp;
+    int i;
+
+    phipp = 2. * (atan (exp (xy.y / Q->kR)) - FORTPI);
+    lampp = xy.x / Q->kR;
+    cp = cos (phipp);
+    phip = aasin (P->ctx, Q->cosp0 * sin (phipp) + Q->sinp0 * cp * cos (lampp));
+    lamp = aasin (P->ctx, cp * sin (lampp) / cos (phip));
+    con = (Q->K - log (tan (FORTPI + 0.5 * phip)))/Q->c;
+    for (i = NITER; i ; --i) {
+        esp = P->e * sin(phip);
+        delp = (con + log(tan(FORTPI + 0.5 * phip)) - Q->hlf_e *
+            log((1. + esp)/(1. - esp)) ) *
+            (1. - esp * esp) * cos(phip) * P->rone_es;
+        phip -= delp;
+        if (fabs(delp) < EPS)
+            break;
+    }
+    if (i) {
+        lp.phi = phip;
+        lp.lam = lamp / Q->c;
+    } else
+        I_ERROR
+    return (lp);
 }
+
+
+#if 0
 FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(somerc)
-	double cp, phip0, sp;
-
-	P->hlf_e = 0.5 * P->e;
-	cp = cos(P->phi0);
-	cp *= cp;
-	P->c = sqrt(1 + P->es * cp * cp * P->rone_es);
-	sp = sin(P->phi0);
-	P->cosp0 = cos( phip0 = aasin(P->ctx, P->sinp0 = sp / P->c) );
-	sp *= P->e;
-	P->K = log(tan(FORTPI + 0.5 * phip0)) - P->c * (
-		log(tan(FORTPI + 0.5 * P->phi0)) - P->hlf_e *
-		log((1. + sp) / (1. - sp)));
-	P->kR = P->k0 * sqrt(P->one_es) / (1. - sp * sp);
-	P->inv = e_inverse;
-	P->fwd = e_forward;
-ENDENTRY(P)
+#endif
+
+
+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);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(somerc) {
+    double cp, phip0, sp;
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+
+    Q->hlf_e = 0.5 * P->e;
+    cp = cos (P->phi0);
+    cp *= cp;
+    Q->c = sqrt (1 + P->es * cp * cp * P->rone_es);
+    sp = sin (P->phi0);
+    Q->cosp0 = cos( phip0 = aasin (P->ctx, Q->sinp0 = sp / Q->c) );
+    sp *= P->e;
+    Q->K = log (tan (FORTPI + 0.5 * phip0)) - Q->c * (
+        log (tan (FORTPI + 0.5 * P->phi0)) - Q->hlf_e *
+        log ((1. + sp) / (1. - sp)));
+    Q->kR = P->k0 * sqrt(P->one_es) / (1. - sp * sp);
+    P->inv = e_inverse;
+    P->fwd = e_forward;
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_somerc_selftest (void) {return 0;}
+#else
+
+int pj_somerc_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=somerc   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=somerc   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {222638.98158654713,  110579.96521824898},
+        {222638.98158654713,  -110579.96521825089},
+        {-222638.98158654713,  110579.96521824898},
+        {-222638.98158654713,  -110579.96521825089},
+    };
+
+    XY s_fwd_expect[] = {
+        {223402.14425527418,  111706.74357494408},
+        {223402.14425527418,  -111706.74357494518},
+        {-223402.14425527418,  111706.74357494408},
+        {-223402.14425527418,  -111706.74357494518},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {0.0017966305682390426,  0.00090436947704129484},
+        {0.0017966305682390426,  -0.00090436947704377105},
+        {-0.0017966305682390426,  0.00090436947704129484},
+        {-0.0017966305682390426,  -0.00090436947704377105},
+    };
+
+    LP s_inv_expect[] = {
+        {0.0017904931097838226,  0.00089524655485801927},
+        {0.0017904931097838226,  -0.00089524655484529714},
+        {-0.0017904931097838226,  0.00089524655485801927},
+        {-0.0017904931097838226,  -0.00089524655484529714},
+    };
+
+    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
diff --git a/src/PJ_stere.c b/src/PJ_stere.c
index a0524627..fbca195f 100644
--- a/src/PJ_stere.c
+++ b/src/PJ_stere.c
@@ -1,239 +1,409 @@
-#define PROJ_PARMS__ \
-	double phits; \
-	double sinX1; \
-	double cosX1; \
-	double akm1; \
-	int	mode;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
+
 PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";
-#define sinph0	P->sinX1
-#define cosph0	P->cosX1
-#define EPS10	1.e-10
-#define TOL	1.e-8
-#define NITER	8
-#define CONV	1.e-10
-#define S_POLE	0
-#define N_POLE	1
-#define OBLIQ	2
-#define EQUIT	3
-	static double
-ssfn_(double phit, double sinphi, double eccen) {
-	sinphi *= eccen;
-	return (tan (.5 * (HALFPI + phit)) *
-	   pow((1. - sinphi) / (1. + sinphi), .5 * eccen));
+
+
+struct pj_opaque {
+    double phits;
+    double sinX1;
+    double cosX1;
+    double akm1;
+    int mode;
+};
+
+#define sinph0  P->opaque->sinX1
+#define cosph0  P->opaque->cosX1
+#define EPS10   1.e-10
+#define TOL 1.e-8
+#define NITER   8
+#define CONV    1.e-10
+#define S_POLE  0
+#define N_POLE  1
+#define OBLIQ   2
+#define EQUIT   3
+
+static double ssfn_ (double phit, double sinphi, double eccen) {
+    sinphi *= eccen;
+    return (tan (.5 * (HALFPI + phit)) *
+       pow ((1. - sinphi) / (1. + sinphi), .5 * eccen));
 }
-FORWARD(e_forward); /* ellipsoid */
-	double coslam, sinlam, sinX=0.0, cosX=0.0, X, A, sinphi;
-
-	coslam = cos(lp.lam);
-	sinlam = sin(lp.lam);
-	sinphi = sin(lp.phi);
-	if (P->mode == OBLIQ || P->mode == EQUIT) {
-		sinX = sin(X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - HALFPI);
-		cosX = cos(X);
-	}
-	switch (P->mode) {
-	case OBLIQ:
-		A = P->akm1 / (P->cosX1 * (1. + P->sinX1 * sinX +
-		   P->cosX1 * cosX * coslam));
-		xy.y = A * (P->cosX1 * sinX - P->sinX1 * cosX * coslam);
-		goto xmul;
-	case EQUIT:
-		A = 2. * P->akm1 / (1. + cosX * coslam);
-		xy.y = A * sinX;
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A, sinphi;
+
+    coslam = cos (lp.lam);
+    sinlam = sin (lp.lam);
+    sinphi = sin (lp.phi);
+    if (Q->mode == OBLIQ || Q->mode == EQUIT) {
+        sinX = sin (X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - HALFPI);
+        cosX = cos (X);
+    }
+
+    switch (Q->mode) {
+    case OBLIQ:
+        A = Q->akm1 / (Q->cosX1 * (1. + Q->sinX1 * sinX +
+           Q->cosX1 * cosX * coslam));
+        xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam);
+        goto xmul; /* but why not just  xy.x = A * cosX; break;  ? */
+
+    case EQUIT:
+        A = 2. * Q->akm1 / (1. + cosX * coslam);
+        xy.y = A * sinX;
 xmul:
-		xy.x = A * cosX;
-		break;
-	case S_POLE:
-		lp.phi = -lp.phi;
-		coslam = - coslam;
-		sinphi = -sinphi;
-	case N_POLE:
-		xy.x = P->akm1 * pj_tsfn(lp.phi, sinphi, P->e);
-		xy.y = - xy.x * coslam;
-		break;
-	}
-	xy.x = xy.x * sinlam;
-	return (xy);
+        xy.x = A * cosX;
+        break;
+
+    case S_POLE:
+        lp.phi = -lp.phi;
+        coslam = - coslam;
+        sinphi = -sinphi;
+    case N_POLE:
+        xy.x = Q->akm1 * pj_tsfn (lp.phi, sinphi, P->e);
+        xy.y = - xy.x * coslam;
+        break;
+    }
+
+    xy.x = xy.x * sinlam;
+    return xy;
 }
-FORWARD(s_forward); /* spheroid */
-	double  sinphi, cosphi, coslam, sinlam;
-
-	sinphi = sin(lp.phi);
-	cosphi = cos(lp.phi);
-	coslam = cos(lp.lam);
-	sinlam = sin(lp.lam);
-	switch (P->mode) {
-	case EQUIT:
-		xy.y = 1. + cosphi * coslam;
-		goto oblcon;
-	case OBLIQ:
-		xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double  sinphi, cosphi, coslam, sinlam;
+
+    sinphi = sin(lp.phi);
+    cosphi = cos(lp.phi);
+    coslam = cos(lp.lam);
+    sinlam = sin(lp.lam);
+
+    switch (Q->mode) {
+    case EQUIT:
+        xy.y = 1. + cosphi * coslam;
+        goto oblcon;
+    case OBLIQ:
+        xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
 oblcon:
-		if (xy.y <= EPS10) F_ERROR;
-		xy.x = (xy.y = P->akm1 / xy.y) * cosphi * sinlam;
-		xy.y *= (P->mode == EQUIT) ? sinphi :
-		   cosph0 * sinphi - sinph0 * cosphi * coslam;
-		break;
-	case N_POLE:
-		coslam = - coslam;
-		lp.phi = - lp.phi;
-	case S_POLE:
-		if (fabs(lp.phi - HALFPI) < TOL) F_ERROR;
-		xy.x = sinlam * ( xy.y = P->akm1 * tan(FORTPI + .5 * lp.phi) );
-		xy.y *= coslam;
-		break;
-	}
-	return (xy);
+        if (xy.y <= EPS10) F_ERROR;
+        xy.x = (xy.y = Q->akm1 / xy.y) * cosphi * sinlam;
+        xy.y *= (Q->mode == EQUIT) ? sinphi :
+           cosph0 * sinphi - sinph0 * cosphi * coslam;
+        break;
+    case N_POLE:
+        coslam = - coslam;
+        lp.phi = - lp.phi;
+    case S_POLE:
+        if (fabs (lp.phi - HALFPI) < TOL) F_ERROR;
+        xy.x = sinlam * ( xy.y = Q->akm1 * tan (FORTPI + .5 * lp.phi) );
+        xy.y *= coslam;
+        break;
+    }
+    return xy;
 }
-INVERSE(e_inverse); /* ellipsoid */
-	double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
-	int i;
-
-	rho = hypot(xy.x, xy.y);
-	switch (P->mode) {
-	case OBLIQ:
-	case EQUIT:
-		cosphi = cos( tp = 2. * atan2(rho * P->cosX1 , P->akm1) );
-		sinphi = sin(tp);
-                if( rho == 0.0 )
-		    phi_l = asin(cosphi * P->sinX1);
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
+    int i;
+
+    rho = hypot (xy.x, xy.y);
+
+    switch (Q->mode) {
+    case OBLIQ:
+    case EQUIT:
+        cosphi = cos ( tp = 2. * atan2 (rho * Q->cosX1 , Q->akm1) );
+        sinphi = sin (tp);
+                if ( rho == 0.0 )
+            phi_l = asin (cosphi * Q->sinX1);
                 else
-		    phi_l = asin(cosphi * P->sinX1 + (xy.y * sinphi * P->cosX1 / rho));
-
-		tp = tan(.5 * (HALFPI + phi_l));
-		xy.x *= sinphi;
-		xy.y = rho * P->cosX1 * cosphi - xy.y * P->sinX1* sinphi;
-		halfpi = HALFPI;
-		halfe = .5 * P->e;
-		break;
-	case N_POLE:
-		xy.y = -xy.y;
-	case S_POLE:
-		phi_l = HALFPI - 2. * atan(tp = - rho / P->akm1);
-		halfpi = -HALFPI;
-		halfe = -.5 * P->e;
-		break;
-	}
-	for (i = NITER; i--; phi_l = lp.phi) {
-		sinphi = P->e * sin(phi_l);
-		lp.phi = 2. * atan(tp * pow((1.+sinphi)/(1.-sinphi),
-		   halfe)) - halfpi;
-		if (fabs(phi_l - lp.phi) < CONV) {
-			if (P->mode == S_POLE)
-				lp.phi = -lp.phi;
-			lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);
-			return (lp);
-		}
-	}
-	I_ERROR;
+            phi_l = asin (cosphi * Q->sinX1 + (xy.y * sinphi * Q->cosX1 / rho));
+
+        tp = tan (.5 * (HALFPI + phi_l));
+        xy.x *= sinphi;
+        xy.y = rho * Q->cosX1 * cosphi - xy.y * Q->sinX1* sinphi;
+        halfpi = HALFPI;
+        halfe = .5 * P->e;
+        break;
+    case N_POLE:
+        xy.y = -xy.y;
+    case S_POLE:
+        phi_l = HALFPI - 2. * atan (tp = - rho / Q->akm1);
+        halfpi = -HALFPI;
+        halfe = -.5 * P->e;
+        break;
+    }
+
+    for (i = NITER;  i--;  phi_l = lp.phi) {
+        sinphi = P->e * sin(phi_l);
+        lp.phi = 2. * atan (tp * pow ((1.+sinphi)/(1.-sinphi), halfe)) - halfpi;
+        if (fabs (phi_l - lp.phi) < CONV) {
+            if (Q->mode == S_POLE)
+                lp.phi = -lp.phi;
+            lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y);
+            return lp;
+        }
+    }
+    I_ERROR;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double  c, rh, sinc, cosc;
+
+    sinc = sin (c = 2. * atan ((rh = hypot (xy.x, xy.y)) / Q->akm1));
+    cosc = cos (c);
+    lp.lam = 0.;
+
+    switch (Q->mode) {
+    case EQUIT:
+        if (fabs (rh) <= EPS10)
+            lp.phi = 0.;
+        else
+            lp.phi = asin (xy.y * sinc / rh);
+        if (cosc != 0. || xy.x != 0.)
+            lp.lam = atan2 (xy.x * sinc, cosc * rh);
+        break;
+    case OBLIQ:
+        if (fabs (rh) <= EPS10)
+            lp.phi = P->phi0;
+        else
+            lp.phi = asin (cosc * sinph0 + xy.y * sinc * cosph0 / rh);
+        if ((c = cosc - sinph0 * sin (lp.phi)) != 0. || xy.x != 0.)
+            lp.lam = atan2 (xy.x * sinc * cosph0, c * rh);
+        break;
+    case N_POLE:
+        xy.y = -xy.y;
+    case S_POLE:
+        if (fabs (rh) <= EPS10)
+            lp.phi = P->phi0;
+        else
+            lp.phi = asin (Q->mode == S_POLE ? - cosc : cosc);
+        lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y);
+        break;
+    }
+    return lp;
 }
-INVERSE(s_inverse); /* spheroid */
-	double  c, rh, sinc, cosc;
-
-	sinc = sin(c = 2. * atan((rh = hypot(xy.x, xy.y)) / P->akm1));
-	cosc = cos(c);
-	lp.lam = 0.;
-	switch (P->mode) {
-	case EQUIT:
-		if (fabs(rh) <= EPS10)
-			lp.phi = 0.;
-		else
-			lp.phi = asin(xy.y * sinc / rh);
-		if (cosc != 0. || xy.x != 0.)
-			lp.lam = atan2(xy.x * sinc, cosc * rh);
-		break;
-	case OBLIQ:
-		if (fabs(rh) <= EPS10)
-			lp.phi = P->phi0;
-		else
-			lp.phi = asin(cosc * sinph0 + xy.y * sinc * cosph0 / rh);
-		if ((c = cosc - sinph0 * sin(lp.phi)) != 0. || xy.x != 0.)
-			lp.lam = atan2(xy.x * sinc * cosph0, c * rh);
-		break;
-	case N_POLE:
-		xy.y = -xy.y;
-	case S_POLE:
-		if (fabs(rh) <= EPS10)
-			lp.phi = P->phi0;
-		else
-			lp.phi = asin(P->mode == S_POLE ? - cosc : cosc);
-		lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y);
-		break;
-	}
-	return (lp);
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+static PJ *setup(PJ *P) {                   /* general initialization */
+    double t;
+    struct pj_opaque *Q = P->opaque;
+
+    if (fabs ((t = fabs (P->phi0)) - HALFPI) < EPS10)
+        Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
+    else
+        Q->mode = t > EPS10 ? OBLIQ : EQUIT;
+    Q->phits = fabs (Q->phits);
+
+    if (P->es) {
+        double X;
+
+        switch (Q->mode) {
+        case N_POLE:
+        case S_POLE:
+            if (fabs (Q->phits - HALFPI) < EPS10)
+                Q->akm1 = 2. * P->k0 /
+                   sqrt (pow (1+P->e,1+P->e) * pow (1-P->e,1-P->e));
+            else {
+                Q->akm1 = cos (Q->phits) /
+                   pj_tsfn (Q->phits, t = sin (Q->phits), P->e);
+                t *= P->e;
+                Q->akm1 /= sqrt(1. - t * t);
+            }
+            break;
+        case EQUIT:
+        case OBLIQ:
+            t = sin (P->phi0);
+            X = 2. * atan (ssfn_(P->phi0, t, P->e)) - HALFPI;
+            t *= P->e;
+            Q->akm1 = 2. * P->k0 * cos (P->phi0) / sqrt(1. - t * t);
+            Q->sinX1 = sin (X);
+            Q->cosX1 = cos (X);
+            break;
+        }
+        P->inv = e_inverse;
+        P->fwd = e_forward;
+    } else {
+        switch (Q->mode) {
+        case OBLIQ:
+            sinph0 = sin (P->phi0);
+            cosph0 = cos (P->phi0);
+        case EQUIT:
+            Q->akm1 = 2. * P->k0;
+            break;
+        case S_POLE:
+        case N_POLE:
+            Q->akm1 = fabs (Q->phits - HALFPI) >= EPS10 ?
+               cos (Q->phits) / tan (FORTPI - .5 * Q->phits) :
+               2. * P->k0 ;
+            break;
+        }
+
+        P->inv = s_inverse;
+        P->fwd = s_forward;
+    }
+    return P;
 }
-FREEUP; if (P) pj_dalloc(P); }
-	static PJ *
-setup(PJ *P) { /* general initialization */
-	double t;
-
-	if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
-		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
-	else
-		P->mode = t > EPS10 ? OBLIQ : EQUIT;
-	P->phits = fabs(P->phits);
-	if (P->es) {
-		double X;
-
-		switch (P->mode) {
-		case N_POLE:
-		case S_POLE:
-			if (fabs(P->phits - HALFPI) < EPS10)
-				P->akm1 = 2. * P->k0 /
-				   sqrt(pow(1+P->e,1+P->e)*pow(1-P->e,1-P->e));
-			else {
-				P->akm1 = cos(P->phits) /
-				   pj_tsfn(P->phits, t = sin(P->phits), P->e);
-				t *= P->e;
-				P->akm1 /= sqrt(1. - t * t);
-			}
-			break;
-		case EQUIT:
-		case OBLIQ:
-			t = sin(P->phi0);
-			X = 2. * atan(ssfn_(P->phi0, t, P->e)) - HALFPI;
-			t *= P->e;
-			P->akm1 = 2. * P->k0 * cos(P->phi0) / sqrt(1. - t * t);
-			P->sinX1 = sin(X);
-			P->cosX1 = cos(X);
-			break;
-		}
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
-		switch (P->mode) {
-		case OBLIQ:
-			sinph0 = sin(P->phi0);
-			cosph0 = cos(P->phi0);
-		case EQUIT:
-			P->akm1 = 2. * P->k0;
-			break;
-		case S_POLE:
-		case N_POLE:
-			P->akm1 = fabs(P->phits - HALFPI) >= EPS10 ?
-			   cos(P->phits) / tan(FORTPI - .5 * P->phits) :
-			   2. * P->k0 ;
-			break;
-		}
-		P->inv = s_inverse;
-		P->fwd = s_forward;
-	}
-	return P;
+
+
+PJ *PROJECTION(stere) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->phits = pj_param (P->ctx, P->params, "tlat_ts").i ?
+               pj_param (P->ctx, P->params, "rlat_ts").f : HALFPI;
+
+    return setup(P);
 }
-ENTRY0(stere)
-	P->phits = pj_param(P->ctx, P->params, "tlat_ts").i ?
-            pj_param(P->ctx, P->params, "rlat_ts").f : HALFPI;
-ENDENTRY(setup(P))
-ENTRY0(ups)
+
+
+PJ *PROJECTION(ups) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
 	/* International Ellipsoid */
 	P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - HALFPI: HALFPI;
 	if (!P->es) E_ERROR(-34);
 	P->k0 = .994;
 	P->x0 = 2000000.;
 	P->y0 = 2000000.;
-	P->phits = HALFPI;
+	Q->phits = HALFPI;
 	P->lam0 = 0.;
-ENDENTRY(setup(P))
+
+    return setup(P);
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_stere_selftest (void) {return 0;}
+#else
+
+int pj_stere_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=stere   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=stere   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 445289.70910023432,  221221.76694834774},
+        { 445289.70910023432, -221221.76694835056},
+        {-445289.70910023432,  221221.76694834774},
+        {-445289.70910023432, -221221.76694835056},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223407.81025950745,  111737.938996443},
+        { 223407.81025950745, -111737.938996443},
+        {-223407.81025950745,  111737.938996443},
+        {-223407.81025950745, -111737.938996443},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017966305682022392,  0.00090436947502443507},
+        { 0.0017966305682022392, -0.00090436947502443507},
+        {-0.0017966305682022392,  0.00090436947502443507},
+        {-0.0017966305682022392, -0.00090436947502443507},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.001790493109747395,  0.00089524655465513144},
+        { 0.001790493109747395, -0.00089524655465513144},
+        {-0.001790493109747395,  0.00089524655465513144},
+        {-0.001790493109747395, -0.00089524655465513144},
+    };
+
+    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
+
+
+
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_ups_selftest (void) {return 0;}
+#else
+
+int pj_ups_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=ups   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {2433455.5634384668,  -10412543.301512826},
+        {2448749.1185681992,  -10850493.419804076},
+        {1566544.4365615332,  -10412543.301512826},
+        {1551250.8814318008,  -10850493.419804076},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {-44.998567498074834,  64.9182362867341},
+        {-44.995702709112308,  64.917020250675748},
+        {-45.004297076028529,  64.915804280954518},
+        {-45.001432287066002,  64.914588377560719},
+    };
+
+    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
diff --git a/src/PJ_sterea.c b/src/PJ_sterea.c
index fff2e952..cbfb8fb4 100644
--- a/src/PJ_sterea.c
+++ b/src/PJ_sterea.c
@@ -23,61 +23,155 @@
 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */
-#define PROJ_PARMS__ \
-	double phic0; \
-	double cosc0, sinc0; \
-	double R2; \
-	void *en;
-
 #define PJ_LIB__
-#include	<projects.h>
-
-PROJ_HEAD(sterea, "Oblique Stereographic Alternative")
-	"\n\tAzimuthal, Sph&Ell";
-# define DEL_TOL	1.e-14
-# define MAX_ITER	10
-
-FORWARD(e_forward); /* ellipsoid */
-	double cosc, sinc, cosl, k;
-
-	lp = pj_gauss(P->ctx, lp, P->en);
-	sinc = sin(lp.phi);
-	cosc = cos(lp.phi);
-	cosl = cos(lp.lam);
-	k = P->k0 * P->R2 / (1. + P->sinc0 * sinc + P->cosc0 * cosc * cosl);
-	xy.x = k * cosc * sin(lp.lam);
-	xy.y = k * (P->cosc0 * sinc - P->sinc0 * cosc * cosl);
-	return (xy);
+#include    <projects.h>
+
+
+struct pj_opaque {
+    double phic0;
+    double cosc0, sinc0;
+    double R2;
+    void *en;
+};
+
+
+PROJ_HEAD(sterea, "Oblique Stereographic Alternative") "\n\tAzimuthal, Sph&Ell";
+# define DEL_TOL    1.e-14
+# define MAX_ITER   10
+
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double cosc, sinc, cosl, k;
+
+    lp = pj_gauss(P->ctx, lp, Q->en);
+    sinc = sin(lp.phi);
+    cosc = cos(lp.phi);
+    cosl = cos(lp.lam);
+    k = P->k0 * Q->R2 / (1. + Q->sinc0 * sinc + Q->cosc0 * cosc * cosl);
+    xy.x = k * cosc * sin(lp.lam);
+    xy.y = k * (Q->cosc0 * sinc - Q->sinc0 * cosc * cosl);
+    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 rho, c, sinc, cosc;
+
+    xy.x /= P->k0;
+    xy.y /= P->k0;
+    if ( (rho = hypot (xy.x, xy.y)) ) {
+        c = 2. * atan2 (rho, Q->R2);
+        sinc = sin (c);
+        cosc = cos (c);
+        lp.phi = asin (cosc * Q->sinc0 + xy.y * sinc * Q->cosc0 / rho);
+        lp.lam = atan2 (xy.x * sinc, rho * Q->cosc0 * cosc - xy.y * Q->sinc0 * sinc);
+    } else {
+        lp.phi = Q->phic0;
+        lp.lam = 0.;
+    }
+    return pj_inv_gauss(P->ctx, lp, Q->en);
+}
+
+
+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);
 }
-INVERSE(e_inverse); /* ellipsoid */
-	double rho, c, sinc, cosc;
-
-	xy.x /= P->k0;
-	xy.y /= P->k0;
-	if((rho = hypot(xy.x, xy.y))) {
-		c = 2. * atan2(rho, P->R2);
-		sinc = sin(c);
-		cosc = cos(c);
-		lp.phi = asin(cosc * P->sinc0 + xy.y * sinc * P->cosc0 / rho);
-		lp.lam = atan2(xy.x * sinc, rho * P->cosc0 * cosc -
-			xy.y * P->sinc0 * sinc);
-	} else {
-		lp.phi = P->phic0;
-		lp.lam = 0.;
-	}
-	return(pj_inv_gauss(P->ctx, lp, P->en));
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-FREEUP; if (P) { if (P->en) free(P->en); free(P); } }
-ENTRYA(sterea)
-
-        P->en=0; 
-ENTRYX
-	double R;
-
-	if (!(P->en = pj_gauss_ini(P->e, P->phi0, &(P->phic0), &R))) E_ERROR_0;
-	P->sinc0 = sin(P->phic0);
-	P->cosc0 = cos(P->phic0);
-	P->R2 = 2. * R;
-	P->inv = e_inverse;
-	P->fwd = e_forward;
-ENDENTRY(P)
+
+
+PJ *PROJECTION(sterea) {
+    double R;
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->en = pj_gauss_ini(P->e, P->phi0, &(Q->phic0), &R);
+    if (0==P)
+        E_ERROR_0;
+
+    Q->sinc0 = sin (Q->phic0);
+    Q->cosc0 = cos (Q->phic0);
+    Q->R2 = 2. * R;
+
+    P->inv = e_inverse;
+    P->fwd = e_forward;
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_sterea_selftest (void) {return 0;}
+#else
+
+int pj_sterea_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=sterea   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=sterea   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222644.89410919772,  110611.09187173686},
+        { 222644.89410919772, -110611.09187173827},
+        {-222644.89410919772,  110611.09187173686},
+        {-222644.89410919772, -110611.09187173827},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223407.81025950745,  111737.93899644315},
+        { 223407.81025950745, -111737.93899644315},
+        {-223407.81025950745,  111737.93899644315},
+        {-223407.81025950745, -111737.93899644315},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017966305682019911,  0.00090436947683099009},
+        { 0.0017966305682019911, -0.00090436947684371233},
+        {-0.0017966305682019911,  0.00090436947683099009},
+        {-0.0017966305682019911, -0.00090436947684371233},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.001790493109747395,  0.00089524655465446378},
+        { 0.001790493109747395, -0.00089524655465446378},
+        {-0.001790493109747395,  0.00089524655465446378},
+        {-0.001790493109747395, -0.00089524655465446378},
+    };
+
+    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
diff --git a/src/PJ_sts.c b/src/PJ_sts.c
index c35d5cfc..689c0489 100644
--- a/src/PJ_sts.c
+++ b/src/PJ_sts.c
@@ -1,54 +1,353 @@
-#define PROJ_PARMS__ \
-	double C_x, C_y, C_p; \
-	int tan_mode;
 #define PJ_LIB__
 # include	<projects.h>
-PROJ_HEAD(kav5, "Kavraisky V") "\n\tPCyl., Sph.";
-PROJ_HEAD(qua_aut, "Quartic Authalic") "\n\tPCyl., Sph.";
-PROJ_HEAD(mbt_s, "McBryde-Thomas Flat-Polar Sine (No. 1)") "\n\tPCyl., Sph.";
-PROJ_HEAD(fouc, "Foucaut") "\n\tPCyl., Sph.";
-FORWARD(s_forward); /* spheroid */
+
+PROJ_HEAD(kav5,    "Kavraisky V")         "\n\tPCyl., Sph.";
+PROJ_HEAD(qua_aut, "Quartic Authalic")    "\n\tPCyl., Sph.";
+PROJ_HEAD(fouc,    "Foucaut")             "\n\tPCyl., Sph.";
+PROJ_HEAD(mbt_s,   "McBryde-Thomas Flat-Polar Sine (No. 1)") "\n\tPCyl., Sph.";
+
+
+struct pj_opaque {
+	double C_x, C_y, C_p; \
+	int tan_mode;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double c;
 
-	xy.x = P->C_x * lp.lam * cos(lp.phi);
-	xy.y = P->C_y;
-	lp.phi *= P->C_p;
+	xy.x = Q->C_x * lp.lam * cos(lp.phi);
+	xy.y = Q->C_y;
+	lp.phi *= Q->C_p;
 	c = cos(lp.phi);
-	if (P->tan_mode) {
+	if (Q->tan_mode) {
 		xy.x *= c * c;
-		xy.y *= tan(lp.phi);
+		xy.y *= tan (lp.phi);
 	} else {
 		xy.x /= c;
-		xy.y *= sin(lp.phi);
+		xy.y *= sin (lp.phi);
 	}
-	return (xy);
+	return xy;
 }
-INVERSE(s_inverse); /* spheroid */
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double c;
-	
-	xy.y /= P->C_y;
-	c = cos(lp.phi = P->tan_mode ? atan(xy.y) : aasin(P->ctx,xy.y));
-	lp.phi /= P->C_p;
-	lp.lam = xy.x / (P->C_x * cos(lp.phi));
-	if (P->tan_mode)
+
+	xy.y /= Q->C_y;
+	c = cos (lp.phi = Q->tan_mode ? atan (xy.y) : aasin (P->ctx, xy.y));
+	lp.phi /= Q->C_p;
+	lp.lam = xy.x / (Q->C_x * cos(lp.phi));
+	if (Q->tan_mode)
 		lp.lam /= c * c;
 	else
 		lp.lam *= c;
-	return (lp);
+	return lp;
+}
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-FREEUP; if (P) pj_dalloc(P); }
-	static PJ *
-setup(PJ *P, double p, double q, int mode) {
-	P->es = 0.;
+
+
+static PJ *setup(PJ *P, double p, double q, int mode) {
+	P->es  = 0.;
 	P->inv = s_inverse;
 	P->fwd = s_forward;
-	P->C_x = q / p;
-	P->C_y = p;
-	P->C_p = 1/ q;
-	P->tan_mode = mode;
+	P->opaque->C_x = q / p;
+	P->opaque->C_y = p;
+	P->opaque->C_p = 1/ q;
+	P->opaque->tan_mode = mode;
 	return P;
 }
-ENTRY0(kav5) ENDENTRY(setup(P, 1.50488, 1.35439, 0))
-ENTRY0(qua_aut) ENDENTRY(setup(P, 2., 2., 0))
-ENTRY0(mbt_s) ENDENTRY(setup(P, 1.48875, 1.36509, 0))
-ENTRY0(fouc) ENDENTRY(setup(P, 2., 2., 1))
+
+
+
+
+
+PJ *PROJECTION(fouc) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+    return setup(P, 2., 2., 1);
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_fouc_selftest (void) {return 0;}
+#else
+int pj_fouc_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=fouc   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=fouc   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {222588.12067589167,  111322.31670069379},
+        {222588.12067589167,  -111322.31670069379},
+        {-222588.12067589167,  111322.31670069379},
+        {-222588.12067589167,  -111322.31670069379},
+    };
+
+    XY s_fwd_expect[] = {
+        {223351.10900341379,  111703.9077217125},
+        {223351.10900341379,  -111703.9077217125},
+        {-223351.10900341379,  111703.9077217125},
+        {-223351.10900341379,  -111703.9077217125},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {0.0017966305685702751,  0.00089831528410111959},
+        {0.0017966305685702751,  -0.00089831528410111959},
+        {-0.0017966305685702751,  0.00089831528410111959},
+        {-0.0017966305685702751,  -0.00089831528410111959},
+    };
+
+    LP s_inv_expect[] = {
+        {0.0017904931101116717,  0.00089524655487369749},
+        {0.0017904931101116717,  -0.00089524655487369749},
+        {-0.0017904931101116717,  0.00089524655487369749},
+        {-0.0017904931101116717,  -0.00089524655487369749},
+    };
+
+    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
+
+
+
+
+
+
+PJ *PROJECTION(kav5) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    P->pfree = freeup;
+    P->descr = des_kav5;
+    return setup(P, 1.50488, 1.35439, 0);
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_kav5_selftest (void) {return 0;}
+#else
+int pj_kav5_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=kav5   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=kav5   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {200360.90530882866,  123685.08247699818},
+        {200360.90530882866,  -123685.08247699818},
+        {-200360.90530882866,  123685.08247699818},
+        {-200360.90530882866,  -123685.08247699818},
+    };
+
+    XY s_fwd_expect[] = {
+        {201047.7031108776,  124109.05062917093},
+        {201047.7031108776,  -124109.05062917093},
+        {-201047.7031108776,  124109.05062917093},
+        {-201047.7031108776,  -124109.05062917093},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {0.0019962591348533314,  0.00080848256185253912},
+        {0.0019962591348533314,  -0.00080848256185253912},
+        {-0.0019962591348533314,  0.00080848256185253912},
+        {-0.0019962591348533314,  -0.00080848256185253912},
+    };
+
+    LP s_inv_expect[] = {
+        {0.0019894397264987643,  0.00080572070962591153},
+        {0.0019894397264987643,  -0.00080572070962591153},
+        {-0.0019894397264987643,  0.00080572070962591153},
+        {-0.0019894397264987643,  -0.00080572070962591153},
+    };
+
+    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
+
+
+
+
+
+PJ *PROJECTION(qua_aut) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+    return setup(P, 2., 2., 0);
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_qua_aut_selftest (void) {return 0;}
+#else
+int pj_qua_aut_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=qua_aut   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=qua_aut   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {222613.54903309655,  111318.07788798446},
+        {222613.54903309655,  -111318.07788798446},
+        {-222613.54903309655,  111318.07788798446},
+        {-222613.54903309655,  -111318.07788798446},
+    };
+
+    XY s_fwd_expect[] = {
+        {223376.62452402918,  111699.65437918637},
+        {223376.62452402918,  -111699.65437918637},
+        {-223376.62452402918,  111699.65437918637},
+        {-223376.62452402918,  -111699.65437918637},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {0.0017966305684046586,  0.00089831528412872229},
+        {0.0017966305684046586,  -0.00089831528412872229},
+        {-0.0017966305684046586,  0.00089831528412872229},
+        {-0.0017966305684046586,  -0.00089831528412872229},
+    };
+
+    LP s_inv_expect[] = {
+        {0.0017904931099477471,  0.00089524655490101819},
+        {0.0017904931099477471,  -0.00089524655490101819},
+        {-0.0017904931099477471,  0.00089524655490101819},
+        {-0.0017904931099477471,  -0.00089524655490101819},
+    };
+
+    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
+
+
+
+
+
+PJ *PROJECTION(mbt_s) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+    return setup(P, 1.48875, 1.36509, 0);
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_mbt_s_selftest (void) {return 0;}
+#else
+int pj_mbt_s_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=mbt_s   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5"};
+    char s_args[] = {"+proj=mbt_s   +a=6400000    +lat_1=0.5 +lat_2=2 +n=0.5"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {204131.51785027285,  121400.33022550763},
+        {204131.51785027285,  -121400.33022550763},
+        {-204131.51785027285,  121400.33022550763},
+        {-204131.51785027285,  -121400.33022550763},
+    };
+
+    XY s_fwd_expect[] = {
+        {204831.24057099217,  121816.46669603503},
+        {204831.24057099217,  -121816.46669603503},
+        {-204831.24057099217,  121816.46669603503},
+        {-204831.24057099217,  -121816.46669603503},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {0.0019593827209883237,  0.00082369854658027549},
+        {0.0019593827209883237,  -0.00082369854658027549},
+        {-0.0019593827209883237,  0.00082369854658027549},
+        {-0.0019593827209883237,  -0.00082369854658027549},
+    };
+
+    LP s_inv_expect[] = {
+        {0.0019526892859206603,  0.00082088471512331508},
+        {0.0019526892859206603,  -0.00082088471512331508},
+        {-0.0019526892859206603,  0.00082088471512331508},
+        {-0.0019526892859206603,  -0.00082088471512331508},
+    };
+
+    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