/*
This is the equation using the normal:
vx * X + vy * Y + vz * Z + D = 0
where D = - px * vx - py * vy - pz * vz
337122
function Plane3D( a , b , c , d ) {
	var self = Object.create( Plane3D.prototype ) ;

	self.normal = new Vector3D( a , b , c ) ;
	self.d = + d ;
	self.normalizeCheck() ;

	return self ;
}
189750
def dot(v1, v2):
    '''
    '''
    x1, y1, z1 = v1
    x2, y2, z2 = v2
    return x1 * x2 + y1 * y2 + z1 * z2
1295840
function(v) {

      // [0, v'] = Q * [0, v] * Q'

      // Q
      var w1 = this['w'];
      var x1 = this['x'];
      var y1 = this['y'];
      var z1 = this['z'];

      // [0, v]
      var w2 = 0;
      var x2 = v[0];
      var y2 = v[1];
      var z2 = v[2];

      // Q * [0, v]
      var w3 = /*w1 * w2*/ -x1 * x2 - y1 * y2 - z1 * z2;
      var x3 = w1 * x2 + /*x1 * w2 +*/ y1 * z2 - z1 * y2;
      var y3 = w1 * y2 + /*y1 * w2 +*/ z1 * x2 - x1 * z2;
      var z3 = w1 * z2 + /*z1 * w2 +*/ x1 * y2 - y1 * x2;

      var w4 = w3 * w1 + x3 * x1 + y3 * y1 + z3 * z1;
      var x4 = x3 * w1 - w3 * x1 - y3 * z1 + z3 * y1;
      var y4 = y3 * w1 - w3 * y1 - z3 * x1 + x3 * z1;
      var z4 = z3 * w1 - w3 * z1 - x3 * y1 + y3 * x1;

      return [x4, y4, z4];
    }
575105
function(x,y,z,normalX,normalY,normalZ,penetration,flip){
        var p=this.points[this.numPoints++];
        p.position.x=x;
        p.position.y=y;
        p.position.z=z;
        var r=this.body1.rotation;
        var rx=x-this.body1.position.x;
        var ry=y-this.body1.position.y;
        var rz=z-this.body1.position.z;

        var tr = r.elements;
        p.localPoint1.x=rx*tr[0]+ry*tr[3]+rz*tr[6];
        p.localPoint1.y=rx*tr[1]+ry*tr[4]+rz*tr[7];
        p.localPoint1.z=rx*tr[2]+ry*tr[5]+rz*tr[8];
        r=this.body2.rotation;
        rx=x-this.body2.position.x;
        ry=y-this.body2.position.y;
        rz=z-this.body2.position.z;
        p.localPoint2.x=rx*tr[0]+ry*tr[3]+rz*tr[6];
        p.localPoint2.y=rx*tr[1]+ry*tr[4]+rz*tr[7];
        p.localPoint2.z=rx*tr[2]+ry*tr[5]+rz*tr[8];

        p.normalImpulse=0;
        if(flip){
            p.normal.x=-normalX;
            p.normal.y=-normalY;
            p.normal.z=-normalZ;
        }else{
            p.normal.x=normalX;
            p.normal.y=normalY;
            p.normal.z=normalZ;
        }
        p.penetration=penetration;
        p.warmStarted=false;
    }
1286122
function phi1z(eccent, qs) {
        var sinphi, cosphi, con, com, dphi;
        var phi = asinz(0.5 * qs);
        if (eccent < EPSLN) {
            return phi;
        }

        var eccnts = eccent * eccent;
        for (var i = 1; i <= 25; i++) {
            sinphi = Math.sin(phi);
            cosphi = Math.cos(phi);
            con = eccent * sinphi;
            com = 1 - con * con;
            dphi = 0.5 * com * com / cosphi * (qs / (1 - eccnts) - sinphi / com + 0.5 / eccent * Math.log((1 - con) / (1 + con)));
            phi = phi + dphi;
            if (Math.abs(dphi) <= 1e-7) {
                return phi;
            }
        }
        return null;
    }
108445
def exp(self):
        
        
        # Init
        vecNorm = self.x**2 + self.y**2 + self.z**2
        wPart = np.exp(self.w)        
        q = Quaternion()
        
        # Calculate
        q.w = wPart * np.cos(vecNorm)
        q.x = wPart * self.x * np.sin(vecNorm) / vecNorm
        q.y = wPart * self.y * np.sin(vecNorm) / vecNorm
        q.z = wPart * self.z * np.sin(vecNorm) / vecNorm
        
        return q
