Python math.acos() Examples

def _get_skew(corners, board):
    """
    Get skew for given checkerboard detection.
    Scaled to [0,1], which 0 = no skew, 1 = high skew
    Skew is proportional to the divergence of three outside corners from 90 degrees.
    """
    # TODO Using three nearby interior corners might be more robust, outside corners occasionally
    # get mis-detected
    up_left, up_right, down_right, _ = _get_outside_corners(corners, board)

    def angle(a, b, c):
        """
        Return angle between lines ab, bc
        """
        ab = a - b
        cb = c - b
        return math.acos(numpy.dot(ab,cb) / (numpy.linalg.norm(ab) * numpy.linalg.norm(cb)))

    skew = min(1.0, 2. * abs((math.pi / 2.) - angle(up_left, up_right, down_right)))
    return skew 

def align_bone_z_axis(obj, bone, vec):
    """ Rolls the bone to align its z-axis as closely as possible to
        the given vector.
        Must be in edit mode.
    """
    bone_e = obj.data.edit_bones[bone]

    vec = bone_e.y_axis.cross(vec)
    vec.normalize()

    dot = max(-1.0, min(1.0, bone_e.x_axis.dot(vec)))
    angle = math.acos(dot)

    bone_e.roll += angle

    dot1 = bone_e.x_axis.dot(vec)

    bone_e.roll -= angle * 2

    dot2 = bone_e.x_axis.dot(vec)

    if dot1 > dot2:
        bone_e.roll += angle * 2 

def interpolate(self, other, this_weight):
        q0, q1 = np.roll(self.q, shift=1), np.roll(other.q, shift=1)
        u = 1 - this_weight
        assert(u >= 0 and u <= 1)
        cos_omega = np.dot(q0, q1)

        if cos_omega < 0:
            result = -q0[:]
            cos_omega = -cos_omega
        else:
            result = q0[:]

        cos_omega = min(cos_omega, 1)

        omega = math.acos(cos_omega)
        sin_omega = math.sin(omega)
        a = math.sin((1-u) * omega)/ sin_omega
        b = math.sin(u * omega) / sin_omega

        if abs(sin_omega) < 1e-6:
            # direct linear interpolation for numerically unstable regions
            result = result * this_weight + q1 * u
            result /= math.sqrt(np.dot(result, result))
        else:
            result = result*a + q1*b
        return Quaternion(np.roll(result, shift=-1))

    # To conversions 

def get_angle_formed_by(p1,p2,p3): # angle formed by three positions in space

    # based on code submitted by Paul Sherwood
    r1 = distance(p1,p2)
    r2 = distance(p2,p3)
    r3 = distance(p1,p3)
    
    small = 1.0e-10
    
    if (r1 + r2 - r3) < small:
        # This seems to happen occasionally for 180 angles 
        theta = math.pi
    else:
        theta = math.acos( (r1*r1 + r2*r2  - r3*r3) / (2.0 * r1*r2) )
    return theta;

#------------------------------------------------------------------------------ 

def setPolygon(self):
        '''Calculate position and rotation of the arc arrow head.'''
        rotDeg = 0            
        xlength = self.pos1.x() - self.pos2.x()
        ylength = self.pos1.y() - self.pos2.y()
        d = math.sqrt( math.pow( xlength , 2) + math.pow( ylength , 2) )        
        if d > 0:
            beta = math.acos( xlength / d )
            rotDeg = math.degrees( beta ) 
        
        self.arrowPolygonObject.setPolygon( QtGui.QPolygonF( [
                                                 QtCore.QPointF( (self.pos2.x() -10),  (self.pos2.y() +5)), 
                                                 QtCore.QPointF( (self.pos2.x() -10) , (self.pos2.y() -5)), 
                                                 QtCore.QPointF(       self.pos2.x() ,      self.pos2.y())
                                                 ] ) ) 
        
        self.arrowPolygonObject.setBrush( QtGui.QBrush(QtCore.Qt.black) )
        
        """ self.angle()!!!!!!!!!"""
#         self.arcLinePolygon.angle()
#         self.arcLinePolygon.rotate(rotDeg)
#         self.arcLinePolygon.setPos( self.pos2 ) 

    #------------------------------------------------------------------------------------------------ 

def setPolygon(self):
        rotDeg = 0            
        xlength = self.pos1.x() - self.pos2.x()
        ylength = self.pos1.y() - self.pos2.y()
        d = math.sqrt( math.pow( xlength , 2) + math.pow( ylength , 2) )        
        if d > 0:
            beta = math.acos( xlength / d )
            rotDeg = math.degrees( beta ) 
        
        self.arcLinePolygon.setPolygon( QtGui.QPolygonF( [
                                                 QtCore.QPointF( (self.pos2.x() -10),  (self.pos2.y() +5)), 
                                                 QtCore.QPointF( (self.pos2.x() -10) , (self.pos2.y() -5)), 
                                                 QtCore.QPointF(       self.pos2.x() ,      self.pos2.y())
                                                 ] ) ) 
        
        self.arcLinePolygon.setBrush( QtGui.QBrush(QtCore.Qt.black) )
        
        """ self.angle()!!!!!!!!!"""
#         self.arcLinePolygon.angle()
#         self.arcLinePolygon.rotate(rotDeg)
#         self.arcLinePolygon.setPos( self.pos2 ) 
    #------------------------------------------------------------------------------------------------ 

