A Study in Applied Spirituality

The following function will generate the animation in Matlab. If you are interested, feel free to use it. Please cut and paste, looking at the HTML source is a waste of time, because of the syntax highlighting directives.

The highlighting is done using highlight.m (© 2003 Guillaume Flandin)

clear ; % clear variables, if any
close ; % close the current figure, if any

h = subplot(1,1,1) ; % get the handle of the plot window
set(h,'xtick',[]) ;  % so that we can disable the X and
set(h,'ytick',[]) ;  % Y tick marks

title('Superluminal motion and our perception of it') ; % title of the plot window
% mov = avifile('grb.avi') % if you want to make the movie from matlab 
fprintf('Hit <enter> if ready to go: ') ;
pause ; % to resize the figure

nStep = 50 ; % number of photons within the window
nY = 25 ;    % height at which the object is flying by
nZm = 4 ;    % a zoom factor, larger numbers make the animation smoother
nMovie = 8 ; % movie frame grabber decimation, kind of inverse of nZm
nMarg = 5 ;  % number of photons outside the window
xMarg = nMarg*nZm ;

beta = 10 ;  % the speed of the fly-by object

% a photon shape
[x,y] = cylinder(1,32) ; % this gives a unit circle, in fact, two rows of the same circle
xG = x(1,:)' ;           % use the first row, x coordinates
yG = y(1,:)' ;           % y coordinates of the circle
% the flying object size
lObj = 3 ;
xObj = xG * lObj ;       % flying object is just a bigger circle
yObj = yG * lObj ;

% max number of photons
maxG = ceil(nStep+nMarg*2) ; % nStep in the window and nMarg on either side
nG = 0 ;                     % init the photon counter

% the order in which the photons, rays, object etc are created is
% important, because the rays should not fall on top of the object or
% the phantoms
hGs = zeros(maxG,1) ; % handles to the photons
hRays = zeros(maxG,1) ; % handles to the light rays
for iG=1:maxG
    hRays(iG) = line([0 -lObj],[0 -lObj]);   % make the photons outside the plotting window
    set(hRays(iG),'color','r') ;             % in red color
    hGs(iG) =  patch(xG-lObj,yG-lObj,'r') ;  % light rays also in red color
end

% make the fly-by object
xStart = 0 ;                 % starting position
xEnd = nStep*nZm ;           % ending position
yStart = nY*nZm ;            % always at the same height
hTrace = line([xStart, xStart],[yStart yStart]) ;  % trace line of the fly-by object
hObj = patch(xObj+xStart,yObj+yStart,'w') ; % object (in white) handle (over the trace)
set(hTrace,'color','k') ;                   % set the trace to black 

% draw two phantoms, which will be moved later
hPhantom1 = patch(xObj-lObj,yObj-lObj,'g') ;
hPhantom2 = patch(xObj-lObj,yObj-lObj,'y') ;

% make the plot window canvas
axis([0 xEnd 0 1.5*yStart]) ;

% observer's position
xO = (xStart+xEnd)/2 ;  % at the center of the plotting window
yO = 0 ;                % at the bottom

% xincrements
dx = zeros(maxG,1) ;    % init to zero
dy = zeros(maxG,1) ;
% store the photon emitted points
xEmit = zeros(maxG,1) ; % init to zero
yEmit = zeros(maxG,1) ;
% to save the frame for the movie
iFrame = 0 ;
% boolean for the first photon crossing the observer
FirstCross = true ;
% move the object and create photons.
% the scaling of the end point is necessary, because if the object is
% moving very fast, it takes quite a while before the photons reach the
% observer.  So we have to keep moving the object
for iStep=-xMarg:xEnd*(1+beta) 
    % move the object
    xNow = xStart+iStep ; % by one step
    yNow = yStart ; % constant
    set(hObj,'xdata',xObj+xNow) ;
    set(hTrace,'xdata',[xStart xNow]) ;
    if (mod(iStep,nZm) == 0 && nG<maxG)
        % add a photon
        nG = nG + 1;
        set(hGs(nG),'xdata',xG+xNow,'ydata',yG+yNow) ; % put the photon (already created outside the loop)
        % photon direction cosine
        ddx = (xO-xNow) ;
        ddy = (yO-yNow) ;
        % normalize the photon speed
        % divide the photon step size by beta because we want the object 
        % to be beta times faster than the photon
        r = sqrt(ddx*ddx+ddy*ddy)*beta ; 
        dx(nG) = ddx/r ;
        dy(nG) = ddy/r ;
        xEmit(nG) = xNow ;
        yEmit(nG) = yNow ;
        % set(hRays(nG),'xdata', [xNow xNow+dx(nG)], 'ydata',[yNow yNow+dy(nG)]) ; 
    end
    % move the photons
    xMin = xEnd ;
    xMax = xStart ;
    nAbove = 0 ;
    for iG=1:nG
        % new position of the photon
        xpos = get(hGs(iG),'xdata')+dx(iG) ;
        ypos = get(hGs(iG),'ydata')+dy(iG) ;
        set(hGs(iG),'xdata',xpos,'ydata',ypos) ; % move the photon to the new position
        % the "center" of the photon
        xcent = mean(xpos) ;
        ycent = mean(ypos) ;
        % redraw the line from the emitted position and the current
        % position
        xpos = [xEmit(iG) ; xcent] ;
        ypos = [yEmit(iG) ; ycent] ;
        set(hRays(iG),'XData',xpos,'YData',ypos) ;
        
        % if the photon has dropped below the observer, draw the phantom
        if (ycent<yO)
            if (FirstCross)  
                % put a black dot ("core") where the two apparent objects are receding from.
                FirstCross = false ;
                patch(xG+xEmit(iG),yG+yEmit(iG),'k')
            end
            if (xEmit(iG) < xO-yEmit(iG)/beta)
                if (xEmit(iG) < xMin) 
                    xMin = xEmit(iG) ;
                end
            else
                if (xEmit(iG) > xMax) 
                    xMax = xEmit(iG) ;
                end
            end
            % color of the phantoms, shows blue shift or redshift
            cF = (xMax-xMin)/(xEnd-xStart) ;
            col = [cF 0 1-cF] ;
            % position of the phantoms
            if (xMin ~= xEnd)
                set(hPhantom1,'xdata',xObj+xMin,'ydata',yObj+yEmit(iG)) ;
            end
            if (xMax ~= xStart)
                set(hPhantom2,'xdata',xObj+xMax,'ydata',yObj+yEmit(iG)) ;
            end
            if (cF>0 && cF<1) 
                col(2) = min(col(1), col(3)) ;
                col = col/max(col) ;
                set(hPhantom1,'facecolor',col) ;
                set(hPhantom2,'facecolor',col) ;
            end
        else
            % count the number of photons above y=0.
            % if none left, we can exit the loop.
            nAbove = nAbove + 1 ; 
        end
    end
    drawnow ; % redraw the window
    if (mod(iStep,nMovie) == 0)
        F = getframe(gca);
        %  mov = addframe(mov,F); % if making the movie
        iFrame = iFrame+1 ;
        M(iFrame) = F ;
    end
    if (xMax > xEnd+xMarg)
        break ;
    end
    % see if the loop is useful
    if (iStep > xEnd && nAbove < 1)
        break ;
    end
end
% mov = close(mov) ; % if making the movie
% movie(M) ; % to review the frames
% movie2avi(M,'grb.avi','compression','none') ; % to save a higer quality movie
return ;