notebook.community

Edit and run 





In [382]:



    


from numpy import roots, random, array, linspace
from numpy import fabs as fb
from math import atan, fabs, pi, tan, cos, sin, sqrt, tan
import sympy
from sympy import Symbol

def cubic (Qx, Qy, Ux, Uy):
    
    a = Uy
    b = (Ux + 2*Qy)
    c = (2*Qx - Uy)
    d = -Ux
    
    det = -4*b*b*b*d + b*b*c*c - 4*a*c*c*c + 18*a*b*c*d - 27*a*a*d*d
    
    if (det < 0):
        return 'c'
    
    if (det > 0):
        
        a = roots([a, b, c, d])
        a = a.real
        
        a = [atan(a[0]), atan(a[1]), atan(a[2])]
            
        U = [Ux*cos(a[0]) + Uy*sin(a[0]), Ux*cos(a[1]) + Uy*sin(a[1]), Ux*cos(a[2]) + Uy*sin(a[2])]
        rightU = [2*sin(a[0])*cos(a[0]), 2*sin(a[1])*cos(a[1]), 2*sin(a[2])*cos(a[2])]
    
        for i in range(0, 3):
            if (U[i] * rightU[i] < 0):
                a[i] = a[i] + pi
        
        print a[0]*(180/pi), a[1]*(180/pi), a[2]*(180/pi)

        a = sorted(a)
        a = [a[0] - a[0], a[1] - a[0], a[2] - a[0]]        
        
        if (a[2] > pi):
            return 'a'
        else:
            return 'b'



    











In [346]:



    


mu, sigma = 0, 1



    











In [524]:



    


Ux = 1.0
Uy = 1.0
Qx = -1.0
Qy = 1.1



    











In [525]:



    


print cubic(Qx, Qy, Ux, Uy)



    


















    






-75.6291843111 44.0609744917 165.171028792
a

















In [526]:



    


U = array([-2*Ux + 2*Uy, -Ux + 2*Uy, 2*Uy, Ux + 2*Uy, 2*Ux + 2*Uy,-2*Ux + Uy, -Ux + Uy, Uy, Ux + Uy, 2*Ux + Uy,-2*Ux, -Ux, 0.0, Ux, 2*Ux,-2*Ux - Uy, -Ux - Uy, -Uy, Ux - Uy, 2*Ux - Uy,-2*Ux - 2*Uy, -Ux - 2*Uy, -2*Uy, Ux - 2*Uy, 2*Ux - 2*Uy])
Q = array([-2*Qx + 2*Qy, -Qx + 2*Qy, 2*Qy, Qx + 2*Qy, 2*Qx + 2*Qy,-2*Qx + Qy, -Qx + Qy, Qy, Qx + Qy, 2*Qx + Qy,-2*Qx, -Qx, 0.0, Qx, 2*Qx,-2*Qx - Qy, -Qx - Qy, -Qy, Qx - Qy, 2*Qx - Qy,-2*Qx - 2*Qy, -Qx - 2*Qy, -2*Qy, Qx - 2*Qy, 2*Qx - 2*Qy])
P_x = Ux * U + Q * Qx
P_y = Uy * U + Q * Qy



    











In [527]:



    


P1_x = array([-2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2])
P1_y = array([2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -2, -2, -2, -2, -2])

P2_x = P1_x + P_x*0.1*np.reciprocal(fb(P_x))
P2_y = P1_y + P_y*0.1*np.reciprocal(fb(P_y))

P3_x = P1_x - P_x*0.1*np.reciprocal(fb(P_x))
P3_y = P1_y - P_y*0.1*np.reciprocal(fb(P_y))



    


















    






/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:4: RuntimeWarning: divide by zero encountered in reciprocal
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:4: RuntimeWarning: invalid value encountered in multiply
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:5: RuntimeWarning: divide by zero encountered in reciprocal
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:5: RuntimeWarning: invalid value encountered in multiply
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:7: RuntimeWarning: divide by zero encountered in reciprocal
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:7: RuntimeWarning: invalid value encountered in multiply
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:8: RuntimeWarning: divide by zero encountered in reciprocal
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:8: RuntimeWarning: invalid value encountered in multiply

