1072072
function( ax, ay,
								   bx, by,
								   cx, cy,
								   px, py ) {

		  var aX, aY, bX, bY;
		  var cX, cY, apx, apy;
		  var bpx, bpy, cpx, cpy;
		  var cCROSSap, bCROSScp, aCROSSbp;

		  aX = cx - bx;  aY = cy - by;
		  bX = ax - cx;  bY = ay - cy;
		  cX = bx - ax;  cY = by - ay;
		  apx= px  -ax;  apy= py - ay;
		  bpx= px - bx;  bpy= py - by;
		  cpx= px - cx;  cpy= py - cy;

		  aCROSSbp = aX*bpy - aY*bpx;
		  cCROSSap = cX*apy - cY*apx;
		  bCROSScp = bX*cpy - bY*cpx;

		  return ( (aCROSSbp >= 0.0) && (bCROSScp >= 0.0) && (cCROSSap >= 0.0) );

	}
662914
function(w, x, y, z) {

      parse(P, w, x, y, z);

      // dot(Q1, Q2) := w1 * w2 + dot(v1, v2)

      return this['w'] * P['w'] + this['x'] * P['x'] + this['y'] * P['y'] + this['z'] * P['z'];
    }
86176
public Vector4d hermite(Vector4dc t0, Vector4dc v1, Vector4dc t1, double t, Vector4d dest) {
        double t2 = t * t;
        double t3 = t2 * t;
        dest.x = (x + x - v1.x() - v1.x() + t1.x() + t0.x()) * t3 + (3.0 * v1.x() - 3.0 * x - t0.x() - t0.x() - t1.x()) * t2 + x * t + x;
        dest.y = (y + y - v1.y() - v1.y() + t1.y() + t0.y()) * t3 + (3.0 * v1.y() - 3.0 * y - t0.y() - t0.y() - t1.y()) * t2 + y * t + y;
        dest.z = (z + z - v1.z() - v1.z() + t1.z() + t0.z()) * t3 + (3.0 * v1.z() - 3.0 * z - t0.z() - t0.z() - t1.z()) * t2 + z * t + z;
        dest.w = (w + w - v1.w() - v1.w() + t1.w() + t0.w()) * t3 + (3.0 * v1.w() - 3.0 * w - t0.w() - t0.w() - t1.w()) * t2 + w * t + w;
        return dest;
    }
304942
function eq (y, c) {
  const J0 = base.horner(y * 0.001, c)
  const T = base.J2000Century(J0)
  const W = 35999.373 * D2R * T - 2.47 * D2R
  const Δλ = 1 + 0.0334 * cos(W) + 0.0007 * cos(2 * W)
  let S = 0
  for (let i = terms.length - 1; i >= 0; i--) {
    const t = terms[i]
    S += t.a * cos((t.b + t.c * T) * D2R)
  }
  return J0 + 0.00001 * S / Δλ
}
151450
def equation(self):
        
        mat = np.empty((self.nunknowns, self.model.neq))
        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        for icp in range(self.ncp):
            istart = icp * self.nlayers
            ieq = 0
            for e in self.model.elementlist:
                if e.nunknowns > 0:
                    qxin, qyin = e.disvecinflayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin)
                    qxout, qyout = e.disvecinflayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout)
                    mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                        (qxin - qxout) * self.cosnorm[icp] + (qyin - qyout) * self.sinnorm[icp]
                    ieq += e.nunknowns
                else:
                    qxin, qyin = e.disveclayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin)
                    qxout, qyout = e.disveclayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout)
                    rhs[istart:istart + self.nlayers] -= (qxin - qxout) * self.cosnorm[icp] + (qyin - qyout) * \
                                                         self.sinnorm[icp]
        return mat, rhs
79851
def _compute_equations(self, x, verbose=False):
        '''
        '''
        # compute the error and the normals.
        normals = []
        values = []
        signs = []
        error = 0.0
        if verbose:
            print()
            print(' '.join('% 10.3e' % val for val in x), end=' ')
            active_str = ''
        for i, (sign, equation) in enumerate(self.equations):
            value, normal = equation(x)
            if (i < len(self.lock) and self.lock[i]) or \
               (sign==-1 and value > -self.threshold) or \
               (sign==0) or (sign==1 and value < self.threshold):
                values.append(value)
                normals.append(normal)
                signs.append(sign)
                error += value**2
                if verbose:
                    active_str += 'X'
                if i < len(self.lock):
                    self.lock[i] = True
            elif verbose:
                active_str += '-'
        error = np.sqrt(error)
        normals = np.array(normals, float)
        values = np.array(values, float)
        signs = np.array(signs, int)
        if verbose:
            print('[%s]' % active_str, end=' ')
            if error < self.threshold:
                print('OK')
            else:
                print('%.5e' % error)
        return normals, values, error, signs
999536
def norm(self):
        
        tmp = self.w**2 + self.x**2 + self.y**2 + self.z**2
        return tmp**0.5
