% 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
%   r,th - vectors of radius and polar angle values
%   r0   - radial position of the concentrated force
%   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    - Radius of membrane 0.6 m

clear
a    = 0.6;             % Radius of membrane
T    = 2000;            % Tension
sigma = 0.26;           % Membrane density
r0   = 0.4;             % Location of concentrated force
w    = 0;               % Frequency of applied force
                        % 0 for one hit, >0 for more
th   = linspace(0,2*pi,100); 
r    = linspace(0,a,30); 
alp  = sqrt(T/sigma);   % Wavespeed
tmax = 10;              % Duration of run

% First 42 Bessel function roots
wsort=alp*[...
  2.4048  3.8317  5.1356  5.5201  6.3801  7.0156
  7.5883 8.4173  8.6537  8.7715  9.7611  9.9362
 10.1735 11.0647 11.0864 11.6199 11.7916 12.2251
 12.3385 13.0152 13.3237 13.3543 13.5893 14.3726
 14.4755 14.7960 14.8213 14.9309 15.5898 15.7002
 16.0378 16.2234 16.4707 16.6983 17.0037 17.2412
 17.6159 17.8014 17.9599 18.0711 18.2876 18.4335];

disp(' '), disp(['The first forty-two ',...
                 'natural frequencies are:'])
disp(wsort)
    
[u,x,y,t] = memcirwv(r,th,r0,alp,w,tmax); % Solution
membanim(u,x,y,t);                        % Animation