Мірошниченко_ГЙМ_ПР№14-16

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks

class ModelParams:
    def __init__(self, tau_hg=1.0e-6):
        self.phi_0g_star = 0.4
        self.g_g = 12.0
        self.M_g = 2.5e5
        self.mu_g = 3.0e5
        self.phi_1g_star = 3.0e-6
        self.e_g = 3.6e-4
        self.phi_2g = 5.6e-13
        self.phi_3g = 3.0e-20
        self.phi_0D_star = 5.0e-9
        self.g_D = 2.0e-8
        self.M_D = 0.0
        self.mu_D = 1.65e-4
        self.phi_1D_star = 1.0e-24
        self.e_D = 6.0e-23
        self.phi_gD = 1.0e-16
        self.psi_gD = 1.0e-23
        self.tau_hg = tau_hg

def _phi0_partial(phi0_star, g_m, M_m, eps):
    return phi0_star + g_m * eps + 0.5 * M_m * (eps**2)

def _phi1(phi1_star, e_m, eps):
    return phi1_star + 2.0 * e_m * eps

def calculate_I2_critical(h_g, epsilon_ii_e, N_D, p):
    phi0g_part = _phi0_partial(p.phi_0g_star, p.g_g, p.M_g, epsilon_ii_e)
    phi0D_part = _phi0_partial(p.phi_0D_star, p.g_D, p.M_D, epsilon_ii_e)

    phi1g = _phi1(p.phi_1g_star, p.e_g, epsilon_ii_e)
    phi1D = _phi1(p.phi_1D_star, p.e_D, epsilon_ii_e)
    if abs(phi1D) < 1e-50:
        phi1D = 1e-50

    h = np.asarray(h_g, dtype=float)

    ratio_phi = p.phi_gD / phi1D
    ratio_psi = p.psi_gD / phi1D

    term_h3 = (2.0 * (p.psi_gD**2) / phi1D - p.phi_3g) * (h**3)
    term_h2 = (p.phi_2g - 3.0 * ratio_psi * p.phi_gD) * (h**2)

    term_h1_coeff = (
        (p.phi_gD**2) / phi1D
        - phi1g
        - 2.0 * ratio_psi * phi0D_part
        - 4.0 * (p.psi_gD**2) / (p.tau_hg * (phi1D**2)) * N_D
    )
    term_h1 = term_h1_coeff * h

    term_h0 = (
        phi0g_part
        + ratio_phi * phi0D_part
        + 2.0 * p.psi_gD * p.phi_gD / (p.tau_hg * (phi1D**2)) * N_D
    )

    numerator = -(term_h3 + term_h2 + term_h1 + term_h0)
    denominator = 2.0 * p.mu_D * (ratio_phi - 2.0 * ratio_psi * h) + 2.0 * p.mu_g

    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        I2 = numerator / denominator

    invalid = (~np.isfinite(I2)) | (np.abs(denominator) < 1e-20) | (np.abs(I2) > 1e10)
    if np.ndim(I2) == 0:
        return float(I2) if not bool(invalid) else float("nan")
    I2 = I2.astype(float, copy=False)
    I2[invalid] = np.nan
    return I2

def _pick_extrema(y, mode):
    if y.size < 3:
        return float("nan"), float("nan")

    y_safe = np.where(np.isfinite(y), y, -1e30)
    peaks, _ = find_peaks(y_safe)
    valleys, _ = find_peaks(-y_safe)

    y_max = float("nan")
    y_min = float("nan")

    if peaks.size:
        idx = int(peaks[0]) if mode == "first" else int(peaks[np.argmax(y_safe[peaks])])
        y_max = float(y[idx]) if np.isfinite(y[idx]) else float("nan")

    if valleys.size:
        idx = int(valleys[0]) if mode == "first" else int(valleys[np.argmin(y_safe[valleys])])
        y_min = float(y[idx]) if np.isfinite(y[idx]) else float("nan")

    return y_max, y_min

def compute_diagram_a(eps_range, h_grid, p, extrema_mode):
    upper = np.full_like(eps_range, np.nan, dtype=float)
    lower = np.full_like(eps_range, np.nan, dtype=float)
    hg0 = np.full_like(eps_range, np.nan, dtype=float)

    for i, eps in enumerate(eps_range):
        y = calculate_I2_critical(h_grid, epsilon_ii_e=float(eps), N_D=0.0, p=p)
        y_max, y_min = _pick_extrema(y, mode=extrema_mode)
        upper[i] = y_max
        lower[i] = y_min
        hg0[i] = calculate_I2_critical(0.0, epsilon_ii_e=float(eps), N_D=0.0, p=p)

    return upper, lower, hg0

def compute_diagram_b(nd_range, epsilon_fixed, h_grid, p, extrema_mode):
    upper = np.full_like(nd_range, np.nan, dtype=float)
    lower = np.full_like(nd_range, np.nan, dtype=float)
    hg0 = np.full_like(nd_range, np.nan, dtype=float)

    for i, nd in enumerate(nd_range):
        y = calculate_I2_critical(h_grid, epsilon_ii_e=epsilon_fixed, N_D=float(nd), p=p)
        y_max, y_min = _pick_extrema(y, mode=extrema_mode)
        upper[i] = y_max
        lower[i] = y_min
        hg0[i] = calculate_I2_critical(0.0, epsilon_ii_e=epsilon_fixed, N_D=float(nd), p=p)

    return upper, lower, hg0

