% plot2osc.m
% plot position v time for two coupled oscillators
%
% Hastily written -- without lots of comments
% Raghuveer Parthasarathy
% Nov. 7, 2007

clear all
close all

w0 = 2*pi;      % =sqrt(g/l) = angular frequency, rad/sec
wc = 0.5*2*pi;  % =sqrt(k/m), coupling
wp = sqrt(w0*w0 + wc*wc);
A = 1;  % amplitude (arbitrary units)

t = 0.0:0.04:20;  % time array, seconds

% Analytic solution to oscillator motions
xa = A*cos(0.5*(wp-w0)*t).*cos(0.5*(wp+w0)*t);  % position of A
xb = A*sin(0.5*(wp-w0)*t).*sin(0.5*(wp+w0)*t);  % position of B



% The rest is just "cosmetic" -- plotting the analytic solution as an "animation"

% Define the equil. position of mass A is x=0, of mass B is +2 units to the right
figure;  hold on
axis([-1.1 3.1 0 3])
L = 2.8;  % pendulum length, for plotting; arbitrary units
ya = sqrt(L*L-xa.*xa);  % y position, for plotting
yb = sqrt(L*L-xb.*xb);
for j=1:length(t)
    if j==1
        % Give the user a chance to make the window bigger
        disp('Resize window; press Enter'); pause
    end
    clf; 
    subplot('Position', [0.1 .45 0.87 0.5]);
    hold on; xlabel('t'); ylabel('x');
    plot([xa(j) 0], [L-ya(j) L], 'b-', 'Linewidth', 2.0);  % pendulum string
    plot([xb(j)+2 2], [L-yb(j) L], 'b-', 'Linewidth', 2.0);  % pendulum string
    % Plot the "spring" as a line of varying thickness
    plot([xa(j) xb(j)+2], [(L-ya(j)) (L-yb(j))], '-', ...
        'Color', [0.7 0.7 0.7], 'Linewidth', 10.0/(abs(xa(j)-xb(j)-2)+0.3));
    plot(xa(j), L-ya(j), 'ko', 'MarkerSize', 15, 'MarkerFaceColor', 'r');
    plot(xb(j)+2, L-yb(j), 'ko', 'MarkerSize', 15, 'MarkerFaceColor', [0.0 0.5 0.0]);
    axis([-1.1 3.1 -0.25 3])
    subplot('Position', [0.1 0.25 0.87 0.13]);  hold on; axis([0 max(t) -1.1 1.1]);
    plot(t(1:j), xa(1:j), 'r.'); ylabel('xa');
    subplot('Position', [0.1 0.07 0.87 0.13]);  hold on; axis([0 max(t) -1.1 1.1]);
    plot(t(1:j), xb(1:j), '.', 'Color', [0.0 0.5 0.0]); xlabel('t'); ylabel('xb');
    pause(0.02)
end