Python numpy.arcsin函数代码示例

def theta_generate(system, plot = False, realonly = False, useapproximations = False):
        global theta_2,theta_3, cos_theta_2, cos_theta_3, tan_theta_2
        (n_1,n_2,n_3) = system.indices
        
        theta_1 = system.theta_1
        freq = system.freq


        cos_theta_2 = (n_2)**-1  * array(n_2**2-((n_1)*sin(theta_1))**2,dtype = complex)
        cos_theta_3 = (n_3)**-1  * array(n_3**2-((n_1)*sin(theta_1))**2,dtype = complex)
        tan_theta_2 = sqrt((1-cos_theta_2**2)/cos_theta_2)

            
        theta_2 = numpy.arcsin(sin(theta_1) * n_1/ n_2,dtype = complex)
        theta_3 = numpy.arcsin(sin(theta_1)* n_1 / n_3,dtype = complex)

        
        
        if plot == True:

                
                py.plot(freq,numpy.real(cos_theta_2[IR]),'r')
                py.plot(freq,numpy.imag(cos_theta_2[IR]),'b')
                py.plot(freq,numpy.real(cos_theta_3[IR]),'g')
                py.plot(freq,numpy.imag(cos_theta_3[IR]),'k')
                py.title("cos Angles of incidence")
                py.legend(('real2','imag2','real3','imag3'))
                py.show()
        
        system.angles = (theta_2,theta_3, cos_theta_2, cos_theta_3, tan_theta_2)
        
        return 0

 def __arc_segment(self,r,a,b):
   """
     Returns length of circle segment in a half-infinite slit.
       r   ... radius of circle which is centered at the origin
       a,b ... borders of slit containing (x,y) with x>a, b>y>0 
   """  
   #print r, a, b
   # scalar version
   if np.isscalar(r):
     if r<a:              # circle outside slit
       return 0;   
     elif r**2>a**2+b**2: # circle does not intersect left border at x=a
       return r*np.arcsin(b/r);
     else:                # circle only intersects left border at x=a
       return r*np.arccos(a/r);
 
   # parallel version
   else:
     r     = np.atleast_1d(r); 
     arc   = np.zeros(r.size);
     index = r>=a;              # select radii which intersect half-slit
     y2_a  = r[index]**2-a**2;  # intersection of circle with left border at x=a
     ymax  = np.where( y2_a > b**2, b, np.sqrt(y2_a) ); 
                                # largest y value on circle segment in half-slit
     arc[index] = r[index] * np.arcsin( ymax / r[index] );
                                # calculate arc length
     # DEBUG
     if self.verbosity>2:
       arc2 = [ self.__arc_segment(ir,a,b) for ir in r ];
       assert(np.allclose( arc, arc2 ));
     return arc;
 
   raise RuntimeError("Should never reach this point");

 def getCurvePoints(self,P0,P1,P2,R,n=20):
     """Generates a curve. P0,P1 and P2 define a corner and R defines the radius of the
     tangential circle element. n gives the number of point to approximate the curve.
     P1 is the corner itself and P0 and P2 give the tangents.
     """
     P0=np.array(P0)
     P1=np.array(P1)
     P2=np.array(P2)
     o1=(P0-P1)/np.sqrt(np.dot(P0-P1,P0-P1))
     o2=(P2-P1)/np.sqrt(np.dot(P2-P1,P2-P1))
     if np.arcsin(np.cross(o1,o2)) > 0:
         a=1.
         b=-1.
     else:
         a=-1.
         b=1.
     
     v1=R*np.dot(np.array([[0.,b],[a,0.]]),o1)
     v2=R*np.dot(np.array([[0.,a],[b,0.]]),o2)
     dv=v2-v1
     a=np.array([[o1[0],-o2[0]],[o1[1],-o2[1]]])
     b=dv
     x=np.linalg.solve(a,b)
     circleCenter= P1+x[0]*o1+v1
     angle = np.arcsin(np.cross(v2/R,v1/R))
     points=[]
     for i in range(n+1):
         x=-i*angle/n
         rot = np.array([[np.cos(x),-np.sin(x)],
                          [np.sin(x),np.cos(x)]])
         points.append(circleCenter+np.dot(rot,-v1))
     return points

    def __remove_particles_randomly(self):

        grid_pairs = np.array([self.grid[0].ravel(), self.grid[1].ravel()]).T

        particles = grid_pairs[
            grid_pairs[:, 1]
            > (
                self.sin_amp * np.sin(self.sin_freq * grid_pairs[:, 0] + np.arcsin(1.0))
                + self.sin_amp
                + self.particle_diameter
            )
        ]
        self.driver = grid_pairs[
            grid_pairs[:, 1]
            <= (
                self.sin_amp * np.sin(self.sin_freq * grid_pairs[:, 0] + np.arcsin(1.0))
                + self.sin_amp
                + self.particle_diameter
            )
        ]

        self.x = particles[:, 0]

        number_of_particles_to_remove = self.__compute_number_of_particles_to_remove()

        np.random.shuffle(particles)

        self.x = particles[:-number_of_particles_to_remove, 0]
        self.y = particles[:-number_of_particles_to_remove, 1]

