Rotation Curve

import numpy as np
import matplotlib.pyplot as plt

# Константы
G = 6.67430e-11 # Гравитационная постоянная
Msun = 1.98847e30 # Масса Солнца

# Создание массива случайных позиций звезд
N = 1000 # Количество звезд
Rmax = 30e3 # Максимальное расстояние от центра галактики в парсеках
R = Rmax*np.sqrt(np.random.rand(N))
theta = 2*np.pi*np.random.rand(N)
x = R*np.cos(theta)
y = R*np.sin(theta)

# Создание массива случайных скоростей звезд
vmax = 1000 # Максимальная скорость в км/с
v = vmax*np.sqrt(np.random.rand(N))

# Моделирование движения галактики
tmax = 50 # Время моделирования в миллион лет
dt = 0.01 # Шаг по времени в миллион лет
t = np.arange(0, tmax, dt)
xgal = np.zeros_like(t)
ygal = np.zeros_like(t)
vgal = np.zeros_like(t)

for i in range(len(t)):
    # Расстояние от каждой звезды до центра галактики
    r = np.sqrt(x**2 + y**2)

    # Суммирование гравитационных сил от всех звезд
    ax = -G*np.sum(Msun*x/r**3)
    ay = -G*np.sum(Msun*y/r**3)

    # Добавление гравитационной силы от центральной массы
    ax -= G*4*Msun*x[0]/r[0]**3
    ay -= G*4*Msun*y[0]/r[0]**3

    # Интегрирование ускорения для получения скорости и позиции
    vx = v*np.cos(theta) + ax*dt
    vy = v*np.sin(theta) + ay*dt
    x += vx*dt
    y += vy*dt

    # Перенос центра координат в центр галактики
    x -= x[0]
    y -= y[0]

    # Вычисление положения центра масс галактики
    xcm = np.sum(Msun*x)/np.sum(Msun)
    ycm = np.sum(Msun*y)/np.sum(Msun)

    # Вычисление скорости и положения центра масс галактики
    vgal[i] = np.sqrt(np.sum((v - np.sqrt(ax**2 + ay**2))**2)/N)

    # Вычисление радиуса и скорости вращения для каждой звезды
    r = np.sqrt(x**2 + y**2)
    vrot = np.sqrt(G*np.sum(Msun)/r)

# Вычисление средней скорости вращения для каждого радиуса
rbin = np.linspace(0, Rmax, 100)
vrot_mean = np.zeros_like(rbin)
for i in range(len(rbin)-1):
    idx = np.where((r >= rbin[i]) & (r < rbin[i+1]))[0]
    vrot_mean[i] = np.mean(vrot[idx])

# Построение графика кривой вращения
plt.plot(rbin[:-1], vrot_mean)
plt.xlabel('Radius (pc)')
plt.ylabel('Rotation Velocity (km/s)')
plt.show()