def build_figure_2_2(tau_hg, eps_min, eps_max, eps_points, nd_min, nd_max, nd_points, epsilon_fixed, hg_max, hg_points, extrema_mode, y_limit_a, y_limit_b):
    p = ModelParams(tau_hg=tau_hg)

    eps_range = np.linspace(eps_min, eps_max, eps_points)
    nd_range = np.linspace(nd_min, nd_max, nd_points)
    h_grid = np.linspace(0.0, hg_max, hg_points)

    upper_a, lower_a, hg0_a = compute_diagram_a(
        eps_range=eps_range,
        h_grid=h_grid,
        p=p,
        extrema_mode=extrema_mode,
    )
    upper_b, lower_b, hg0_b = compute_diagram_b(
        nd_range=nd_range,
        epsilon_fixed=epsilon_fixed,
        h_grid=h_grid,
        p=p,
        extrema_mode=extrema_mode,
    )

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 5.5), constrained_layout=True)

    x_a = eps_range * 1e3
    ax1.plot(x_a, upper_a, color="tab:blue", linewidth=2.0)
    ax1.plot(x_a, lower_a, color="tab:green", linewidth=2.0)
    ax1.plot(x_a, hg0_a, color="tab:blue", linestyle=":", linewidth=2.0)
    ax1.set_title("а)")
    ax1.set_xlabel(r"$\epsilon_{ii}^e,\ \times 10^{-3}$")
    ax1.set_ylabel(r"$I_2$")
    ax1.set_ylim(-y_limit_a, y_limit_a)
    ax1.grid(False)
    ax1.ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
    ax1.ticklabel_format(style="sci", axis="x", scilimits=(0, 0))

    ax1.text(float(x_a.min() + 0.15 * (x_a.max() - x_a.min())), 0.55 * y_limit_a, "B", fontsize=14)
    ax1.text(float(x_a.min() + 0.20 * (x_a.max() - x_a.min())), 0.15 * y_limit_a, "A", fontsize=14)
    ax1.text(float(x_a.min() + 0.10 * (x_a.max() - x_a.min())), -0.45 * y_limit_a, "A*", fontsize=14)
    ax1.text(float(x_a.min() + 0.25 * (x_a.max() - x_a.min())), -0.75 * y_limit_a, "B*", fontsize=14)

    x_b = nd_range * 1e15
    ax2.plot(x_b, upper_b, color="tab:blue", linewidth=2.0)
    ax2.plot(x_b, lower_b, color="tab:green", linewidth=2.0)
    ax2.plot(x_b, hg0_b, color="tab:blue", linestyle=":", linewidth=2.0)
    ax2.set_title("б)")
    ax2.set_xlabel(r"$N_D,\ \mathrm{Дж}^2\cdot\mathrm{с}\cdot\mathrm{м}^{-2}\ \times 10^{-15}$")
    ax2.set_ylabel(r"$I_2$")
    ax2.set_ylim(-y_limit_b, y_limit_b)
    ax2.grid(False)
    ax2.ticklabel_format(style="sci", axis="y", scilimits=(0, 0))

    ax2.text(float(x_b.min() + 0.28 * (x_b.max() - x_b.min())), 0.55 * y_limit_b, "B", fontsize=14)
    ax2.text(float(x_b.min() + 0.45 * (x_b.max() - x_b.min())), -0.01 * y_limit_b, "A", fontsize=14)
    ax2.text(float(x_b.min() + 0.02 * (x_b.max() - x_b.min())), -0.35 * y_limit_b, "A*", fontsize=14)
    ax2.text(float(x_b.min() + 0.45 * (x_b.max() - x_b.min())), -0.85 * y_limit_b, "B*", fontsize=14)

    plt.show()

def main():
    CONFIG = {
        "tau_hg": 1.0e-6,
        "eps_min": -1.5e-3,
        "eps_max": 1.0e-3,
        "eps_points": 220,
        "nd_min": 0.0,
        "nd_max": 3.0e-15,
        "nd_points": 180,
        "eps_fixed": -1.0e-3,
        "hg_max": 3.0e7,
        "hg_points": 6000,
        "extrema_mode": "best",
        "y_limit_a": 1.5e-5,
        "y_limit_b": 1.0e-5,
    }

    build_figure_2_2(
        tau_hg=CONFIG["tau_hg"],
        eps_min=CONFIG["eps_min"],
        eps_max=CONFIG["eps_max"],
        eps_points=CONFIG["eps_points"],
        nd_min=CONFIG["nd_min"],
        nd_max=CONFIG["nd_max"],
        nd_points=CONFIG["nd_points"],
        epsilon_fixed=CONFIG["eps_fixed"],
        hg_max=CONFIG["hg_max"],
        hg_points=CONFIG["hg_points"],
        extrema_mode=CONFIG["extrema_mode"],
        y_limit_a=CONFIG["y_limit_a"],
        y_limit_b=CONFIG["y_limit_b"],
    )

if __name__ == "__main__":
    main()