In [528]:



    


P_x*0.1*np.reciprocal(fb(P_x))



    


















    






/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in reciprocal
  if __name__ == '__main__':
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: invalid value encountered in multiply
  if __name__ == '__main__':










    
Out[528]:








array([-0.1, -0.1, -0.1,  0.1,  0.1, -0.1, -0.1, -0.1,  0.1,  0.1, -0.1,
       -0.1,  nan,  0.1,  0.1, -0.1, -0.1,  0.1,  0.1,  0.1, -0.1, -0.1,
        0.1,  0.1,  0.1])

















In [529]:



    


print P2_x
print P2_y



    


















    






[-2.1 -1.1 -0.1  1.1  2.1 -2.1 -1.1 -0.1  1.1  2.1 -2.1 -1.1  nan  1.1  2.1
 -2.1 -1.1  0.1  1.1  2.1 -2.1 -1.1  0.1  1.1  2.1]
[ 2.1  2.1  2.1  2.1  2.1  1.1  1.1  1.1  1.1  1.1  0.1  0.1  nan -0.1 -0.1
 -1.1 -1.1 -1.1 -1.1 -1.1 -2.1 -2.1 -2.1 -2.1 -2.1]

















In [530]:



    


print P_x*np.reciprocal(P_y)



    


















    






[ -0.90909091  -0.48672566  -0.04524887   0.41666667   0.90047393
  -1.70124481  -0.90909091  -0.04524887   0.90047393   1.94029851 -20.         -20.
          nan -20.         -20.           1.94029851   0.90047393
  -0.04524887  -0.90909091  -1.70124481   0.90047393   0.41666667
  -0.04524887  -0.48672566  -0.90909091]










    






/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in reciprocal
  if __name__ == '__main__':
/Users/heyfaraday/anaconda/envs/new2/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: invalid value encountered in multiply
  if __name__ == '__main__':

















In [531]:



    


P2_x



    


















    
Out[531]:








array([-2.1, -1.1, -0.1,  1.1,  2.1, -2.1, -1.1, -0.1,  1.1,  2.1, -2.1,
       -1.1,  nan,  1.1,  2.1, -2.1, -1.1,  0.1,  1.1,  2.1, -2.1, -1.1,
        0.1,  1.1,  2.1])

















In [532]:



    


import matplotlib.pyplot as plt
%matplotlib inline



    











In [533]:



    


plt.plot(P1_x, P1_y, 'ro')
plt.plot(P2_x, P2_y, 'go')
plt.plot(P3_x, P3_y, 'bo')

plt.show()



    


















    
























In [339]:



    


Ux = random.normal(mu, sigma, 1000000)
Uy = random.normal(mu, sigma, 1000000)
Qx = random.normal(mu, sigma, 1000000)
Qy = random.normal(mu, sigma, 1000000)



    











In [340]:



    


a, b, c = 0, 0, 0

for i in range(0, 100000):
    ch = ''
    ch = cubic(Qx[i], Qy[i], Ux[i], Uy[i])
    if (ch == 'a'):
        a = a + 1
    if (ch == 'b'):
        b = b + 1
    if (ch == 'c'):
        c = c + 1
        
print a/100000.0, b/100000.0, c/100000.0



    


















    






0.52896 0.09177 0.37927

















In [205]:



    


P_x



    


















    
Out[205]:








array([  6.,  11.,  16.,  21.,  26.,  -2.,   3.,   8.,  13.,  18., -10.,
        -5.,   0.,   5.,  10., -18., -13.,  -8.,  -3.,   2., -26., -21.,
       -16., -11.,  -6.])

















In [210]:



    


P_y



    


















    
Out[210]:








array([ 10.,  18.,  26.,  34.,  42.,  -3.,   5.,  13.,  21.,  29., -16.,
        -8.,   0.,   8.,  16., -29., -21., -13.,  -5.,   3., -42., -34.,
       -26., -18., -10.])

















In [ ]:

Content source: heyfaraday/CMB_test

Similar notebooks:

notebook.community | gallery | about