Get best of the week

Modeling the universe: celestial mechanics clearly



Let's imagine that we need to launch a soccer ball into the orbit of the Earth. No rockets needed! Enough mountains, a height of 100 kilometers and remarkable strength. But how much do you need to kick the ball so that it never returns to Earth? How to send the ball on a journey to the stars, having only brute force and knowledge of celestial mechanics?

Today in the program:

  • Infinite possibilities of one formula
  • How to take energy from Jupiter
  • Where do planets come from
  • How mathematics helped discover Neptune

Fortunately, we live in the age of computer technology. We do not need to climb a high mountain and kick the ball with all our strength, everything can be modeled! Let's get started.

One formula


The same one, known from the lessons of physics and astronomy:



It shows how strongly the bodies will interact, depending on their masses, the distance between them and the gravitational constant G.

I wrote a program in which you can place balls interacting with each other by gravitational forces, with Each ball has its own mass, speed and coordinates. For clarity, the balls leave a trail.

Let's put a large and massive blue ball (Earth) and a small red ball near it. Run the simulation:



He fell!

To enter orbit, speed is needed so that the ball falls and misses the earth all the time. But WHAT speed? And again, school knowledge comes to the rescue:

The minimum speed required to enter the Earth’s orbit is calledfirst cosmic speed .

For the Earth, it is 7.91 km / s. And for simulation it can be easily calculated: We



disperse the ball and see the result:



Normal flight!

The ball describes a circle with the Earth in the center. What will happen if you give it a little more speed? Now



let's check: Now the shape of the orbit is elliptical, we can distinguish 2 very important points - apogee and perigee .

Apogee is the point at which the ball is as far from the Earth as possible.

Perigee - on the contrary, the point closest to the Earth.

With an increase in the initial velocity, the perigee does not change, but the apogee is getting farther, and in the end has an infinite distance to the Earth. Here we come close to the conceptsecond space velocity . This is the speed that must be given to the ball so that it overcomes the gravity of the Earth and flies away to plow the expanses of the universe. For land, it is 11.2 km / s.

An interesting trick: if we multiply the first cosmic velocity by √2, we get the second cosmic velocity.

Multiplied. Launched. Received:



He flew away irrevocably! By the way, now it has a parabolic orbit. And if you run the ball even stronger, we get a hyperbola. It turns out interestingly, mathematics haunts us everywhere.

In this case, the formula remains the same. The circle turns into an ellipse, an ellipse into a parabola, and a parabola into a hyperbola due to the elongation of the orbit (increased eccentricity ).

How to take energy from Jupiter?


Let's expand our model, add the Sun, make the Earth revolve around it.



Imagine that the ball needs to be given such speed that it flies outside the solar system - the third cosmic speed . In the real world, it is 16.7 km / s. Unfortunately, this speed is too high, I'm afraid we will not have enough strength ...

Wait! But what if you take a little speed from some massive body, for example, Jupiter. We can fly up to something very massive and make a gravitational maneuver . When flying past Jupiter, the gravitational forces mutually attract the ball and the gas giant, but the mass of the ball is so small that it has almost no effect on the movement of Jupiter, and Jupiter itself accelerates a body passing by at high speeds.

Talk less, work more:



The moment of the gravitational maneuver - the ball flew up to Jupiter.



Hooray! We got a speed sufficient to exit the solar system, while not spending anything. True, Jupiter began to move a little slower, but we definitely will not notice it.

All spacecraft launched by man beyond the limits of the solar system (Voyagers 1 and 2, Pioneers 10 and 11, New Horizons) used just such a method for acceleration.

Zoom in!


I added the friction of the particles so that when they collide, they transfer part of the energy to each other. I also introduced the power of a normal reaction, now the particles respect their personal space, pushing others away from themselves.

We put the random generation of balls and give them a random direction and speed. Let them be, say, 100 pieces.



Complete chaos, each particle moves wherever it wants, but nevertheless, the gravitational forces take their toll and clusters of balls begin to form:



And after a while a large body is formed, consisting of 99 balls and one single ball orbiting around it:



With another launch, the following :



Two massive bodies orbiting a common center of mass. If we imagine that these two objects are stars, then we get a double star. Interestingly, about half of the stars in our galaxy are binary. If our Sun had a companion star, then in the sky we could observe the following picture:



Where do planets come from?


The main reason for the appearance of rings is the destruction of satellites flying too close to the planet, or rather, crossing the Roche limit . In this case, the tidal forces caused by the gravity of the planet become larger than the forces holding the satellite intact, and it breaks into many parts, leaving behind a ring that encircles the planet. Let's simulate this situation:



The satellite is a little further than the Roche limit, it rotates around the planet in a stable circular orbit. But what happens if you generate it a little closer to the planet?



