Finds strike and dip of plane given normal vector having components n, e,
    and u.

    Adapted from MATLAB script
    `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_
    written by Andy Michael and Oliver Boyd.
325777
def StrikeDip(n, e, u):
    
    r2d = 180 / np.pi
    if u < 0:
        n = -n
        e = -e
        u = -u

    strike = np.arctan2(e, n) * r2d
    strike = strike - 90
    while strike >= 360:
        strike = strike - 360
    while strike < 0:
        strike = strike + 360
    x = np.sqrt(np.power(n, 2) + np.power(e, 2))
    dip = np.arctan2(x, u) * r2d
    return (strike, dip)
540338
private void computeAngleOfRotation(PointTrack t, Vector3D_F64 pointingPlane) {
		VoTrack p = t.getCookie();

		// Compute ground pointing vector
		groundCurr.x = pointingPlane.z;
		groundCurr.y = -pointingPlane.x;
		double norm = groundCurr.norm();
		groundCurr.x /= norm;
		groundCurr.y /= norm;

		// dot product.  vectors are normalized to 1 already
		double dot = groundCurr.x * p.ground.x + groundCurr.y * p.ground.y;
		// floating point round off error some times knocks it above 1.0
		if (dot > 1.0)
			dot = 1.0;
		double angle = Math.acos(dot);

		// cross product to figure out direction
		if (groundCurr.x * p.ground.y - groundCurr.y * p.ground.x > 0)
			angle = -angle;

		farAngles.add(angle);
	}
1274661
def sdot( U, V ):
    '''
    
    '''
    nu = U.ndim
    #nv = V.ndim
    return np.tensordot( U, V, axes=(nu-1,0) )
330752
function intersectTriangle (out, pt, dir, tri) {
  var EPSILON = 0.000001
  var edge1 = [0, 0, 0]
  var edge2 = [0, 0, 0]
  var tvec = [0, 0, 0]
  var pvec = [0, 0, 0]
  var qvec = [0, 0, 0]

  vec3.subtract(edge1, tri[1], tri[0])
  vec3.subtract(edge2, tri[2], tri[0])

  vec3.cross(pvec, dir, edge2)
  var det = vec3.dot(edge1, pvec)

  if (det < EPSILON) return null
  vec3.subtract(tvec, pt, tri[0])
  var u = vec3.dot(tvec, pvec)
  if (u < 0 || u > det) return null
  vec3.cross(qvec, tvec, edge1)
  var v = vec3.dot(dir, qvec)
  if (v < 0 || u + v > det) return null

  var t = vec3.dot(edge2, qvec) / det
  out[0] = pt[0] + t * dir[0]
  out[1] = pt[1] + t * dir[1]
  out[2] = pt[2] + t * dir[2]
  return t
}
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 ;
}
443097
def fitPlaneSVD(XYZ):
    
    [rows,cols] = XYZ.shape
    # Set up constraint equations of the form  AB = 0,
    # where B is a column vector of the plane coefficients
    # in the form b(1)*X + b(2)*Y +b(3)*Z + b(4) = 0.
    p = (np.ones((rows,1)))
    AB = np.hstack([XYZ,p])
    [u, d, v] = np.linalg.svd(AB,0)        
    # Solution is last column of v.
    B = np.array(v[3,:])
    coeff = -B[[0, 1, 3]]/B[2]
    return coeff
448580
function(matrix, x0, y0)
    {
        var epsilon = 1e-8;
        var i;
        // Vertices of the quad with transform matrix applied
        var points = [];
        for (i = 0; i < 4; ++i)
            points[i] = WebInspector.Geometry.multiplyVectorByMatrixAndNormalize(new WebInspector.Geometry.Vector(this.vertices[i * 3], this.vertices[i * 3 + 1], this.vertices[i * 3 + 2]), matrix);
        // Calculating quad plane normal
        var normal = WebInspector.Geometry.crossProduct(WebInspector.Geometry.subtract(points[1], points[0]), WebInspector.Geometry.subtract(points[2], points[1]));
        // General form of the equation of the quad plane: A * x + B * y + C * z + D = 0
        var A = normal.x;
        var B = normal.y;
        var C = normal.z;
        var D = -(A * points[0].x + B * points[0].y + C * points[0].z);
        // Finding t from the equation
        var t = -(D + A * x0 + B * y0) / C;
        // Point of the intersection
        var pt = new WebInspector.Geometry.Vector(x0, y0, t);
        // Vectors from the intersection point to vertices of the quad
        var tVects = points.map(WebInspector.Geometry.subtract.bind(null, pt));
        // Intersection point lies inside of the polygon if scalar products of normal of the plane and
        // cross products of successive tVects are all nonstrictly above or all nonstrictly below zero
        for (i = 0; i < tVects.length; ++i) {
            var product = WebInspector.Geometry.scalarProduct(normal, WebInspector.Geometry.crossProduct(tVects[i], tVects[(i + 1) % tVects.length]));
            if (product < 0)
                return undefined;
        }
        return t;
    }
1319661
def project_and_to2dim(self, pps, plane_center):
        
        proj = self.projectionpoints(pps)
        [u1, u2, u3] = self.orthonormal_vectors()
        PP = np.array([[u1[0], u2[0], u3[0]],
                       [u1[1], u2[1], u3[1]],
                       [u1[2], u2[2], u3[2]]])
        xypps = list()
        for pp in proj:
            xyzpp = np.dot(pp, PP)
            xypps.append(xyzpp[0:2])
        if str(plane_center) == str('mean'):
            mean = np.zeros(2, np.float)
            for pp in xypps:
                mean += pp
            mean /= len(xypps)
            xypps = [pp - mean for pp in xypps]
        elif plane_center is not None:
            projected_plane_center = self.projectionpoints([plane_center])[0]
            xy_projected_plane_center = np.dot(projected_plane_center, PP)[0:2]
            xypps = [pp - xy_projected_plane_center for pp in xypps]
        return xypps
572434
function(vector$$1) {
                var V = vector$$1.elements || vector$$1;
                return Plane$1.create([
                  this.anchor.elements[0] + V[0],
                  this.anchor.elements[1] + V[1],
                  this.anchor.elements[2] + (V[2] || 0)
                ], this.normal);
              }
707636
function (plane, rightvector) {
  if (arguments.length < 2) {
    // choose an arbitrary right hand vector, making sure it is somewhat orthogonal to the plane normal:
    rightvector = vec3.random(plane)
  } else {
    rightvector = rightvector
  }
  this.v = vec3.unit(vec3.cross(plane, rightvector))
  this.u = vec3.cross(this.v, plane)
  this.plane = plane
  this.planeorigin = vec3.scale(plane[3], plane)
}
736779
def get_nodal_planes(self):
        
        # Convert reference frame to NED
        self.tensor, self.tensor_sigma = self._to_ned()
        self.ref_frame = 'NED'
        # Eigenvalue decomposition
        # Tensor
        _, evect = utils.eigendecompose(self.tensor)
        # Rotation matrix
        _, rot_vec = utils.eigendecompose(np.matrix([[0., 0., -1],
                                                     [0., 0., 0.],
                                                     [-1., 0., 0.]]))
        rotation_matrix = (np.matrix(evect * rot_vec.T)).T
        if np.linalg.det(rotation_matrix) < 0.:
            rotation_matrix *= -1.
        flip_dc = np.matrix([[0., 0., -1.],
                             [0., -1., 0.],
                             [-1., 0., 0.]])
        rotation_matrices = sorted(
            [rotation_matrix, flip_dc * rotation_matrix],
            cmp=cmp_mat)
        nodal_planes = GCMTNodalPlanes()
        dip, strike, rake = [(180. / pi) * angle
                             for angle in utils.matrix_to_euler(rotation_matrices[0])]
        # 1st Nodal Plane
        nodal_planes.nodal_plane_1 = {'strike': strike % 360,
                                      'dip': dip,
                                      'rake': -rake}

        # 2nd Nodal Plane
        dip, strike, rake = [(180. / pi) * angle
                             for angle in utils.matrix_to_euler(rotation_matrices[1])]
        nodal_planes.nodal_plane_2 = {'strike': strike % 360.,
                                      'dip': dip,
                                      'rake': -rake}
        return nodal_planes
228397
def _principal_angle(a, B):
  ''''''
  # TODO: check case for t = d-1
  if a.shape[1] == 1:
    return a.T.dot(B)[0,:,0]

  # find normals that maximize distance when projected
  x1 = np.einsum('abc,adc->abd', B, B).dot(a) - a   # b.dot(b.T).dot(a) - a
  x2 = np.einsum('ab,cad->cbd', a.dot(a.T), B) - B  # a.dot(a.T).dot(b) - b
  xx = np.vstack((x1, x2))

  # batch PCA (1st comp. only)
  xx -= xx.mean(axis=1)[:,None]
  c = np.einsum('abc,abd->acd', xx, xx)
  _, vecs = np.linalg.eigh(c)
  fpc = vecs[:,:,-1]
  fpc1 = fpc[:len(x1)]
  fpc2 = fpc[len(x1):]

  # a.dot(fpc1).dot(b.dot(fpc2))
  lhs = a.dot(fpc1.T).T
  rhs = np.einsum('abc,ac->ab', B, fpc2)
  return np.einsum('ij,ij->i', lhs, rhs)
570022
def calc_ABDE_from_lamination_parameters(self):
        
        # dummies used to unpack vector results
        du1, du2, du3, du4, du5, du6 = 0, 0, 0, 0, 0, 0
        # A matrix terms
        A11,A22,A12, du1,du2,du3, A66,A16,A26 =\
            (self.h       ) * np.dot(self.matobj.u, self.xiA)
        # B matrix terms
        B11,B22,B12, du1,du2,du3, B66,B16,B26 =\
            (self.h**2/4. ) * np.dot(self.matobj.u, self.xiB)
        # D matrix terms
        D11,D22,D12, du1,du2,du3, D66,D16,D26 =\
            (self.h**3/12.) * np.dot(self.matobj.u, self.xiD)
        # E matrix terms
        du1,du2,du3, E44,E55,E45, du4,du5,du6 =\
            (self.h       ) * np.dot(self.matobj.u, self.xiE)

        self.A = np.array([[A11, A12, A16],
                           [A12, A22, A26],
                           [A16, A26, A66]], dtype=np.float64)

        self.B = np.array([[B11, B12, B16],
                           [B12, B22, B26],
                           [B16, B26, B66]], dtype=np.float64)

        self.D = np.array([[D11, D12, D16],
                           [D12, D22, D26],
                           [D16, D26, D66]], dtype=np.float64)

        # printing E acoordingly to Reddy definition for E44, E45 and E55
        self.E = np.array([[E55, E45],
                           [E45, E44]], dtype=np.float64)

        self.ABD = np.array([[A11, A12, A16, B11, B12, B16],
                             [A12, A22, A26, B12, B22, B26],
                             [A16, A26, A66, B16, B26, B66],
                             [B11, B12, B16, D11, D12, D16],
                             [B12, B22, B26, D12, D22, D26],
                             [B16, B26, B66, D16, D26, D66]], dtype=np.float64)

        # printing ABDE acoordingly to Reddy definition for E44, E45 and E55
        self.ABDE = np.array([[A11, A12, A16, B11, B12, B16, 0, 0],
                              [A12, A22, A26, B12, B22, B26, 0, 0],
                              [A16, A26, A66, B16, B26, B66, 0, 0],
                              [B11, B12, B16, D11, D12, D16, 0, 0],
                              [B12, B22, B26, D12, D22, D26, 0, 0],
                              [B16, B26, B66, D16, D26, D66, 0, 0],
                              [0, 0, 0, 0, 0, 0, E55, E45],
                              [0, 0, 0, 0, 0, 0, E45, E44]],
                               dtype=np.float64)
352938
def xstep(self):
        
        """

        ngsit = 0
        gsrrs = np.inf
        YU = self.Y - self.U
        SYU = self.S + YU[..., -1]
        YU[..., -1] = 0.0
        ATYU = self.cnst_AT(YU)
        while gsrrs > self.opt['GSTol'] and ngsit < self.opt['MaxGSIter']:
            self.X = self.GaussSeidelStep(
                SYU, self.X, ATYU, 1.0, self.lcw, 1.0)
            gsrrs = sl.rrs(
                self.cnst_AT(self.cnst_A(self.X)),
                self.cnst_AT(self.cnst_c() - self.cnst_B(self.Y) - self.U)
            )
            ngsit += 1

        self.xs = (ngsit, gsrrs)
131394
def integrator(integrand,xmin,xmax,n_points,factor=2):
    '''
    
    '''
    integral_vector = np.empty([n_points+1])
    dx  = (xmax-xmin)/n_points
    # integrate
    for i in xrange(n_points+1):
        xnow = xmin + i * dx
        integral, error = integrate.quad(rotate(integrand,
                xnow), xmin*factor, xmax*factor)
        integral_vector[i] = integral
    # normalize
    normalization = np.average(integral_vector)*(xmax-xmin)
    normalized_vector = integral_vector/normalization
    return normalized_vector
1360213
function computeSplineControlPoints(K)
{
	p1=[];
	p2=[];
	n = K.length-1;
	
	/*rhs vector*/
	a=[];
	b=[];
	c=[];
	r=[];
	
	/*left most segment*/
	a[0]=0;
	b[0]=2;
	c[0]=1;
	r[0] = K[0]+2*K[1];
	
	/*internal segments*/
	for (i = 1; i < n - 1; i++)
	{
		a[i]=1;
		b[i]=4;
		c[i]=1;
		r[i] = 4 * K[i] + 2 * K[i+1];
	}
			
	/*right segment*/
	a[n-1]=2;
	b[n-1]=7;
	c[n-1]=0;
	r[n-1] = 8*K[n-1]+K[n];
	
	/*solves Ax=b with the Thomas algorithm (from Wikipedia)*/
	for (i = 1; i < n; i++)
	{
		m = a[i]/b[i-1];
		b[i] = b[i] - m * c[i - 1];
		r[i] = r[i] - m*r[i-1];
	}
 
	p1[n-1] = r[n-1]/b[n-1];
	for (i = n - 2; i >= 0; --i)
		p1[i] = (r[i] - c[i] * p1[i+1]) / b[i];
		
	/*we have p1, now compute p2*/
	for (i=0;i<n-1;i++)
		p2[i]=2*K[i+1]-p1[i+1];
	
	p2[n-1]=0.5*(K[n]+p1[n-1]);
	
	return {p1:p1, p2:p2};
}
533951
public Vector4d frustumPlane(int plane, Vector4d dest) {
        switch (plane) {
        case PLANE_NX:
            dest.set(m03 + m00, m13 + m10, m23 + m20, m33 + m30).normalize3(dest);
            break;
        case PLANE_PX:
            dest.set(m03 - m00, m13 - m10, m23 - m20, m33 - m30).normalize3(dest);
            break;
        case PLANE_NY:
            dest.set(m03 + m01, m13 + m11, m23 + m21, m33 + m31).normalize3(dest);
            break;
        case PLANE_PY:
            dest.set(m03 - m01, m13 - m11, m23 - m21, m33 - m31).normalize3(dest);
            break;
        case PLANE_NZ:
            dest.set(m03 + m02, m13 + m12, m23 + m22, m33 + m32).normalize3(dest);
            break;
        case PLANE_PZ:
            dest.set(m03 - m02, m13 - m12, m23 - m22, m33 - m32).normalize3(dest);
            break;
        default:
            throw new IllegalArgumentException("plane"); //$NON-NLS-1$
        }
        return dest;
    }