checkpoint (still doesn't work)

scripts/string_model_de.py

@@ -101,68 +101,147 @@ import math
 # y(-1, n) = -y(1, n)
 # y(N + 1, n) = -y(N - 1, n)
 
+# simulation parameters
+f1 = 50 # fundamental frequency
+f_e = 44100 # sampling frequency
+N = 50 # number of string segments
+delta_t = 1/f_e # time step
+L = 0.5 # string length, meters
+delta_x = L/N # spatial step
+H = int(f_e * 0.1) # length of simulation in time
+
 # string parameters
-E = 1 # youngs modulus
+E = 200 * 10**9 # youngs modulus, steel = 200GPa
-mu = 1 # linear mass density
+rho = 8000 # density, steel, kg/m^3
-kappa = 1 # radius of gyration
+radius = 0.002 # meters
-L = 1 # string length
+kappa = radius/2 # radius of gyration, r/2 for a circular string
+S = math.pi*radius**2 # string cross sectional area, assuming circular
+mu = S*rho # linear mass density
 M_S = mu*L # string mass
-S = 1 # string cross sectional area
+T = 10000 # string tension, N
-T = 1 # string tension
 c = math.sqrt(T/mu) # transverse wave velocity
-stiffness = 1 # string stiffness parameter
+stiffness = kappa**2 * delta_t**2 / delta_x**4 # string stiffness parameter
 sigma = 1 # decay rate
 tau = 1/sigma # decay time
-omega = 1 # angular frequency
+omega = 1/f1 # angular frequency
 
 # hammer parameters
-M_H = 1 # hammer mass
+M_H = 0.5 # hammer mass, kg
 HSMR = M_H/M_S # hammer-mass string ratio
-V_H_0 = 1 # initial hammer velocity at t=0
+V_H_0 = 10 # initial hammer velocity at t=0, m/s
-x_0 = 1 # distance of hammer from agraffe
+x_0 = L/2 # 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 
+b_3 = 0.001 # some constant 
-K = 1 # hammer stiffness
+K = 10 # 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_1 = (2 - 2*r**2 + b_3/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
 
-
+# string
 x = [0] * N # current string position
+x_initial = [0] * N # assuming string at rest
 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_sample = int(alpha*N) # location where we sample the string position for signal
 
-x_out = np.zeros(H) # taking this as the sound output at some artibraty point along the string
+# hammer
+F_H_current = 0 # current force exterted by hammer
+eta_current = 0 # current hammer position
+eta_last = 0 # last hammer position
+i_H = 16 # hammer width
+x_hammer = int(alpha * N) # location where hammer strikes
+eta_0 = 0 # initial hammer displacement
+
+x_out = np.zeros(H) # taking this as the sound output at some arbitrary point along the string
 t = np.arange(0, H/f_e, delta_t)
 
+def plot_current():
+    plt.plot(np.arange(0, N, 1), x)
+    plt.title("Response")
+    plt.xlabel("t")
+    plt.ylabel("x(t, x_sample)")
+    plt.grid()
+    plt.show()
+
+def spatial_window(i):
+    if(i < x_hammer+i_H/2 and i > x_hammer-i_H/2):
+        return 1
+    else:
+        return 0
+
+# first iteration, n = 1
+last_x1 = x_initial.copy()
+eta_current = V_H_0 * delta_t + eta_0
+
+for i in range(1, N-1):
+    x[i] = (last_x1[i+1] + last_x1[i-1])/2 # irrelavant if string is at rest at hammer strike, but relevant if a repeat hit
+last_x1 = x.copy()
+F_H_current = K * abs(eta_current - x[x_hammer])**p
+
+# second iteration, n = 2
+eta_next = 2*eta_current - eta_0 - (delta_t**2*N*F_H_current)/M_S
+eta_last = eta_current
+eta_current = eta_next
+
+for i in range(1, N-1):
+    x[i] = last_x1[i+1] + last_x1[i-1] - x_initial[i] + ((delta_t**2)*N*F_H_current*spatial_window(i))/M_S
+
+plot_current()
+
+last_x2 = last_x1.copy()
+last_x1 = x.copy()
+F_H_current = K * abs(eta_current - x[x_hammer])**p
+
+# rest of the iterations
 for n in range(H): # time
 
-    for i in range(N): # space
+    if(n <= 2):
-        x_out[n] = math.sin(n/10000)
+        continue
+
+    eta_next = 2*eta_current - eta_last - (delta_t**2*N*F_H_current)/M_S
+    eta_last = eta_current
+    eta_current = eta_next
+
+    push = eta_current - x[x_hammer]
+    if(push > 0):
+        F_H_current = K * abs(push)**p
+    else:
+        F_H_current = 0
+
+    next_x = [0] * N
+
+    for i in range(2, N-2): # space
+
+        term_1 = a_1*x[i] + a_2*last_x1[i]
+        term_2 = a_3*(x[i+1] + x[i-1]) + a_4*(x[i+2] + x[i-2])
+        term_3 = a_5*(last_x1[i+1] + last_x1[i-1] + last_x2[i])
+        term_4 = (delta_t**2 * N*F_H_current * spatial_window(i))/M_S
+
+        next_x[i] = term_1 + term_2 + term_3 + term_4
+
+    last_x2 = last_x1.copy()
+    last_x1 = x.copy()
+    x = next_x.copy()
+
+    if(n == 30):
+        plot_current()
+
+    x_out[n] = x[x_sample]
 
 # plotting
 plt.plot(t, x_out)
-plt.title("Step Response")
+plt.title("Response")
 plt.xlabel("t")
-plt.ylabel("y")
+plt.ylabel("x(t, x_sample)")
 plt.grid()
 plt.show()