Simulation of the Solar system.
In this article, we will look at how to create a full—fledged model of the Solar System that not only demonstrates gravitational interactions, but also allows us to study the consequences of changes in key parameters of the system, such as the mass of the central star. The presented implementation expands the classical problems of celestial mechanics, including thermodynamic aspects, which makes it possible to study the complex effect of changes in the mass of the Sun on both the orbital dynamics and the climatic conditions of the planets.
Key aspects of implementation:
-
Implementation of Newton's law of universal gravitation for a system of 9 interacting bodies (the Sun + 8 planets), where each body affects all the others, not just the central one.
-
The model includes a fundamental mass-luminosity ratio for stars (L~ M3·⁵), which allows us to estimate how a change in the mass of the Sun affects its energy release. Based on this relationship, the equilibrium temperatures of the planets are calculated in degrees Celsius.
-
A separate calculation and visualization system has been created, showing how the temperature of the planets changes with varying solar masses in the range from 1 to 10 solar masses. The results are presented in the form of graphs and tables with temperature markers (🔥 for temperatures above 1000°C, ❄ — below 0°C).
-
Using the Euler method to solve equations of motion with adaptive pitch, providing a balance between accuracy and performance when modeling on time scales up to 10 years.
This project serves as an excellent starting point for studying Engee, demonstrating how a programming language handles tasks that require computational efficiency, clear code, and visual visualization. In the following sections, we will examine in detail the architecture of the code, the physical foundations, the implementation of temperature calculations and the educational aspects of this model, as well as consider how such simulations can be used in the educational process and popular science projects.
Now let's move on to analyzing the code blocks and start with the constants used in the work.
Gravitational constant (G = 6.67430×10-11 m3·s⁻2·kg⁻1) is a fundamental physical constant that determines the force of gravitational attraction between bodies with mass. In classical Newtonian mechanics, it is included in the law of universal gravitation, according to which the force of attraction between two bodies is proportional to the product of their masses and inversely proportional to the square of the distance between their centers. It is this constant that connects the geometry of space-time with the distribution of matter, determining the intensity of gravitational interaction in the non-relativistic limit. In our simulation, the gravitational constant allows us to calculate the acceleration with which the Sun attracts the planets, and vice versa, forming realistic orbital dynamics.
Astronomical unit (AU = 1.496×1011 meters) is a historically established unit of measurement for distances in astronomy, approximately equal to the average distance from the Earth to the Sun. The introduction of this constant into the model serves two important purposes. First, it provides a natural scale for representing the orbits of planets, since the semi-major axes of their orbits are traditionally expressed in astronomical units. Secondly, using AU simplifies the perception and adjustment of model parameters, since researchers are used to operating with these values in an astronomical context, rather than in absolute metric units.
Stefan-Boltzmann constant (σ = 5.670374419×10⁻⁸ W·m⁻2·K⁻⁴) is a fundamental physical constant describing the Stefan-Boltzmann law for blackbody radiation. According to this law, the total energy radiated by a unit of the surface of a blackbody per unit of time is proportional to the fourth power of its thermodynamic temperature. In the context of our model, this constant is a key element for calculating the equilibrium temperature of planets. It allows you to switch from the luminosity of the Sun (determined through mass) to the heat flow reaching the planet, and further to the calculated temperature of its surface, taking into account the albedo. It is thanks to this constant that the model acquires a thermodynamic component that connects gravitational dynamics with the climatic conditions on the planets.
const G = 6.67430e-11 # The gravitational constant
const AU = 1.496e11 # Astronomical unit (m)
const YEAR = 365.25 * 86400.0 # Year in seconds
const SOLAR_MASS = 1.989e30 # Mass of the Sun (kg)
const SOLAR_RADIUS = 6.957e8 # Radius of the Sun (m)
const STEFAN_BOLTZMANN = 5.670374419e-8
Structure Body It is a fundamental element of the model that encapsulates the complete physical and visual profile of a celestial body. Being mutable, it allows you to dynamically update the state of objects at each step of the simulation, reflecting their evolution over time.
It includes three categories of parameters:
- Dynamic characteristics (position
x, y, speedvx, vy, accelerationax, ay) - physical properties (mass
mass, physical and orbital radii) - Attributes for visualization and analysis (name, color, temperature).
The field is especially important temperature, which is calculated based on the luminosity of the Sun and the distance to the body.
mutable struct Body
name::String
mass::Float64 # Weight (kg)
radius::Float64 # Physical radius (m)
orbit_radius::Float64 # Initial orbital radius (m)
color::Symbol
x::Float64 # Position X (m)
y::Float64 # Y position (m)
vx::Float64 # Speed X (m/s)
vy::Float64 # Speed Y (m/s)
ax::Float64 # Acceleration X (m/s2)
ay::Float64 # Acceleration Y (m/s2)
temperature::Float64 # Equilibrium temperature (°C)
end
Array planets_data It contains the initial astrophysical parameters and visual attributes of the eight planets of the Solar System, which serve as the initial conditions for creating objects of the type Body in the simulation.
planets_data = [
("Mercury", 3.285e23, 2.4397e6, 0.39*AU, :orange),
("Venus", 4.867e24, 6.0518e6, 0.72*AU, :green),
("Earth", 5.972e24, 6.371e6, 1.0*AU, :blue),
("Mars", 6.39e23, 3.3895e6, 1.52*AU, :red),
("Jupiter", 1.898e27, 6.9911e7, 5.2*AU, :brown),
("Saturn", 5.683e26, 5.8232e7, 9.58*AU, :gold),
("Uranus", 8.681e25, 2.5362e7, 19.22*AU, :cyan),
("Neptune", 1.024e26, 2.4622e7, 30.05*AU, :darkblue)
]
Function solar_luminosity implements the astrophysical law "mass-luminosity" (L ∼ M3·⁵), converting the relative change in the mass of the Sun into an appropriate multiplier for its energy release. Function equilibrium_temperature_c Based on this luminosity and the distance to the planet, it calculates its equilibrium temperature in degrees Celsius using a model of energy absorption and re-emission, taking into account albedo and the Stefan-Boltzmann law.
function solar_luminosity(mass_factor)
return mass_factor^3.5 # L ~ M^3.5
end
function equilibrium_temperature_c(luminosity_factor, distance_au)
# Simplified formula for temperature in degrees Celsius
albedo = 0.3
solar_constant = 1361.0 # W/m2 for Land
# Temperature in Kelvin
temp_k = (solar_constant * luminosity_factor * (1 - albedo) /
(4 * STEFAN_BOLTZMANN))^0.25
# Distance adjustment and conversion to Celsius
return (temp_k / sqrt(distance_au)) - 273.15
end
Function create_temperature_stats performs a comprehensive thermodynamic analysis, calculating the equilibrium temperatures of all planets for a wide range of solar masses (from 1 to 10 solar masses), creates two types of graphs — the temperature dependence on the mass of the star and on the distance to it, displays a table of results and saves both graphs to files for future use, acting as an autonomous module for visualizing climate scenarios.
function create_temperature_stats()
sun_mass_factors = [1, 10, 100, 1000, 1e4, 1e5, 1e6, 1e7]
planets_info = [
("Mercury", 0.39, :orange),
("Venus", 0.72, :green),
("Earth", 1.0, :blue),
("Mars", 1.52, :red),
("Jupiter", 5.2, :brown),
("Saturn", 9.58, :gold),
("Uranus", 19.22, :cyan),
("Neptune", 30.05, :darkblue)
]
# Calculation of temperatures
temperatures = Dict()
for (name, distance_au, color) in planets_info
temps = Float64[]
for mass_factor in sun_mass_factors
lum_factor = solar_luminosity(mass_factor)
temp_c = equilibrium_temperature_c(lum_factor, distance_au)
push!(temps, temp_c)
end
temperatures[name] = temps
end
p1 = plot(size=(1000, 600),
title="The dependence of the temperature of the planets on the mass of the Sun",
xlabel="The mass of the Sun (in solar masses)",
ylabel="Temperature (°C)",
xscale=:log10,
yscale=:linear,
grid=true, gridstyle=:dash,
legend=:topleft,
framestyle=:box)
for (name, _, color) in planets_info
plot!(p1, sun_mass_factors, temperatures[name],
label=name, color=color, linewidth=2, marker=:circle)
end
hline!(p1, [0], label="Freezing point of water", color=:blue, linestyle=:dash)
hline!(p1, [100], label="Boiling point of water", color=:red, linestyle=:dash)
annotate!(p1, 1e5, 1500, text("Extreme conditions", 10, :gray))
p2 = plot(size=(1000, 600),
title="The temperature of the planets depends on the distance",
xlabel="Distance from the Sun (AU)",
ylabel="Temperature (°C)",
xscale=:log10,
yscale=:linear,
grid=true,
legend=:topright)
distances = [d for (_, d, _) in planets_info]
planet_names = [n for (n, _, _) in planets_info]
for (i, mass_factor) in enumerate([1, 10, 100, 1000])
temps = [temperatures[name][i] for name in planet_names]
plot!(p2, distances, temps,
label="M☉ = $(Int(mass_factor))×",
linewidth=2,
marker=:circle)
end
display(p1)
display(p2)
println("\n" * "="^70)
println("TEMPERATURES OF PLANETS AT DIFFERENT SUN MASSES (°C)")
println("="^70)
print("The mass of the Sun: ")
for mf in sun_mass_factors
print(lpad("$(Int(mf))×", 10))
end
println()
for (name, _, _) in planets_info
print(rpad(name, 12))
for temp in temperatures[name]
if temp > 1000
print(lpad("🔥$(round(Int, temp))", 10))
elseif temp > 100
print(lpad("$(round(Int, temp))", 10))
elseif temp > 0
print(lpad("$(round(Int, temp))", 10))
else
print(lpad("❄$(round(Int, temp))", 10))
end
end
println()
end
println("="^70)
savefig(p1, "temperature_vs_mass.png")
savefig(p2, "temperature_vs_distance.png")
return temperatures
end
Function calculate_initial_velocity calculates the initial orbital velocity of a celestial body, ensuring its movement along a circular trajectory around the Sun, based on the given position of the body and the mass of the central star in accordance with classical Newtonian mechanics.
function calculate_initial_velocity(sun_mass, body_x, body_y)
r = sqrt(body_x^2 + body_y^2)
v_circ = sqrt(G * sun_mass / r)
return (-body_y / r * v_circ, body_x / r * v_circ)
end
The following is the main function simulate_solar_system it initializes and performs a complete simulation cycle of the gravitational dynamics and thermal state of the system. It starts with creating an array of bodies bodies, which includes the Sun with a mass adjusted by a given coefficient sun_mass_factor, and all planets with calculated initial velocities for circular orbits. Integration parameters such as time step dt and the total number of steps n_steps, convert the input astronomical values (years, days) into computational intervals, preparing the system for numerical solution.
Within the main loop, the function recalculates the gravitational accelerations for all pairs of bodies at each step using Newton's law and updates their velocities and coordinates using the Euler method. At the same time, the current equilibrium temperature for the planets is calculated based on the luminosity of the Sun and the distance to it, which adds a thermodynamic context to a purely mechanical simulation. With the specified frequency, the system state is saved for subsequent visualization.
Upon completion of the calculations, the function creates an animated GIF file in which the frames sequentially display the evolution of the system: the Sun as a minimalistic dot in the center, planets colored according to their temperature (from cold blue to red-hot red), and a barely noticeable transparent zone for scale. The animation title dynamically reflects key parameters, including the temperature of the Earth, providing a visual representation of the results of a computational experiment on changing the mass of the central star.
function simulate_solar_system(sun_mass_factor=1.0; years=10, step_days=5, filename="solar_system.gif")
sun_mass = SOLAR_MASS * sun_mass_factor
sun_luminosity = solar_luminosity(sun_mass_factor)
bodies = [Body("Sun", sun_mass, SOLAR_RADIUS, 0.0, :yellow,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5778.0 - 273.15)]
for (name, mass, radius, orbit_r, color) in planets_data
vx, vy = calculate_initial_velocity(sun_mass, orbit_r, 0.0)
body = Body(name, mass, radius, orbit_r, color,
orbit_r, 0.0, vx, vy, 0.0, 0.0, 0.0)
push!(bodies, body)
end
dt = step_days * 86400.0
n_steps = Int(round(years * YEAR / dt))
frames_per_gif = 150
step_skip = max(1, n_steps ÷ frames_per_gif)
positions_history = []
function compute_accelerations!(bodies)
for body in bodies
body.ax = 0.0
body.ay = 0.0
end
n = length(bodies)
for i in 1:n
for j in i+1:n
dx = bodies[j].x - bodies[i].x
dy = bodies[j].y - bodies[i].y
r_sq = dx^2 + dy^2
r = sqrt(r_sq)
if r > 1e-6
force = G * bodies[i].mass * bodies[j].mass / r_sq
ax_i = force * dx / (r * bodies[i].mass)
ay_i = force * dy / (r * bodies[i].mass)
ax_j = -force * dx / (r * bodies[j].mass)
ay_j = -force * dy / (r * bodies[j].mass)
bodies[i].ax += ax_i
bodies[i].ay += ay_i
bodies[j].ax += ax_j
bodies[j].ay += ay_j
end
end
end
end
println("Simulation: M☉ = $(sun_mass_factor)×, $(years) years, step $(step_days) days")
for step in 0:n_steps
if step % step_skip == 0
pos = []
for body in bodies
if body.mass > 0.0
if body.name != "Sun"
dist_au = sqrt(body.x^2 + body.y^2) / AU
body.temperature = equilibrium_temperature_c(sun_luminosity, dist_au)
end
push!(pos, (body.x, body.y, body.color, body.name, body.temperature))
end
end
push!(positions_history, pos)
end
step == n_steps && break
compute_accelerations!(bodies)
for body in bodies
body.vx += body.ax * dt
body.vy += body.ay * dt
body.x += body.vx * dt
body.y += body.vy * dt
end
end
max_dist = 1.2 * maximum([b.orbit_radius for b in bodies if b.name != "Sun"])
anim = @animate for (frame, positions) in enumerate(positions_history)
time_years = frame * step_skip * dt / YEAR
earth_temp = 0.0
for (x, y, color, name, temp) in positions
if name == "Earth"
earth_temp = temp
break
end
end
p = plot(size=(1000, 800), aspect_ratio=:equal,
xlim=(-max_dist, max_dist), ylim=(-max_dist, max_dist),
title="M☉ = $(sun_mass_factor)× | L☉ = $(round(sun_luminosity, digits=1))× | $(round(time_years, digits=1)) years | T♁ = $(round(earth_temp, digits=1))°C",
legend=:none, framestyle=:box,
background_color=:black,
titlefontsize=12)
scatter!(p, [0.0], [0.0],
markersize=2, color=:yellow, marker=:circle,
markerstrokewidth=0, markercolor=:yellow, alpha=0.8)
theta = range(0, 2π, length=100)
zone_radius = 0.05 * AU
circle_x = zone_radius .* cos.(theta)
circle_y = zone_radius .* sin.(theta)
plot!(p, circle_x, circle_y, color=:lightgray, linewidth=0.5, alpha=0.15, label="")
for (x, y, color, name, temp) in positions
if name != "Sun"
body = first([b for b in bodies if b.name == name])
size_factor = 5 + 2 * log10(body.mass / 5.972e24)
marker_color = if temp > 1000
:red
elseif temp > 100
:orange
elseif temp > 0
:green
else
:lightblue
end
scatter!(p, [x], [y],
markersize=max(3, size_factor),
color=marker_color, marker=:circle,
markerstrokewidth=1)
end
end
plot!(p, [max_dist*0.7, max_dist*0.7 + AU], [max_dist*0.9, max_dist*0.9],
color=:white, linewidth=2)
annotate!(p, max_dist*0.7 + AU/2, max_dist*0.88,
text("1 AU", 10, :white))
end
gif(anim, filename, fps=20)
println("Animation saved: $filename")
return bodies
end
Let's run the test and look at the results.
create_temperature_stats()
simulate_solar_system(1.0, years=10, step_days=3, filename="normal_sun.gif");
simulate_solar_system(10.0, years=5, step_days=2, filename="massive_sun_10x.gif");
simulate_solar_system(10000.0, years=2, step_days=1, filename="extreme_sun_10000x.gif");
Conclusion
The tests performed demonstrate two key physical consequences of changes in the mass of the Sun.
Temperature analysis (create_temperature_stats()) clearly shows that the dependence of luminosity on mass leads to an extreme increase in the temperature of the planets. Even at 10 times the mass of the Sun, temperatures on the inner planets exceed the boiling point of water, and when the mass increases to 10–-10⁷ solar, all planets heat up above 1000 °C, which illustrates the fundamental relationship between the structure of the star and the climatic conditions in the system.
Simulations of orbital dynamics confirm that an increase in the mass of the Sun accelerates the orbital motion of the planets (according to Kepler's law) and changes the stability of their trajectories due to the dominance of gravity of the central body. Animations visualize these effects: at normal mass (normal_sun.gif) planets exhibit classical Keplerian orbits, and with its growth (massive_sun_10x.gif, extreme_sun_10000x.gif their movement becomes faster, and the color changes from cold shades to red-hot.