2 files changed, 14 insertions, 5 deletions
diff --git a/src/projections/lsat.cpp b/src/projections/lsat.cpp
index 5b7520d3..f9eec1b9 100644
--- a/src/projections/lsat.cpp
+++ b/src/projections/lsat.cpp
@@ -112,7 +112,6 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
     int nn;
     double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
-
     lamdp = xy.x / Q->b;
     nn = 50;
     do {
@@ -135,10 +134,14 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
         lamdp -= TOL;
     spp = sin(phidp);
     sppsq = spp * spp;
+    const double denom = 1. - sppsq * (1. + Q->u);
+    if( denom == 0.0 ) {
+        proj_errno_set(P, PJD_ERR_INVALID_X_OR_Y);
+        return proj_coord_error().lp;
+    }
     lamt = atan(((1. - sppsq * P->rone_es) * tan(lamdp) *
         Q->ca - spp * Q->sa * sqrt((1. + Q->q * dd) * (
-        1. - sppsq) - sppsq * Q->u) / cos(lamdp)) / (1. - sppsq
-        * (1. + Q->u)));
+        1. - sppsq) - sppsq * Q->u) / cos(lamdp)) / denom);
     sl = lamt >= 0. ? 1. : -1.;
     scl = cos(lamdp) >= 0. ? 1. : -1;
     lamt -= M_HALFPI * (1. - scl) * sl;
diff --git a/src/projections/omerc.cpp b/src/projections/omerc.cpp
index c33f6489..c0278043 100644
--- a/src/projections/omerc.cpp
+++ b/src/projections/omerc.cpp
@@ -191,6 +191,9 @@ PJ *PROJECTION(omerc) {
                 gamma = alpha_c;
         } else
             alpha_c = aasin(P->ctx, D*sin(gamma0 = gamma));
+        if( fabs(fabs(P->phi0) - M_HALFPI) <= TOL ) {
+            return pj_default_destructor(P, PJD_ERR_LAT_0_OR_ALPHA_EQ_90);
+        }
         P->lam0 = lamc - aasin(P->ctx, .5 * (F - 1. / F) *
            tan(gamma0)) / Q->B;
     } else {
@@ -206,8 +209,11 @@ PJ *PROJECTION(omerc) {
             lam2 += M_TWOPI;
         P->lam0 = adjlon(.5 * (lam1 + lam2) - atan(
            J * tan(.5 * Q->B * (lam1 - lam2)) / p) / Q->B);
-        gamma0 = atan(2. * sin(Q->B * adjlon(lam1 - P->lam0)) /
-           (F - 1. / F));
+        const double denom = F - 1. / F;
+        if( denom == 0 ) {
+            return pj_default_destructor(P, PJD_ERR_INVALID_ECCENTRICITY);
+        }
+        gamma0 = atan(2. * sin(Q->B * adjlon(lam1 - P->lam0)) / denom);
         gamma = alpha_c = aasin(P->ctx, D * sin(gamma0));
     }
     Q->singam = sin(gamma0);