#!/usr/bin/env python3.10

"""Solve the fourth-order beam-deflection boundary-value problem."""

from pathlib import Path
import shutil

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np


def build_linear_system(interior_points, beam_length, stiffness_parameter, load_parameter):
    """
    Build the finite-difference linear system for the beam problem.

    Returns
    -------
    tuple[numpy.ndarray, numpy.ndarray, float]
        System matrix, right-hand side, and grid spacing.
    """
    step_size = beam_length / (interior_points + 1)
    coefficient = 6.0 + stiffness_parameter * step_size**4
    system_matrix = np.zeros((interior_points, interior_points), dtype=float)
    right_hand_side = -load_parameter * step_size**4 * np.ones(interior_points, dtype=float)

    for row_index in range(interior_points):
        system_matrix[row_index, row_index] = coefficient
        if row_index - 1 >= 0:
            system_matrix[row_index, row_index - 1] = -4.0
        if row_index + 1 < interior_points:
            system_matrix[row_index, row_index + 1] = -4.0
        if row_index - 2 >= 0:
            system_matrix[row_index, row_index - 2] = 1.0
        if row_index + 2 < interior_points:
            system_matrix[row_index, row_index + 2] = 1.0

    system_matrix[0, 0] = 7.0 + stiffness_parameter * step_size**4
    system_matrix[-1, -1] = 5.0 + stiffness_parameter * step_size**4
    return system_matrix, right_hand_side, step_size


def configure_plot_style():
    """Configure a consistent LaTeX-like plotting style."""
    mpl.rcParams.update(
        {
            "text.usetex": shutil.which("latex") is not None,
            "font.family": "serif",
            "font.serif": ["Computer Modern Roman", "DejaVu Serif"],
            "mathtext.fontset": "cm",
            "axes.titlesize": 15,
            "axes.labelsize": 13,
            "xtick.labelsize": 11,
            "ytick.labelsize": 11,
            "legend.fontsize": 11,
        }
    )


def main():
    """Run the beam-deflection simulation from the project."""
    configure_plot_style()

    beam_length = 1.0
    interior_points = 2**9
    load_parameter = 10.0
    stiffness_parameter = 0.5

    system_matrix, right_hand_side, step_size = build_linear_system(
        interior_points, beam_length, stiffness_parameter, load_parameter
    )
    solution_vector = np.linalg.solve(system_matrix, right_hand_side)

    grid_points = np.linspace(0.0, beam_length, interior_points + 2)
    full_solution = np.zeros(interior_points + 2, dtype=float)
    full_solution[1:-1] = solution_vector

    output_directory = Path(__file__).resolve().parent.parent / "Figuras"
    output_directory.mkdir(exist_ok=True)

    plt.figure(figsize=(7.0, 5.0))
    plt.plot(grid_points, full_solution, linewidth=2.0)
    plt.xlabel("x")
    plt.ylabel("u(x)")
    plt.grid(True, linewidth=0.4, alpha=0.5)
    plt.tight_layout()
    plt.savefig(output_directory / "proj_sol.pdf", bbox_inches="tight")
    plt.close()

    print(f"step_size = {step_size:.10f}")
    print(f"minimum_deflection = {np.min(full_solution):.10f}")


if __name__ == "__main__":
    main()