The satellite scattered into many small parts that formed rings around the planet. So it is in the real world. Triton (satellite of Neptune) is gradually approaching the planet, and in 2 billion years it will be torn, and Neptune will have more rings than Saturn.

How was Neptune discovered and what does mathematics have to do with it?


Since we are talking about Neptune, let's talk about its discovery. “A planet open at the tip of a pen” has a mass, which means it acts on objects around. Astronomers of the 19th century noticed changes in the orbit of Uranus, its orbit was different from the calculated one, apparently, something affected it. The orbit of Uranus had disturbances:



This exaggerated model shows how an unknown body beyond Uranus affected its orbit. Astronomers could only calculate the position of a secret planet and look through a telescope. Indeed, the planet Neptune was exactly where it was predicted!



Conclusion


Of course, this simulation does not generalize all the laws and phenomena occurring in space, for example, Einstein's theory of relativity is not taken into account here, since the speed of particles is far from the speed of light. But there are many more interesting things that can be implemented in this simulation. Try it yourself! All you need is Python3 and the Pygame library.

Source
# 
Track = True

# 
Track_time = 5

# 
G = 5

#  (  -
#,  - )
MagnConst = 0

# 
count = 100

#  
kv = 6

#  
RANDOM = True

#  
r = 3

#   
WIN_WIDTH, WIN_HEIGHT = 900, 650


''',  ,   '''

#   
zg = 2

#   
zm = 2

#. ,    -   
k = 40

# 
antiG = 0.1

max_speed = 3

ResDist = 1

# 
EarthG = 0

#   
Mirror = True

import pygame
from math import hypot, ceil, sqrt
from random import randint, random


def custom_pos():
    '''      '''
    '''   RANDOM = FALSE'''
    B.append(Ball(200, 300, YELLOW, r = 10, mass = 200, vx = 0.151))    #x, y, col, r, vx, vy, mass
    B.append(Ball(200, 50, GREEN, r = 6, mass = 10, vx = -(200 * G / 250)**0.5))
    

class Ball:
    def __init__(self, x, y, col, r = 4, vx = 0, vy = 0, mass = 4):
        self.x = x
        self.y = y
        self.r = r
        self.col = col
        self.vx = vx
        self.vy = vy
        self.mass = mass
        
    def move(self, Walls, WIN_WIDTH, WIN_HEIGHT, ofs_x, ofs_y):
        if Walls:
            x = self.x - ofs_x
            y = self.y - ofs_y
            if x <= 0 and self.vx < 0:
                if Mirror:
                    self.vx = -self.vx
                else:
                    self.x += WIN_WIDTH
                self.vx, self.vy = self.v_norm(self.vx, self.vy)
            if x >= WIN_WIDTH and self.vx > 0:
                if Mirror:
                    self.vx = -self.vx
                else:
                    self.x -= WIN_WIDTH
                self.vx, self.vy = self.v_norm(self.vx, self.vy)
            if y <= 0 and self.vy < 0:
                if Mirror:
                    self.vy = -self.vy
                else:
                    self.y += WIN_HEIGHT
                self.vx, self.vy = self.v_norm(self.vx, self.vy)
            if y >= WIN_HEIGHT and self.vy > 0:
                if Mirror:
                    self.vy = -self.vy
                else:
                    self.y -= WIN_HEIGHT
                self.vx, self.vy = self.v_norm(self.vx, self.vy)
            
        self.x += self.vx
        self.y += self.vy

        
    def force(self, ind, selfind):
        ox = B[ind].x
        oy = B[ind].y
        m = B[ind].mass
        if m < 0.01 and self.mass < 0.01:
            return 0
        r = B[ind].r
        vx = B[ind].vx
        vy = B[ind].vy
        dist = hypot(self.x - ox, self.y - oy)
        min_dist = (self.r + B[ind].r) * ResDist
        f = 0
        m_relative = self.mass / B[ind].mass
        if dist <= min_dist:
            newVx = (vx * m + self.vx * self.mass) / (m + self.mass)
            newVy = (vy * m + self.vy * self.mass) / (m + self.mass)
            self.vx = (newVx + k * self.vx) / (k + 1)
            B[ind].vx = (newVx + k * B[ind].vx) / (k + 1)
            self.vy = (newVy + k * self.vy) / (k + 1)
            B[ind].vy = (newVy + k * B[ind].vy) / (k + 1)
            f -= antiG * min(abs(min_dist - dist), min(m, self.mass) * 3)
        else:
            f += min(self.mass * B[ind].mass * G  / (dist ** zg), G / 10)
            mf = MagnConst * self.mass / (dist ** zm)
            if B[ind].col == B[selfind].col:
                mf = - mf
            f += mf
        fx = f * ((ox - self.x) / dist)
        fy = f * ((oy - self.y) / dist)
        ax = fx / self.mass
        ay = fy / self.mass
        self.vx += ax
        self.vy += ay + EarthG
        B[ind].vx -= ax * m_relative
        B[ind].vy -= ay * m_relative - EarthG

    @staticmethod
    def v_norm(vx, vy):
        v = hypot(vx, vy)
        if v > max_speed:
            vx = max_speed * (vx / v)
            vy = max_speed * (vy / v)
        return vx, vy


