% animate_2D_vib.m
%
% This MATLAB m-file "animates" a vibrational mode
% for the 2-D lattice.
% K. Beers. MIT ChE. 10/4/2002
function iflag_main = animate_2D_vib();
iflag_main = 0;

% read in results of normal mode analysis
load lattice_2D_vib.mat;

% display normal modes and their frequencies
disp(' ');
disp('Normal mode frequencies (rad/sec) : ');
disp('-----------------------------------');
for j=1:length(freq)
    disp([int2str(j), '  ', num2str(freq(j))]);
end
disp(' ');
    
% select normal mode to animate
k_show = input('Enter number of mode to visualize : ');

% Create vector of theta values for a single
% oscillation.
num_frames = 50;
theta_vect = linspace(0,2*pi,num_frames);
[sin_max,iframe_max] = max(sin(theta_vect));
if(length(iframe_max) > 1)
    iframe_max = iframe_max(1);
end

% Set amplitude of oscillation so that the
% maximum displacement is equal to 1/3 times
% the equilibirum bond length.
displace_max = max(abs(V_eig(:,k_show)));
a = l_b/3/displace_max;

% Create movie by plotting 2-D lattice positions
title_phrase = ['2-D lattice, N = ', int2str(N), ...
        ', mode # = ', int2str(k_show), ...
        ', \omega = ', num2str(freq(k_show))];

movie_fig = figure;
for iframe = 1:num_frames
    hold off;
    % compute departure vector
    delta = a*sin(theta_vect(iframe))*V_eig(:,k_show);
    q = q_min + delta;
    % plot positions of atoms
    for i=1:N
        for j=1:N
            [k,kx,ky] = get_master_index(i,j,N);
            plot(q(kx),q(ky),'o');
            hold on;
        end
    end
    axis([ (-l_b) (N+l_b) (-l_b) (N+l_b) ]);
    title(title_phrase);
    % get frame for movie
    Movie_of_mode(iframe) = getframe;
    % if maximum amplitude, make separate plot
    if(iframe == iframe_max)
        figure;
        for i=1:N
        for j=1:N
            [k,kx,ky] = get_master_index(i,j,N);
            plot(q(kx),q(ky),'o');
            hold on;
            plot(q_min(kx),q_min(ky),'x');
            x_line = [q_min(kx), q(kx)];
            y_line = [q_min(ky), q(ky)];
            plot(x_line,y_line);
        end
        end
        axis([ (-l_b) (N+l_b) (-l_b) (N+l_b) ]);
        title([title_phrase, ', max. amplitude']);
        figure(movie_fig);
    end
    
end

% Show several periods of oscillation, repeat
% again and again until stop.
iquit = 1;
times_repeat = 10;
FPS = num_frames/2;  % frames per second
while(iquit)
    movie(Movie_of_mode,times_repeat,FPS);
    iquit = input('Quit (0) or Show again (1) : ');
end

iflag_main = 1;
return;



% ==============================================================
% get_master_index.m
%
% This function calculates the master index label of
% a site in the 2D lattice of NxN sites.
function [k,kx,ky] = get_master_index(i,j,N);

k = (i-1)*N+j;
kx = 2*k-1;
ky = 2*k;

return;