---

kill(all)$
numer : true$

/* --- БЛОК 1: МЕТОД БИСЕКЦИИ --- */
sign_change(expr, int) := block(
    [var, a, b],
    [var, a, b]: int,
    is(subst(a, var, expr) * subst(b, var, expr) < 0)
)$

bisect(expr, int, eps, imax) := block(
    [var, beg, end, xold, xnew, cnv, fxnew, i],
    if not(sign_change(expr, int)) then return(false),
    [var, beg, end]: int, 
    i: 0, [xold, xnew]: [beg, end],
    cnv: abs(xnew - xold),
    while cnv > eps do (
        i: i + 1, xold: xnew, xnew: (beg + end) / 2.0,
        cnv: abs(xnew - xold),
        fxnew: subst(xnew, var, expr),
        if fxnew = 0.0 then return([xnew, i, cnv]),
        if i > imax then return([xnew, i, cnv]),
        if sign_change(expr, [var, beg, xnew]) then end: xnew else beg: xnew
    ),
    [xnew, i, cnv]
)$

sign_change_ints(expr, int, dvar) := block(
    [var, vmin, vmax, a, b, out],
    [var, vmin, vmax]: int,
    [a, b]: [vmin, vmin + dvar], out: [],
    while b < vmax do (
        if sign_change(expr, [var, a, b]) then out: cons([a, b], out),
        a: b, b: b + dvar
    ),
    if sign_change(expr, [var, a, vmax]) then out: cons([a, vmax], out),
    reverse(out)
)$

roots_bisect(expr, int, dvar, eps, imax) := block(
    [len, ints, res, tmp],
    ints: sign_change_ints(expr, int, dvar),
    len: length(ints),
    if len = 0 then return([]),
    res: [],
    for i:1 thru len do (
        tmp: bisect(expr, cons(first(int), ints[i]), eps, imax),
        if tmp # false then res: endcons(first(tmp), res)
    ),
    res
)$

/* --- БЛОК 2: ДИСПЕРСИОННОЕ УРАВНЕНИЕ (Вар. 6) --- */
F(G, V, prms) := block(
    [m, ep1, mu1, ep2, mu2, B1, B2, c1_sq, c2_sq, ry1_sq, ry2_sq, k1, k2, sn1, cs1, sn2, cs2, res],
    [m, ep1, mu1, ep2, mu2, B1, B2] : prms,
    
    c1_sq : 1 - (m/V)^2,
    ry1_sq : %pi^2 * c1_sq * (1 - G),
    c2_sq : 1 - (ep1*mu1/(ep2*mu2)) * (m/V)^2,
    ry2_sq : %pi^2 * (ep2*mu2/(ep1*mu1)) * c2_sq - %pi^2 * c1_sq * G,
    
    k1 : sqrt(abs(ry1_sq)), k2 : sqrt(abs(ry2_sq)),
    
    if ry1_sq >= 0 then (sn1: sin(B1*V*k1), cs1: cos(B1*V*k1))
    else (sn1: sinh(B1*V*k1), cs1: cosh(B1*V*k1)),
    
    if ry2_sq >= 0 then (sn2: sin(B2*V*k2), cs2: cos(B2*V*k2))
    else (sn2: sinh(B2*V*k2), cs2: cosh(B2*V*k2)),
    
    res : k1 * mu2 * cs1 * sn2 + k2 * mu1 * cs2 * sn1,
    return(res)
)$

