Вращающаяся черная дыра

using LinearAlgebra

function kerr_newman_metric(M, a, Q, r, θ)
    ρ² = r^2 + a^2*cos(θ)^2     # Функция горизонта и сингулярности
    Δ = r^2 - 2M*r + a^2 + Q^2  # Дискриминант для горизонтов событий
    
    g = zeros(4,4)
    g[1,1] = -(Δ - a^2*sin(θ)^2)/ρ²                      # g_tt - временная компонента
    g[1,4] = g[4,1] = -a*sin(θ)^2*(2M*r - Q^2)/ρ²        # g_tϕ - компонента вращения
    g[2,2] = ρ²/Δ                                        # g_rr - радиальная компонента  
    g[3,3] = ρ²                                          # g_θθ - угловая компонента
    g[4,4] = ((r^2+a^2)^2 - Δ*a^2*sin(θ)^2)*sin(θ)^2/ρ²  # g_ϕϕ - азимутальная компонента
    
    return g
end

function compute_everything(M, a, Q, r, θ)
    # Всё в одной функции
    g = kerr_newman_metric(M, a, Q, r, θ)  # Метрический тензор
    
    # Четырехвектор скорости для статического наблюдателя
    u = zeros(4); u[4] = 1.0  # Начальное условие: dt/dτ = 1
    
    # Нормировка: g_μν u^μ u^ν = -1 (для времениподобного вектора)
    norm_squared = dot(u, g, u)  # Это должно быть отрицательным для времениподобного вектора
    if norm_squared >= 0
        # Если вектор пространственноподобный, это проблема - используем безопасную нормировку
        u ./= sqrt(abs(norm_squared))
    else
        u ./= sqrt(-norm_squared)  # Корень из положительного числа
    end
    
    # Простой расчет кривизны через вторую производную
    h = 1e-6  # Шаг для численного дифференцирования
    g_rp = kerr_newman_metric(M, a, Q, r+h, θ)  # Метрика в точке r+h
    g_rm = kerr_newman_metric(M, a, Q, r-h, θ)  # Метрика в точке r-h
    R_simple = (g_rp[1,1] - 2g[1,1] + g_rm[1,1]) / h^2  # Конечная разность второй производной g_tt
    
    return g, u, R_simple
end

# Параметры черной дыры и координаты точки наблюдения
M = 1.0    # Масса черной дыры в натуральных единицах
a = 0.5    # Параметр вращения (a = J/M)
Q = 0.3    # Электрический заряд черной дыры
r = 2.0    # Радиальная координата (расстояние от центра)
θ = π/4    # Угловая координата (π/4 = 45°)

# Вычисление всех величин в заданной точке
g, u, R = compute_everything(M, a, Q, r, θ)

# Вывод результатов
println("Метрика Керра-Ньюмена в точке (r=$r, θ=$θ):")
display(g)
println("\nЧетырехвектор скорости статического наблюдателя:")
println("u = $u")  # Компоненты: [u^r, u^θ, u^ϕ, u^t]
println("\nУпрощенная мера кривизны (вторая производная g_tt):")
println("R ≈ $R")

using Plots.PlotMeasures

# Создаем сетку для визуализации
r_range = range(1.5, 5.0, length=50)  # Радиус от 1.5 до 5 (избегаем сингулярность)
θ_range = range(0.1, π-0.1, length=50)  # Угол от 0.1 до π-0.1 (избегаем полюса)

# Вычисляем компоненты метрики и кривизны на сетке
g_tt_values = [kerr_newman_metric(M, a, Q, r, θ)[1,1] for r in r_range, θ in θ_range]
g_rr_values = [kerr_newman_metric(M, a, Q, r, θ)[2,2] for r in r_range, θ in θ_range]
curvature_values = [compute_everything(M, a, Q, r, θ)[3] for r in r_range, θ in θ_range]

p1 = heatmap(r_range, θ_range, g_tt_values', 
            title="Временная компонента метрики g_tt",
            ylabel="Угол θ",
            color=:viridis, clim=(-2, 0), titlefont=font(12), topmargin=2mm, titlelocation=:left)

p2 = heatmap(r_range, θ_range, g_rr_values', 
            title="Радиальная компонента метрики g_rr", 
            ylabel="Угол θ",
            color=:plasma, clim=(0, 10), titlefont=font(12), topmargin=2mm, titlelocation=:left)

p3 = heatmap(r_range, θ_range, curvature_values',
            title="Упрощенная мера кривизны",
            xlabel="Радиус r", ylabel="Угол θ", 
            color=:hot, clim=(-1, 1), titlefont=font(12), topmargin=2mm, titlelocation=:left)

# Отображаем все графики вместе
plot(p1, p2, p3, layout=(3,1), size=(800, 900))