codigos_python/Blackhole1 at main · carlosroxo1/codigos_python · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import matplotlib.pyplot as plt
import math
import numpy as np
from mpl_toolkits import mplot3d
N=200
theta=1.7
X=[]
Y=[]
Z=[]
Xi=[]
Yi=[]
Zi=[]
Vx=[]
Vy=[]
Vz=[]
C=[]
m=1

xmin = -8
xmax = 8
ymin = -8
ymax = 8
for x in range(N):
  for y in range(N):
    X.append(x*(xmax-xmin)/N+xmin)
    Y.append(y*(ymax-ymin)/N+ymin)
    Z.append(-10)
    Vx.append(0)
    Vy.append(0)
    Vz.append(1)
    C.append([0,0,0])
    Xi.append(x*(xmax-xmin)/N+xmin)
    Yi.append(y*(ymax-ymin)/N+ymin)
    Zi.append(-10)

N = len(X)
dt = 0.1
for i in range(N):
  x = X[i]
  y = Y[i]
  z = Z[i]
  vx = Vx[i]
  vy = Vy[i]
  vz = Vz[i]
  for t in range(250):
    r = math.sqrt(x*x+y*y+z*z)
    a=m/(r*r)
    ax = -a*x/r
    ay = -a*y/r
    az = -a*z/r
    vx=vx+ax*dt
    vy=vy+ay*dt
    vz=vz+az*dt

    v = math.sqrt(vx*vx+vy*vy+vz*vz)
    if (v!=0 and v!=1):
      vx=vx/v
      vy=vy/v
      vz=vz/v

    x=x+vx*dt
    y=y+vy*dt
    z=z+vz*dt

    x1 = x
    y1 = math.cos(theta)*y-math.sin(theta)*z
    z1 = math.sin(theta)*y+math.cos(theta)*z
    if math.pow(5-math.sqrt(x1*x1+y1*y1),2)+math.pow(z1*5,2)<4:
      C[i]=[0,1,0]
      break

plt.scatter(Xi,Yi,color=C)