% Vibrating Membranes from Advanced Mathematics and Mechanics
% Applications Using MATLAB by Howard B. Wilson, Louis H. Turcotte,
% and David Halpern: http://www.crcpress.com/product/isbn/9781584882626
%
% Functions needed
%      - memcirwv, membanim, besjroot

% Input Data
%   x,y,t -  vectors of position and time values
%   x0,y0 -  coordinates of the point where the force acts
%   alp   -  wave propagation velocity in the membrane
%   w     - frequency of the applied force
%   tmax  - maximum time for computing the solution

% Physical dimensions - defaults from 
%       http://hyperphysics.phy-astr.gsu.edu/hbase/music/cirmem.html
%   T    - Tension 2000 N/m
%   sigma- Density 0.26 kg/m^2
%   a    - Length of membrane 0.6 m
%   b    - Height of membrane 0.6 m

clear
a    = 0.6;             % Length of membrane
b    = 0.6;             % Height of membrane
T    = 2000;            % Tension
sigma = 0.26;           % Membrane density
x0   = a/4;             % Location of concentrated force
y0   = b/4;
w    = 0;               % Frequency of applied force
                        % 0 for one hit, >0 for more
alp  = sqrt(T/sigma);   % Wavespeed
tmax = 5;              % Duration of run

N=40; M=40;             
pan=pi/a*(1:N)'; 
pbm=pi/b*(1:M);
W=alp*sqrt(repmat(pan.^2,1,M)+repmat(pbm.^2,N,1));
    wsort=sort(W(:)); wsort=reshape(wsort(1:42),6,7)';
    disp(' ') 
    disp(['The first forty-two natural ',...
          'frequencies are:'])
    disp(wsort)
   
[u,x,y,t] = memrecwv(a,b,alp,w,x0,y0,tmax); % Solution
membanim(u,x,y,t);                          % Animation