run.jl

# ===================================================================
# Full Animation Script for Bifurcation Project
# (Version with 9 GIFs, including nested cycles and FHN-C)
# ===================================================================

# --- 0. Package Setup ---
# (Uncomment the following lines if you don't have the packages)
# using Pkg
# Pkg.add("DifferentialEquations")
# Pkg.add("Plots")
# Pkg.add("Printf")
# Pkg.add("Statistics")

# --- 1. Load All Packages ---
using DifferentialEquations
using Plots
using Plots.PlotMeasures # Added for mm padding
using Printf             # For formatting text in the HUD
using Statistics         # For mean() in period calculation

# Set a default plotting backend (optional, but good for consistency)
gr() 

# ===================================================================
# --- 2. ODE Model Definitions ---
# ===================================================================

# --- PART A: Van der Pol ---
# p = [eps, L]
function vdp!(du, u, p, t)
    x, y = u
    eps, L = p
    du[1] = y
    du[2] = -x - y*(eps + L*x^2)
end

# --- PART B: Normal Form ---
# p = [eps, L1]
function normal_form!(du, u, p, t)
    x, y = u
    eps, L1 = p
    r_sq = x^2 + y^2
    du[1] = y - x*(eps + L1*r_sq)
    du[2] = -x - y*(eps + L1*r_sq)
end

# --- PART B (Extended): Normal Form with r^4 term (for GIF 8 & 9) ---
# p = [eps, L1, L2]
function normal_form_r4!(du, u, p, t)
    x, y = u
    eps, L1, L2 = p
    r_sq = x^2 + y^2
    r_4 = r_sq^2
    du[1] = y - x*(eps + L1*r_sq + L2*r_4)
    du[2] = -x - y*(eps + L1*r_sq + L2*r_4)
end


# --- PART C: Fitzhugh-Nagumo (SNIC model for GIF 10) ---
# p = [eps, a]
function fhn_C!(du, u, p, t)
    x, y = u
    eps, a = p
    du[1] = x - x^3 - y
    du[2] = eps*(x - a) # The y'-nullcline is x = a
end

# --- PART D: Fitzhugh-Nagumo (Hopf model for GIF 5, 6, 7) ---
# p = [eps, a]
function fhn_D!(du, u, p, t)
    x, y = u
    eps, a = p
    du[1] = x - x^3 - y
    du[2] = eps*(x - a*y) # The y'-nullcline is y = x/a
end

# FHN-D model for GIF 7 (Subcritical)
# p = [eps, a]
function fhn_D_subcritical!(du, u, p, t)
    x, y = u
    eps, a = p
    du[1] = x - (x^3 / 3) - y  # The only change
    du[2] = eps*(x - a*y)
end

# --- FHN Nullcline Helper Functions ---
fhn_v_nullcline(x) = x - x^3
fhn_v_nullcline_sub(x) = x - (x^3 / 3)
fhn_w_nullcline_D(x, a) = x / a # For Part D
# (Part C nullcline is a vertical line, plotted directly)


# ===================================================================
# --- 2.5 Vector Field Helper Functions ---
# ===================================================================

function get_vector_field_vdp(p, X_grid, Y_grid)
    eps, L = p
    U = Y_grid
    V = -X_grid .- Y_grid .* (eps .+ L .* X_grid.^2)
    return (U, V)
end

function get_vector_field_normal_form(p, X_grid, Y_grid)
    eps, L1 = p
    r_sq = X_grid.^2 .+ Y_grid.^2
    U = Y_grid .- X_grid .* (eps .+ L1 .* r_sq)
    V = -X_grid .- Y_grid .* (eps .+ L1 .* r_sq)
    return (U, V)
end

function get_vector_field_normal_form_r4(p, X_grid, Y_grid)
    eps, L1, L2 = p
    r_sq = X_grid.^2 .+ Y_grid.^2
    r_4 = r_sq.^2
    U = Y_grid .- X_grid .* (eps .+ L1 .* r_sq .+ L2 .* r_4)
    V = -X_grid .- Y_grid .* (eps .+ L1 .* r_sq .+ L2 .* r_4)
    return (U, V)
end

function get_vector_field_fhn_C(p, X_grid, Y_grid)
    eps, a = p
    U = X_grid .- X_grid.^3 .- Y_grid
    V = eps .* (X_grid .- a)
    return (U, V)
end

function get_vector_field_fhn_D(p, X_grid, Y_grid, v_nullcline_func)
    eps, a = p
    U = v_nullcline_func.(X_grid) .- Y_grid
    V = eps .* (X_grid .- a .* Y_grid)
    return (U, V)
end

# ===================================================================
# --- 3. Analysis Helper Functions ---
# ===================================================================

