text: import {abs, asin, cos, epsilon, pi, sin, sqrt} from "./math.js";
  
  export default function(a, b, c, d, e, f, g, h) {
    if (arguments.length < 8) h = 0;
  
    function forward(lambda, phi) {
      if (!phi) return [a * lambda / pi, 0];
      var phi2 = phi * phi,
          xB = a + phi2 * (b + phi2 * (c + phi2 * d)),
          yB = phi * (e - 1 + phi2 * (f - h + phi2 * g)),
          m = (xB * xB + yB * yB) / (2 * yB),
          alpha = lambda * asin(xB / m) / pi;
      return [m * sin(alpha), phi * (1 + phi2 * h) + m * (1 - cos(alpha))];
    }
  
    forward.invert = function(x, y) {
      var lambda = pi * x / a,
          phi = y,
          deltaLambda, deltaPhi, i = 50;
      do {
        var phi2 = phi * phi,
            xB = a + phi2 * (b + phi2 * (c + phi2 * d)),
            yB = phi * (e - 1 + phi2 * (f - h + phi2 * g)),
            p = xB * xB + yB * yB,
            q = 2 * yB,
            m = p / q,
            m2 = m * m,
            dAlphadLambda = asin(xB / m) / pi,
            alpha = lambda * dAlphadLambda,
            xB2 = xB * xB,
            dxBdPhi = (2 * b + phi2 * (4 * c + phi2 * 6 * d)) * phi,
            dyBdPhi = e + phi2 * (3 * f + phi2 * 5 * g),
            dpdPhi = 2 * (xB * dxBdPhi + yB * (dyBdPhi - 1)),
            dqdPhi = 2 * (dyBdPhi - 1),
            dmdPhi = (dpdPhi * q - p * dqdPhi) / (q * q),
            cosAlpha = cos(alpha),
            sinAlpha = sin(alpha),
            mcosAlpha = m * cosAlpha,
            msinAlpha = m * sinAlpha,
            dAlphadPhi = ((lambda / pi) * (1 / sqrt(1 - xB2 / m2)) * (dxBdPhi * m - xB * dmdPhi)) / m2,
            fx = msinAlpha - x,
            fy = phi * (1 + phi2 * h) + m - mcosAlpha - y,
            deltaxDeltaPhi = dmdPhi * sinAlpha + mcosAlpha * dAlphadPhi,
            deltaxDeltaLambda = mcosAlpha * dAlphadLambda,
            deltayDeltaPhi = 1 + dmdPhi - (dmdPhi * cosAlpha - msinAlpha * dAlphadPhi),
            deltayDeltaLambda = msinAlpha * dAlphadLambda,
            denominator = deltaxDeltaPhi * deltayDeltaLambda - deltayDeltaPhi * deltaxDeltaLambda;
        if (!denominator) break;
        lambda -= deltaLambda = (fy * deltaxDeltaPhi - fx * deltayDeltaPhi) / denominator;
        phi -= deltaPhi = (fx * deltayDeltaLambda - fy * deltaxDeltaLambda) / denominator;
      } while ((abs(deltaLambda) > epsilon || abs(deltaPhi) > epsilon) && --i > 0);
      return [lambda, phi];
    };
  
    return forward;
  }
  