def convert_galactic_equa(l, b):
   '''l and b in degrees'''
   r = pi / 180.0
   dec = arcsin(cos(b * r) * sin((l - 33) * r) * sin(62.6 * r) + sin(b * r) * cos(62.6 * r))
   ra  = 282.25 + arcsin((cos(b * r) * sin((l - 33) * r) * cos(62.6 * r) - sin(b * r) * sin(62.6 * r)) / cos(dec)) / r
   dec = dec / r
   return ra, dec

  def test_cartesianToSpherical(self):
    """
    Check correct working of the Cartesian to Spherical coordinates transformation.
    """
    r, phi, theta = cartesianToSpherical(array([1,0,0]), array([0,1,0]), array([0,0,1]))
    assert_array_almost_equal(r,array([1.0,1.0,1.0]))
    assert_array_almost_equal(phi,array([0.0,pi/2.0,0.0]))
    assert_array_almost_equal(theta,array([0.0,0.0,pi/2.0]))

    r, phi, theta = cartesianToSpherical(1.0, 1.0, 0.0)
    assert_almost_equal(r, sqrt(2.0))
    assert_almost_equal(phi, pi/4.0)
    assert_almost_equal(theta, 0.0)

    r, phi, theta = cartesianToSpherical(1.0, 1.0, 1.0)
    assert_almost_equal(r, sqrt(3.0))
    assert_almost_equal(phi, pi/4.0)
    assert_almost_equal(theta, arcsin(1.0/sqrt(3.0)))

    r, phi, theta = cartesianToSpherical(1.0, 1.0, -1.0)
    assert_almost_equal(r, sqrt(3.0))
    assert_almost_equal(phi, pi/4.0)
    assert_almost_equal(theta, -arcsin(1.0/sqrt(3.0)))

    r, phi, theta = cartesianToSpherical(-1.0, -1.0, -1.0)
    assert_almost_equal(r, sqrt(3.0))
    assert_almost_equal(phi, pi/4.0+pi)
    assert_almost_equal(theta, -arcsin(1.0/sqrt(3.0)))

    assert_raises(Exception, cartesianToSpherical, 0.0, 0.0, 0.0)
    assert_raises(Exception, cartesianToSpherical, array([1,0,0,0]), array([0,1,0,0]), array([0,0,0,1]))

