phase angle from cos elongation and distances
@private
@param {number} Δ
@param {number} R
@param {number} cψ
@returns {number}
169492
function pa (Δ, R, cψ) {
  const sψ = Math.sin(Math.acos(cψ))
  let i = Math.atan(R * sψ / (Δ - R * cψ))
  if (i < 0) {
    i += Math.PI
  }
  return i
}
247028
public float angleCos(Vector3fc v) {
        double length1Squared = x * x + y * y + z * z;
        double length2Squared = v.x() * v.x() + v.y() * v.y() + v.z() * v.z();
        double dot = x * v.x() + y * v.y() + z * v.z();
        return (float) (dot / (Math.sqrt(length1Squared * length2Squared)));
    }
784751
def beta_r(self, r, kwargs):
        
        if self._type == 'const':
            return self.const_beta(kwargs)
        elif self._type == 'OsipkovMerritt':
            return self.ospikov_meritt(r, kwargs)
        elif self._type == 'Colin':
            return self.colin(r, kwargs)
        elif self._type == 'isotropic':
            return self.isotropic()
        elif self._type == 'radial':
            return self.radial()
        else:
            raise ValueError('anisotropy type %s not supported!' % self._type)
443947
def C(w, Xs):
    '''
    '''
    n = len(Xs)
    P = projection_matrix(w)
    Ys = [np.dot(P, X) for X in Xs]
    A = calc_A(Ys)
    A_hat = calc_A_hat(A, skew_matrix(w))

    return np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))
846994
function(angle, origin) {
        angle *= Math.PI / 180;
        var radius = this.distanceTo(origin);
        var theta = angle + Math.atan2(this.y - origin.y, this.x - origin.x);
        this.x = origin.x + (radius * Math.cos(theta));
        this.y = origin.y + (radius * Math.sin(theta));
        this.clearBounds();
    }
539829
function angle(a, b) {
    var tempA = fromValues(a[0], a[1], a[2])
    var tempB = fromValues(b[0], b[1], b[2])
 
    normalize(tempA, tempA)
    normalize(tempB, tempB)
 
    var cosine = dot(tempA, tempB)

    if(cosine > 1.0){
        return 0
    } else {
        return Math.acos(cosine)
    }     
}
1088823
def ZernikeSequence(x,y,maxRadial=5,epsilon=1e-10):
    
    if maxRadial>5:
        raise ValueError('Code for higher radial orders not implemented')
    # Derive radius and exp(i*theta) 
    temp = x + 1j*y
    r1 = np.abs(temp)
    e1 = temp/(r1+epsilon)
     # Generate powers of r recursively
    r2 = r1*r1
    r3 = r2*r1
    r4 = r3*r1
    r5 = r4*r1
    # Generate cos and sin terms recursively from exp(i*theta)
    e2 = e1*e1
    e3 = e2*e1
    e4 = e3*e1
    e5 = e4*e1
    ctheta = e1.real
    stheta = e1.imag
    c2theta = e2.real
    s2theta = e2.imag
    c3theta = e3.real
    s3theta = e3.imag
    c4theta = e4.real
    s4theta = e4.imag
    c5theta = e5.real
    s5theta = e5.imag
    # Generate all the zernikes
    zernike = np.zeros((21,)+x.shape )
    zernike[0] =  1.0
    zernike[1] =  2.0*r1*ctheta
    zernike[2] =  2.0*r1*stheta
    zernike[3] =  sqrt(3.0)*(2.0*r2 - 1.0)
    zernike[4] =  sqrt(6.0)*r2*s2theta
    zernike[5] =  sqrt(6.0)*r2*c2theta
    zernike[6] =  sqrt(8.0)*(3.0*r3 - 2.0*r1)*stheta
    zernike[7] =  sqrt(8.0)*(3.0*r3 - 2.0*r1)*ctheta
    zernike[8] =  sqrt(8.0)*r3*s3theta
    zernike[9] = sqrt(8.0)*r3*c3theta
    zernike[10] = sqrt(5.0)*(6.*r4 - 6.*r2 + 1.)
    zernike[11] = sqrt(10.)*(4.*r4 - 3.*r2)*c2theta
    zernike[12] = sqrt(10.)*(4.*r4 - 3.*r2)*s2theta
    zernike[13] = sqrt(10.)*r4*c4theta
    zernike[14] = sqrt(10.)*r4*s4theta
    zernike[15] = sqrt(12.)*(10*r5-12*r3+3*r1)*ctheta
    zernike[16] = sqrt(12.)*(10*r5-12*r3+3*r1)*stheta
    zernike[17] = sqrt(12.)*(5*r5-4*r3)*c3theta
    zernike[18] = sqrt(12.)*(5*r5-4*r3)*s3theta
    zernike[19] = sqrt(12.)*r5*c5theta
    zernike[20] = sqrt(12.)*r5*s5theta
    # Make zernike zero outside unit circle (useful for dot product)
    zernike = zernike*np.less_equal(r1, 1.0)
    return(zernike[:NumZernike(maxRadial)])