def eLowDotOperator(stack, z, mode):
	if mode == 1:   # num
		stack.append(utilities.formatNum(math.acos(z)))
	#elif mode == 2:
	elif mode == 3: # str or list
		if len(z) == 0:
			stack.append([])
		else:
			result = ""
			for i in z:
				i = utilities.castToList(i)
				if len(i) >= 2:
					if type(i[1]) == str and type(result) == str:
						result += (i[1] * utilities.castToNumber(i[0]))
					else:
						result = list(result)
						result += [i[1]] * utilities.castToNumber(i[0])
			stack.append(result)
	else:
		monadNotImplemented(mode, '')

# ? 

def get_max_width_for_circles(self, rad1, rad2, max_centering_proportion):
        r"""
        (the "r" in the above line is to keep pylint happy)
           __
          /__\ <- compute the line width which is drawable between 2 circles.
         /  _ \   max_centering_proportion : 0, touching the circle1, 1,
        |  |_| |   touching the circle2, 0.5 : middle between the 2 circles
        |      |
         \    /
          \__/
                 basically, max_centering_proportion is
                 max_centering_proportion/nb_lines
        """
        # radius at the center of the 2 circles
        rmid = rad2 - (rad2-rad1)*max_centering_proportion
        return sin(acos(rmid/rad2)) * rad2 * 2 

def _Angle(u, v):
    """Return angle between two vectors.

    Args:
      u: (float, float)
      v: (float, float)
    Returns:
      float - angle in radians between u and v, where
        it is +/- depending on sign of ux * vy - uy * vx
    """

    (ux, uy) = u
    (vx, vy) = v
    costheta = (ux * vx + uy * vy) / \
        (math.sqrt(ux ** 2 + uy ** 2) * math.sqrt(vx ** 2 + vy ** 2))
    if costheta > 1.0:
        costheta = 1.0
    if costheta < -1.0:
        costheta = -1.0
    theta = math.acos(costheta)
    if ux * vy - uy * vx < 0.0:
        theta = -theta
    return theta 

def Angle(a, b, c, points):
    """Return Angle abc in degrees, in range [0,180),
    where a,b,c are indices into points."""

    u = Sub2(points.pos[c], points.pos[b])
    v = Sub2(points.pos[a], points.pos[b])
    n1 = Length2(u)
    n2 = Length2(v)
    if n1 == 0.0 or n2 == 0.0:
        return 0.0
    else:
        costheta = Dot2(u, v) / (n1 * n2)
        if costheta > 1.0:
            costheta = 1.0
        if costheta < - 1.0:
            costheta = - 1.0
        return math.acos(costheta) * 180.0 / math.pi 

def proj_xy(self, t, next=None):
        """
            length of projection of sections at crossing line / circle intersections
            deformation unit vector for profil in xy axis
            so f(x_profile) = position of point in xy plane
        """
        if next is None:
            return self.normal(t).v.normalized(), 1
        v0 = self.normal(1).v.normalized()
        v1 = next.normal(0).v.normalized()
        direction = v0 + v1
        adj = (v0 * self.length) * (v1 * next.length)
        hyp = (self.length * next.length)
        c = min(1, max(-1, adj / hyp))
        size = 1 / cos(0.5 * acos(c))
        return direction.normalized(), min(3, size) 

def Encode(self, normal):
        x = normal[0]
        y = normal[1]
        z = normal[2]
        # normalize
        l = math.sqrt((x*x) + (y*y) + (z*z))
        if l == 0:
            return 0
        x = x/l
        y = y/l
        z = z/l

        if (x == 0.0) & (y == 0.0) :
            if z > 0.0:
                return 0
            else:
                return (128 << 8)

        lng = math.acos(z) * 255 / (2 * math.pi)
        lat = math.atan2(y, x) * 255 / (2 * math.pi)
        retval = ((int(lat) & 0xFF) << 8) | (int(lng) & 0xFF)
        return retval 

def Angle(a, b, c, points):
    """Return Angle abc in degrees, in range [0,180),
    where a,b,c are indices into points."""

    u = Sub2(points.pos[c], points.pos[b])
    v = Sub2(points.pos[a], points.pos[b])
    n1 = Length2(u)
    n2 = Length2(v)
    if n1 == 0.0 or n2 == 0.0:
        return 0.0
    else:
        costheta = Dot2(u, v) / (n1 * n2)
        if costheta > 1.0:
            costheta = 1.0
        if costheta < - 1.0:
            costheta = - 1.0
        return math.acos(costheta) * 180.0 / math.pi 

def set_mat(obj, bone_name, matrix):
        """ Sets the bone to have the given transform matrix.
        """
        a = obj.data.edit_bones[bone_name]

        a.head = (0, 0, 0)
        a.tail = (0, 1, 0)

        a.transform(matrix)

        d = acos(a.matrix.to_quaternion().dot(matrix.to_quaternion())) * 2.0

        roll_1 = a.roll + d
        roll_2 = a.roll - d

        a.roll = roll_1
        d1 = a.matrix.to_quaternion().dot(matrix.to_quaternion())
        a.roll = roll_2
        d2 = a.matrix.to_quaternion().dot(matrix.to_quaternion())

        if d1 > d2:
            a.roll = roll_1
        else:
            a.roll = roll_2 

def angle_on_plane(plane, vec1, vec2):
    """ Return the angle between two vectors projected onto a plane.
    """
    plane.normalize()
    vec1 = vec1 - (plane * (vec1.dot(plane)))
    vec2 = vec2 - (plane * (vec2.dot(plane)))
    vec1.normalize()
    vec2.normalize()

    # Determine the angle
    angle = math.acos(max(-1.0, min(1.0, vec1.dot(vec2))))

    if angle < 0.00001:  # close enough to zero that sign doesn't matter
        return angle

    # Determine the sign of the angle
    vec3 = vec2.cross(vec1)
    vec3.normalize()
    sign = vec3.dot(plane)
    if sign >= 0:
        sign = 1
    else:
        sign = -1

    return angle * sign 

