function f=modesavi(aviname)
% movie of oscillation modes

if nargin==0, aviname='modes.avi'; end;
n=40;
for i=1:n
    phi=2*pi*i/n;
    mode(0,phi);
    mode(1,phi);
    mode(2,phi);
    mode(3,phi);
    axis equal; axis off; light; lighting gouraud;
    view(70,15);
    hold off;
    f(i)=getframe;
end;
if nargout==0, 
    movie2avi(f, aviname, 'quality',95);
end;

function mode(type,phi,pos)
% oscillation mode of cubic cell

if nargin<1, type=0; end;
if nargin<2, phi=pi/2; end;
if nargin<3, pos=type; end;

u = 5;
delta = 24;
s=2;
a = [0.25 0.25 0.25];
if type>0
    a = 0.8*a;
    a(1,type) = -0.3;
end;
x0 = 0; y0 = 0; z0 = 0;
if pos<2, z0 = delta; end;
if mod(pos,2)==1, y0 = delta; end;
U = u*eye(3);
col(1,:) = [1 0.2 0.2];
col(2,:) = [0.2 1 0.2];
col(3,:) = [0.2 0.2 1];
for i=1:6
    k=ceil(i/2);
    r(i,:)=(-1)^(i-1)*U(k,:)*(1+a(k)*cos(phi));
    Atom(r(i,:)+[x0 y0 z0],s,col(k,:));
    hold on;
end;
[X,Y]=meshgrid(2*u*[-1,1]);
Z=X*0;
h1=surf(X+x0,Y+y0,Z+z0);
h2=surf(Z+x0,Y+y0,X+z0);
h3=surf(X+x0,Z+y0,Y+z0);
set([h1 h2 h3],'Linewidth',2,'Edgealpha',0.4,'Facealpha',0);
for i=1:3
    Connect(2*U(i,:)+[x0 y0 z0],-2*U(i,:)+[x0 y0 z0],0.2*s,0.8*[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);

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