1072450
def rpole_to_pot_aligned(rpole, sma, q, F, d, component=1):
    
    q = q_for_component(q, component=component)
    rpole_ = np.array([0, 0, rpole/sma])
    logger.debug("libphoebe.roche_Omega(q={}, F={}, d={}, rpole={})".format(q, F, d, rpole_))
    pot = libphoebe.roche_Omega(q, F, d, rpole_)
    return pot_for_component(pot, component, reverse=True)
463554
def cos_r(self, N=None):  # percent=0.9
		""""""

		if not hasattr(self, 'F') or self.F.shape[1] < self.rank:
			self.fs_r(N=self.rank)  # generate F
		self.dr = norm(self.F, axis=1)**2
		# cheaper than diag(self.F.dot(self.F.T))?

		return apply_along_axis(lambda _: _/self.dr, 0, self.F[:, :N]**2)
365773
def R_r(self, r):
        
        phi = np.random.uniform(0, 2*np.pi)
        theta = np.random.uniform(0, np.pi)
        x = r * np.sin(theta) * np.cos(phi)
        y = r * np.sin(theta) * np.sin(phi)
        R = np.sqrt(x**2 + y**2)
        return R, x, y
44087
def derivative(self, t, n=1):
        """"""
        angle = radians(self.theta + t*self.delta)
        phi = radians(self.rotation)
        rx = self.radius.real
        ry = self.radius.imag
        k = (self.delta*2*pi/360)**n  # ((d/dt)angle)**n

        if n % 4 == 0 and n > 0:
            return rx*cos(phi)*cos(angle) - ry*sin(phi)*sin(angle) + 1j*(
                rx*sin(phi)*cos(angle) + ry*cos(phi)*sin(angle))
        elif n % 4 == 1:
            return k*(-rx*cos(phi)*sin(angle) - ry*sin(phi)*cos(angle) + 1j*(
                -rx*sin(phi)*sin(angle) + ry*cos(phi)*cos(angle)))
        elif n % 4 == 2:
            return k*(-rx*cos(phi)*cos(angle) + ry*sin(phi)*sin(angle) + 1j*(
                -rx*sin(phi)*cos(angle) - ry*cos(phi)*sin(angle)))
        elif n % 4 == 3:
            return k*(rx*cos(phi)*sin(angle) + ry*sin(phi)*cos(angle) + 1j*(
                rx*sin(phi)*sin(angle) - ry*cos(phi)*cos(angle)))
        else:
            raise ValueError("n should be a positive integer.")
190706
def coordinate(i, j)
      x = @pos.x + (i * @delta_col * @cosines[0]) + (j * @delta_row * @cosines[3])
      y = @pos.y + (i * @delta_col * @cosines[1]) + (j * @delta_row * @cosines[4])
      z = @pos.z + (i * @delta_col * @cosines[2]) + (j * @delta_row * @cosines[5])
      Coordinate.new(x, y, z)
    end
984007
function(inPositions) {
			var p = inPositions;
			// yay math!, rad -> deg
			var a = Math.asin(p.y / p.h) * (180 / Math.PI);
			// fix for range limits of asin (-90 to 90)
			// Quadrants II and III
			if (p.x < 0) {
				a = 180 - a;
			}
			// Quadrant IV
			if (p.x > 0 && p.y < 0) {
				a += 360;
			}
			return a;
		}