def align_bone_x_axis(obj, bone, vec):
    """ Rolls the bone to align its x-axis as closely as possible to
        the given vector.
        Must be in edit mode.
    """
    bone_e = obj.data.edit_bones[bone]

    vec = vec.cross(bone_e.y_axis)
    vec.normalize()

    dot = max(-1.0, min(1.0, bone_e.z_axis.dot(vec)))
    angle = math.acos(dot)

    bone_e.roll += angle

    dot1 = bone_e.z_axis.dot(vec)

    bone_e.roll -= angle * 2

    dot2 = bone_e.z_axis.dot(vec)

    if dot1 > dot2:
        bone_e.roll += angle * 2 

def set_mat(obj, bone_name, matrix):
        """ Sets the bone to have the given transform matrix.
        """
        a = obj.data.edit_bones[bone_name]

        a.head = (0, 0, 0)
        a.tail = (0, 1, 0)

        a.transform(matrix)

        d = acos(a.matrix.to_quaternion().dot(matrix.to_quaternion())) * 2.0

        roll_1 = a.roll + d
        roll_2 = a.roll - d

        a.roll = roll_1
        d1 = a.matrix.to_quaternion().dot(matrix.to_quaternion())
        a.roll = roll_2
        d2 = a.matrix.to_quaternion().dot(matrix.to_quaternion())

        if d1 > d2:
            a.roll = roll_1
        else:
            a.roll = roll_2 

def angle_on_plane(plane, vec1, vec2):
    """ Return the angle between two vectors projected onto a plane.
    """
    plane.normalize()
    vec1 = vec1 - (plane * (vec1.dot(plane)))
    vec2 = vec2 - (plane * (vec2.dot(plane)))
    vec1.normalize()
    vec2.normalize()

    # Determine the angle
    angle = math.acos(max(-1.0, min(1.0, vec1.dot(vec2))))

    if angle < 0.00001:  # close enough to zero that sign doesn't matter
        return angle

    # Determine the sign of the angle
    vec3 = vec2.cross(vec1)
    vec3.normalize()
    sign = vec3.dot(plane)
    if sign >= 0:
        sign = 1
    else:
        sign = -1

    return angle * sign 

def get_sph_cor(x,y,z):    #using radian
    lon=0
    if(x==0 and y>0):
        lon=PI/2
    elif(x==0 and y<0):
        lon=3*PI/2
    elif(x==0 and y==0):
        print ("error")
        return
    elif(x>0 and y==0):
        lon=0
    elif(x>0 and y>0):
        lon=math.atan(float(y)/float(x))
    elif(x>0 and y<0):
        lon=2*PI+math.atan(float(y)/float(x))
    elif(x<0 and y==0):
        lon=PI
    elif(x<0 and y>0):
        lon=PI+math.atan(float(y)/float(x))
    elif(x<0 and y<0):
        lon=PI+math.atan(float(y)/float(x))

    lat=PI/2-math.acos(z/1.0)

    return lon,lat 

def calc_rotation_matrix(q1, q2, ref_q1, ref_q2):
  ref_nv = np.cross(ref_q1, ref_q2) 
  q_nv = np.cross(q1, q2)
  if min(norm(ref_nv), norm(q_nv)) == 0.:  # avoid 0 degree including angle
    return np.identity(3)
  axis = np.cross(ref_nv, q_nv)
  angle = rad2deg(acos(ref_nv.dot(q_nv) / (norm(ref_nv) * norm(q_nv))))
  R1 = axis_angle_to_rotation_matrix(axis, angle)
  rot_ref_q1, rot_ref_q2 = R1.dot(ref_q1), R1.dot(ref_q2)  # rotate ref_q1,2 plane to q1,2 plane

  cos1 = max(min(q1.dot(rot_ref_q1) / (norm(rot_ref_q1) * norm(q1)), 1.), -1.)  # avoid math domain error
  cos2 = max(min(q2.dot(rot_ref_q2) / (norm(rot_ref_q2) * norm(q2)), 1.), -1.)
  angle1 = rad2deg(acos(cos1))
  angle2 = rad2deg(acos(cos2))
  angle = (angle1 + angle2) / 2.
  axis = np.cross(rot_ref_q1, q1)
  R2 = axis_angle_to_rotation_matrix(axis, angle)

  R = R2.dot(R1)
  return R 

def calc_angle(ac1, ac2, ac3, fk, th0):
    th0 = math.pi/180.0 * th0   # degrees to radians
    rji1 = ac1[0] - ac2[0]
    rji2 = ac1[1] - ac2[1]
    rji3 = ac1[2] - ac2[2]
    rjk1 = ac3[0] - ac2[0]
    rjk2 = ac3[1] - ac2[1]
    rjk3 = ac3[2] - ac2[2]
    
    # get the angle from the dot product equation ( A*B = |A|*|B|*cos(theta) )
    # where A and B are vectors bewteen atoms (1->2 and 2->3)
    bji2inv = 1./(rji1**2 + rji2**2 + rji3**2)
    bjk2inv = 1./(rjk1**2 + rjk2**2 + rjk3**2)
    bjiinv = math.sqrt(bji2inv)
    bjkinv = math.sqrt(bjk2inv)
    scp = (rji1*rjk1 + rji2*rjk2 + rji3*rjk3)
    scp = scp * bjiinv* bjkinv
    if scp > 1.0:
        scp =  1.0
    elif scp < -1.0:
        scp = -1.0
    theta = math.acos(scp)   # in radians
    dtheta = (theta-th0)  
    en = 0.5*fk*dtheta**2 
    return en,theta*180.0/math.pi    # kcal/mol and degrees 