class Point:
    def __init__(self, x, y, col, r = 0, max_age = Track_time):
        self.age = 0
        self.x = x
        self.y = y
        self.col = col
        self.r = r
        self.max_age = max_age
    def vis(self, ofs_x, ofs_y):
        pygame.draw.circle(sc, self.col, (round(self.x - ofs_x),
                                          round(self.y - ofs_y)), self.r, 0)
        self.age += 1
        if self.age > self.max_age:
            T.remove(self)
        
def rand(count, WIN_WIDTH, WIN_HEIGHT):
    global kv
    B = []
    for i in range(count):
        m = r ** 2
        x = randint(0, WIN_WIDTH) + random()
        y = randint(0, WIN_HEIGHT) + random()
        vx = kv * randint(-100, 100) / 100
        vy = kv * randint(-100, 100) / 100
        col = Colors[randint(0, len(Colors) - 1)]
        B.append(Ball(x, y, col, r = r, vx = vx, vy = vy, mass = m))
    return B

def createBall(col, x, y, r = r, m = r):
    m = r
    B.append(Ball(x, y, col))

def get_offset(B):
    sum_x, sum_y = 0, 0
    m = 0
    for i in range(len(B)):
        sum_x += B[i].x * B[i].mass
        sum_y += B[i].y * B[i].mass
        m += B[i].mass
    if len(B) == 0:
        return 0, 0
    return sum_x / m, sum_y / m

def visBalls(B):
    for i in range(len(B)):
        pygame.draw.circle(sc, B[i].col, (round(B[i].x - ofs_x),
                                          round(B[i].y - ofs_y)), B[i].r, 0)
        T.append(Point(B[i].x, B[i].y, B[i].col))
        
FPS = 60
darkblue = (0, 2, 25)
ORANGE = (255, 200, 150)
RED = (255, 150, 150)
GREEN = (150, 255, 150)
BLUE = (150, 150, 255)
YELLOW = (255, 255, 0)
Colors = [RED, BLUE]#, GREEN]#, ORANGE]                       
pygame.init() 
clock = pygame.time.Clock()
sc = pygame.display.set_mode((WIN_WIDTH, WIN_HEIGHT))
sc.fill(darkblue)

maxF = 0.3
minv = 0.01
Walls = True
Collisions = True
Same = True
Check = False
tt = []

B = []
if RANDOM:
    B = rand(count, WIN_WIDTH, WIN_HEIGHT)
else:
    custom_pos()
    
Pause = False
delay = 0

if Track:
    T = []
for z in range(100000):
    sc = pygame.display.set_mode((WIN_WIDTH, WIN_HEIGHT))
    sc.fill(darkblue)
    ofs_x, ofs_y = get_offset(B)
    ofs_x -= WIN_WIDTH // 2
    ofs_y -= WIN_HEIGHT // 2
    for i in pygame.event.get():
        if i.type == pygame.QUIT:
            pygame.quit()
            quit()
        if i.type == pygame.KEYDOWN:
            if i.key == pygame.K_SPACE:
                Pause = not Pause
            elif i.key == pygame.K_w:
                WIN_HEIGHT += 10
                WIN_WIDTH += 10
            elif i.key == pygame.K_s:
                WIN_HEIGHT -= 10
                WIN_WIDTH -= 10
                
    pressed = pygame.mouse.get_pressed()
    pos = pygame.mouse.get_pos()
    x = pos[0]
    y = pos[1]
    if pressed[0] and delay < 0:
        delay = 20
        createBall(RED, x + ofs_x, y + ofs_y)
    if pressed[2] and delay < 0:
        delay = 20
        createBall(BLUE, x + ofs_x, y + ofs_y )
    delay -= 1
    
    if not Pause:
        for i in range(len(B)):
            for j in range(i + 1, len(B)):
                B[i].force(j, i)
        for i in range(len(B)):
            B[i].move(Walls, WIN_WIDTH, WIN_HEIGHT, ofs_x, ofs_y)
    for i in range(len(T)):
        try:
            T[i].vis(ofs_x, ofs_y)
        except IndexError:
            pass
    visBalls(B)
    
    pygame.display.update()
    clock.tick(FPS)


I hope this article has been informative for you. Thank you for the attention!

More articles:


All Articles