Converted qsc

2 files changed, 383 insertions, 294 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 931a853c..e3fd8ae9 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -366,7 +366,6 @@ int pj_longlat_selftest (void) {return 10000;}
 
 int pj_ob_tran_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;}
diff --git a/src/PJ_qsc.c b/src/PJ_qsc.c
index 12cb9d63..b9e9eb0d 100644
--- a/src/PJ_qsc.c
+++ b/src/PJ_qsc.c
@@ -28,7 +28,7 @@
  * Other projection centers will not work!
  *
  * In the projection code below, each cube face is handled differently.
- * See the computation of the face parameter in the ENTRY0(qsc) function
+ * See the computation of the face parameter in the PROJECTION(qsc) function
  * and the handling of different face values (FACE_*) in the forward and
  * inverse projections.
  *
@@ -38,15 +38,18 @@
  * three areas of a cube face are handled by rotation of AREA_0.
  */
 
-#define PROJ_PARMS__ \
-        int face; \
-        double a_squared; \
-        double b; \
-        double one_minus_f; \
-        double one_minus_f_squared;
 #define PJ_LIB__
-#include        <projects.h>
+#include <projects.h>
+
+struct pj_opaque {
+        int face;
+        double a_squared;
+        double b;
+        double one_minus_f;
+        double one_minus_f_squared;
+};
 PROJ_HEAD(qsc, "Quadrilateralized Spherical Cube") "\n\tAzi, Sph.";
+
 #define EPS10 1.e-10
 
 /* The six cube faces. */
@@ -66,316 +69,403 @@ PROJ_HEAD(qsc, "Quadrilateralized Spherical Cube") "\n\tAzi, Sph.";
 
 /* Helper function for forward projection: compute the theta angle
  * and determine the area number. */