def solve_acos_plus_bsin_plus_c(a, b, c):
    """
    Solves a*cos(x) + b*sin(x) + c = 0 for x in [-pi/2, pi/2].

    The method used is to change variable by setting t = sin(x).

    Args:
        a, b, c: coefficients of the equation

    Returns:
        A solution x if it exists, None in the other cases.
    """
    out = None
    if a * a + b * b - c * c > 0:
        poly = np.poly1d([a * a + b * b, 2 * b * c, c * c - a * a])
        if not any([-1 <= t <= 1 for t in poly.r]):
            print("no solutions between -1 and 1: ", poly.r)
        elif all([-1 <= t <= 1 for t in poly.r]):
            i = np.argmin(np.abs(np.array([a * np.sqrt(1 - t * t) + b * t + c for t in poly.r])))
            out = np.arcsin(poly.r[i])
        else:
            if -1 <= poly.r[0] <= 1:
                out = np.arcsin(poly.r[0])
            else:
                np.arcsin(poly.r[1])

    # check that the solution is in [-pi/2, pi/2] and satisfies the equation
    assert -np.pi / 2 < out < np.pi / 2
    np.testing.assert_allclose(a * np.cos(out) + b * np.sin(out) + c, 0, atol=1e-12)
    return out

def tthToRad(twoTheta, unit, wavelength=None, directDist=None):
    """
    Convert a two theta angle from original `unit` to radian.

    `directDist = ai.getFit2D()["directDist"]`
    """
    if isinstance(twoTheta, numpy.ndarray):
        pass
    elif isinstance(twoTheta, collections.Iterable):
        twoTheta = numpy.array(twoTheta)

    if unit == units.TTH_RAD:
        return twoTheta
    elif unit == units.TTH_DEG:
        return numpy.deg2rad(twoTheta)
    elif unit == units.Q_A:
        if wavelength is None:
            raise AttributeError("wavelength have to be specified")
        return numpy.arcsin((twoTheta * wavelength) / (4.e-10 * numpy.pi)) * 2.0
    elif unit == units.Q_NM:
        if wavelength is None:
            raise AttributeError("wavelength have to be specified")
        return numpy.arcsin((twoTheta * wavelength) / (4.e-9 * numpy.pi)) * 2.0
    elif unit == units.R_MM:
        if directDist is None:
            raise AttributeError("directDist have to be specified")
        # GF: correct formula?
        return numpy.arctan(twoTheta / directDist)
    elif unit == units.R_M:
        if directDist is None:
            raise AttributeError("directDist have to be specified")
        # GF: correct formula?
        return numpy.arctan(twoTheta / (directDist * 0.001))
    else:
        raise ValueError("Converting from 2th to unit %s is not supported", unit)

    def reconstruct_angle(self, event, offsets=None):
        """Reconstruct angles from a single event"""

        c = 3.00e+8

        if offsets is not None:
            self._correct_offsets(event, offsets)

        dt1 = event['t1'] - event['t3']
        dt2 = event['t1'] - event['t4']

        station = self.cluster.stations[event['station_id']]
        r1, phi1 = station.calc_r_and_phi_for_detectors(1, 3)
        r2, phi2 = station.calc_r_and_phi_for_detectors(1, 4)

        phi = arctan2((dt2 * r1 * cos(phi1) - dt1 * r2 * cos(phi2)),
                      (dt2 * r1 * sin(phi1) - dt1 * r2 * sin(phi2)) * -1)
        theta1 = arcsin(c * dt1 * 1e-9 / (r1 * cos(phi - phi1)))
        theta2 = arcsin(c * dt2 * 1e-9 / (r2 * cos(phi - phi2)))

        e1 = sqrt(self.rel_theta1_errorsq(theta1, phi, phi1, phi2, r1, r2))
        e2 = sqrt(self.rel_theta2_errorsq(theta2, phi, phi1, phi2, r1, r2))

        theta_wgt = (1 / e1 * theta1 + 1 / e2 * theta2) / (1 / e1 + 1 / e2)

        if theta_wgt < 0:
            theta_wgt *= -1
            phi += pi
            phi = (phi + pi) % (2 * pi) - pi

        return theta_wgt, phi