def get_spherical_rotatation(p1, p2, width, height, theta_multiplier):
    v1 = get_sphere_mapping(p1[0], p1[1], width, height)
    v2 = get_sphere_mapping(p2[0], p2[1], width, height)

    d = min(max([dot(v1, v2), -1]), 1)

    if abs(d - 1.0) < 0.000001:
        return None

    raxis = norm( cross(v1, v2) )
    rtheta = theta_multiplier * rad2deg * _acos(d)
    #rtheta = 2.0 * rad2deg * _acos(d)

    glPushMatrix()
    glLoadIdentity()
    glRotatef(rtheta, *raxis)
    mat = (c_float*16)()
    glGetFloatv(GL_MODELVIEW_MATRIX, mat)
    glPopMatrix()

    return mat 

def _get_angle_axis(self):
        lim = 1e-12
        norm = np.linalg.norm(self.q)
        if norm < lim:
            angle = 0
            axis = [0, 0, 0]
        else:
            rnorm = 1.0 / norm
            angle = acos(max(-1, min(1, rnorm * self.q[3])));
            sangle = sin(angle)
            if sangle < lim:
                axis = [0, 0, 0]
            else:
                axis = (rnorm / sangle) * np.array(self.q[0:3])

            angle *= 2

        return (angle, axis) 

def slerp(qa, qb, t):
        # Calculate angle between them.
        #qa.w * qb.w + qa.x * qb.x + qa.y * qb.y + qa.z * qb.z;
        cosHalfTheta = np.dot(qa.q,qb.q)
        #if qa=qb or qa=-qb then theta = 0 and we can return qa
        if (abs(cosHalfTheta) >= 1.0):
            return Quat(np.copy(qa.q));
        #Calculate temporary values.
        halfTheta = acos(cosHalfTheta);
        sinHalfTheta = sqrt(1.0 - cosHalfTheta*cosHalfTheta);
        
        #if theta = 180 degrees then result is not fully defined
        #we could rotate around any axis normal to qa or qb
        if(abs(sinHalfTheta) < 0.001):
            return Quat(qa.q * 0.5 + qb.q * 0.5);
        
        ratioA = sin((1 - t) * halfTheta) / sinHalfTheta;
        ratioB = sin(t * halfTheta) / sinHalfTheta; 
        #calculate Quaternion for general case.
        return Quat(qa.q * ratioA + qb.q * ratioB); 

def __init__(self, lattice, m, n):
        self.lattice = lattice
        self.m = m
        self.n = n

        # Chiral vector
        self.c = lattice.pos(m,n)
        self.magC = mag(self.c)
        # Translation vector 
        d = gcd(2*n+m,2*m+n)
        self.t = lattice.pos((2*n+m)/d, -(2*m+n)/d);
        self.magT = mag(self.t)
        # Chiral rotation matrix (rotate a1 along x-axis)
        self.theta = acos(norm(self.c)[0]*copysign(1, self.c[1]))
        self.rotM = np.array([
                [cos(self.theta), sin(self.theta)],
                [-sin(self.theta), cos(self.theta)]]).T

        # Calculate atoms and bonds in unit cell
        self._boundsErr = mag(lattice.pos(0,0,0) - lattice.pos(0,0,1))
        self.indices = self._calcIndices(m, n)
        self.atoms = self._calcAtoms(self.indices)
        self.bonds = self._calcBonds(self.indices) 

def compute_distance_to_wall(a, r1, r2):
    """ Given a sighting that determines a distance of r1 to the wall, and then a rotation by angle *a*
        to the left and a sighting of distance r2, returns the angle that we are currently rotated by
        from the perpendicular to the wall and the current distance to the wall as a pair tuple. Angle
        is given in degrees before the rotation from r1 to r2. Rotating to the right by this angle
        should cause us to face the wall directly.
    """
    try:
        if r1 < r2:
            r = r1/r2
            i = 1.0
        elif r1 > r2:
            r = r2/r1
            i = -1.0
        else:
            return 0, r2
        d = radians(a)
        c = cos(d)
        s = sin(d)
        dt = sqrt(1 - c * c * r * r)
        x = c * c * r + s * dt
        return degrees(acos(x)) * i, r2 * x
    except ValueError:
        traceback.print_exc()
        return None, None 

def calc_eff_area( self, v ) :

		# Note. #nvn is a vector similar to inflow particle bulk 
		# velocity and ndir is the look direction. 


		# Normalize the particle velocity.

		vn  = calc_arr_norm( v )
		nvn = tuple( [ -c for c in vn ] )

		# Calculate the particle inflow angle (in degrees) relative to
		# the cup normal (i.e., the cup pointing direction).

		psi = acos( calc_arr_dot( self['dir'], nvn ) )*pi/180.
		if ( psi > 90. ) :
			return 0. 
 		
		# Return the effective collecting area corresponding to "psi".

		return interp( psi, self._spec._eff_deg, self._spec._eff_area )

	#-----------------------------------------------------------------------
	# DEFINE THE FUNCTION TO CALCULATE EXPECTED MAXWELLIAN CURRENT.
	#----------------------------------------------------------------------- 

