% 这个octave/matlab脚本基于基本的小球和轻弹簧模型，演示了机械波的传播
% 小球的驱动力源于两侧弹簧变形的不均匀。
% 可以理解为固体物理中的格波，也可以理解为波动方程的离散形式
% CC-BY

clc
clear

dt = 0.1;

k = 1;
m = 1;
N = 50;
_x = 1:N;

u2 = zeros(N,1); %displacement from equilibrium 
u1 = zeros(N,1);
u0 = zeros(N,1);

figure();
hold on
axis equal

for tick = 0:1:floor(10/dt)
  t = tick*dt;
  for i = 2:N-1
    fr = k*(u1(i+1)-u1(i)); % force between i and i+1 particles.
    fl = k*(u1(i)-u1(i-1)); % force between i-1 and i particles.
    f = fr-fl; %total force on i particle
    u2(i) = 2*u1(i) - u0(i) + (dt)^2/m*f;
  end

  u2(1) = 0.2*sin(2*pi/10*t); %boundary conditions
  u2(N) = 0;

  u0 = u1;
  u1 = u2;

  if mod(tick,5) == 0
    clf;
    hold on
    axis equal
    axis([0 10 -2 2])

    line([_x(1)+u2(1) _x(N)+u2(N)],[0 0])
    for i = 1:N
      h0 = scatter(_x(i), 0); #initial location
      h = scatter(_x(i)+u2(i), 0);
      set(h0, 'MarkerEdgeColor',[0.8 0.8 0.8]);
    end
    drawnow;
    pause(0.01)
  end

end