def tpi2k(id, p, dI):
    p2brho = 10./2.99792458
    deg = pi/180.0
    if id not in table.keys():
        print( "type identifier is unknown")
        print( "known types:")
        print( list(table.keys())[:10])
        print( list(table.keys())[10:20])
        print( list(table.keys())[20:30])
        print( list(table.keys())[35:])
        return None
    magnet_type = table[id]["type"]
    a = table[id]["A"]
    el = table[id]["EL"]
    km = table[id]["KM"]
    dk0 = p3(a, dI ) * km/(p*p2brho)
    dk = 0.
    if magnet_type == "S" or magnet_type == "Q":
        #case kQ:
        dk = dk0
    elif magnet_type == "B":
        phi = 2.0*np.arcsin(0.5*el * dk0)
        dk = phi/deg
    elif magnet_type ==  "C":
        phi = 2.0*np.arcsin(0.5*el * dk0)
        dk = phi*1.e3
    elif magnet_type == "T":
        phi = 2.0*np.arcsin(0.5*el * dk0)
        dk = phi*1.e3
    return dk

def test_rphi_to_dl_2d():
    #This is a tangent point
    r,phi= 6., numpy.arccos(0.75)
    d,l= bovy_coords.rphi_to_dl_2d(r,phi,degree=False,ro=8.,phio=0.)
    l= numpy.arcsin(0.75)
    d= 6./numpy.tan(l)
    assert numpy.fabs(d-6./numpy.tan(numpy.arcsin(0.75))) < 10.**-10., 'dl_to_rphi_2d conversion did not work as expected'
    assert numpy.fabs(l-numpy.arcsin(0.75)) < 10.**-10., 'rphi_to_dl_2d conversion did not work as expected'
    #This is another point
    r,phi= 2., 55.
    d,l= bovy_coords.rphi_to_dl_2d(r,phi,degree=True,ro=2.*numpy.sqrt(2.),
                                   phio=10.)
    assert numpy.fabs(d-2.) < 10.**-10., 'rphi_to_dl_2d conversion did not work as expected'
    assert numpy.fabs(l-45.) < 10.**-10., 'rphi_to_dl_2d conversion did not work as expected'
    #This is another point, for arrays
    r,phi= 2., 45.
    os= numpy.ones(2)
    d,l= bovy_coords.rphi_to_dl_2d(os*r,os*phi,
                                   degree=True,ro=2.*numpy.sqrt(2.),
                                   phio=0.)
    assert numpy.all(numpy.fabs(d-2.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    assert numpy.all(numpy.fabs(l-45.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    #This is another point, for lists, which for some reason I support
    r,phi= 2., 45.
    d,l= bovy_coords.rphi_to_dl_2d([r,r],[phi,phi],
                                   degree=True,ro=2.*numpy.sqrt(2.),
                                   phio=0.)
    d= numpy.array(d)
    l= numpy.array(l)
    assert numpy.all(numpy.fabs(d-2.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    assert numpy.all(numpy.fabs(l-45.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    return None

 def _cal_center(p0,p1,p2):
     #p0 is of type up-down 
     origin=(p1+p2)/2
     y_v=f3(np.zeros(3),p1-origin)
     x_v=f3(np.zeros(3),p0-origin)
     z_v=np.cross(x_v,y_v)
     T=f1(x0_v,y0_v,z0_v,x_v,y_v,z_v)
     #print T
     r=f2(p1,p2)/2*np.tan(np.pi/2-self.open_angle/2.)
     phi=0.
     L1=f2(p1,p2)/2./np.tan(self.open_angle/2.)
     L2=f2(p0,origin)
     #look at document#1 in binder for detail
     a,b,c=1+np.tan(self.theta_top_down)**2,2*L1/L2*np.tan(self.theta_top_down),(L1/L2)**2-1
     sin_list=[(b+(b**2-4*a*c)**0.5)/2./a,(b-(b**2-4*a*c)**0.5)/2./a]
     theta=0
     if (sin_list[0]<1)&(sin_list[0]>0):
         if self.switch==False:
             theta=np.pi/2+np.arcsin(sin_list[0])
         elif self.switch==True:
             theta=np.pi/2-np.arcsin(sin_list[0])
     else:
         if self.switch==False:
             theta=np.pi/2+np.arcsin(sin_list[1])
         elif self.switch==True:
             theta=np.pi/2-np.arcsin(sin_list[1])
         
     center_point_new=np.array([r*np.cos(phi)*np.sin(theta),r*np.sin(phi)*np.sin(theta),r*np.cos(theta)])
     center_point_org=np.dot(inv(T),center_point_new)+origin
     return center_point_org

def R2axis(M):

    m00 = M[0, 0]
    m01 = M[0, 1]
    m02 = M[0, 2]
    m10 = M[1, 0]
    m11 = M[1, 1]
    m12 = M[1, 2]
    m20 = M[2, 0]
    m21 = M[2, 1]
    m22 = M[2, 2]
    # symmetric matrix K
    K = np.array(
        [
            [m00 - m11 - m22, 0.0, 0.0, 0.0],
            [m01 + m10, m11 - m00 - m22, 0.0, 0.0],
            [m02 + m20, m12 + m21, m22 - m00 - m11, 0.0],
            [m21 - m12, m02 - m20, m10 - m01, m00 + m11 + m22],
        ]
    )
    K /= 3.0
    # quaternion is eigenvector of K that corresponds to largest eigenvalue
    w, V = np.linalg.eigh(K)
    q = V[[3, 0, 1, 2], np.argmax(w)]

    # quaternion to Axis
    nrm = np.linalg.norm(q[1:4])
    if nrm <= 0.0:
        a = np.array([0.0, 0.0, 0.0])
    elif q[0] < 0.0:
        a = (np.pi - np.arcsin(nrm)) * 2 * q[1:4] / nrm
    else:
        a = np.arcsin(nrm) * 2 * q[1:4] / nrm

    return a

def funDR1R2(x, v1, v2, v3, z1, z2):
# def funDR1R2(x, v1, v2, v3, z1, z2):
	"""
		Computes arrival time of direct, and two critically refracted waves from three layer models

		x: offset (array)
		v1: velocity of 1st layer (float)
		v2: velocity of 2nd layer (float)
		v3: velocity of 3rd layer (float)
		z1: thickness of 1st layer (float)
		z2: thickness of 2nd layer (float)

	"""
	direct = 1./v1*x
	theta1 = np.arcsin(v1/v2)    
	theta2 = np.arcsin(v2/v3)
	ti1 = 2*z1*np.cos(theta1)/v1
	ti2 = 2*z2*np.cos(theta2)/v2
	xc1 = 2*z1*np.tan(theta1)
	xc2 = 2*z2*np.tan(theta2)
	act1 = x > xc1
	act2 = x > xc2   
	ref1 = np.ones_like(x)*np.nan
	ref1[act1] = 1./v2*x[act1] + ti1    
	ref2 = np.ones_like(x)*np.nan
	ref2[act2] = 1./v3*x[act2] + ti2 + ti1
	t0 = 2.*z1/v1
	refl1 = np.sqrt(t0**2+x**2/v1**2)

	return np.c_[direct, ref1, ref2, refl1]

def evalExactφ(sqp,sqpnext,sqpdesired):

	a=sqp/sqpdesired
	b=sqpnext/sqpdesired
	if(abs(a+1) > float("1e-12") and (a*a + b*b - 1)>0):

		φ1=2*np.arctan(     ( b -  (a*a + b*b - 1)**0.5 )/(a+1)    )
		φ2=2*np.arctan(     ( b +  (a*a + b*b - 1)**0.5 )/(a+1)    )

		if(np.isnan(φ1)):
			print("Ok nan found, locating butter chicken")
			print( (a*a + b*b -1) )
			print( (b + (a*a + b*b -1)**0.5 ) )

		φ1=np.arcsin(np.sin(φ1))
		φ2=np.arcsin(np.sin(φ2))
		if(abs(φ1)<abs(φ2)):
			print(φ1,φ2)
			return φ1
		else:
			print(φ2,φ1)
			return φ2
	else:
		print("Glaba, something went globular :P")
		return 0