def __ComputeCurved(vpercent, w, vec, via, pts, segs):
    """Compute the curves part points"""

    radius = via[1]/2.0

    # Compute the bezier middle points
    req_angle = asin(vpercent/100.0)
    oppside = tan(req_angle)*(radius-(w/sin(req_angle)))
    length = sqrt(radius*radius + oppside*oppside)
    d = req_angle - acos(radius/length)
    vecBC = [vec[0]*cos(d)+vec[1]*sin(d), -vec[0]*sin(d)+vec[1]*cos(d)]
    pointBC = via[0] + wxPoint(int(vecBC[0] * length), int(vecBC[1] * length))
    d = -d
    vecAE = [vec[0]*cos(d)+vec[1]*sin(d), -vec[0]*sin(d)+vec[1]*cos(d)]
    pointAE = via[0] + wxPoint(int(vecAE[0] * length), int(vecAE[1] * length))

    curve1 = __Bezier(pts[1], pointBC, pts[2], n=segs)
    curve2 = __Bezier(pts[4], pointAE, pts[0], n=segs)

    return curve1 + [pts[3]] + curve2 

def getAngle(self, center, p1, p2):
        dx1 = p1.x() - center.x();
        dy1 = p1.y() - center.y();

        dx2 = p2.x() - center.x();
        dy2 = p2.y() - center.y();

        c = math.sqrt(dx1*dx1 + dy1*dy1) * math.sqrt(dx2*dx2 + dy2*dy2)
        if c == 0: return 0
        y = (dx1*dx2+dy1*dy2)/c
        if y>1: return 0
        angle = math.acos(y)

        if (dx1*dy2-dx2*dy1)>0:   
            return angle
        else:
            return -angle 

def earth_distance(lat1, lon1, lat2, lon2):
    """ Distance in meters between two points specified in degrees. """
    x1 = calc_rad(lat1) * math.cos(degree_to_radian(lon1)) * math.sin(degree_to_radian(90-lat1))
    x2 = calc_rad(lat2) * math.cos(degree_to_radian(lon2)) * math.sin(degree_to_radian(90-lat2))
    y1 = calc_rad(lat1) * math.sin(degree_to_radian(lon1)) * math.sin(degree_to_radian(90-lat1))
    y2 = calc_rad(lat2) * math.sin(degree_to_radian(lon2)) * math.sin(degree_to_radian(90-lat2))
    z1 = calc_rad(lat1) * math.cos(degree_to_radian(90-lat1))
    z2 = calc_rad(lat2) * math.cos(degree_to_radian(90-lat2))
    a = (x1*x2 + y1*y2 + z1*z2)/pow(calc_rad((lat1+lat2)/2), 2)
    # a should be in [1, -1] but can sometimes fall outside it by
    # a very small amount due to rounding errors in the preceding
    # calculations (this is prone to happen when the argument points
    # are very close together).  Thus we constrain it here.
    if abs(a) > 1:
        a = 1
    elif a < -1:
        a = -1
    return calc_rad((lat1+lat2) / 2) * math.acos(a) 

def get_zenith(Latitude, Longitude, d, hour, minute, timezone):
    gamma_val = ((2 * math.pi) / 365) * ((d - 1) + (hour - 12) / 24)
    decl_angle = 0.006918 - (0.399912 * math.cos(gamma_val)) + 0.070257 * math.sin(gamma_val) - 0.006758 * math.cos(
        2 * gamma_val) \
                 + 0.000907 * math.sin(2 * gamma_val) - 0.002697 * math.cos(3 * gamma_val) + 0.00148 * math.sin(
        3 * gamma_val)
    eq_time = 229.18 * (
        0.000075 + 0.001868 * math.cos(gamma_val) - 0.032077 * math.sin(gamma_val) - 0.014615 * math.cos(2 * gamma_val)
        - 0.040849 * math.sin(2 * gamma_val))
    time_offset = eq_time - 4 * Longitude + 60*timezone
    true_solar_time = hour * 60 + minute + time_offset
    solar_hour_angle = true_solar_time / 4 - 180

    Z_deg = (180/math.pi)*math.acos((math.sin(Latitude * (math.pi / 180)) * math.sin(decl_angle)) + (
        math.cos(Latitude * (math.pi / 180)) * math.cos(decl_angle) * math.cos(solar_hour_angle * (math.pi / 180))))
    return Z_deg 

def get_zenith(Latitude, Longitude, d, hour, minute, timezone):
    gamma_val = ((2 * math.pi) / 365) * ((d - 1) + (hour - 12) / 24)
    decl_angle = 0.006918 - (0.399912 * math.cos(gamma_val)) + 0.070257 * math.sin(gamma_val) - 0.006758 * math.cos(
        2 * gamma_val) \
                 + 0.000907 * math.sin(2 * gamma_val) - 0.002697 * math.cos(3 * gamma_val) + 0.00148 * math.sin(
        3 * gamma_val)
    eq_time = 229.18 * (
        0.000075 + 0.001868 * math.cos(gamma_val) - 0.032077 * math.sin(gamma_val) - 0.014615 * math.cos(2 * gamma_val)
        - 0.040849 * math.sin(2 * gamma_val))
    time_offset = eq_time - 4 * Longitude + 60*timezone
    true_solar_time = hour * 60 + minute + time_offset
    solar_hour_angle = true_solar_time / 4 - 180

    Z_deg = (180/math.pi)*math.acos((math.sin(Latitude * (math.pi / 180)) * math.sin(decl_angle)) + (
        math.cos(Latitude * (math.pi / 180)) * math.cos(decl_angle) * math.cos(solar_hour_angle * (math.pi / 180))))
    return Z_deg 

