Physics Beyond Our Senses
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 ;