"""
Normalizes a vector field for clean plotting with `quiver`.
Returns (U_norm, V_norm)
"""
function normalize_vectors(U, V, scale_factor)
    M = sqrt.(U.^2 .+ V.^2) .+ 1e-9 # Magnitude, with epsilon
    U_norm = U .* scale_factor ./ M
    V_norm = V .* scale_factor ./ M
    return (U_norm, V_norm)
end

"""
Finds the final "settled" radius of a limit cycle.
We find the max(x) in the last 25% of the simulation.
"""
function calculate_radius(sol)
    # Ensure simulation ran long enough
    if length(sol.t) < 10
        return 0.0
    end
    t_start_analysis = sol.t[1] + 0.75 * (sol.t[end] - sol.t[1])
    indices = findall(t -> t > t_start_analysis, sol.t)
    
    # Check if we have any points in the analysis window
    if isempty(indices)
        # Fallback: just use the last 10% of points
        indices = max(1, length(sol) - 10):length(sol)
    end
    
    return maximum(abs.(sol[1, indices]))
end

"""
Finds the period of the limit cycle.
We find the time between peaks in the last half of the data.
"""
function calculate_period(sol)
    if length(sol.t) < 20
        return 0.0
    end
    
    t_start_analysis = sol.t[1] + 0.5 * (sol.t[end] - sol.t[1])
    indices = findall(t -> t > t_start_analysis, sol.t)
    
    if length(indices) < 3
        return 0.0 # Not enough data
    end
    
    sol_data = sol[1, indices]
    peak_indices = []
    for i in 2:(length(sol_data)-1)
        # Find local maxima
        if sol_data[i] > sol_data[i-1] && sol_data[i] > sol_data[i+1]
            push!(peak_indices, i)
        end
    end
    
    if length(peak_indices) < 2
        return 0.0 # Not enough peaks to find a period
    end
    
    peak_times = sol.t[indices[peak_indices]]
    return mean(diff(peak_times)) # Average period
end

"""
Creates a non-linear range of parameter values that clusters
points around the bifurcation at 0.0.
"""
function create_bifurcation_range(p_min, p_max, n_frames)
    # Handle case where one side is zero
    if p_min == 0.0
        return p_max .* (range(0, 1, length=n_frames).^2)
    elseif p_max == 0.0
        return p_min .* (range(1, 0, length=n_frames).^2)
    end

    n_pos = Int(floor(n_frames * abs(p_max) / (abs(p_min) + abs(p_max))))
    n_neg = n_frames - n_pos
    
    # Use a quadratic sweep to cluster points near 0
    pos_range = p_max .* (range(0, 1, length=n_pos).^2)
    neg_range = p_min .* (range(1, 0, length=n_neg).^2)
    
    # Combine and remove duplicate 0
    full_range = [neg_range; pos_range[2:end]]
    return full_range
end

"""
Creates and returns a new plot object showing
the eigenvalues in the complex plane.
"""
function create_eigenvalue_plot(eps_i)
    # 1. Calculate Eigenvalue Data
    # For our models, λ = -eps ± i
    re_lambda = -eps_i
    im_lambda = [1.0, -1.0]
    
    # 2. Determine Stability Color
    dot_color = (re_lambda < 0) ? :green : :red

    # 3. Define Plot Limits
    # Our Re(λ) sweeps from ~ -0.5 (from eps=0.5) to ~ +1.5 (from eps=-1.5)
    plot_xlims = (-0.6, 1.6) # (min, max)
    plot_ylims = (-1.5, 1.5)
    
    # 4. Create the base plot
    p_eig = plot(
        xlims = plot_xlims,
        ylims = plot_ylims,
        xlabel = "Re(λ)",
        ylabel = "Im(λ)",
        title = "Eigenvalues (λ)",
        legend = false,
        margin = 3mm # Added padding
    )
    
    # 5. Draw the Imaginary Axis (Stability Boundary)
    # FIX: Pass y-coordinates as a Vector, not a Tuple
    plot!(p_eig, [0, 0], [plot_ylims[1], plot_ylims[2]],
        color = :black,
        linestyle = :dash
    )
    
    # 6. Plot the two eigenvalue dots
    scatter!(p_eig,
        [re_lambda, re_lambda], # Real parts
        im_lambda,             # Imaginary parts
        color = dot_color,
        markersize = 5
    )
    
    return p_eig
end

