Elliptic integrals with JavaScript - depends on math.js

ellipticIntegrals.js

/*
Carlson elliptic integrals and incomplete elliptic integrals.
Author: Stéphane Laurent.
Depends on math.js <https://mathjs.org/>.
*/

var ellipticIntegrals = (function () {

  let CarlsonRF = function(x, y, z, err) {
    if(err <= 0) {
      throw "The value of `err` must be nonnegative.";
    }
    let nzeros = math.equal(x, 0) + math.equal(y, 0) + math.equal(z, 0);
    if(nzeros > 1) {
      throw "At most one of `x`, `y`, `z` can be 0.";
    }
    let dx = 2 * err;
    let dy = dx;
    let dz = dx;
    let A;
    while(dx > err || dy > err || dz > err) {
      let srx = math.sqrt(x);
      let sry = math.sqrt(y);
      let srz = math.sqrt(z);
      let lambda = math.add(
        math.multiply(srx, sry),
        math.add(
          math.multiply(sry, srz),
          math.multiply(srz, srx)
        )
      );
      x = math.divide(math.add(x, lambda), 4);
      y = math.divide(math.add(y, lambda), 4);
      z = math.divide(math.add(z, lambda), 4);
      A = math.divide(math.add(x, math.add(y, z)), 3);
      dx = math.abs(math.divide(math.subtract(A, x), A));
      dy = math.abs(math.divide(math.subtract(A, y), A));
      dz = math.abs(math.divide(math.subtract(A, z), A));
    }
    let E2 = dx*dy + dy*dz + dz*dx;
    let E3 = dy*dx*dz;
    let h = 1 - E2/10 + E3/14 + E2*E2/24 - 3*E2*E3/44 - 5*E2*E2*E2/208
            + 3*E3*E3/104 + E2*E2*E3/16;
    return math.divide(h, math.sqrt(A));
  };

  let CarlsonRD = function(x, y, z, err) {
    if(err <= 0) {
      throw "The value of `err` must be nonnegative.";
    }
    let nzeros = math.equal(x, 0) + math.equal(y, 0) + math.equal(z, 0);
    if(nzeros > 1) {
      throw "At most one of `x`, `y`, `z` can be 0.";
    }
    let dx = 2 * err;
    let dy = dx;
    let dz = dx;
    let s = math.complex(0, 0);
    let fac = 1;
    let A;
    while(dx > err || dy > err || dz > err) {
      let srx = math.sqrt(x);
      let sry = math.sqrt(y);
      let srz = math.sqrt(z);
      let lambda = math.add(
        math.multiply(srx, sry),
        math.add(
          math.multiply(sry, srz),
          math.multiply(srz, srx)
        )
      );
      s = math.add(
        s, math.divide(fac, math.multiply(srz, math.add(z, lambda)))
      );
      fac /= 4;
      x = math.divide(math.add(x, lambda), 4);
      y = math.divide(math.add(y, lambda), 4);
      z = math.divide(math.add(z, lambda), 4);
      A = math.divide(math.add(x, math.add(y, math.multiply(3, z))), 5);
      dx = math.abs(math.divide(math.subtract(A, x), A));
      dy = math.abs(math.divide(math.subtract(A, y), A));
      dz = math.abs(math.divide(math.subtract(A, z), A));
    }
    let E2 = dx*dy + 3*dy*dz + 3*dz*dz + 3*dx*dz;
    let E3 = dz*dz*dz + 3*dx*dz*dz + 3*dx*dy*dz + 3*dy*dz*dz;
    let E4 = dy*dz*dz*dz + dx*dz*dz*dz + 3*dx*dy*dz*dz;
    let E5 = dx*dy*dz*dz*dz;
    let h = fac * (1 - 3*E2/14 + E3/6 + 9*E2*E2/88 - 3*E4/22
               - 9*E2*E3/52 + 3*E5/26 - E2*E2*E2/16 + 3*E3*E3/40
               + 3*E2*E4/20 + 45*E2*E2*E3/272 - 9*(E3*E4 + E2*E5)/68);
    return math.add(
      math.multiply(3, s), math.divide(h, math.multiply(A, math.sqrt(A)))
    );
  };

  let CarlsonRJ = function(x, y, z, p, err) {
    if(err <= 0) {
      throw "The value of `err` must be nonnegative.";
    }
    let nzeros = math.equal(x, 0) + math.equal(y, 0) + 
      math.equal(z, 0) + math.equal(p, 0);
    if(nzeros > 1) {
      throw "At most one of `x`, `y`, `z`, `p` can be 0.";
    }
    let A0 = math.divide(
      math.add(math.add(math.add(math.add(x, y), z), p), p), 5
    );
    let A = math.clone(A0);
    let delta = 
      math.multiply(
        math.multiply(
          math.subtract(p, x), math.subtract(p, y)
        ), 
        math.subtract(p, z)
      );
    let M = Math.max(
              math.abs(math.subtract(A, x)),
              math.abs(math.subtract(A, y)),
              math.abs(math.subtract(A, z))
            );
    let Q = Math.pow(4/err, 1/6) * M;
    let d = [];
    let e = [];
    let f = 1;
    let fac = 1;
    while(math.abs(A) <= Q) {
      let srx = math.sqrt(x);
      let sry = math.sqrt(y);
      let srz = math.sqrt(z);
      let srp = math.sqrt(p);
      let dnew = math.multiply(
        math.add(srp, srx), 
        math.multiply(math.add(srp, sry), math.add(srp, srz))
      );
      d.push(math.multiply(f, dnew));
      e.push(math.divide(math.multiply(fac, delta), math.square(dnew)));
      f *= 4;
      fac /= 64;
      let lambda = math.add(
        math.multiply(srx, sry),
        math.add(
          math.multiply(sry, srz),
          math.multiply(srz, srx)
        )
      );
      x = math.divide(math.add(x, lambda), 4);
      y = math.divide(math.add(y, lambda), 4);
      z = math.divide(math.add(z, lambda), 4);
      p = math.divide(math.add(p, lambda), 4);
      A = math.divide(math.add(A, lambda), 4);
      Q /= 4;
    }
    let fA = math.multiply(f, A);
    let X = math.divide(math.subtract(A0, x), fA);
    let Y = math.divide(math.subtract(A0, y), fA);
    let Z = math.divide(math.subtract(A0, z), fA);
    let P = math.divide(math.unaryMinus(math.add(math.add(X, Y), Z)), 2);
    let P2 = math.square(P);
    let E2 = math.subtract(
      math.add(
        math.add(math.multiply(X, Y), math.multiply(Y, Z)), math.multiply(Z, X)
      ), 
      math.multiply(3, P2)
    );
    let E3 = math.add(
      math.multiply(math.multiply(X, Y), Z), 
      math.multiply(math.add(math.multiply(2, E2), math.multiply(4, P2)), P)
    );
    let E4 = math.multiply(
      P, 
      math.add(
        math.multiply(2, math.multiply(math.multiply(X, Y), Z)), 
        math.multiply(P, math.add(E2, math.multiply(3, P2)))
      )
    );
    let E5 = math.multiply(math.multiply(math.multiply(X, Y), Z), P2);
    let h = math.divide(
      math.add(
        math.subtract(
          math.subtract(
            math.add(
              math.add(
                math.subtract(1, math.divide(math.multiply(3, E2), 14)), 
                math.divide(E3, 6)
              ), 
              math.divide(math.multiply(math.multiply(9, E2), E2), 88)
            ), 
            math.divide(math.multiply(3, E4), 22)
          ), 
          math.divide(math.multiply(math.multiply(9, E2), E3), 52)
        ), 
        math.divide(math.multiply(3, E5), 26)
      ), 
      f
    );
    let s = math.complex(0, 0);
    let one = math.complex(1, 0);
    let n = e.length;
    for(let i = 0; i < n; i++) {
      let srei = math.sqrt(e[i]);
      let a = math.equal(srei, 0) ? one : math.divide(math.atan(srei), srei);
      s = math.add(s, math.divide(a, d[i]));
    }
    return math.add(
      math.divide(h, math.multiply(A, math.sqrt(A))), 
      math.multiply(6, s)
    );
  };

  let ellipticE = function(phi, m, err) {
    let out;
    let phi_r = math.re(phi);
    let PI = Math.PI;
    let PI_2 = PI / 2;
    if(math.equal(phi, 0)) {
      out = math.complex(0, 0);
    } else if(phi_r >= -PI_2 && phi_r <= PI_2) {
      if(math.equal(m, 0)) {
        out = phi;
      } else if(math.equal(m, 1)) {
        out = math.sin(phi);
      } else {
        let sine = math.sin(phi);
        let sine2 = math.square(sine);
        let cosine2 = math.subtract(1, sine2);
        let oneminusmsine2 = math.subtract(1, math.multiply(m, sine2));
        out = math.multiply(
          sine, 
          math.subtract(
            CarlsonRF(cosine2, oneminusmsine2, 1, err),
            math.divide(
              math.multiply(
                m, 
                math.multiply(
                  sine2, 
                  CarlsonRD(cosine2, oneminusmsine2, 1, err)
                )
              ), 
              3
            )
          )
        );
      }
    } else {
      let k = phi_r > PI_2 ?
        Math.ceil(phi_r/PI - 0.5) : -Math.floor(0.5 - phi_r/PI);
      phi = math.subtract(phi, k*PI);
      out = math.add(
        math.multiply(2*k, ellipticE(PI_2, m, err)), 
        ellipticE(phi, m, err)
      );
    }
    return out;
  };

  let isInf = function(x) {
    return x === Infinity || x === -Infinity;
  };

  let ellipticF = function(phi, m, err) {
    let out;
    let phi_r = math.re(phi);
    let phi_i = math.im(phi);
    let m_r = math.re(m);
    let m_i = math.im(m);
    let PI = Math.PI;
    let PI_2 = PI / 2;
    if(math.equal(phi, 0) || isInf(m_r) || isInf(m_i)) {
      out = math.complex(0,0);
    } else if(
        phi_r == 0 && isInf(phi_i) && m_i === 0 && m_r > 0 && m_r < 1
    ) {
      let minv = 1 / m_r;
      let msqrt = Math.sqrt(m_r);
      out = math.subtract(
        ellipticF(PI_2, m, err), 
        math.divide(ellipticF(PI_2, minv, err), msqrt)
      );
      if(phi_i < 0) {
        out = math.unaryMinus(out);
      }
    } else if((phi_r === PI_2 || phi_r === -PI_2) && math.equal(m, 1)) {
      out = math.complex(NaN, NaN);
    } else if(phi_r >= -PI_2 && phi_r <= PI_2) {
      if(math.equal(m, 1) && Math.abs(phi_r) < PI_2) {
        out = math.asinh(math.tan(phi));
      } else if(math.equal(m, 0)) {
        out = phi;
      } else {
        let sine = math.sin(phi);
        let sine2 = math.square(sine);
        let cosine2 = math.subtract(1, sine2);
        let oneminusmsine2 = math.subtract(1, math.multiply(m, sine2));
        out = math.multiply(sine, CarlsonRF(cosine2, oneminusmsine2, 1, err));
      }
    } else {
      let k = phi_r > PI_2 ?
        Math.ceil(phi_r/PI - 0.5) : -Math.floor(0.5 - phi_r/PI);
      phi = math.subtract(phi, k*PI);
      out = math.add(
        math.multiply(2*k, ellipticF(PI_2, m, err)), 
        ellipticF(phi, m, err)
      );
    }
    return out;
  };
  
  let ellipticZ = function(phi, m, err) {
    let out;
    let phi_r = math.re(phi);
    let m_r = math.re(m);
    let m_i = math.im(m);
    let PI = Math.PI;
    let PI_2 = PI / 2;
    if(isInf(m_r) && m_i === 0) {
      out = math.complex(NaN, NaN);
    } else if(math.equal(m, 1)) {
      if(Math.abs(phi_r) <= PI_2) {
        out = math.sin(phi);
      } else {
        let k = phi_r > PI_2 ?
          Math.ceil(phi_r/PI - 0.5) : -Math.floor(0.5 - phi_r/PI);
        phi = math.subtract(phi, k*PI);
        out = math.sin(phi);
      }
    } else {
      out = math.subtract(
        ellipticE(phi, m, err),
        math.multiply(
          math.divide(
            ellipticE(PI_2, m, err), 
            ellipticF(PI_2, m, err)
          ), 
          ellipticF(phi, m, err)
        )
      );
    }
    return out;
  };
  
  let ellipticPI = function(phi, n, m, err) {
    let out;
    let phi_r = math.re(phi);
    let phi_i = math.im(phi);
    let m_r = math.re(m);
    let m_i = math.im(m);
    let n_r = math.re(n);
    let n_i = math.im(n);
    let PI = Math.PI;
    let PI_2 = PI / 2;
    let complete = math.equal(phi, PI_2);
    if(
        math.equal(phi, 0) ||
          (isInf(n_r) && n_i === 0) ||
          (isInf(m_r) && m_i === 0)
    ) {
      out = math.complex(0, 0);
    } else if(complete && math.equal(m, 1) && n_r !== 1 && n_i === 0) {
      out = math.complex(n_r > 1.0 ? -Infinity : Infinity, 0);
    } else if(complete && math.equal(n, 1)) {
      out = math.complex(NaN, NaN);
    } else if(complete && math.equal(m, 0)) {
      out = math.divide(PI_2, math.sqrt(math.subtract(1, n)));
    } else if(complete && math.equal(n, m)) {
      out = math.divide(ellipticE(PI_2, m, err), math.subtract(1, m));
    } else if(complete && math.equal(n, 0)) {
      out = ellipticF(PI_2, m, err);
    } else if(phi_r >= -PI_2 && phi_r <= PI_2) {
      let sine = math.sin(phi);
      let sine2 = math.square(sine);
      let cosine2 = math.subtract(1, sine2);
      let oneminusmsine2 = math.subtract(1, math.multiply(m, sine2));
      let oneminusnsine2 = math.subtract(1, math.multiply(n, sine2));
      out = math.multiply(
        sine, 
        math.add(
          CarlsonRF(cosine2, oneminusmsine2, 1, err),
          math.divide(
            math.multiply(
              math.multiply(n, sine2),
              CarlsonRJ(cosine2, oneminusmsine2, 1, oneminusnsine2, err)
            ),
            3
          )
        )
      );
    } else {
      let k = phi_r > PI_2 ?
        Math.ceil(phi_r/PI - 0.5) : -Math.floor(0.5 - phi_r/PI);
      phi = math.subtract(phi, k*PI);
      out = math.add(
        math.multiply(2*k, ellipticPI(PI_2, n, m, err)), 
        ellipticPI(phi, n, m, err)
      );
    }
    return out;
  };


  return {
    CarlsonRF: CarlsonRF,
    CarlsonRD: CarlsonRD,
    CarlsonRJ: CarlsonRJ,
    ellipticE: ellipticE,
    ellipticF: ellipticF,
    ellipticZ: ellipticZ,
    ellipticPI: ellipticPI
  }
})();