N = 32 # number of masses
beta =.3 # cubic string spring
A = 10 # amplitude
modes = 3 # number of modes to plot


function calculate_x_i(mass_num, node_num, num_masses, amplitutude)
    return A *
        sqrt(2 / (num_masses + 1)) * 
        sin((mass_num * node_num * π) / (num_masses + 1))
end

function calculate_x_s_for_mode(mass_num, node_num, num_masses, amplitutude)
    x = 0
    for i in 1:num_masses
        x += calculate_x_i(mass_num, node_num, num_masses, amplitutude)
    end
    return x
end

# calculate the x_s for the first three modes

x_s = zeros(modes, N+2) # add two for zeros on ends
# add the first and last 0s
for i in 1:modes
    for j in 1:N
        x_s[i, j+1] = calculate_x_s_for_mode(j, i, N, A)
    end
end

# plot the first three modes
using Plots
plot(x_s[1, :], label="Mode 1", marker=:circle, xlabel="Mass Number", ylabel="Displacement", title="First Three Modes")
plot!(x_s[2, :], label="Mode 2", marker=:circle)
plot!(x_s[3, :], label="Mode 3", marker=:circle)