Code covered by the BSD License  

Highlights from
Quaternion belt-trick

image thumbnail

Quaternion belt-trick

by

 

This program demonstrates the quaternion belt-trick.

quaternion_belt_trick
% This program creates a Matlab AVI movie of the quaternion belt-trick, also known as 
% the "Dirac belt-trick" and the "waiter trick".
%
% Bill Davidson, quellen@yahoo.com
% created 14 November 2008
%
% Notes:
%
% (1) If you want to save time and see the end-result, just watch the included AVI movie.  
% It opens in Windows Explorer with Windows Media Player.
%
% (2) The movie is stereoscopic.  If you sit back at least a meter from the screen and 
% cross your eyes slightly to fuse the two images in your "mind's eye", you should get 
% a nice stereoscopic effect.
%
% (3) The program takes a few minutes to generate the AVI movie - about 10 minutes on 
% a 3GHz machine.  It produces one figure for each frame of video and takes a little 
% software snapshot of each one. There are 96 frames in the movie, the msximum allowed.  
% Any screensavers, pop-up windows, or moving about on the desktop while the program is 
% running can mess up the resulting movie.  This problem is due to the way MatLab does 
% its movies and it seems to be unavoidable.
%
% (4) The program utilizes a quaternion data class for the quaternion calculations.  
% I wrote the data class and it is very compact but also full-featured. You can separately 
% download it for free from the same place in Matlab Central where you got this program.  
% Matlab data classes, while a somewhat obscure topic, are very powerful, easy, and 
% interesting and are well worth learning about.  With data classes you can create a 
% working algebra with nearly any properties you can imagine.  My motivation for writing 
% the quaternion data class was to fully understand something mathematicians call the 
% "Cayley-Dickson Construction" which defines quaternions in terms of nested complex numbers.
%
% (5) The AVI file included is a compressed version of the original AVI output, which was 
% quite large, about 140Mb.
%
% Credits: The general idea for this program came from the book, Knots and Physics, by 
% Louis H Kauffman, World Scientific. Dr Kauffman included a clever BASIC program in that 
% book that inspired this program.

function quaternion_belt_trick

clear all;
close all;
clear classes;
clc;
addpath('.\my_classes'); % This line assumes that the program resides in the directory just above "my_classes"
format compact
format short

methods(quaternion); % This line loads the quaternion data class

rad1=100; % outer circle
rad2=20; % inner circle

center1=[rad1,rad1]; % left
center2=[rad1+200,rad1]; % right

u1234=quaternion(0,0,1,0); % gets rotated
vl234=quaternion(0,0,0,1); % axis of rotation

max_frames=96;
%max_frames=5;
dphi=2*pi/(max_frames-1); % include both ends
for k=1:max_frames

    [s,err]=sprintf('k=%d of %d',k,max_frames);disp(s);

    current_figure=figure;
    resize_figure(current_figure);

    hold on;
    %axis equal;

    plot_circles(current_figure,center1,center2,rad1,rad2);

    [s,err]=sprintf('%2d of %d',k,max_frames);
    text(center1(1)-rad1*0.9,center1(2)+rad1*0.9,s);
    text(center2(1)-rad1*0.9,center2(2)+rad1*0.9,s);
    
    plot_belts(current_figure,u1234,center1,center2,rad1,rad2);

    u1234=rotation(dphi,vl234,u1234);

    F(k)=getframe;

    close(current_figure);

end

save_movie(F);

%------------------------------------

function resize_figure(current_figure)

screen_size=get(0,'ScreenSize');
screen_xwidth=screen_size(3);
screen_ywidth=screen_size(4);
screen_center(1)=screen_size(3)/2;
screen_center(2)=screen_size(4)/2;

xfactor=1.0; % width
yfactor=0.6; % height

w=get(current_figure,'Position');
xwidth=screen_xwidth*xfactor;
ywidth=screen_ywidth*yfactor;
w(1)=screen_center(1)-(xwidth/2);
w(2)=screen_center(2)-(ywidth/2);
w(3)=xwidth;
w(4)=ywidth;
set(current_figure,'Position',w);

%------------------------------------

function plot_circles(current_figure,center1,center2,rad1,rad2)

mi=100;
iis=1:(mi+1);
thetas=2*pi*(iis-1)/mi;
als=cos(thetas);
bls=sin(thetas);
for ii=2:(mi+1)
    plot(center1(1)+rad1*[als(ii-1),als(ii)],center1(2)+rad1*[bls(ii-1),bls(ii)]); % left one
    plot(center1(1)+rad2*[als(ii-1),als(ii)],center1(2)+rad2*[bls(ii-1),bls(ii)]);
    plot(center2(1)+rad1*[als(ii-1),als(ii)],center2(2)+rad1*[bls(ii-1),bls(ii)]); % right one
    plot(center2(1)+rad2*[als(ii-1),als(ii)],center2(2)+rad2*[bls(ii-1),bls(ii)]);
end
drawnow;

%------------------------------------

function plot_belts(current_figure,u1234,center1,center2,rad1,rad2)

mi=80;
mj=10;
rxs=zeros(1,mj+1);
rys=zeros(1,mj+1);
iis=1:(mi+1);
rs=rad1-((rad1-rad2)*(iis-1)/mi);
for ii=1:(mi+1)
    jjs=1:(mj+1);
    xs=(20/mj)*((jjs-1)-(mj/2));
    ys=sqrt(rs(ii)*rs(ii)-xs.*xs);
    theta=2*pi*(ii-1)/mi;
    for jj=1:(mj+1)
        efgh=quaternion(0,xs(jj),ys(jj),0);
        wxyz=rotation(theta,u1234,efgh);
        dwxyz=double(wxyz);
        rxs(jj)=dwxyz{2};
        rys(jj)=dwxyz{3};
        rzs(jj)=dwxyz{4};
    end
    plot1(center1(1)+rxs-0.1*rzs,center1(2)+rys,rzs,center1(1),center1(2),rad2); % left one
    plot1(center2(1)+rxs+0.1*rzs,center2(2)+rys,rzs,center2(1),center2(2),rad2); % right one
    drawnow;

end

%------------------------------------

function plot1(sxs,sys,szs,xc,yc,rad2)
hidz=(szs<0);
txs=sxs-xc;
tys=sys-yc;
hidr=(((txs.*txs)+(tys.*tys))<(rad2*rad2));
hid=hidr&hidz;
sxs(hid)=NaN;
sys(hid)=NaN;
plot(sxs,sys); % left one

%------------------------------------

function save_movie(F)
for i=0:100
    [movie_fname,err]=sprintf('quaternion_belt_trick_movie_%d.avi',i)
    if 0==exist(movie_fname),break,end
end
movie2avi(F,movie_fname,'compression','None'); % there are a lot of chices here but this one seems to work well. 1204k per frame

Contact us