/* --- БЛОК 3: СБОР ТОЧЕК ДХ --- */
eval_dsp(prms, NG, GMap, Vmin, Vmax, dV, eps, imax) := block(
    [V, %V, %%G, %G, g_val, n, num_modes],
    num_modes : 4, 
    %%G : makelist([], n, 1, num_modes),
    
    %G : float(makelist(i / NG, i, 0, NG - 1)),
    %G : map(lambda([val], subst(val, 'x, GMap)), %G),
    
    for i:1 thru length(%G) do (
        g_val : %G[i],
        %V : roots_bisect(F(g_val, V, prms), [V, Vmin, Vmax], dV, eps, imax),
        if length(%V) > 0 then (
            %V : sort(%V), 
            for n:1 thru min(length(%V), num_modes) do (
                %%G[n] : endcons([%V[n], g_val], %%G[n])
            )
        )
    ),
    %%G
)$

/* --- БЛОК 4: РАСЧЕТ И ВЫВОД РЕЗУЛЬТАТОВ --- */
prms : [0, 15.0, 1.0, 9.0, 1.0, 1/3, 2/3]$

%%G0 : block([%p],
    %p: cons(0, rest(prms)),
    eval_dsp(%p, 40, 1 - (x - 1)^2, 0.01, 14.0, 0.15, 1E-6, 100)
)$

%%G1 : block([%p],
    %p: cons(1, rest(prms)),
    eval_dsp(%p, 40, 1 - (x - 1)^2, 1.01, 14.0, 0.15, 1E-6, 100)
)$

/* Печать точных значений частот отсечки */
print("--------------------------------------------------")$
print("ЗНАЧЕНИЯ НОРМИРОВАННЫХ ЧАСТОТ ОТСЕЧКИ (при G -> 0):")$
for idx:1 thru 4 do (
    if length(%%G0[idx]) > 0 then print("Мода LE0", idx, ": V =", %%G0[idx][1][1]),
    if length(%%G1[idx]) > 0 then print("Мода LE1", idx, ": V =", %%G1[idx][1][1])
                             else print("Мода LE1", idx, ": Вне диапазона расчёта (V > 14)")
)$
print("--------------------------------------------------")$

/* Рисуем Окно 1: Дисперсионные характеристики */
load(draw)$
draw2d(
    terminal = wxt,
    title = "Дисперсионные характеристики LE-волн (Вар. 6)",
    points_joined = true, grid = true,
    xlabel = "V (Нормированная частота)",
    ylabel = "G (Норм. пост. распространения)",
    
    /* m=0 */
    key = "LE01", color = red, points(%%G0[1]),
    key = "LE02", color = orange, points(%%G0[2]),
    key = "LE03", color = dark_yellow, points(%%G0[3]),
    key = "LE04", color = dark_green, points(%%G0[4]),
    
    /* m=1 */
    points_joined = false, point_size = 1,
    key = "LE11", color = red, point_type = 3, points(%%G1[1]),
    key = "LE12", color = orange, point_type = 4, points(%%G1[2]),
    key = "LE13", color = dark_yellow, point_type = 5, points(%%G1[3]),
    key = "LE14", color = dark_green, point_type = 6, points(%%G1[4])
)$

/* --- БЛОК 5: СТРОИТЕЛЬСТВО ЭПЮРЫ ПОЛЯ Ez --- */
ezy(m, Y, G, V, prms) := block(
    [ep1, mu1, ep2, mu2, B1, B2, c1_sq, c2_sq, ry1_sq, ry2_sq, k1, k2, A, B],
    [m, ep1, mu1, ep2, mu2, B1, B2]: prms,
    c1_sq : 1 - (m/V)^2,
    ry1_sq : %pi^2 * c1_sq * (1 - G),
    c2_sq : 1 - (ep1*mu1/(ep2*mu2)) * (m/V)^2,
    ry2_sq : %pi^2 * (ep2*mu2/(ep1*mu1)) * c2_sq - %pi^2 * c1_sq * G,
    
    k1 : sqrt(abs(ry1_sq)), k2 : sqrt(abs(ry2_sq)),
    
    if Y >= 0 then (
        if ry1_sq >= 0 then A : sin(V * k1 * (B1 - Y)) else A : sinh(V * k1 * (B1 - Y)),
        return(A)
    ) else (
        if ry2_sq >= 0 then B : (sin(V * k1 * B1) / sin(V * k2 * B2)) * sin(V * k2 * (B2 + Y))
        else B : (sinh(V * k1 * B1) / sinh(V * k2 * B2)) * sinh(V * k2 * (B2 + Y)),
        return(B)
    )
)$

/* Берем стабильную точку из середины графика моды LE01 */
_pt : %%G0[1][fix(length(%%G0[1])/2)]$
_V : _pt[1]$
_G : _pt[2]$
[ B1, B2 ] : rest(prms, 5)$

/* Рисуем Окно 2: Эпюра поля Ez */
draw2d(
    terminal = wxt,
    title = concat("Распределение поля Ez вдоль слоев (LE01) при V=", string(_V)),
    xlabel = "Y (Ось поперек слоев)",
    ylabel = "Ez (относ. единицы)",
    grid = true,
    color = blue,
    explicit(ezy(0, Y, _G, _V, prms), Y, -B2, B1)
)$