# Pkg.add(["Plots", "DifferentialEquations"])
using DifferentialEquations, Printf, DataFrames
gr()

G = 1.0      # 引力常数（归一化）
c = 1.0      # 光速（归一化）
M = 1.0      # 黑洞的质量（单位为GM/c2=1->r_s=2）
r_s = 2*G*M/c^2   # Schwarzschild半径=2
println("📊Schwarzschild半径：r_s=$r_s")

r_newton = range(0.1, 10, length=300)
v_escape_newton = @. sqrt(2*G*M/r_newton)          # 经典配方
r_einstein = range(r_s, 10, length=300)
v_escape_einstein = @. c * sqrt(r_s / r_einstein)  # 光子的相对论

p1 = plot(r_newton, v_escape_newton,
     label="牛顿",
     linewidth=3, color=:blue,
     xlabel="距离中心r",
     ylabel="逃逸速度v/c",
     title="预测黑洞：牛顿vs爱因斯坦",
     legend=:topright, grid=true, ylim=(0, 2.5))
plot!(p1, r_einstein, v_escape_einstein,
      label="光子",
      linewidth=3, color=:red, ls=:dash)
hline!(p1, [c], ls=:dot, color=:black, label="光速", linewidth=2)
vline!(p1, [r_s], ls=:dash, color=:green, label="地平线", linewidth=2)
annotate!(p1, r_s+0.5, 2.2, text("事件视界", :green, :bold, 11))
annotate!(p1, 3.0, 1.8, text("随着r<r_s：v_esc>c→没有什么可以离开！", :red, 10))
display(p1)

r_time = range(r_s*1.01, 10, length=300)
time_dilation = @. sqrt(1 - r_s / r_time)           # dτ/dt
redshift = @. 1 / sqrt(1 - r_s / r_time) - 1;       # z = Δλ/λ

p2 = plot(r_time, time_dilation,
     linewidth=3, color=:purple,
     xlabel="距离r", ylabel="Dt/dt比率",
     title="引力时间扩张",
     legend=:bottomright, grid=true, ylim=(0, 1.1))
hline!(p2, [1.0], ls=:dot, color=:green, label="远离BH（时间正常）", linewidth=1)
vline!(p2, [r_s], ls=:dash, color=:red, label="地平线（时间停止）", linewidth=2)
annotate!(p2, r_s+0.3, 0.1, text("在地平线上：时间停止！", :red, :bold, 9))
annotate!(p2, 5.0, 0.9, text("远离黑洞：时间正常运行", :green, 9))
annotate!(p2, 2.0, 0.7, text("在R=1.5r_s:1小时近Bh=1.7小时在地球上", :purple, 8))
display(p2)

p3 = plot(r_time, redshift,
     linewidth=3, color=:orange,
     xlabel="距离r", ylabel="红移z",
     title="引力红移",
     legend=:topright, grid=true, ylim=(0, 10))
vline!(p3, [r_s], ls=:dash, color=:red, label="地平线(z→←)", linewidth=2)
display(p3)

ΔL = 2.0   # 人体身高~2米
r_tidal = range(r_s*1.01, 20, length=300)
tidal_accel = @. (2 * G * M / r_tidal^3) * ΔL;

p4 = plot(r_tidal, tidal_accel,
     linewidth=3, color=:darkred,
     xlabel="距离r", ylabel="潮汐加速度（相对单位）",
     title="潮汐力—"意大利面条"",
     legend=:topright, grid=true)

r_spaghetti = (2*G*M*ΔL / 100)^(1/3) # 有条件的破坏阈值
vline!(p4, [r_spaghetti], ls=:dot, color=:blue, label="破坏的门槛")
annotate!(p4, r_spaghetti+0.5, maximum(tidal_accel)*0.8,
          text("对于r<$(round(r_spaghetti,digits=2))，对象正在被撕开", :blue, 8))
display(p4)

function radial_fall_full!(du, u, p, τ)
    r = u[1]
    E = p[1]
    r_s = p[2]
    radicand = E^2 - (1 - r_s/r)
    if radicand < 0
        du[1] = 0.0
    else
        du[1] = -sqrt(radicand)   
    end
end
E = 1.0                    
u0 = [10.0]                
τspan = (0.0, 40.0)        
prob_rad = ODEProblem(radial_fall_full!, u0, τspan, (E, r_s))
sol_rad = solve(prob_rad, Tsit5(), saveat=0.05, dtmax=0.1)
r_vals = sol_rad[1,:]
τ_vals = sol_rad.t
t_vals = zeros(length(τ_vals))
last_t = 0.0
for i in 2:length(τ_vals)
    r_mid = (r_vals[i-1] + r_vals[i]) / 2
    if r_mid > r_s
        dt_dτ = E / (1 - r_s / r_mid)
        last_t += dt_dτ * (τ_vals[i] - τ_vals[i-1])
        t_vals[i] = last_t
    else
        t_vals[i] = last_t
    end
end
println("相对论的堕落:")
println("   到达地平线的适当时间:$(round(t_vals[end],digits=2))")
println("   坐标时间到地平线:$(round(t_vals[end],digits=2))(diverges)")

dilation_vals = [r > r_s ? sqrt(1 - r_s/r) : 0.0 for r in r_vals]

