module GrayScottSolver

include("utils/laplacian.jl")

using Observables

using .Laplacian

"""
    step_fhn!(U, V, params_obs::Observable; dx=1, dt=0.01)

    calculation of new "matrix" for each iteration (step) with GrayScott

    # Arguments:
    `U`: activator matrix
    `V`: inhibitor matrix
    `param_obs`: used parameters from CombinedPDEParams
    `dx`: dx
"""
function step_gray_scott!(U, V, params_obs::Observable; dx=1)
    # Extract parameters
    p = params_obs[]
    Du, Dv = p.Du, p.Dv
    F, k = p.F, p.k

    # Compute Laplacians on the interior
    lap_u = laplacian(U, dx)  # shape (N-2, N-2)
    lap_v = laplacian(V, dx)

    # Extract interior values
    u = U[2:end-1, 2:end-1]
    v = V[2:end-1, 2:end-1]

    # Gray-Scott update equations (Euler-style)
    uvv = u .* v .* v
    u_new = u .+ (Du .* lap_u .- uvv .+ F .* (1 .- u))
    v_new = v .+ (Dv .* lap_v .+ uvv .- (F .+ k) .* v)

    # Update interior
    U[2:end-1, 2:end-1] .= u_new
    V[2:end-1, 2:end-1] .= v_new

    # Apply periodic boundary conditions
    U[1, :] .= U[end-1, :]
    U[end, :] .= U[2, :]
    U[:, 1] .= U[:, end-1]
    U[:, end] .= U[:, 2]

    V[1, :] .= V[end-1, :]
    V[end, :] .= V[2, :]
    V[:, 1] .= V[:, end-1]
    V[:, end] .= V[:, 2]

    return U  # for visualization
end

end # Module end