% 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;