947326
def _G(self, x, p):
        

        prefactor = (p + p ** 3) ** -1 * p

        if isinstance(x, np.ndarray):

            inds0 = np.where(x * p == 1)
            inds1 = np.where(x * p < 1)
            inds2 = np.where(x * p > 1)

            func = np.ones_like(x)

            func[inds0] = np.log(0.25 * x[inds0] ** 2 * p ** 2) + np.pi * p * (self._u(x[inds0]) - 1) + \
                          2 * p ** 2 * (self._u(x[inds0]) * np.arctanh(self._u(x[inds0]) ** -1) +
                                        np.log(0.5 * x[inds0]))

            func[inds1] = np.log(0.25 * x[inds1] ** 2 * p ** 2) + np.pi * p * (self._u(x[inds1]) - 1) + \
                          2 * p ** 2 * (self._u(x[inds1]) * np.arctanh(self._u(x[inds1]) ** -1) +
                                        np.log(0.5 * x[inds1])) + 2 * self._g(x[inds1], p) * np.arctanh(
                self._g(x[inds1], p))

            func[inds2] = np.log(0.25 * x[inds2] ** 2 * p ** 2) + np.pi * p * (self._u(x[inds2]) - 1) + \
                          2 * p ** 2 * (self._u(x[inds2]) * np.arctanh(self._u(x[inds2]) ** -1) +
                                        np.log(0.5 * x[inds2])) - 2 * self._f(x[inds2], p) * np.arctan(
                self._f(x[inds2], p))


        else:

            if x * p == 1:

                func = np.log(0.25 * x ** 2 * p ** 2) + np.pi * p * (self._u(x) - 1) + \
                       2 * p ** 2 * (self._u(x) * np.arctanh(self._u(x) ** -1) +
                                     np.log(0.5 * x))

            elif x * p < 1:

                func = np.log(0.25 * x ** 2 * p ** 2) + np.pi * p * (self._u(x) - 1) + \
                       2 * p ** 2 * (self._u(x) * np.arctanh(self._u(x) ** -1) +
                                     np.log(0.5 * x)) + 2 * self._g(x, p) * np.arctanh(self._g(x, p))

            else:

                func = np.log(0.25 * x ** 2 * p ** 2) + np.pi * p * (self._u(x) - 1) + \
                       2 * p ** 2 * (self._u(x) * np.arctanh(self._u(x) ** -1) +
                                     np.log(0.5 * x)) - 2 * self._f(x, p) * np.arctan(self._f(x, p))

        return func * prefactor
1006809
function QuadSD_ak1(NN, u, v, p, q, iPar){
                   // Divides p by the quadratic 1, u, v placing the quotient in q and the remainder in a, b
                   // iPar is a dummy variable for passing in the two parameters--a and b--by reference
                   q[0] = iPar.b = p[0];
                   q[1] = iPar.a = -(u*iPar.b) + p[1];

                   for (var i = 2; i < NN; i++){
                       q[i] = -(u*iPar.a + v*iPar.b) + p[i];
                       iPar.b = iPar.a;
                       iPar.a = q[i];
                   } 
                   return;
               }
1062243
def cal_pel(self, v, temp):
        
        if (self.eqn_el is None) or (self.params_el is None):
            return np.zeros_like(v)
        params = self._set_params(self.params_el)
        return func_el[self.eqn_el](v, temp, *params,
                                    self.n, self.z, t_ref=self.t_ref,
                                    three_r=self.three_r)
372405
function calculateNormalZProbability(z) {
    if (z == 0.0) {
        var x = 0.0;
    }
    else {
        var y = 0.5 * Math.abs(z);

        // Here comes the magic.

        if (y >= Z_MAX * 0.5) {
            x = 1.0;
        }
        else if (y < 1.0) {
            var w = y * y;

            x = 0.000124818987 * w -0.001075204047;
            x *= w +0.005198775019;
            x *= w -0.019198292004;
            x *= w +0.059054035642;
            x *= w -0.151968751364;
            x *= w +0.319152932694;
            x *= w -0.531923007300;
            x *= w +0.797884560593;
            x *= y * 2.0;
        }
        else {
            y -= 2.0;

            x = -0.000045255659 * y +0.000152529290;
            x *= y -0.000019538132;
            x *= y -0.000676904986;
            x *= y +0.001390604284;
            x *= y -0.000794620820;
            x *= y -0.002034254874;
            x *= y +0.006549791214;
            x *= y -0.010557625006;
            x *= y -0.010557625006;
            x *= y +0.011630447319;
            x *= y -0.009279453341;
            x *= y +0.005353579108;
            x *= y -0.002141268741;
            x *= y +0.000535310849;
            x *= y +0.999936657524;
        }
    }

    return z > 0.0 ? (x + 1.0) * 0.5 : (1.0 - x) * 0.5;
}