# L-39 MCS 471 Mon 21 Nov 2022 : wavediff.jl

# Application of the finite difference method to an example of the wave
# equation, taken from section 8.2 of Numerical Analysis by Timothy Sauer.

using Printf
using LinearAlgebra

Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

"""
    eigwavemat(c::Float64,a::Float64,b::Float64,T::Float64,M::Int,N::Int)

shows the eigenvalues of the matrix in the stability analysis.
"""
function eigwavemat(c::Float64,a::Float64,b::Float64,T::Float64,
                    M::Int,N::Int)
    h = (b-a)/(M+1)
    k = T/(N-1)
    sigma = c*k/h
    Adiag = (2-2*sigma^2)*ones(M)
    Asubd = (sigma^2)*ones(M-1)
    A = diagm(Adiag) + diagm(+1 => Asubd) + diagm(-1 => Asubd)
    A = diagm(Adiag) + diagm(+1 => Asubd) + diagm(-1 => Asubd)
    Id = Matrix{Float64}(I,M,M)
    B = zeros(M,M)
    AA = [A -Id; Id B]
    println("The eigenvalues  of AA :")
    v = eigvals(AA)
    show(stdout, "text/plain", v); println("")
    absv = [abs(v[i]) for i=1:length(v)]
    println("The absolute values :")
    show(stdout, "text/plain", absv); println("")
end

"""
    wavediff(c::Float64,f0::Function,g0::Function,lb::Function,rb::Function,
             a::Float64,b::Float64,T::Float64,M::Int,N::Int)

returns the M-by-N matrix of approximated values over [a,b]
with time running from 0 to T.

The speed of the wave equations is determined by c.

The function f0 is a function in the space variable x
and gives the initial position over [a,b].

The function g0 is a function in the space variable x
and gives the initial velocity over [a,b].

The function lb is a function in the time variable
and defines the values at the left bound of [a,b].

The function rb is a function in the time variable 
and defines the values at the right bound of [a,b].

EXAMPLE :
    f(x) = sin(pi*x)
    g(x) = 0
    lb(t) = 0
    rb(t) = 0
    sol = wavediff(2.0,f,g,lb,rb,0.0,1.0,1.0,20,40)
"""
function wavediff(c::Float64,f0::Function,g0::Function,lb::Function,
                  rb::Function,a::Float64,b::Float64,T::Float64,M::Int,N::Int)
    result = zeros(M,N)
    offset = zeros(M)
    h = (b-a)/(M+1)
    k = T/(N-1)
    sigma = c*k/h
    Adiag = (2-2*sigma^2)*ones(M)
    Asubd = (sigma^2)*ones(M-1)
    A = diagm(Adiag) + diagm(+1 => Asubd) + diagm(-1 => Asubd)
    result[:,1] = Array([f0(i*h) for i=1:M])
    gvals = Array([g0(i*h) for i=1:M])
    offset[1] = 0.5*(sigma^2)*lb(0)
    offset[M] = 0.5*(sigma^2)*rb(0)
    result[:,2] = 0.5*A*result[:,1] + offset + k*gvals
    for j=2:N-1
        offset[1] = (sigma^2)*lb((j-1)*k)
        offset[M] = (sigma^2)*rb((j-1)*k)
        result[:,j+1] = A*result[:,j] - result[:,j-1] + offset
    end
    return result
end

"""
    Runs a particular example of the wave equation.
"""
function main()
    f(x) = sin(pi*x)
    g(x) = 0
    lb(t) = 0
    rb(t) = 0
    sol = wavediff(2.0,f,g,lb,rb,0.0,1.0,1.0,20,25)
    println("The solution :")
    show(stdout, "text/plain", sol); println("")
    eigwavemat(2.0,0.0,1.0,1.0,20,25)
end

main()