function F=coupledavi(mode)
% movie of coupled oscillations
% mode may be 'Symmetric', 'Antisymmetric', or 'Beat' 

if nargin<1, mode='Sym'; end;
if lower(mode(1))=='s',
    n = 40; A = 0; B = pi/18;
elseif lower(mode(1))=='a',
    n = 48; A = pi/18; B = 0;
else
    n = 240; A = pi/36; B = A;
end;

for k=1:n
    phi1=A*sin(k*pi/24)+B*sin(k*pi/20);
    phi2=A*sin(k*pi/24)-B*sin(k*pi/20);
    coupled(phi1,phi2);
    F(k)=getframe;
end;

if nargout==0,
    movie2avi(F,'coupled.avi','quality',95);
end;
end

function coupled(phi1,phi2)
% coupled oscillations

if nargin<1, phi1=pi/12; end;
if nargin<2, phi2=-phi1; end;

Connect([-4 -1 0],[10 -1 0],0.2,[0 0.3 0]);
hold on;
Connect([-4 0 -12],[10 0 -12],0.5,[0 0.3 0]);
A0=[0 0 0];
B0=[6 0 0];
A1=A0-4*[sin(phi1) 0 cos(phi1)];
B1=B0-4*[sin(phi2) 0 cos(phi2)];
Pendulum(A0,phi1,[1 0.2 0.2]);
Pendulum(B0,phi2,[0.2 0.2 1]);
Spring(A1,B1,0.5,0.5,2);
view(0,0);
hold off;
set(gcf,'Color',0.95*[1 1 1]);
end

function h=Atom(pos,r,color)
% atom at position pos, radius r, and color

if nargin<3, color=0.8*[1 1 1]; end;
if nargin<2, r=1; end;
if nargin<1, pos=[0 0 0]; end;

[x,y,z]=sphere;
h=surf(pos(1)+r*x,pos(2)+r*y,pos(3)+r*z);
set(h,'EdgeAlpha',0,'FaceColor',color);
end

function h=Connect(a,b,r,color)
% connection cylinder between atoms

if nargin<4, color=0.9*[1 1 1]; end;
if nargin<3, r=0.1; end;
if nargin<2, return; end;

d = b-a;
l = sqrt(sum(d.*d));
if l==0, return; end;
[x,y,z] = cylinder([r r]);
z = l*z;
h=surf(a(1)+x, a(2)+y, a(3)+z);
set(h,'EdgeAlpha',0,'FaceColor',color);
phirot = 180/pi*acos(d(3)/l);
ax = [-d(2) d(1) 0];
if (sum(ax.*ax)==0)
    ax = [1 0 0];
end;
rotate(h,ax,phirot,a);
end

function Pendulum(a, phi, col)
% pendulum for coupled oscillations

if nargin<1, a=[0 0 0]; end;
if nargin<2, phi=pi/10; end;
if nargin<3, col=0.4*[1 1 1]; end;
L=[0 0 -10];
b=a+L;
hc=Connect(a,b,0.1);
hold on;
hs=Atom(b,0.8,col);
hb=Connect(a,a+0.05*L,0.3);
hspr=Connect(a+0.38*L,a+0.42*L,0.3);
set([hb hc hspr],'FaceColor',0.8*[1 1 1]);
rotate([hb hs hc hspr],[0 1 0],phi*180/pi,a);
axis equal; axis off; camlight; lighting gouraud;
end

function hc=Spring(a, b, e, r0, LW)
% plot 3D spring

if nargin<1, a=[0 0 0]; end;
if nargin<2, b=a+[0 0 1]; end;
if nargin<3, e=sqrt(dot(b-a,b-a))/6; end;
if nargin<4, r0=e/2; end;
if nargin<5, LW=1; end;
PHI = 16*pi;     % 8 windings
d = b-a;
l = sqrt(sum(d.*d));
L = l-2*e;
r = sqrt(r0*r0-L*L/PHI/PHI);
angle = linspace(0,PHI,400);
x = [0 0 r*sin(angle) 0 0]+a(1);
y = [0 0 r*cos(angle) 0 0]+a(2);
z = [0 e linspace(e,l-e,400) l-e l]+a(3);
hc = plot3(x,y,z,'-k','Linewidth',LW);
phirot = 180/pi*acos(d(3)/l);
ax = [-d(2) d(1) 0];
if (sum(ax.*ax)==0)
    ax = [1 0 0];
end;
rotate(hc,ax,phirot,a);
end