"""
Creates a static plot object to be used as an information panel.
"""
function create_info_panel(equation_text::String, ic_text::String)
    # Create a blank canvas
    p_info = plot(grid=false, showaxis=false, xticks=false, yticks=false, 
                  border_stroke=nothing, ylims=(0,1), xlims=(0,1),
                  margin = 3mm) # Added padding
    
    # Add the left-aligned equation text
    # FIX: Replaced Unicode with ASCII equivalents to prevent font errors
    annotate!(p_info, 0.05, 0.5, 
        text(equation_text, :left, 9, :black, "Courier"))
    
    # Add the right-aligned IC text
    # FIX: Replaced Unicode with ASCII equivalents
    annotate!(p_info, 0.95, 0.5, 
        text(ic_text, :right, 9, :black, "Courier"))
        
    return p_info
end


# ===================================================================
# --- 4. GIF Generator Functions ---
# ===================================================================

"""
GIF 1: Supercritical Bifurcation (L1 = +1)
"""
function generate_gif_1_supercritical(L1_val, filename)
    println("  Generating $filename...")
    prob = ODEProblem(normal_form!, [0.1, 0.1], (0.0, 50.0), [0.5, L1_val])
    prob_outer = ODEProblem(normal_form!, [2.0, 2.0], (0.0, 50.0), [0.5, L1_val])
    
    # Use non-linear range to cluster frames near eps = 0
    eps_range = create_bifurcation_range(-0.5, 0.5, 120)
    
    plot_lims = ((-1.5, 1.5), (-1.5, 1.5))
    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]
    
    # --- Create Static Info Panel ---
    eq_string = "Normal Form (Supercritical):\n" *
                "x' = y - x(eps + L1*r^2)\n" *
                "y' = -x - y(eps + L1*r^2)\n" *
                "L1 = +1.0"
                
    ic_string = "Initial Conditions (u0):\n" *
                "Inner: [0.1, 0.1]\n" *
                "Outer: [2.0, 2.0]"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------

    anim = @animate for eps_i in eps_range
        p_new = [eps_i, L1_val]
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_normal_form(p_new, X_grid, Y_grid)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        sol = solve(remake(prob, p=p_new), Tsit5(), saveat=0.1)
        sol_outer = solve(remake(prob_outer, p=p_new), Tsit5(), saveat=0.1)
        
        hud_text = @sprintf("eps = %.3f", eps_i)
        title_str = "Supercritical Hopf Bifurcation"
        
        plot!(p1, sol, vars=(1,2), label="Inner Trajectory", linewidth=2, color=:blue)
        plot!(p1, sol_outer, vars=(1,2), label="Outer Trajectory", linewidth=2, color=:red)
        scatter!(p1, [0], [0], label="Origin", color=:black, markersize=3)
        plot!(p1, title="$title_str\n$hud_text", legend=:topright)
        
        # --- Plot 2: Eigenvalue Plane ---
        p_eig = create_eigenvalue_plot(eps_i)
        
        # --- Combine ---
        p_main = plot(p1, p_eig, layout = @layout([A B]))
        plot(p_main, p_info, layout=@layout([A{0.85h}; B{0.15h}]), size=(800, 500))
    end
    gif(anim, filename, fps=25)
end

"""
GIF 3: Bifurcation Diagram Tracer
"""
function generate_gif_3()
    println("  Generating gif_3_bifurcation_tracer.gif...")
    prob = ODEProblem(normal_form!, [2.0, 2.0], (0.0, 100.0), [0.5, 1.0]) # p=[eps, L1]
    
    # Use non-linear range
    eps_range = create_bifurcation_range(-1.0, 0.5, 120)
    
    bifurcation_data = Tuple{Float64, Float64}[] # To store (eps, r)
    plot_lims = ((-1.5, 1.5), (-1.5, 1.5))
    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]
    
    # --- Create Static Info Panel ---
    eq_string = "Normal Form (Supercritical):\n" *
                "x' = y - x(eps + L1*r^2)\n" *
                "y' = -x - y(eps + L1*r^2)\n" *
                "L1 = +1.0"
                
    ic_string = "Initial Condition (u0):\n" *
                "[2.0, 2.0]\n" *
                "r = sqrt(-eps / L1)"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------
    
    # We need to pre-generate all data to set plot limits correctly
    all_radii = []
    for eps_i in eps_range
        # FIX: Added saveat=0.1 for robust radius calculation
        sol = solve(remake(prob, p=[eps_i, 1.0]), Tsit5(), saveat=0.1)
        r = (eps_i < 0) ? calculate_radius(sol) : 0.0
        push!(all_radii, r)
        push!(bifurcation_data, (eps_i, r))
    end
    
    max_radius = maximum(all_radii) * 1.1

    anim = @animate for i in 1:length(eps_range)
        eps_i = eps_range[i]
        p_new = [eps_i, 1.0]
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_normal_form(p_new, X_grid, Y_grid)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        sol = solve(remake(prob, p=p_new), Tsit5(), saveat=0.1)
        plot!(p1, sol, vars=(1,2), title="Phase Plane\neps = $(@sprintf("%.3f", eps_i))", label="Trajectory", color=:blue)
        scatter!(p1, [0], [0], label="Origin", color=:black, markersize=3)

        # --- Plot 2: Bifurcation Diagram ---
        # FIX: Corrected xlims from (0.5, -1.0) to (-1.0, 0.5)
        p2 = plot(xlims=(-1.0, 0.5), ylims=(0, max_radius), title="Bifurcation Diagram", xlabel="eps", ylabel="Radius (r)",
                  margin = 3mm) # Added padding
        plot_data_so_far = bifurcation_data[1:i]
        scatter!(p2, [d[1] for d in plot_data_so_far], [d[2] for d in plot_data_so_far], label="", color=:red, markersize=2)

        # --- Plot 3: Eigenvalue Plane ---
        p_eig = create_eigenvalue_plot(eps_i)
        
        # --- Combine ---
        p_main = plot(p1, p2, p_eig, layout = @layout([A B C]))
        plot(p_main, p_info, layout=@layout([A{0.85h}; B{0.15h}]), size=(1200, 500))
    end
    gif(anim, "gif_3_bifurcation_tracer.gif", fps=20)
