% membraneWave.m
% Simulates a wave in a membrane
% JMS, Nov.2019

% Set parameters
clear all
nx = 200;     % number of beads along x
ny = 100;     % number of beads along y
k = 1;        % spring constant (N/m)
m1 = 1;       % beads mass in section 1 (kg)
m2 = 3;       % beads mass in section 2 (kg)
x12 = 100;    % frontier between sections 1 and 2 (m)
dx = 1;       % distance between beads (m)
T = 1;        % string tension (N)
ww = 15;      % wavepacket width (m)
lambda = 5;   % initial wavelength (m)
x0 = 50;      % initial coordinates of wave packet center (m)
y0 = 30;
k0x = 1;      % relative components of initial wavevector
k0y = 0.3;
dt = 0.1;     % time step (s)
tmax = 150;   % final time (s)

% Set initial displacements
c = sqrt(T*dx/m1);                         % wave velocity in section 1
k0 = sqrt(k0x^2+k0y^2);
cx = c*k0x/k0;                             % wave packet velocity
cy = c*k0y/k0;
k0x = (2*pi/lambda)*k0x/k0;                % initial wave vector, with correct norm
k0y = (2*pi/lambda)*k0y/k0;
w = c*norm(k);                             % wave frequency
dy = dx;
x = (1:nx)*dx;                             % bead positions along x
y = (1:ny)*dy;                             % bead positions along y
[x,y] = ndgrid(x,y);                       % all-bead coordinates
mass = m1*(x<=x12) + m2*(x>x12);           % individual bead masses
t0 = 0;                                    % initial time
wa = exp(-0.5*((x-x0-cx*t0).^2+ ...        % wave packet amplitude
               (y-y0-cy*t0).^2)/ww^2);
dwadt = ( cx*(x-x0-cx*t0) + ...            % derivative Dwa/Dt
          cy*(y-y0-cy*t0) )/ww^2.*wa;
z0 =    wa.*cos(k0x*x+k0y*y-w*t0);         % displacements at t0
v0 = dwadt.*cos(k0x*x+k0y*y-w*t0) + ...    % velocities at t0
      w*wa.*sin(k0x*x+k0y*y-w*t0);
z0 = z0 - sum(mass.*z0)/sum(mass);         % subtract center of mass displacement
v0 = v0 - sum(mass.*v0)/sum(mass);         % subtract center of mass velocity

% Integrate Newton's equation
myforce = @(z)membraneForces(z,nx,ny,T,dx);
[z,v] = verlet( myforce, mass(:)', z0(:)', v0(:)', dt, tmax );

% Reshape trajectory arrays
t = (0:dt:tmax);
nt = numel(t);
z = reshape(z,nx,ny,nt);
v = reshape(v,nx,ny,nt);

% Plot trajectory
figure(1)
xmax = (nx+1)*dx;
ymax = (ny+1)*dy;
for it = 1:5:nt
  surf(x',y',z(:,:,it)')
  grid on
  axis([0,xmax,0,ymax,-1.2,+1.2])
  daspect([xmax,xmax,2.4])
  view(30,60)
  pause(0.1)
end