979174
private void calculateRotationalAxis(Matrix rotation, Atom translation,
			double c) {
		// Calculate magnitude of rotationAxis components, but not signs
		double sum=0;
		double[] rotAx = new double[3];
		for(int i=0;i<3;i++) {
			rotAx[i] = Math.sqrt(rotation.get(i, i)-c);
			sum+=rotAx[i]*rotAx[i];
		}
		for(int i=0;i<3;i++) {
			rotAx[i] /= Math.sqrt(sum);
		}

		// Now determine signs
		double d0 = rotation.get(2,1)-rotation.get(1,2); //=2u[0]*sin(theta)
		double d1 = rotation.get(0,2)-rotation.get(2,0); //=2u[1]*sin(theta)
		double d2 = rotation.get(1,0)-rotation.get(0,1); //=2u[2]*sin(theta)

		double s12 = rotation.get(2,1)+rotation.get(1,2); //=2*u[1]*u[2]*(1-cos(theta))
		double s02 = rotation.get(0,2)+rotation.get(2,0); //=2*u[0]*u[2]*(1-cos(theta))
		double s01 = rotation.get(1,0)+rotation.get(0,1); //=2*u[0]*u[1]*(1-cos(theta))

		//Only use biggest d for the sign directly, for numerical stability around 180deg
		if( Math.abs(d0) < Math.abs(d1) ) { // not d0
			if( Math.abs(d1) < Math.abs(d2) ) { //d2
				if(d2>=0){ //u[2] positive
					if( s02 < 0 ) rotAx[0] = -rotAx[0];
					if( s12 < 0 ) rotAx[1] = -rotAx[1];
				} else { //u[2] negative
					rotAx[2] = -rotAx[2];
					if( s02 >= 0 ) rotAx[0] = -rotAx[0];
					if( s12 >= 0 ) rotAx[1] = -rotAx[1];
				}
			} else { //d1
				if(d1>=0) {//u[1] positive
					if( s01 < 0) rotAx[0] = -rotAx[0];
					if( s12 < 0) rotAx[2] = -rotAx[2];
				} else { //u[1] negative
					rotAx[1] = -rotAx[1];
					if( s01 >= 0) rotAx[0] = -rotAx[0];
					if( s12 >= 0) rotAx[2] = -rotAx[2];
				}
			}
		} else { // not d1
			if( Math.abs(d0) < Math.abs(d2) ) { //d2
				if(d2>=0){ //u[2] positive
					if( s02 < 0 ) rotAx[0] = -rotAx[0];
					if( s12 < 0 ) rotAx[1] = -rotAx[1];
				} else { //u[2] negative
					rotAx[2] = -rotAx[2];
					if( s02 >= 0 ) rotAx[0] = -rotAx[0];
					if( s12 >= 0 ) rotAx[1] = -rotAx[1];
				}
			} else { //d0
				if(d0>=0) { //u[0] positive
					if( s01 < 0 ) rotAx[1] = -rotAx[1];
					if( s02 < 0 ) rotAx[2] = -rotAx[2];
				} else { //u[0] negative
					rotAx[0] = -rotAx[0];
					if( s01 >= 0 ) rotAx[1] = -rotAx[1];
					if( s02 >= 0 ) rotAx[2] = -rotAx[2];
				}
			}
		}

		rotationAxis = new AtomImpl();
		rotationAxis.setCoords(rotAx);

		// Calculate screw = (rotationAxis dot translation)*u
		double dotProduct = Calc.scalarProduct(rotationAxis, translation);

		screwTranslation = Calc.scale(rotationAxis, dotProduct);
		otherTranslation = Calc.subtract(translation, screwTranslation);

		Atom hypot = Calc.vectorProduct(otherTranslation,rotationAxis);
		Calc.scaleEquals(hypot,.5/Math.tan(theta/2.0));

		// Calculate rotation axis position
		rotationPos = Calc.scaleAdd(.5,otherTranslation, hypot);
	}
687079
def _dK_dR(self, R):
        """"""
        return -self._ns * self._N / R**2 / self._sin_alpha
1293319
def K(self, r, R, kwargs):
        
        u = r / R
        if np.min(u) < 1:
            raise ValueError("3d radius is smaller than projected radius! Does not make sense.")

        if self._type == 'const_wrong':
            beta = kwargs['beta']
            k = 1./2. * u**(2*beta - 1.) * ((3./2 - beta) * np.sqrt(np.pi) * special.gamma(beta - 1./2)/special.gamma(beta)
                        + beta * self._B(x=1./u**2, a=beta+1./2, b=1./2) - self._B(x=1./u**2, a=beta-1./2, b=1./2))
        elif self._type == 'const':
            beta = kwargs['beta']
            k = np.sqrt(1 - 1./u**2) / (1. - 2*beta) + np.sqrt(np.pi)/2 * special.gamma(beta - 1./2)/special.gamma(beta)\
                * (3./2 - beta) * u**(2*beta - 1.) * (1 - special.betainc(beta+1./2, 1./2, 1./u**2))
        elif self._type == 'isotropic':
            k = np.sqrt(1 - 1./u**2)
        elif self._type == 'radial':
            k = np.pi/4 * u - 1./2 * np.sqrt(1 - 1./u**2) - u/2. * np.arcsin(1./u)
        elif self._type == 'OsipkovMerritt':
            r_ani = kwargs['r_ani']
            ua = r_ani / R
            k = (ua**2 + 1./2) / (ua**2 + 1)**(3./2) * (u**2 + ua**2) / u * np.arctan(np.sqrt((u**2 - 1) / (ua**2 + 1))) \
                - 1./2/(ua**2 + 1) * np.sqrt(1 -1./u**2)
        elif self._type == 'Colin':
            r_ani = kwargs['r_ani']
            ua = r_ani / R
            if ua == 1:
                k = (1 + 1./u) * np.arccosh(u) - 1./6 * (8./u + 7) * np.sqrt((u-1.)/(u+1.))
            else:
                k = 0.5 / (ua**2 - 1) * np.sqrt(1 - 1./u**2) + (1. + ua/u) * np.cosh(u) - np.sign(ua - 1) * ua * \
                            (ua**2 - 0.5) / np.abs(ua**2-1)**(3./2) * (1. + ua/u) * np.arccosh((ua*u + 1)/(u + ua))
        else:
            raise ValueError('anisotropy type %s not supported!' % self._type)
        return k
529664
function angleFromThreePoints (ax, ay, bx, by, cx, cy) {
  var c = distance(ax, ay, bx, by)
  var a = distance(bx, by, cx, cy)
  var b = distance(ax, ay, cx, cy)
  return Math.acos((a * a + c * c - b * b) / (2 * a * c))
}