end

"""
GIF 4: Period Tracker
"""
function generate_gif_4()
    println("  Generating gif_4_period_tracker.gif...")
    t_span = (0.0, 150.0) 
    prob_vdp = ODEProblem(vdp!, [1.0, 1.0], t_span, [0.0, 1.0]) # p=[eps, L]
    prob_nf = ODEProblem(normal_form!, [1.0, 1.0], t_span, [0.0, 1.0]) # p=[eps, L1]
    
    # Use non-linear range
    eps_range = create_bifurcation_range(-1.5, 0.0, 120) # Only negative eps
    
    plot_lims_vdp = ((-3,3), (-4,4))
    plot_lims_nf = ((-3,3), (-3,3))
    
    # Grid for VDP
    # --- Finer Grid ---
    x_range_vdp = range(plot_lims_vdp[1][1], plot_lims_vdp[1][2], length=30)
    y_range_vdp = range(plot_lims_vdp[2][1], plot_lims_vdp[2][2], length=30)
    X_grid_vdp = [x for x in x_range_vdp, y in y_range_vdp]
    Y_grid_vdp = [y for x in x_range_vdp, y in y_range_vdp]
    
    # Grid for NF
    # --- Finer Grid ---
    x_range_nf = range(plot_lims_nf[1][1], plot_lims_nf[1][2], length=30)
    y_range_nf = range(plot_lims_nf[2][1], plot_lims_nf[2][2], length=30)
    X_grid_nf = [x for x in x_range_nf, y in y_range_nf]
    Y_grid_nf = [y for x in x_range_nf, y in y_range_nf]
    
    # --- Create Static Info Panel ---
    eq_string = "Left: Van der Pol (L=1.0)\n" *
                "x' = y\n" *
                "y' = -x - y(eps + L*x^2)\n" *
                "\nRight: Normal Form (L1=1.0)\n" *
                "r' = -r(eps + L1*r^2)"
                
    ic_string = "Natural Period (T0 ~ 2*pi/Im(lambda)): 6.28\n" *
                "\nInitial Condition (u0):\n" *
                "[1.0, 1.0]"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------

    anim = @animate for eps_i in eps_range
        # --- Plot 1: VDP ---
        p_vdp = [eps_i, 1.0]
        (U_vdp, V_vdp) = get_vector_field_vdp(p_vdp, X_grid_vdp, Y_grid_vdp)
        # --- Finer Grid Aesthetics ---
        (U_vdp_norm, V_vdp_norm) = normalize_vectors(U_vdp, V_vdp, 0.15)
        p1 = quiver(X_grid_vdp, Y_grid_vdp, quiver=(U_vdp_norm, V_vdp_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims_vdp[1], ylims=plot_lims_vdp[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        # FIX: Added saveat=0.1 for robust period calculation
        sol_vdp = solve(remake(prob_vdp, p=p_vdp), Tsit5(), saveat=0.1)
        period_vdp = calculate_period(sol_vdp)
        plot!(p1, sol_vdp, vars=(1,2), title="Van der Pol", label="", color=:blue)
        annotate!(p1, 0, 3.5, text(@sprintf("Period: %.2f", period_vdp), 10))

        # --- Plot 2: Normal Form ---
        p_nf = [eps_i, 1.0]
        (U_nf, V_nf) = get_vector_field_normal_form(p_nf, X_grid_nf, Y_grid_nf)
        # --- Finer Grid Aesthetics ---
        (U_nf_norm, V_nf_norm) = normalize_vectors(U_nf, V_nf, 0.15)
        p2 = quiver(X_grid_nf, Y_grid_nf, quiver=(U_nf_norm, V_nf_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims_nf[1], ylims=plot_lims_nf[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        # FIX: Added saveat=0.1 for robust period calculation
        sol_nf = solve(remake(prob_nf, p=p_nf), Tsit5(), saveat=0.1)
        period_nf = calculate_period(sol_nf)
        plot!(p2, sol_nf, vars=(1,2), title="Normal Form", label="", color=:red)
        annotate!(p2, 0, 2.8, text(@sprintf("Period: %.2f", period_nf), 10))

        # --- Plot 3: Eigenvalue Plane ---
        # (Eigenvalues are the same for both)
        p_eig = create_eigenvalue_plot(eps_i)

        # --- Combine ---
        p_main = plot(p1, p2, p_eig, layout = @layout([A B C]), 
                      plot_title=@sprintf("eps = %.3f", eps_i))
        plot(p_main, p_info, layout=@layout([A{0.85h}; B{0.15h}]), size=(1200, 500))
    end
    gif(anim, "gif_4_period_tracker.gif", fps=20)
end

"""
GIF 5 & 7: FHN `a` Sweep (Activation & Super/Sub-critical)
Now with Oscilloscope view!
"""
function generate_gif_5_7(model, v_nullcline_func, filename)
    println("  Generating $filename...")
    prob = ODEProblem(model, [-0.5, -0.5], (0.0, 150.0), [0.1, 2.0]) # p = [eps, a]
    
    # WIDENED RANGE: More change
    a_range = range(2.0, 0.2, length=120)
    
    x_null_range = range(-2, 2, length=100)
    plot_lims = ((-2, 2), (-1, 1)) # (xlims, ylims)
    
    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]
    
    # --- Create Static Info Panel ---
    is_subcritical = occursin("subcritical", filename)
    eq_x = is_subcritical ? "x' = x - x^3/3 - y\n" : "x' = x - x^3 - y\n"
    title_str = is_subcritical ? "FHN-D (Subcritical)" : "FHN-D (Supercritical)"
    
    eq_string = "$title_str (Part D):\n" *
                eq_x *
                "y' = eps(x - ay)\n" *
                "eps = 0.1"
                
    ic_string = "Initial Condition (u0):\n" *
                "[-0.5, -0.5]"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------

    anim = @animate for a_i in a_range
        p_new = [0.1, a_i]
        sol = solve(remake(prob, p=p_new), Tsit5(), saveat=0.1)
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_fhn_D(p_new, X_grid, Y_grid, v_nullcline_func)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        plot!(p1, sol, vars=(1,2), label="Trajectory", color=:blue, linewidth=2)
        
        # Plot nullclines
        plot!(p1, x_null_range, x -> v_nullcline_func(x), label="v-nullcline (x' = 0)", color=:red, linestyle=:dash)
        plot!(p1, x_null_range, x -> fhn_w_nullcline_D(x, a_i), label="w-nullcline (y' = 0)", color=:green, linestyle=:dash)
        
        plot!(p1, title="Phase Plane (a = $(@sprintf("%.2f", a_i)))", legend=:bottomright)
        
        # --- Plot 2: Oscilloscope ---
        t_span = (sol.t[1], sol.t[end])
        p2 = plot(sol.t, sol[1, :], label="Voltage (x)", xlabel="Time (t)", ylabel="Voltage (x)",
                  xlims=t_span, ylims=plot_lims[1], legend=false, color=:blue,
                  margin = 3mm) # Added padding
        
        # --- Combine ---
        plot(p1, p2, p_info, layout=@layout([A{0.50h}; B{0.35h}; C{0.15h}]), size=(700, 900))
    end
    gif(anim, filename, fps=20)
end

"""
GIF 6: FHN Relaxation Oscillator
Now with Oscilloscope view!
"""
function generate_gif_6()
    println("  Generating gif_6_relaxation_oscillator.gif...")
    prob = ODEProblem(fhn_D!, [-0.5, -0.5], (0.0, 200.0), [1.0, 0.5]) # p = [eps, a]
    
    # WIDENED RANGE: More stiffness
    eps_range = exp.(range(log(1.0), log(0.005), length=120)) # Log sweep
    
    x_null_range = range(-2, 2, length=100)
    plot_lims = ((-2, 2), (-1, 1))

    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]
    
    # --- Create Static Info Panel ---
    eq_string = "FHN-D (Relaxation Oscillator):\n" *
                "x' = x - x^3 - y\n" *
                "y' = eps(x - ay)\n" *
                "a = 0.5"
                
    ic_string = "Initial Condition (u0):\n" *
                "[-0.5, -0.5]"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------

    anim = @animate for eps_i in eps_range
        p_new = [eps_i, 0.5]
        
        # Robustness: Check for stiffness.
        alg = (eps_i < 0.1) ? Rosenbrock23() : Tsit5()
        
        # FIX: Added saveat=0.1 to get a smooth plot even with stiff solvers
        sol = solve(remake(prob, p=p_new), alg, saveat=0.1)
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_fhn_D(p_new, X_grid, Y_grid, fhn_v_nullcline)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        plot!(p1, sol, vars=(1,2), label="Trajectory", color=:blue, linewidth=2)
        
        # Static nullclines
        plot!(p1, x_null_range, fhn_v_nullcline, label="v-nullcline (x' = 0)", color=:red, linestyle=:dash)
        plot!(p1, x_null_range, x -> fhn_w_nullcline_D(x, 0.5), label="w-nullcline (y' = 0)", color=:green, linestyle=:dash)
        
        plot!(p1, title="Phase Plane (eps = $(@sprintf("%.3f", eps_i)))", legend=:bottomright)
        
        # --- Plot 2: Oscilloscope ---
        t_span = (sol.t[1], sol.t[end])
        p2 = plot(sol.t, sol[1, :], label="Voltage (x)", xlabel="Time (t)", ylabel="Voltage (x)",
                  xlims=t_span, ylims=plot_lims[1], legend=false, color=:blue,
                  margin = 3mm) # Added padding
        
        # --- Combine ---
        plot(p1, p2, p_info, layout=@layout([A{0.50h}; B{0.35h}; C{0.15h}]), size=(700, 900))
    end
    gif(anim, "gif_6_relaxation_oscillator.gif", fps=20)
end

# ===================================================================
# --- 4.5 NEW GIF Generator Functions (8, 9, 10) ---
# ===================================================================

"""
GIF 8: True Subcritical & Nested Cycles (Replaces old GIF 2)
Answers Part B: "Examine... with (x^2+y^2)^2 terms..."
"""
function generate_gif_8_nested_cycles()
    filename = "gif_8_nested_cycles.gif"
    println("  Generating $filename...")
    
    # p = [eps, L1, L2]
    p_initial = [0.3, -1.0, 1.0] # eps=0.3, L1=-1, L2=1
    
    # Trajectories
    prob_inner = ODEProblem(normal_form_r4!, [0.1, 0.1], (0.0, 100.0), p_initial)
    prob_middle = ODEProblem(normal_form_r4!, [0.8, 0.0], (0.0, 100.0), p_initial)
    prob_outer = ODEProblem(normal_form_r4!, [2.0, 2.0], (0.0, 100.0), p_initial)
    
    # This range (0.3 to -0.3) is specifically chosen to show all events:
    # eps > 0.25 (Stable origin)
    # 0 < eps < 0.25 (Stable origin, unstable cycle, stable cycle)
    # eps < 0 (Unstable origin, stable cycle)
    eps_range = range(0.3, -0.3, length=120)
    
    plot_lims = ((-2.5, 2.5), (-2.5, 2.5))
    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]
    
    # --- Create Static Info Panel ---
    eq_string = "Normal Form (Subcritical + r^4):\n" *
                "r' = -r(eps + L1*r^2 + L2*r^4)\n" *
                "L1 = -1.0, L2 = +1.0"
                
    ic_string = "Initial Conditions (u0):\n" *
                "Inner: [0.1, 0.1]\n" *
                "Middle: [0.8, 0.0]\n" *
                "Outer: [2.0, 2.0]"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------

    anim = @animate for eps_i in eps_range
        p_new = [eps_i, -1.0, 1.0]
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_normal_form_r4(p_new, X_grid, Y_grid)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        # --- Trajectories ---
        sol_inner = solve(remake(prob_inner, p=p_new), Tsit5(), saveat=0.1)
        sol_middle = solve(remake(prob_middle, p=p_new), Tsit5(), saveat=0.1)
        sol_outer = solve(remake(prob_outer, p=p_new), Tsit5(), saveat=0.1)
        
        plot!(p1, sol_inner, vars=(1,2), label="Inner (u0 = 0.1)", linewidth=2, color=:blue)
        plot!(p1, sol_middle, vars=(1,2), label="Middle (u0 = 0.8)", linewidth=2, color=:orange)
        plot!(p1, sol_outer, vars=(1,2), label="Outer (u0 = 2.0)", linewidth=2, color=:red)
        scatter!(p1, [0], [0], label="Origin", color=:black, markersize=3)
        
        hud_text = @sprintf("eps = %.3f", eps_i)
        
        plot!(p1, title="Subcritical Bifurcation with r^4 term\n$hud_text", legend=:topright)
        
        # --- Plot 2: Eigenvalue Plane ---
        p_eig = create_eigenvalue_plot(eps_i)

        # --- Combine ---
        p_main = plot(p1, p_eig, layout = @layout([A B]))
        plot(p_main, p_info, layout=@layout([A{0.85h}; B{0.15h}]), size=(800, 500))
    end
    gif(anim, filename, fps=25)
end

"""
GIF 9: The Hysteresis Loop
Demonstrates memory in the subcritical (r^4) system.
"""
function generate_gif_9_hysteresis()
    filename = "gif_9_hysteresis.gif"
    println("  Generating $filename...")
    
    # p = [eps, L1, L2]
    p_initial = [-0.5, -1.0, 1.0]
    prob = ODEProblem(normal_form_r4!, [2.0, 0.0], (0.0, 100.0), p_initial)
    
    # Create a sweep up and back down
    # This range is tuned to show the jumps at eps=0 and eps=0.25
    eps_down = range(0.3, -0.3, length=100)
    eps_up = range(-0.3, 0.3, length=100)
    eps_range = [eps_down; eps_up]
    
    bifurcation_data = Tuple{Float64, Float64}[]
    
    # --- Create Static Info Panel ---
    eq_string = "Normal Form (Subcritical + r^4):\n" *
                "r' = -r(eps + L1*r^2 + L2*r^4)\n" *
                "L1 = -1.0, L2 = +1.0"
                
    ic_string = "Initial Condition (u0):\n" *
                "Starts at [2.0, 0.0] and follows\n" *
                "the stable attractor."

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------
    
    # Pre-generate data for plot limits
    all_radii = []
    u0 = [2.0, 0.0] # Start on the outer cycle
    for eps_i in eps_range
        sol = solve(remake(prob, p=[eps_i, -1.0, 1.0], u0=u0), Tsit5(), saveat=0.1)
        r = calculate_radius(sol)
        push!(all_radii, r)
        push!(bifurcation_data, (eps_i, r))
        u0 = sol.u[end] # Use end of last sim as start of next
    end
    max_radius = maximum(all_radii) * 1.1
    
    plot_lims = ((-2.5, 2.5), (-2.5, 2.5))
    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]

    anim = @animate for i in 1:length(eps_range)
        eps_i = eps_range[i]
        p_new = [eps_i, -1.0, 1.0]
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_normal_form_r4(p_new, X_grid, Y_grid)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        # Get the radius for this frame
        r_current = bifurcation_data[i][2]
        # Draw a circle representing the cycle
        theta = 0:0.1:(2π+0.1)
        plot!(p1, r_current .* cos.(theta), r_current .* sin.(theta), label="Limit Cycle", color=:blue, linewidth=3)
        scatter!(p1, [0], [0], label="Origin", color=:black, markersize=3)
        
        plot!(p1, title="Phase Plane\neps = $(@sprintf("%.3f", eps_i))")

        # --- Plot 2: Bifurcation Diagram ---
        p2 = plot(xlims=(-0.3, 0.3), ylims=(0, max_radius), title="Hysteresis Loop", xlabel="eps", ylabel="Radius (r)",
                  margin = 3mm) # Added padding
        plot_data_so_far = bifurcation_data[1:i]
        scatter!(p2, [d[1] for d in plot_data_so_far], [d[2] for d in plot_data_so_far], label="", color=:red, markersize=2)
        
        # Add text for sweep direction
        direction = (i <= 100) ? "Sweeping Down" : "Sweeping Up"
        annotate!(p2, 0, max_radius * 0.9, text(direction, 10))

        # --- Plot 3: Eigenvalue Plane ---
        p_eig = create_eigenvalue_plot(eps_i)
        
        # --- Combine ---
        p_main = plot(p1, p2, p_eig, layout = @layout([A B C]))
        plot(p_main, p_info, layout=@layout([A{0.85h}; B{0.15h}]), size=(1200, 500))
    end
    gif(anim, filename, fps=25)
end

"""
GIF 10: FHN Part C (SNIC Bifurcation)
Answers Part C: "y' = eps(x-a)"
"""
function generate_gif_10_fhn_c()
    filename = "gif_10_fhn_c_snic.gif"
    println("  Generating $filename...")
    
    # p = [eps, a]
    prob = ODEProblem(fhn_C!, [-0.5, -0.5], (0.0, 150.0), [0.1, 1.5]) 
    
    # This range is tuned to show the bifurcation at a_crit ≈ 0.577
    a_range = range(0.7, 0.5, length=120) 
    
    x_null_range = range(-2, 2, length=100)
    plot_lims = ((-2, 2), (-1, 1)) # (xlims, ylims)
    
    # --- Finer Grid ---
    x_range = range(plot_lims[1][1], plot_lims[1][2], length=30)
    y_range = range(plot_lims[2][1], plot_lims[2][2], length=30)
    X_grid = [x for x in x_range, y in y_range]
    Y_grid = [y for x in x_range, y in y_range]
    
    # --- Create Static Info Panel ---
    eq_string = "FHN - Part C:\n" *
                "x' = x - x^3 - y\n" *
                "y' = eps(x - a)\n" *
                "eps = 0.1"
                
    ic_string = "Initial Condition (u0):\n" *
                "[-0.5, -0.5]"

    p_info = create_info_panel(eq_string, ic_string)
    # -------------------------------

    anim = @animate for a_i in a_range
        p_new = [0.1, a_i]
        sol = solve(remake(prob, p=p_new), Tsit5(), saveat=0.1)
        
        # --- Plot 1: Phase Plane ---
        (U, V) = get_vector_field_fhn_C(p_new, X_grid, Y_grid)
        # --- Finer Grid Aesthetics ---
        (U_norm, V_norm) = normalize_vectors(U, V, 0.15)
        p1 = quiver(X_grid, Y_grid, quiver=(U_norm, V_norm),
                   color=:grey, alpha=0.25, arrowhead=0.2, label="", 
                   xlims=plot_lims[1], ylims=plot_lims[2], xlabel="x", ylabel="y",
                   margin = 3mm) # Added padding

        plot!(p1, sol, vars=(1,2), label="Trajectory", color=:blue, linewidth=2)
        
        # Plot nullclines
        plot!(p1, x_null_range, fhn_v_nullcline, label="v-nullcline (x' = 0)", color=:red, linestyle=:dash)
        # Plot the vertical w-nullcline
        plot!(p1, [a_i, a_i], [plot_lims[2][1], plot_lims[2][2]], label="w-nullcline (y' = 0)", color=:green, linestyle=:dash, linewidth=2)
        
        plot!(p1, title="Phase Plane (a = $(@sprintf("%.2f", a_i)))", legend=:bottomright)
        
        # --- Plot 2: Oscilloscope ---
        t_span = (sol.t[1], sol.t[end])
        p2 = plot(sol.t, sol[1, :], label="Voltage (x)", xlabel="Time (t)", ylabel="Voltage (x)",
                  xlims=t_span, ylims=plot_lims[1], legend=false, color=:blue,
                  margin = 3mm) # Added padding
        
        # --- Combine ---
        plot(p1, p2, p_info, layout=@layout([A{0.50h}; B{0.35h}; C{0.15h}]), size=(700, 900))
    end
    gif(anim, filename, fps=20)
end


# ===================================================================
# --- 5. Main Execution ---
# ===================================================================

function main()
    println("Starting animation generation... (this will take several minutes)")
    
    println("\n[1/9] Generating GIF 1 (Supercritical)...")
    generate_gif_1_supercritical(1.0, "gif_1_supercritical.gif")
    
    # Note: GIF 2 (old subcritical) is now replaced by GIF 8.
    
    println("\n[2/9] Generating GIF 3 (Bifurcation Tracer)...")
    generate_gif_3()

    println("\n[3/9] Generating GIF 4 (Period Tracker)...")
    generate_gif_4()
    
    println("\n[4.a/9] Generating GIF 5 (FHN-D Activation, Supercritical)...")
    generate_gif_5_7(fhn_D!, fhn_v_nullcline, "gif_5_fhn_D_activation.gif")

    println("\n[5/9] Generating GIF 6 (Relaxation Oscillator)...")
    generate_gif_6() # FIX: Corrected typo from generate_fig_6

    println("\n[6/9] Generating GIF 7 (FHN-D Activation, Subcritical)...")
    generate_gif_5_7(fhn_D_subcritical!, fhn_v_nullcline_sub, "gif_7_fhn_D_subcritical.gif")
    
    println("\n[7/9] Generating GIF 8 (True Subcritical & Nested Cycles)...")
    generate_gif_8_nested_cycles()
    
    println("\n[8/9] Generating GIF 9 (Hysteresis Loop)...")
    generate_gif_9_hysteresis()
    
    println("\n[9/9] Generating GIF 10 (FHN-C SNIC Bifurcation)...")
    generate_gif_10_fhn_c()
    
    println("\n✅ All 9 animations complete. Check your script directory for the .gif files.")
end

# --- Run the main function ---
# FIX: Wrap in `invokelatest` to prevent World Age errors
Base.invokelatest(main)