def camPosToQuaternion(cx, cy, cz):
    camDist = math.sqrt(cx * cx + cy * cy + cz * cz)
    cx = cx / camDist
    cy = cy / camDist
    cz = cz / camDist
    axis = (-cz, 0, cx)
    angle = math.acos(cy)
    a = math.sqrt(2) / 2
    b = math.sqrt(2) / 2
    w1 = axis[0]
    w2 = axis[1]
    w3 = axis[2]
    c = math.cos(angle / 2)
    d = math.sin(angle / 2)
    q1 = a * c - b * d * w1
    q2 = b * c + a * d * w1
    q3 = a * d * w2 + b * d * w3
    q4 = -b * d * w2 + a * d * w3
    return (q1, q2, q3, q4) 

def SAM(s1, s2):
    """
    Computes the spectral angle mapper between two vectors (in radians).

    Parameters:
        s1: `numpy array`
            The first vector.

        s2: `numpy array`
            The second vector.

    Returns: `float`
            The angle between vectors s1 and s2 in radians.
    """
    try:
        s1_norm = math.sqrt(np.dot(s1, s1))
        s2_norm = math.sqrt(np.dot(s2, s2))
        sum_s1_s2 = np.dot(s1, s2)
        angle = math.acos(sum_s1_s2 / (s1_norm * s2_norm))
    except ValueError:
        # python math don't like when acos is called with
        # a value very near to 1
        return 0.0
    return angle 

def distGEO(x1,y1,x2,y2):
    print("Implementation is wrong")
    assert False
    PI = 3.141592
    deg = int(x1 + .5)
    min_ = x1 - deg
    lat1 = PI * (deg + 5.*min_/3)/180.
    deg = int(y1 + .5)
    min_ = y1 - deg
    long1 = PI * (deg + 5.*min_/3)/180.
    deg = int(x2 + .5)
    min_ = x2 - deg
    lat2 = PI * (deg + 5.*min_/3)/180.
    deg = int(y2 + .5)
    min_ = y2 - deg
    long2 = PI * (deg + 5.*min_/3)/180.


    RRR = 6378.388
    q1 = math.cos( long1 - long2 );
    q2 = math.cos( lat1 - lat2 );
    q3 = math.cos( lat1 + lat2 );
    return int(RRR * math.acos(.5*((1.+q1)*q2 - (1.-q1)*q3)) + 1.) 

def distGEO(x1,y1,x2,y2):
    print("Implementation is wrong")
    assert False
    PI = 3.141592
    deg = int(x1 + .5)
    min_ = x1 - deg
    lat1 = PI * (deg + 5.*min_/3)/180.
    deg = int(y1 + .5)
    min_ = y1 - deg
    long1 = PI * (deg + 5.*min_/3)/180.
    deg = int(x2 + .5)
    min_ = x2 - deg
    lat2 = PI * (deg + 5.*min_/3)/180.
    deg = int(y2 + .5)
    min_ = y2 - deg
    long2 = PI * (deg + 5.*min_/3)/180.


    RRR = 6378.388
    q1 = math.cos( long1 - long2 );
    q2 = math.cos( lat1 - lat2 );
    q3 = math.cos( lat1 + lat2 );
    return int(RRR * math.acos(.5*((1.+q1)*q2 - (1.-q1)*q3)) + 1.) 

def angleOfVector(p1, v, p2):
    """
    Accepts p1,  v,  and p2 of class Point or coordinate pairs
    Returns the angle in radians between the vector v-to-p1 and vector v-to-p2
    """
    #L1=distance(p1, p2)
    #L2=distance(p1, p3)
    #L3=distance(p2, p3)
    #return math.acos((L1**2+L2**2-L3**2)/(2*L1*L2))
    p1 = dPnt(p1)
    p2 = dPnt(p2)
    angle1 = angleOfLine(v, p1)
    angle2 = angleOfLine(v, p2)
    #get the absolute angle
    angle = abs(angle1 - angle2)
    #get the smallest angle of the vector, should not be greater than a straight line
    if angle > pi:
        angle = 2*pi - angle
    return angle 

def angleOfVector(p1, v, p2):
    """
    Accepts p1,  v,  and p2 of class Point or coordinate pairs
    Returns the angle in radians between the vector v-to-p1 and vector v-to-p2
    """
    #L1=distance(p1, p2)
    #L2=distance(p1, p3)
    #L3=distance(p2, p3)
    #return math.acos((L1**2+L2**2-L3**2)/(2*L1*L2))
    p1 = dPnt(p1)
    p2 = dPnt(p2)
    angle1 = angleOfLine(v, p1)
    angle2 = angleOfLine(v, p2)
    #get the absolute angle
    angle = abs(angle1 - angle2)
    #get the smallest angle of the vector, should not be greater than a straight line
    if angle > pi:
        angle = 2*pi - angle
    return angle 

def dist_sqr_at_l(self, p):
        """return squared distance and arc length at nearest point to p on biarc"""
        p1, t1, cha1 = biarc.point_on_circle(self.p0, self.J, self.t0, p)
        p2, t2, cha2 = biarc.point_on_circle(self.J, self.p1, self.Jt, p)
        c1, k1, a1, l1, c2, k2, a2, l2 = self.circleparameters()
        mind, minl = None, None
        if 0 < t1 < 1:
            mind = la.norm2(p1 - p)
            aa1 = 2 * m.acos(cha1)
            minl = aa1 / abs(k1) if k1 != 0 else t1 * l1
        if 0 < t2 < 1 and (not mind or la.norm2(p2 - p) < mind):
            mind = la.norm2(p2 - p)
            aa2 = 2 * m.acos(cha2)
            minl = l1 + (aa2 / abs(k2) if k2 != 0 else t2 * l2)
        return mind, minl
