Распространение напряженности магнитного поля

try
    using Plots
catch
    error("Пакет Plots не установлен. В терминале Julia выполните:\nimport Pkg; Pkg.add(\"Plots\")")
end

# 1. Физические параметры и константы

f  = 2.45e9          # Частота волны (Гц)
E0 = 500             # Амплитуда электрического поля (В/м)
c  = 3e8             # Скорость света в вакууме (м/с)
ε0 = 8.854e-12       # Диэлектрическая постоянная вакуума (Ф/м)
μ0 = 4pi * 1e-7      # Магнитная постоянная вакуума (Гн/м)
ω  = 2 * pi * f      # Круговая частота (рад/с)

# Волновое сопротивление вакуума
Z0 = sqrt(μ0/ε0)     # ≈ 377 Ом

# Амплитуда магнитного поля из соотношения E/H = Z0
H0 = E0 / Z0         

# Границы зон
z0 = 0.0
z1 = 0.3             # Конец вакуума
z2 = 0.50            # Начало стены
z3 = 1.0             # Конец расчетной области

# Зона 1 (Вакуум)
k1 = ω / c

# Зона 2 (Среда с потерями)
εr2 = 5.0
σ2  = 0.2
εa2 = εr2 * ε0
tg_δ2 = σ2 / (ω * εa2)
δ2 = atan(tg_δ2)              # угол диэлектрических потерь
φ_c2  = δ2 / 2                # Сдвиг фазы H относительно E (отставание)

α2 = ω * sqrt((εa2 * μ0 / 2) * (sqrt(1 + tg_δ2^2) - 1)) 
k2 = ω * sqrt((εa2 * μ0 / 2) * (sqrt(1 + tg_δ2^2) + 1)) 

# Зона 3 (Стена-поглотитель)
εr3 = 8.0
σ3  = 2.0
εa3 = εr3 * ε0
tg_δ3 = σ3 / (ω * εa3)
δ3 = atan(tg_δ3)
φ_c3  = δ3 / 2

α3 = ω * sqrt((εa3 * μ0 / 2) * (sqrt(1 + tg_δ3^2) - 1)) 
k3 = k2  #для простоты

# 2. Вычисление поля H(z) с учётом фазового сдвига 

z_arr = z0:0.001:z3

function calculate_H_and_env(z)
    if z <= z1
        # Зона 1: без потерь, сдвига нет
        env = H0
        phase = k1 * z
        return env * cos(phase), env
        
    elseif z <= z2
        # Зона 2: затухание + МАГНИТНОЕ ПОЛЕ ОТСТАЁТ НА φ_c2
        env = H0 * exp(-α2 * (z - z1))
        phase = k2 * (z - z1) + k1 * z1 - φ_c2
        return env * cos(phase), env
        
    else
        # Зона 3: сильное поглощение + отставание на φ_c3
        env_at_z2 = H0 * exp(-α2 * (z2 - z1))
        env = env_at_z2 * exp(-α3 * (z - z2))
        phase = k3 * (z - z2) + k2 * (z2 - z1) + k1 * z1 - φ_c3
        return env * cos(phase), env
    end
end

results = calculate_H_and_env.(z_arr)
H_arr   = first.(results)
env_arr = last.(results)


# 3. Визуализация графика

p = plot(
    size=(900, 500), 
    dpi=300, 
    grid=true, 
    gridalpha=0.3,
    title="Распространение напряжённости магнитного поля H(z)",
    xlabel="Координата z (м)",
    ylabel="Напряжённость поля H (А/м)",
    margin=5Plots.mm,
    legend=:bottomright
)

# Построение огибающих
plot!(p, z_arr, env_arr,  color=:red, linestyle=:dash, linewidth=1.5, label="Огибающая затухания")
plot!(p, z_arr, -env_arr, color=:red, linestyle=:dash, linewidth=1.5, label="")

# Построение самого поля
plot!(p, z_arr, H_arr, color=:blue, linewidth=1.2, label="H(z)")

# Обозначение границ сред
vline!(p, [z1], color=:black, linestyle=:dot, linewidth=2, label="Начало среды с потерями")
vline!(p, [z2], color=:orange, linestyle=:dot, linewidth=2, label="Начало стены")

# Текстовые аннотации 
max_H = maximum(abs.(H_arr))
y_annot = max_H * 0.9  

annotate!(p, 0.15, y_annot, text("Вакуум", 10, :black))
annotate!(p, 0.40, y_annot, text("Среда с потерями", 10, :black))
annotate!(p, 0.75, y_annot, text("Стена", 10, :black))



# Отображение
display(p)