anim = @animate for i in 1:5:length(τ_vals)
    r_cur = r_vals[i]
    τ_cur = τ_vals[i]
    t_cur = t_vals[i]
    dil_cur = dilation_vals[i]
    title_str = @sprintf "Bh下降：t=％。2f，τ=%。2f，dtt/dt=%。3f" t_cur τ_cur dil_cur
    
    p1 = plot(τ_vals[1:i], r_vals[1:i],
              color=:blue, linewidth=2,
              label="r(τ)",
              xlabel="Τ的适当时间",
              ylabel="半径r",
              legend=:bottomleft,
              title=title_str,
              ylim=(0, maximum(r_vals)+0.5))
    hline!(p1, [r_s], ls=:dash, color=:red, label="地平线r_s")
    scatter!(p1, [τ_cur], [r_cur], color=:orange, ms=6, label="粒子;粒子")

    p2 = plot(τ_vals[1:i], dilation_vals[1:i],
              color=:purple, linewidth=2,
              label="dτ/dt",
              xlabel="Τ的适当时间",
              ylabel="dτ/dt",
              legend=:bottomleft,
              ylim=(0, 1.1))
    hline!(p2, [0], ls=:dot, color=:red, label="地平线和内部")
    scatter!(p2, [τ_cur], [dil_cur], color=:red, ms=6, label="当前值")

    p3 = plot(t_vals[1:i], τ_vals[1:i],
              color=:green, linewidth=2,
              label="τ(t)",
              xlabel="坐标时间t",
              ylabel="Τ的适当时间",
              legend=:bottomleft)
    plot!(p3, t_vals[1:i], t_vals[1:i], ls=:dot, color=:gray, label="τ=t（无BH）")
    scatter!(p3, [t_cur], [τ_cur], color=:orange, ms=6, label="粒子;粒子")
    x_min = minimum(t_vals[1:i])
    y_min = minimum(τ_vals[1:i])
    annotate!(p3, x_min + 0.05*(maximum(t_vals[1:i])-x_min),
                   y_min + 0.05*(maximum(τ_vals[1:i])-y_min),
              text("地平线在τ≈13.55处交叉", :black, :left, 8))

    plot(p1, p2, p3, layout=(3,1), size=(800, 900))
end

gif(anim, "relativistic_fall_full.gif", fps=20)

r_photon = 1.5 * r_s
r_ISCO = 3.0 * r_s

println("附加轨道:")
println("   光子球：r=△r_photon")
println("   最后稳定轨道（ISCO）：r=$r_ISCO")

L_isco = 2.0 * sqrt(3) * M
r_pot = range(r_s+0.01, 20, length=500)
V_eff = @. (1 - r_s/r_pot) * (1 + L_isco^2 / r_pot^2)

p5 = plot(r_pot, V_eff,
     label="V_eff(r) = (1 - r_s/r) * (1 + L^2 / r^2)",
     linewidth=2, color=:darkgreen,
     xlabel="距离r", ylabel="V_eff的有效电位",
     title="有效电位：ISCO和光子球",
     legend=:topright, grid=true)
vline!(p5, [r_photon], ls=:dash, color=:cyan, label="光子球(1.5r_s)")
vline!(p5, [r_ISCO], ls=:dash, color=:magenta, label="ISCO (3 r_s)")
hline!(p5, [1.0], ls=:dot, color=:gray, label="E=1（在无穷远处休息）")
display(p5)
savefig("effective_potential.png")
println(" 图5保存并显示：efficie_potential。巴布亚新几内亚")

G_si = 6.67430e-11
c_si = 299792458
M_sun = 1.989e30
ΔL = 2.0          # 人的身高，m
g_earth = 9.8     # m/s2

r_s_km(mass_solar) = 2 * G_si * (mass_solar * M_sun) / c_si^2 / 1000 # Schwarzschild半径（公里）
function r_tidal_km(mass_solar)
    M_kg = mass_solar * M_sun
    r_m = (2 * G_si * M_kg * ΔL / g_earth)^(1/3)
    return r_m / 1000
end
function dilation_factor(r_km, r_s_km)
    if r_km <= r_s_km
        return Inf
    else
        return 1 / sqrt(1 - r_s_km / r_km)
    end
end

objects = [
    ("恒星质量的BH", 10),
    ("射手座A*", 4.3e6),
    ("M87*", 6.5e9),
    ("TON 618", 6.6e10)
]

rows = []
for (name, mass) in objects
    rs = r_s_km(mass)
    r_tidal = r_tidal_km(mass)
    r_min = max(r_tidal, rs * 1.0001)
    dilation = dilation_factor(r_min, rs)
    push!(rows, (name, mass, rs, r_tidal, r_min, dilation))
end

df = DataFrame(
    Symbol("一个对象") => [row[1] for row in rows],
    Symbol("M (M☉)") => [row[2] for row in rows],
    Symbol("r_s(公里)") => [round(row[3], sigdigits=4) for row in rows],
    Symbol("潮(公里)") => [round(row[4], sigdigits=4) for row in rows],
    Symbol("r_min(公里)") => [round(row[5], sigdigits=4) for row in rows],
    Symbol("减速(1小时→小时)") => [row[6] for row in rows]
)

println("表：对于每个对象，表示可能不被潮汐力撕裂的最小距离（但不超过1.0001r_s），以及相应的时间膨胀。")
println(df)