diff --git a/scripts/string_model_de.py b/scripts/string_model_de.py index 2d91cda..4660552 100644 --- a/scripts/string_model_de.py +++ b/scripts/string_model_de.py @@ -100,3 +100,69 @@ import math # eq 23: temporal boundary conditions # y(-1, n) = -y(1, n) # y(N + 1, n) = -y(N - 1, n) + +# string parameters +E = 1 # youngs modulus +mu = 1 # linear mass density +kappa = 1 # radius of gyration +L = 1 # string length +M_S = mu*L # string mass +S = 1 # string cross sectional area +T = 1 # string tension +c = math.sqrt(T/mu) # transverse wave velocity +stiffness = 1 # string stiffness parameter +sigma = 1 # decay rate +tau = 1/sigma # decay time +omega = 1 # angular frequency + +# hammer parameters +M_H = 1 # hammer mass +HSMR = M_H/M_S # hammer-mass string ratio +V_H_0 = 1 # initial hammer velocity at t=0 +x_0 = 1 # distance of hammer from agraffe +alpha = x_0 / L # relative hammer striking position + +# simulation parameters +f1 = 440 # fundamental frequency +f_e = 44100 # sampling frequency +N = 100 # number of string segments +delta_t = 1/f_e # time step +delta_x = L/N # spatial step +H = f_e * 10 # length of simulation in time + +# empirical constants +b_1 = 1 # some constant +b_3 = 1 # some constant +K = 1 # hammer stiffness +p = 1 # stiffness nonlinear exponent + +# derived components +D = 1 + b_1*delta_t + 2*b_3/delta_t +r = c*delta_t/delta_x +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 + + +x = [0] * N # current string position +last_x1 = [0] * N # string position from last timestep +last_x2 = [0] * N # string position from two timesteps ago +x_next = [0] * N # buffer for next string position + +x_out = np.zeros(H) # taking this as the sound output at some artibraty point along the string +t = np.arange(0, H/f_e, delta_t) + +for n in range(H): # time + + for i in range(N): # space + x_out[n] = math.sin(n/10000) + +# plotting +plt.plot(t, x_out) +plt.title("Step Response") +plt.xlabel("t") +plt.ylabel("y") +plt.grid() +plt.show()