########################################################################### 

def proj_xy(self, t, next=None):
        """
            length of projection of sections at crossing line / circle intersections
            deformation unit vector for profil in xy axis
            so f(x_profile) = position of point in xy plane
        """
        if next is None:
            return self.normal(t).v.normalized(), 1
        v0 = self.normal(1).v.normalized()
        v1 = next.normal(0).v.normalized()
        direction = v0 + v1
        adj = (v0 * self.length) * (v1 * next.length)
        hyp = (self.length * next.length)
        c = min(1, max(-1, adj / hyp))
        size = 1 / cos(0.5 * acos(c))
        return direction.normalized(), min(3, size) 

def _get_defects_count(array, contour, defects, verbose = False):
    ndefects = 0
    
    for i in range(defects.shape[0]):
        s,e,f,_ = defects[i,0]
        beg     = tuple(contour[s][0])
        end     = tuple(contour[e][0])
        far     = tuple(contour[f][0])
        a       = _get_eucledian_distance(beg, end)
        b       = _get_eucledian_distance(beg, far)
        c       = _get_eucledian_distance(end, far)
        angle   = math.acos((b ** 2 + c ** 2 - a ** 2) / (2 * b * c)) * 57
            
        if angle <= 90:
            ndefects = ndefects + 1

            if verbose:
                cv2.circle(array, far, 3, _COLOR_RED, -1)

        if verbose:
            cv2.line(array, beg, end, _COLOR_RED, 1)

    return array, ndefects 

def orien(current,initial,others):
    temp=[]
    #print others
    f=0
    for i in others:
        if  i != current:
            #print i,current,initial
            a=length(current,initial)
            b=length(i,initial)
            c=length(current,i)
            d=math.acos((a*a + b*b - c*c) / (2*a*b))
            temp.append((d,f))
        else:
            temp.append((0,f))
        f+=1
    return temp 

def orien(current,initial,others):
    temp=[]
    #print others
    f=0
    for i in others:
        if  i != current:
            #print i,current,initial
            a=length(current,initial)
            b=length(i,initial)
            c=length(current,i)
            d=math.acos((a*a + b*b - c*c) / (2*a*b))
            temp.append((d,f))
        else:
            temp.append((0,f))
        f+=1
    return temp 

def to_angle_axis(self):
        """ Return axis-angle representation """
        q = np.roll(self.q, shift=1)
        halftheta = math.acos(q[0])
        if abs(halftheta) < 1e-12:
            return 0, np.array((0, 0, 1))
        else:
            theta = halftheta * 2
            axis = np.array(q[1:4]) / math.sin(halftheta)
            return theta, axis 

def quaternion_slerp(quat0, quat1, fraction, spin=0, shortestpath=True):
    """Return spherical linear interpolation between two quaternions.

    >>> q0 = random_quaternion()
    >>> q1 = random_quaternion()
    >>> q = quaternion_slerp(q0, q1, 0.0)
    >>> numpy.allclose(q, q0)
    True
    >>> q = quaternion_slerp(q0, q1, 1.0, 1)
    >>> numpy.allclose(q, q1)
    True
    >>> q = quaternion_slerp(q0, q1, 0.5)
    >>> angle = math.acos(numpy.dot(q0, q))
    >>> numpy.allclose(2.0, math.acos(numpy.dot(q0, q1)) / angle) or \
        numpy.allclose(2.0, math.acos(-numpy.dot(q0, q1)) / angle)
    True

    """
    q0 = unit_vector(quat0[:4])
    q1 = unit_vector(quat1[:4])
    if fraction == 0.0:
        return q0
    elif fraction == 1.0:
        return q1
    d = numpy.dot(q0, q1)
    if abs(abs(d) - 1.0) < _EPS:
        return q0
    if shortestpath and d < 0.0:
        # invert rotation
        d = -d
        q1 *= -1.0
    angle = math.acos(d) + spin * math.pi
    if abs(angle) < _EPS:
        return q0
    isin = 1.0 / math.sin(angle)
    q0 *= math.sin((1.0 - fraction) * angle) * isin
    q1 *= math.sin(fraction * angle) * isin
    q0 += q1
    return q0 

def quaternion_slerp(quat0, quat1, fraction, spin=0, shortestpath=True):
    """Return spherical linear interpolation between two quaternions.

    >>> q0 = random_quaternion()
    >>> q1 = random_quaternion()
    >>> q = quaternion_slerp(q0, q1, 0)
    >>> numpy.allclose(q, q0)
    True
    >>> q = quaternion_slerp(q0, q1, 1, 1)
    >>> numpy.allclose(q, q1)
    True
    >>> q = quaternion_slerp(q0, q1, 0.5)
    >>> angle = math.acos(numpy.dot(q0, q))
    >>> numpy.allclose(2, math.acos(numpy.dot(q0, q1)) / angle) or \
        numpy.allclose(2, math.acos(-numpy.dot(q0, q1)) / angle)
    True

    """
    q0 = unit_vector(quat0[:4])
    q1 = unit_vector(quat1[:4])
    if fraction == 0.0:
        return q0
    elif fraction == 1.0:
        return q1
    d = numpy.dot(q0, q1)
    if abs(abs(d) - 1.0) < _EPS:
        return q0
    if shortestpath and d < 0.0:
        # invert rotation
        d = -d
        numpy.negative(q1, q1)
    angle = math.acos(d) + spin * math.pi
    if abs(angle) < _EPS:
        return q0
    isin = 1.0 / math.sin(angle)
    q0 *= math.sin((1.0 - fraction) * angle) * isin
    q1 *= math.sin(fraction * angle) * isin
    q0 += q1
    return q0 

