function F=chainavi(n, m)
% make movie using chain(n,m,phi)
%   n chain atoms, eigenmode m, 
%   n>0 for acoustic, n<0 for optic mode
%   m>0 for longitudinal, m<0 for transverse

if nargin<2, m=-2; end;
if nargin<1, n=16; end;

for k=1:40
    chain(n,m,k*pi/20);
    F(k)=getframe;
end;

if nargout==0,
    movie2avi(F,'chain.avi','quality',100);
end;
end

function hc=chain(n, m, phi)
% chain(n,m,phi) linear chain with n members
%   eigenmode m, at phase phi (in radians). 
%   n>0 for acoustic, n<0 for optic mode
%   m>0 for longitudinal, m<0 for transverse

if nargin<3, phi=pi/2; end;
if nargin<2; m=-2; end;
if nargin<1; n=-16; end;

if m*n==0 return; end;
if n<0
    n = -n;
    ao = -0.55;
else
    ao = 1; 
end;
if m<0
    m = -m;
    a = [0 0 4];
    d = [3.5 0 0];
else
    a = [1 0 0]; 
    d = [3.5 0 0];
end;

lstick = [0 0 8];
Connect(-lstick,lstick,1,0.8*[1 1 1]);
hold on;
oldpos = [0 0 0];
for i=1:n+1
    phase = m*pi*i/(n+1);
    if mod(i,2)==0
        col = [0.1 0.1 1];
        pos = i*d+a*sin(phase)*sin(phi);
        r = 1;
    else
        col = [1 0.1 0.1];
        pos = i*d+ao*a*sin(phase)*sin(phi);
        r = 1.15;
    end;
    if i<n+1, Atom(pos,r,col); end;
    Spring(oldpos, pos, 1.5, 0.4, 0.2);
    oldpos = pos;
end;
Connect(pos-lstick,pos+lstick,1,0.8*[1 1 1]);

view(-60,2); 
axis equal; axis off; camlight; lighting gouraud;
hold off;
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 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
