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

# https://people.cs.uchicago.edu/~ridg/stabil/pianostring.pdf
# a differential equations model for a piano string

# eq 1: governing wave equation
# d2y/dt2 = c^2*d2y/dx2 - stiffness*c^2*L^2*d4y/dx4 - 2*b_1*dy/dt + 2*b_3*d3y/dt3 + f(x, x_0, t)
# where y = string's transverse displacement, b_1, b_3 = damping coefficients, f = force density
# c = sqrt(T/mu), transverse wave velocity , T = string tension, mu = string linear mass density)

# eq 2: string stiffness
# stiffness = K^2*(E*S/(T*L^2))
# where K = string's radius of gyration (r/2 for a circular string), E = string's Young's modulus,
# S = cross sectional string area, T = string tension, L = string length

# eq 3: decay
# sigma = 1/tau = b_1 + b_3*omega^2
# where sigma = decay rate, tau = decay time, omega = angular frequency

# eq 4: excitation
# f(x, x_0, t) = f_H(t) * g(x, x_0)
# where f_H(t) = hammer force, g(x, x_0) = hammer dimensional effect

# eq 5: hammer time history
# read the paper if you care, useful only for derivation

# eq 6: power law
# F_H(t) = K*|eta(t) - y(x_0, t)|^p
# where eta(t) is the transverse displacement of the hammer head
# and p = stiffness nonlinear exponent

# eq 7: hammer displacement
# M_H*d2eta/dt2 = -F_H(t)
# where M_H is the mass of the hamemr head

# eq 8: boundary conditions
# y(0, t) = y(L, t) = 0, fixed ends dont move
# d2y/dx2(0, t) = d2y/dx2(L, t) = 0, displacement along the string approaching the ends is continuous

# discrete time implementation

# eq 9: continuous to discrete
# y(x, t) -> y(x_i, t_n) -> y(i, n)
# divide the string into i segments and iterate over n timesteps
# where x_i = delta_x * i and t_n = delta_t * n

# eq 10: recurrence derivation
# y(i, n+1) = a_1*y(i, n) + a_2*y(i, n-1) + a_3*[y(i+1, n) + y(i-1, n)] + a_4*[y(i+2,n) + y(i-2,n)]
#               + a_5*[y(i+1, n-1) + y(i-1, n-1) + y(i, n-2)]
#               + [delta_t^2 * N*F_H(n) * g(i, i_0)]/M_S
# where a_1 through a_5 are defined as follows:
# a_1 = [2 - 2*r^2 + b_1/delta_t - 6*stiffness*N^2*r^2]/D
# a_2 = [-1 + b_1*delta_t + 2*b_3/delta_t]/D
# a_3 = [r^2*(1 + 4*stiffness*N^2)]/D
# a_4 = [b_3/delta_t - stiffness*N^2*r^2]/D
# a_5 = [-b_3/delta_t]/D
# where D = 1 + b_1*delta_t + 2*b_3/delta_t
# and r = c*delta_t/delta_x

# eq 11: stability condition
# N_max = sqrt{[-1 + sqrt(1+16*stiffness*gamma^2)]/(8*stiffness)}
# where
# eq 12: idk what gamma represents
# gamma = f_e/(2*f_1), f_e = sampling frequency and f_1 = fundamental frequency
# or
# eq 13: if neglecting stiffness
# N_max = gamma

# eq 14: rest condition
# y(i, 0) = 0

# eq 15: discrete hammer displacement
# at t = delta_t (n = 1)
# eta(1) = V_H_0 * delta_t

# eq 16: discrete truncated taylor series
# y(i, 1) = [y(i + 1, 0) + y(i - 1, 0)]/2

# eq 17: hammer force exertion
# F_H(1) = K*|eta(1) - y(i_0, 1)|^p

# eq 18: string displacement iteration
# y(i, 2) = y(i-1, 1)] + y(i + 1, 1) - y(i, 0) + [delta_t^2 * N*F_H(1) * g(i, i_0)]/M_S

# eq 19: hammer displacement iteration
# eta(2) = 2*eta(1) - eta(0) - [delta_t^2 * F_H(1)]/M_H

# eq 20: hammer force iteration
# F_H(2) = K*|eta(2) - y(i_0, 2)|^p

# eq 21: hammer rest condition
# eta(n + 1) < y(i_0, n + 1)

# eq 22: spacial boundary conditions
# y(0, n) = y(N, n) = 0

# eq 23: temporal boundary conditions
# y(-1, n) = -y(1, n)
# y(N + 1, n) = -y(N - 1, n)