def __init__(self, arg):
        if isinstance(arg, numbers.Real):
            # We precompute sin and cos for rotations
            self.angle = arg
            self.cos = math.cos(self.angle)
            self.sin = math.sin(self.angle)
        elif isinstance(arg, Point):
            # Point angle is the trigonometric angle of the vector
            # [origin, Point]
            pt = arg
            try:
                self.cos = pt.x / pt.length()
                self.sin = pt.y / pt.length()
            except ZeroDivisionError:
                self.cos = 1
                self.sin = 0

            self.angle = math.acos(self.cos)
            if self.sin < 0:
                self.angle = -self.angle
        else:
            raise TypeError("Angle is defined by a number or a Point") 

def __calc(self):
  """
  Perform the actual calculations for sunrise, sunset and
  a number of related quantities.
  
  The results are stored in the instance variables
  sunrise_t, sunset_t and solarnoon_t
  """
  timezone = self.timezone # in hours, east is positive
  longitude= self.long     # in decimal degrees, east is positive
  latitude = self.lat      # in decimal degrees, north is positive

  time  = self.time # percentage past midnight, i.e. noon  is 0.5
  day      = self.day     # daynumber 1=1/1/1900
 
  Jday     =day+2415018.5+time-timezone/24 # Julian day
  Jcent    =(Jday-2451545)/36525    # Julian century

  Manom    = 357.52911+Jcent*(35999.05029-0.0001537*Jcent)
  Mlong    = 280.46646+Jcent*(36000.76983+Jcent*0.0003032)%360
  Eccent   = 0.016708634-Jcent*(0.000042037+0.0001537*Jcent)
  Mobliq   = 23+(26+((21.448-Jcent*(46.815+Jcent*(0.00059-Jcent*0.001813))))/60)/60
  obliq    = Mobliq+0.00256*cos(rad(125.04-1934.136*Jcent))
  vary     = tan(rad(obliq/2))*tan(rad(obliq/2))
  Seqcent  = sin(rad(Manom))*(1.914602-Jcent*(0.004817+0.000014*Jcent))+sin(rad(2*Manom))*(0.019993-0.000101*Jcent)+sin(rad(3*Manom))*0.000289
  Struelong= Mlong+Seqcent
  Sapplong = Struelong-0.00569-0.00478*sin(rad(125.04-1934.136*Jcent))
  declination = deg(asin(sin(rad(obliq))*sin(rad(Sapplong))))
  
  eqtime   = 4*deg(vary*sin(2*rad(Mlong))-2*Eccent*sin(rad(Manom))+4*Eccent*vary*sin(rad(Manom))*cos(2*rad(Mlong))-0.5*vary*vary*sin(4*rad(Mlong))-1.25*Eccent*Eccent*sin(2*rad(Manom)))

  hourangle= deg(acos(cos(rad(90.833))/(cos(rad(latitude))*cos(rad(declination)))-tan(rad(latitude))*tan(rad(declination))))

  self.solarnoon_t=(720-4*longitude-eqtime+timezone*60)/1440
  self.sunrise_t  =self.solarnoon_t-hourangle*4/1440
  self.sunset_t   =self.solarnoon_t+hourangle*4/1440 

def get_angle(v1,v2): # v1,v2 must be unit vectors
    denom = (math.sqrt(((v1[0]*v1[0]) + (v1[1]*v1[1]) + (v1[2]*v1[2]))) *
                math.sqrt(((v2[0]*v2[0]) + (v2[1]*v2[1]) + (v2[2]*v2[2]))))
    if denom>1e-10:
        result = ( (v1[0]*v2[0]) + (v1[1]*v2[1]) + (v1[2]*v2[2]) ) / denom
    else:
        result = 0.0
    result = math.acos(result)
    return result

#------------------------------------------------------------------------------ 

def quaternion_slerp(quat0, quat1, fraction, spin=0, shortestpath=True):
    """Return spherical linear interpolation between two quaternions.

    >>> q0 = random_quaternion()
    >>> q1 = random_quaternion()
    >>> q = quaternion_slerp(q0, q1, 0.0)
    >>> numpy.allclose(q, q0)
    True
    >>> q = quaternion_slerp(q0, q1, 1.0, 1)
    >>> numpy.allclose(q, q1)
    True
    >>> q = quaternion_slerp(q0, q1, 0.5)
    >>> angle = math.acos(numpy.dot(q0, q))
    >>> numpy.allclose(2.0, math.acos(numpy.dot(q0, q1)) / angle) or \
        numpy.allclose(2.0, math.acos(-numpy.dot(q0, q1)) / angle)
    True

    """
    q0 = unit_vector(quat0[:4])
    q1 = unit_vector(quat1[:4])
    if fraction == 0.0:
        return q0
    elif fraction == 1.0:
        return q1
    d = numpy.dot(q0, q1)
    if abs(abs(d) - 1.0) < _EPS:
        return q0
    if shortestpath and d < 0.0:
        # invert rotation
        d = -d
        q1 *= -1.0
    angle = math.acos(d) + spin * math.pi
    if abs(angle) < _EPS:
        return q0
    isin = 1.0 / math.sin(angle)
    q0 *= math.sin((1.0 - fraction) * angle) * isin
    q1 *= math.sin(fraction * angle) * isin
    q0 += q1
    return q0