-static double
-qsc_fwd_equat_face_theta(double phi, double y, double x, int *area) {
-        double theta;
-        if (phi < EPS10) {
+static double qsc_fwd_equat_face_theta(double phi, double y, double x, int *area) {
+    double theta;
+    if (phi < EPS10) {
+        *area = AREA_0;
+        theta = 0.0;
+    } else {
+        theta = atan2(y, x);
+        if (fabs(theta) <= FORTPI) {
             *area = AREA_0;
-            theta = 0.0;
+        } else if (theta > FORTPI && theta <= HALFPI + FORTPI) {
+            *area = AREA_1;
+            theta -= HALFPI;
+        } else if (theta > HALFPI + FORTPI || theta <= -(HALFPI + FORTPI)) {
+            *area = AREA_2;
+            theta = (theta >= 0.0 ? theta - PI : theta + PI);
         } else {
-            theta = atan2(y, x);
-            if (fabs(theta) <= FORTPI) {
-                *area = AREA_0;
-            } else if (theta > FORTPI && theta <= HALFPI + FORTPI) {
-                *area = AREA_1;
-                theta -= HALFPI;
-            } else if (theta > HALFPI + FORTPI || theta <= -(HALFPI + FORTPI)) {
-                *area = AREA_2;
-                theta = (theta >= 0.0 ? theta - PI : theta + PI);
-            } else {
-                *area = AREA_3;
-                theta += HALFPI;
-            }
+            *area = AREA_3;
+            theta += HALFPI;
         }
-        return (theta);
+    }
+    return theta;
 }
 
 /* Helper function: shift the longitude. */
-static double
-qsc_shift_lon_origin(double lon, double offset) {
-        double slon = lon + offset;
-        if (slon < -PI) {
-            slon += TWOPI;
-        } else if (slon > +PI) {
-            slon -= TWOPI;
-        }
-        return slon;
+static double qsc_shift_lon_origin(double lon, double offset) {
+    double slon = lon + offset;
+    if (slon < -PI) {
+        slon += TWOPI;
+    } else if (slon > +PI) {
+        slon -= TWOPI;
+    }
+    return slon;
 }
 
-/* Forward projection, ellipsoid */
-FORWARD(e_forward);
-        double lat, lon;
-        double theta, phi;
-        double t, mu; /* nu; */
-        int area;
-
-        /* Convert the geodetic latitude to a geocentric latitude.
-         * This corresponds to the shift from the ellipsoid to the sphere
-         * described in [LK12]. */
-        if (P->es) {
-            lat = atan(P->one_minus_f_squared * tan(lp.phi));
-        } else {
-            lat = lp.phi;
-        }
 
-        /* Convert the input lat, lon into theta, phi as used by QSC.
-         * This depends on the cube face and the area on it.
-         * For the top and bottom face, we can compute theta and phi
-         * directly from phi, lam. For the other faces, we must use
-         * unit sphere cartesian coordinates as an intermediate step. */
-        lon = lp.lam;
-        if (P->face == FACE_TOP) {
-            phi = HALFPI - lat;
-            if (lon >= FORTPI && lon <= HALFPI + FORTPI) {
-                area = AREA_0;
-                theta = lon - HALFPI;
-            } else if (lon > HALFPI + FORTPI || lon <= -(HALFPI + FORTPI)) {
-                area = AREA_1;
-                theta = (lon > 0.0 ? lon - PI : lon + PI);
-            } else if (lon > -(HALFPI + FORTPI) && lon <= -FORTPI) {
-                area = AREA_2;
-                theta = lon + HALFPI;
-            } else {
-                area = AREA_3;
-                theta = lon;
-            }
-        } else if (P->face == FACE_BOTTOM) {
-            phi = HALFPI + lat;
-            if (lon >= FORTPI && lon <= HALFPI + FORTPI) {
-                area = AREA_0;
-                theta = -lon + HALFPI;
-            } else if (lon < FORTPI && lon >= -FORTPI) {
-                area = AREA_1;
-                theta = -lon;
-            } else if (lon < -FORTPI && lon >= -(HALFPI + FORTPI)) {
-                area = AREA_2;
-                theta = -lon - HALFPI;
-            } else {
-                area = AREA_3;
-                theta = (lon > 0.0 ? -lon + PI : -lon - PI);
-            }
-        } else {
-            double q, r, s;
-            double sinlat, coslat;
-            double sinlon, coslon;
-
-            if (P->face == FACE_RIGHT) {
-                lon = qsc_shift_lon_origin(lon, +HALFPI);
-            } else if (P->face == FACE_BACK) {
-                lon = qsc_shift_lon_origin(lon, +PI);
-            } else if (P->face == FACE_LEFT) {
-                lon = qsc_shift_lon_origin(lon, -HALFPI);
-            }
-            sinlat = sin(lat);
-            coslat = cos(lat);
-            sinlon = sin(lon);
-            coslon = cos(lon);
-            q = coslat * coslon;
-            r = coslat * sinlon;
-            s = sinlat;
-
-            if (P->face == FACE_FRONT) {
-                phi = acos(q);
-                theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
-            } else if (P->face == FACE_RIGHT) {
-                phi = acos(r);
-                theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
-            } else if (P->face == FACE_BACK) {
-                phi = acos(-q);
-                theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
-            } else if (P->face == FACE_LEFT) {
-                phi = acos(-r);
-                theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
-            } else {
-                /* Impossible */
-                phi = theta = 0.0;
-                area = AREA_0;
-            }
-        }
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double lat, lon;
+    double theta, phi;
+    double t, mu; /* nu; */
+    int area;
 
-        /* Compute mu and nu for the area of definition.
-         * For mu, see Eq. (3-21) in [OL76], but note the typos:
-         * compare with Eq. (3-14). For nu, see Eq. (3-38). */
-        mu = atan((12.0 / PI) * (theta + acos(sin(theta) * cos(FORTPI)) - HALFPI));
-        t = sqrt((1.0 - cos(phi)) / (cos(mu) * cos(mu)) / (1.0 - cos(atan(1.0 / cos(theta)))));
-        /* nu = atan(t);        We don't really need nu, just t, see below. */
+    /* Convert the geodetic latitude to a geocentric latitude.
+     * This corresponds to the shift from the ellipsoid to the sphere
+     * described in [LK12]. */
+    if (P->es) {
+        lat = atan(Q->one_minus_f_squared * tan(lp.phi));
+    } else {
+        lat = lp.phi;
+    }
 
-        /* Apply the result to the real area. */
-        if (area == AREA_1) {
-            mu += HALFPI;
-        } else if (area == AREA_2) {
-            mu += PI;
-        } else if (area == AREA_3) {
-            mu += HALFPI + PI;
+    /* Convert the input lat, lon into theta, phi as used by QSC.
+     * This depends on the cube face and the area on it.
+     * For the top and bottom face, we can compute theta and phi
+     * directly from phi, lam. For the other faces, we must use
+     * unit sphere cartesian coordinates as an intermediate step. */
+    lon = lp.lam;
+    if (Q->face == FACE_TOP) {
+        phi = HALFPI - lat;
+        if (lon >= FORTPI && lon <= HALFPI + FORTPI) {
+            area = AREA_0;
+            theta = lon - HALFPI;
+        } else if (lon > HALFPI + FORTPI || lon <= -(HALFPI + FORTPI)) {
+            area = AREA_1;
+            theta = (lon > 0.0 ? lon - PI : lon + PI);
+        } else if (lon > -(HALFPI + FORTPI) && lon <= -FORTPI) {
+            area = AREA_2;
+            theta = lon + HALFPI;
+        } else {
+            area = AREA_3;
+            theta = lon;
         }
-
-        /* Now compute x, y from mu and nu */
-        /* t = tan(nu); */
-        xy.x = t * cos(mu);
-        xy.y = t * sin(mu);
-        return (xy);
-}
-
-/* Inverse projection, ellipsoid */
-INVERSE(e_inverse);
-        double mu, nu, cosmu, tannu;
-        double tantheta, theta, cosphi, phi;
-        double t;
-        int area;
-
-        /* Convert the input x, y to the mu and nu angles as used by QSC.
-         * This depends on the area of the cube face. */
-        nu = atan(sqrt(xy.x * xy.x + xy.y * xy.y));
-        mu = atan2(xy.y, xy.x);
-        if (xy.x >= 0.0 && xy.x >= fabs(xy.y)) {
+    } else if (Q->face == FACE_BOTTOM) {
+        phi = HALFPI + lat;
+        if (lon >= FORTPI && lon <= HALFPI + FORTPI) {
             area = AREA_0;
-        } else if (xy.y >= 0.0 && xy.y >= fabs(xy.x)) {
+            theta = -lon + HALFPI;
+        } else if (lon < FORTPI && lon >= -FORTPI) {
             area = AREA_1;
-            mu -= HALFPI;
-        } else if (xy.x < 0.0 && -xy.x >= fabs(xy.y)) {
+            theta = -lon;
+        } else if (lon < -FORTPI && lon >= -(HALFPI + FORTPI)) {
             area = AREA_2;
-            mu = (mu < 0.0 ? mu + PI : mu - PI);
+            theta = -lon - HALFPI;
         } else {
             area = AREA_3;
-            mu += HALFPI;
+            theta = (lon > 0.0 ? -lon + PI : -lon - PI);
         }
+    } else {
+        double q, r, s;
+        double sinlat, coslat;
+        double sinlon, coslon;
 
-        /* Compute phi and theta for the area of definition.
-         * The inverse projection is not described in the original paper, but some
-         * good hints can be found here (as of 2011-12-14):
-         * http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
-         * (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
-        t = (PI / 12.0) * tan(mu);
-        tantheta = sin(t) / (cos(t) - (1.0 / sqrt(2.0)));
-        theta = atan(tantheta);
-        cosmu = cos(mu);
-        tannu = tan(nu);
-        cosphi = 1.0 - cosmu * cosmu * tannu * tannu * (1.0 - cos(atan(1.0 / cos(theta))));
-        if (cosphi < -1.0) {
-            cosphi = -1.0;
-        } else if (cosphi > +1.0) {
-            cosphi = +1.0;
+        if (Q->face == FACE_RIGHT) {
+            lon = qsc_shift_lon_origin(lon, +HALFPI);
+        } else if (Q->face == FACE_BACK) {
+            lon = qsc_shift_lon_origin(lon, +PI);
+        } else if (Q->face == FACE_LEFT) {
+            lon = qsc_shift_lon_origin(lon, -HALFPI);
         }
+        sinlat = sin(lat);
+        coslat = cos(lat);
+        sinlon = sin(lon);
+        coslon = cos(lon);
+        q = coslat * coslon;
+        r = coslat * sinlon;
+        s = sinlat;
 
-        /* Apply the result to the real area on the cube face.
-         * For the top and bottom face, we can compute phi and lam directly.
-         * For the other faces, we must use unit sphere cartesian coordinates
-         * as an intermediate step. */
-        if (P->face == FACE_TOP) {
-            phi = acos(cosphi);
-            lp.phi = HALFPI - phi;
-            if (area == AREA_0) {
-                lp.lam = theta + HALFPI;
-            } else if (area == AREA_1) {
-                lp.lam = (theta < 0.0 ? theta + PI : theta - PI);
-            } else if (area == AREA_2) {
-                lp.lam = theta - HALFPI;
-            } else /* area == AREA_3 */ {
-                lp.lam = theta;
-            }
-        } else if (P->face == FACE_BOTTOM) {
-            phi = acos(cosphi);
-            lp.phi = phi - HALFPI;
-            if (area == AREA_0) {
-                lp.lam = -theta + HALFPI;
-            } else if (area == AREA_1) {
-                lp.lam = -theta;
-            } else if (area == AREA_2) {
-                lp.lam = -theta - HALFPI;
-            } else /* area == AREA_3 */ {
-                lp.lam = (theta < 0.0 ? -theta - PI : -theta + PI);
-            }
+        if (Q->face == FACE_FRONT) {
+            phi = acos(q);
+            theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
+        } else if (Q->face == FACE_RIGHT) {
+            phi = acos(r);
+            theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
+        } else if (Q->face == FACE_BACK) {
+            phi = acos(-q);
+            theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
+        } else if (Q->face == FACE_LEFT) {
+            phi = acos(-r);
+            theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
         } else {
-            /* Compute phi and lam via cartesian unit sphere coordinates. */
-            double q, r, s, t;
-            q = cosphi;
-            t = q * q;
-            if (t >= 1.0) {
-                s = 0.0;
-            } else {
-                s = sqrt(1.0 - t) * sin(theta);
-            }
-            t += s * s;
-            if (t >= 1.0) {
-                r = 0.0;
-            } else {
-                r = sqrt(1.0 - t);
-            }
-            /* Rotate q,r,s into the correct area. */
-            if (area == AREA_1) {
-                t = r;
-                r = -s;
-                s = t;
-            } else if (area == AREA_2) {
-                r = -r;
-                s = -s;
-            } else if (area == AREA_3) {
-                t = r;
-                r = s;
-                s = -t;
-            }
-            /* Rotate q,r,s into the correct cube face. */
-            if (P->face == FACE_RIGHT) {
-                t = q;
-                q = -r;
-                r = t;
-            } else if (P->face == FACE_BACK) {
-                q = -q;
-                r = -r;
-            } else if (P->face == FACE_LEFT) {
-                t = q;
-                q = r;
-                r = -t;
-            }
-            /* Now compute phi and lam from the unit sphere coordinates. */
-            lp.phi = acos(-s) - HALFPI;
-            lp.lam = atan2(r, q);
-            if (P->face == FACE_RIGHT) {
-                lp.lam = qsc_shift_lon_origin(lp.lam, -HALFPI);
-            } else if (P->face == FACE_BACK) {
-                lp.lam = qsc_shift_lon_origin(lp.lam, -PI);
-            } else if (P->face == FACE_LEFT) {
-                lp.lam = qsc_shift_lon_origin(lp.lam, +HALFPI);
-            }
+            /* Impossible */
+            phi = theta = 0.0;
+            area = AREA_0;
         }
+    }
 
-        /* Apply the shift from the sphere to the ellipsoid as described
-         * in [LK12]. */
-        if (P->es) {
-            int invert_sign;
-            double tanphi, xa;
-            invert_sign = (lp.phi < 0.0 ? 1 : 0);
-            tanphi = tan(lp.phi);
-            xa = P->b / sqrt(tanphi * tanphi + P->one_minus_f_squared);
-            lp.phi = atan(sqrt(P->a * P->a - xa * xa) / (P->one_minus_f * xa));
-            if (invert_sign) {
-                lp.phi = -lp.phi;
-            }
-        }
-        return (lp);
+    /* Compute mu and nu for the area of definition.
+     * For mu, see Eq. (3-21) in [OL76], but note the typos:
+     * compare with Eq. (3-14). For nu, see Eq. (3-38). */
+    mu = atan((12.0 / PI) * (theta + acos(sin(theta) * cos(FORTPI)) - HALFPI));
+    t = sqrt((1.0 - cos(phi)) / (cos(mu) * cos(mu)) / (1.0 - cos(atan(1.0 / cos(theta)))));
+    /* nu = atan(t);        We don't really need nu, just t, see below. */
+
+    /* Apply the result to the real area. */
+    if (area == AREA_1) {
+        mu += HALFPI;
+    } else if (area == AREA_2) {
+        mu += PI;
+    } else if (area == AREA_3) {
+        mu += HALFPI + PI;
+    }
+
+    /* Now compute x, y from mu and nu */
+    /* t = tan(nu); */
+    xy.x = t * cos(mu);
+    xy.y = t * sin(mu);
+    return xy;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(qsc)
-        P->inv = e_inverse;
-        P->fwd = e_forward;
-        /* Determine the cube face from the center of projection. */
-        if (P->phi0 >= HALFPI - FORTPI / 2.0) {
-            P->face = FACE_TOP;
-        } else if (P->phi0 <= -(HALFPI - FORTPI / 2.0)) {
-            P->face = FACE_BOTTOM;
-        } else if (fabs(P->lam0) <= FORTPI) {
-            P->face = FACE_FRONT;
-        } else if (fabs(P->lam0) <= HALFPI + FORTPI) {
-            P->face = (P->lam0 > 0.0 ? FACE_RIGHT : FACE_LEFT);
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double mu, nu, cosmu, tannu;
+    double tantheta, theta, cosphi, phi;
+    double t;
+    int area;
+
+    /* Convert the input x, y to the mu and nu angles as used by QSC.
+     * This depends on the area of the cube face. */
+    nu = atan(sqrt(xy.x * xy.x + xy.y * xy.y));
+    mu = atan2(xy.y, xy.x);
+    if (xy.x >= 0.0 && xy.x >= fabs(xy.y)) {
+        area = AREA_0;
+    } else if (xy.y >= 0.0 && xy.y >= fabs(xy.x)) {
+        area = AREA_1;
+        mu -= HALFPI;
+    } else if (xy.x < 0.0 && -xy.x >= fabs(xy.y)) {
+        area = AREA_2;
+        mu = (mu < 0.0 ? mu + PI : mu - PI);
+    } else {
+        area = AREA_3;
+        mu += HALFPI;
+    }
+
+    /* Compute phi and theta for the area of definition.
+     * The inverse projection is not described in the original paper, but some
+     * good hints can be found here (as of 2011-12-14):
+     * http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
+     * (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
+    t = (PI / 12.0) * tan(mu);
+    tantheta = sin(t) / (cos(t) - (1.0 / sqrt(2.0)));
+    theta = atan(tantheta);
+    cosmu = cos(mu);
+    tannu = tan(nu);
+    cosphi = 1.0 - cosmu * cosmu * tannu * tannu * (1.0 - cos(atan(1.0 / cos(theta))));
+    if (cosphi < -1.0) {
+        cosphi = -1.0;
+    } else if (cosphi > +1.0) {
+        cosphi = +1.0;
+    }
+
+    /* Apply the result to the real area on the cube face.
+     * For the top and bottom face, we can compute phi and lam directly.
+     * For the other faces, we must use unit sphere cartesian coordinates
+     * as an intermediate step. */
+    if (Q->face == FACE_TOP) {
+        phi = acos(cosphi);
+        lp.phi = HALFPI - phi;
+        if (area == AREA_0) {
+            lp.lam = theta + HALFPI;
+        } else if (area == AREA_1) {
+            lp.lam = (theta < 0.0 ? theta + PI : theta - PI);
+        } else if (area == AREA_2) {
+            lp.lam = theta - HALFPI;
+        } else /* area == AREA_3 */ {
+            lp.lam = theta;
+        }
+    } else if (Q->face == FACE_BOTTOM) {
+        phi = acos(cosphi);
+        lp.phi = phi - HALFPI;
+        if (area == AREA_0) {
+            lp.lam = -theta + HALFPI;
+        } else if (area == AREA_1) {
+            lp.lam = -theta;
+        } else if (area == AREA_2) {
+            lp.lam = -theta - HALFPI;
+        } else /* area == AREA_3 */ {
+            lp.lam = (theta < 0.0 ? -theta - PI : -theta + PI);
+        }
+    } else {
+        /* Compute phi and lam via cartesian unit sphere coordinates. */
+        double q, r, s, t;
+        q = cosphi;
+        t = q * q;
+        if (t >= 1.0) {
+            s = 0.0;
         } else {
-            P->face = FACE_BACK;
+            s = sqrt(1.0 - t) * sin(theta);
         }
-        /* Fill in useful values for the ellipsoid <-> sphere shift
-         * described in [LK12]. */
-        if (P->es) {
-            P->a_squared = P->a * P->a;
-            P->b = P->a * sqrt(1.0 - P->es);
-            P->one_minus_f = 1.0 - (P->a - P->b) / P->a;
-            P->one_minus_f_squared = P->one_minus_f * P->one_minus_f;
+        t += s * s;
+        if (t >= 1.0) {
+            r = 0.0;
+        } else {
+            r = sqrt(1.0 - t);
         }
-ENDENTRY(P)
+        /* Rotate q,r,s into the correct area. */
+        if (area == AREA_1) {
+            t = r;
+            r = -s;
+            s = t;
+        } else if (area == AREA_2) {
+            r = -r;
+            s = -s;
+        } else if (area == AREA_3) {
+            t = r;
+            r = s;
+            s = -t;
+        }
+        /* Rotate q,r,s into the correct cube face. */
+        if (Q->face == FACE_RIGHT) {
+            t = q;
+            q = -r;
+            r = t;
+        } else if (Q->face == FACE_BACK) {
+            q = -q;
+            r = -r;
+        } else if (Q->face == FACE_LEFT) {
+            t = q;
+            q = r;
+            r = -t;
+        }
+        /* Now compute phi and lam from the unit sphere coordinates. */
+        lp.phi = acos(-s) - HALFPI;
+        lp.lam = atan2(r, q);
+        if (Q->face == FACE_RIGHT) {
+            lp.lam = qsc_shift_lon_origin(lp.lam, -HALFPI);
+        } else if (Q->face == FACE_BACK) {
+            lp.lam = qsc_shift_lon_origin(lp.lam, -PI);
+        } else if (Q->face == FACE_LEFT) {
+            lp.lam = qsc_shift_lon_origin(lp.lam, +HALFPI);
+        }
+    }
+
+    /* Apply the shift from the sphere to the ellipsoid as described
+     * in [LK12]. */
+    if (P->es) {
+        int invert_sign;
+        double tanphi, xa;
+        invert_sign = (lp.phi < 0.0 ? 1 : 0);
+        tanphi = tan(lp.phi);
+        xa = Q->b / sqrt(tanphi * tanphi + Q->one_minus_f_squared);
+        lp.phi = atan(sqrt(P->a * P->a - xa * xa) / (Q->one_minus_f * xa));
+        if (invert_sign) {
+            lp.phi = -lp.phi;
+        }
+    }
+    return lp;
+}
+
+
+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(qsc) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    P->inv = e_inverse;
+    P->fwd = e_forward;
+    /* Determine the cube face from the center of projection. */
+    if (P->phi0 >= HALFPI - FORTPI / 2.0) {
+        Q->face = FACE_TOP;
+    } else if (P->phi0 <= -(HALFPI - FORTPI / 2.0)) {
+        Q->face = FACE_BOTTOM;
+    } else if (fabs(P->lam0) <= FORTPI) {
+        Q->face = FACE_FRONT;
+    } else if (fabs(P->lam0) <= HALFPI + FORTPI) {
+        Q->face = (P->lam0 > 0.0 ? FACE_RIGHT : FACE_LEFT);
+    } else {
+        Q->face = FACE_BACK;
+    }
+    /* Fill in useful values for the ellipsoid <-> sphere shift
+     * described in [LK12]. */
+    if (P->es) {
+        Q->a_squared = P->a * P->a;
+        Q->b = P->a * sqrt(1.0 - P->es);
+        Q->one_minus_f = 1.0 - (P->a - Q->b) / P->a;
+        Q->one_minus_f_squared = Q->one_minus_f * Q->one_minus_f;
+    }
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_qsc_selftest (void) {return 0;}
+#else
+
+int pj_qsc_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=qsc   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=qsc   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 304638.450843852363,  164123.870923793991},
+        { 304638.450843852363, -164123.870923793991},
+        {-304638.450843852363,  164123.870923793962},
+        {-304638.450843852421, -164123.870923793904},
+    };
+
+    XY s_fwd_expect[] = {
+        { 305863.792402890511,  165827.722754715243},
+        { 305863.792402890511, -165827.722754715243},
+        {-305863.792402890511,  165827.722754715243},
+        {-305863.792402890569, -165827.722754715156},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.00132134098471627126,  0.000610652900922527926},
+        { 0.00132134098471627126, -0.000610652900922527926},
+        {-0.00132134098471627126,  0.000610652900922527926},
+        {-0.00132134098471627126, -0.000610652900922527926},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.00131682718763827234,  0.000604493198178676161},
+        { 0.00131682718763827234, -0.000604493198178676161},
+        {-0.00131682718763827234,  0.000604493198178676161},
+        {-0.00131682718763827234, -0.000604493198178676161},
+    };
+
+    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