% Runs membrane codes with specific parameters
clear

% Rectangular (1) vs Circular (2)
type = 2;

% Dimensions
a    = 0.6;             % Length of membrane
b    = 0.6;             % Height of membrane
R    = 0.6;             % Radius of membrane

% Parameters of membrane
T    = 2000;            % Tension
sigma = 0.26;           % Membrane density
alp  = sqrt(T/sigma);   % Wavespeed

% Applied Impulse
f0    = -1.0;           % Strength of applied force                  
x0   = a/2;             % Location of concentrated force
y0   = b/2;
r0   = R/2;           % Location of concentrated force

% Time steps (Need to adjust for more accurate spectral analysis)
fmax = 500;            % Max frequency desired
fs   = 2*fmax;          % sampling frequency
dt   = 1/fs;            % time step
tmax = 1;               % Duration of run
Nt   = ceil(tmax/dt);   % Number of timesteps   
fram = 30*tmax;
Dt   = ceil(Nt/fram);
Dt=1;

% Grid of independent variables
figure(1)
if type == 1
    [u,x,y,t] = recmem(a,b,alp,f0,x0,y0,tmax,Nt); % Solution
else
    th   = linspace(0,2*pi,100); 
    r    = linspace(0,R,30);
    [u,x,y,t] = circmem(r,th,r0,alp,f0,tmax,Nt); % Solution
end

% Animation
xmin=min(x(:)); xmax=max(x(:));
ymin=min(y(:)); ymax=max(y(:));
xmid=(xmin+xmax)/2; ymid=(ymin+ymax)/2;
d=max(xmax-xmin,ymax-ymin)/2; Nt=length(t);
range=[xmid-d,xmid+d,ymid-d,ymid+d,3*min(u(:)),3*max(u(:))];
for j=1:Dt:Nt
    surf(x,y,u(:,:,j)), axis(range)
    xlabel('x axis'), ylabel('y axis'),  zlabel('u axis')
    titl=sprintf('MEMBRANE POSITION AT T=%6.4f',t(j));
    title(titl), 
    shading interp
    colormap copper
    drawnow, shg, pause(.1)
  end

% Frequency analysis
s=size(u);
N=s(3);
[c,j0]=max(max(abs(u(:,:,2))));
[c,i0]=max(abs(u(:,j0,2)));

for j=1:N
    Y(j)=u(i0,j0,j);
end
figure(2)
plot(t,Y)
figure(3)
fanalf(Y-mean(Y),max(t));

% Add movies?