lab1 / lab1

import matplotlib.pyplot as plt  # библиотека графиков
import scipy.constants as const  # библиотека констант
import numpy as np  # математическая библиотека
import scipy.special as sp  # библиотека специальных функций

# исходные данные
eps = 5.7
Rw = 0.01  # внешний радиус
Rc = 0.004  # внутренний радиус
c = const.c  # скорость света
nu = 0  # значение ню

# данные диапазона частот
start_omega = 5e9
end_omega = 1.15e11
step_omega = 1e9

# данные диапазона beta
step_beta = 2 * 10 ** -4
start_beta = 1 / np.sqrt(eps) + step_beta
end_beta = 1.5 * start_beta


# метод бисекции
def bisect(F, k, left, right, eps):
    E = eps * 2;
    FLeft = F(k, left)
    FRight = F(k, right)
    X = (left + right) / 2

    if (FLeft * FRight > 0):
        print("Ошибка интервала")

    if (FLeft == 0):
        return left
    if (FRight == 0):
        return right

    while (right - left) >= E:
        X = (left + right) / 2
        Y = F(k, X)
        if (Y == 0):
            return X
        if (Y * FLeft < 0):
            right = X
        else:
            left = X
            FLeft = Y
    return X


# вспомогательные функции для решения дисперсионного уравнения
def calc_kz(omega):
    return 2 * np.pi * omega / c



def calc_hi1(beta):
    return (eps - 1 / beta ** 2) ** 0.5


def calc_hi0(beta):
    return ((1 - beta ** 2) / beta ** 2) ** 0.5


def InD1(z):
    return (sp.iv(nu, z) + sp.iv(nu + 2, z)) / 2


def JnD1(z):
    return sp.jv(nu, z) - sp.jv(nu + 2, z) / 2


def KnD1(z):
    return - (sp.kn(nu, z) + sp.kn(nu + 2, z)) / 2


def YnD1(z):
    return (sp.yv(nu, z) - sp.yv(nu + 2, z)) / 2


def ez1D(k, beta):
    hi1 = calc_hi1(beta)
    X = JnD1(hi1 * k * Rc) * sp.yv(nu, hi1 * k * Rw)
    Y = - YnD1(hi1 * k * Rc) * sp.jv(nu, hi1 * k * Rw)
    return X + Y


def ez1(k, beta):
    hi1 = calc_hi1(beta)
    X = sp.jv(nu, hi1 * k * Rc) * sp.yv(nu, hi1 * k * Rw)
    Y = - sp.yv(nu, hi1 * k * Rc) * sp.jv(nu, hi1 * k * Rw)
    return X + Y


# фенкция дисперсионного уравнения
def Disp(k, beta):
    hi1 = calc_hi1(beta)
    hi0 = calc_hi0(beta)
    X = 1 / hi0 * InD1(hi0 * k * Rc) * ez1(k, beta)
    Y = eps / hi1 * ez1D(k, beta) * sp.iv(nu, hi0 * k * Rc)
    return X + Y


data = {}

range_omega = np.arange(step_omega, end_omega, step_omega)  # диапазон частот
range_beta = np.arange(start_beta, end_beta, step_beta)  # дипапазон beta
indexs = range(1, len(range_beta))  # дипазон индексов
points = [[] for i in range(5)]  # линии 5 графиков


# функция поиска корней дисперсионного уравнения
def search_root(omega_range):
    for i, omega in enumerate(omega_range):

        temp = Disp(calc_kz(omega), range_beta[0])  # начальное приближение диспрс. уравнения
        count = 0  # счётчик точек

        for index in indexs:
            next_temp = Disp(calc_kz(omega), range_beta[index])
            if count >= 5:
                break
            if temp * next_temp < 0:
                x = bisect(Disp, calc_kz(omega), range_beta[index - 1], range_beta[index], 10 ** -6)
                points[count].append([omega, x])
                count += 1
            temp = next_temp


search_root(range_omega)

for data in points:
    array_omegas = [d[0] for d in data]
    array_velocities = [d[1] for d in data]
    array_group_velocity = [(array_omegas[i + 1] - array_omegas[i]) / (array_omegas[i + 1] / array_velocities[i + 1]
                            - array_omegas[i] / array_velocities[i]) for i in range(len(array_velocities) - 1)]

    plt.plot(array_omegas, array_velocities, color="red")
    plt.plot(array_omegas[:-1], array_group_velocity, color="blue")

plt.xlabel("w")
plt.ylabel("V/c")
plt.show()