% spheregridmovie.m
% Produces a movie illustrating figure 1(b) of 
%  Wave Propagation Software, Computational Science, and Reproducible Research
%  by R. J. LeVeque, http://www.amath.washington.edu/~rjl/pubs/icm06

n = 20;
[xc,yc] = meshgrid(linspace(-1,1,n+1),linspace(-1,1,n+1));
[xp1,yp1,zp1]=mapc2m(xc,yc);  % upper hemisphere

xp0 = xc;
yp0 = yc;
zp0 = zeros(size(xc));

nframes = 11;   % number of frames in the movie
for Frame=1:nframes
  % for this frame, interpolate between rectangle and sphere:
  alf = (Frame-1)/(nframes-1);
  xp = (1-alf)*xp0 + alf*xp1;
  yp = (1-alf)*yp0 + alf*yp1;
  zp = (1-alf)*zp0 + alf*zp1;

  % now find corresponding points on lower hemisphere, and add rotation:
  theta = -alf*pi;
  xpn1 = -2 - xp;
  xpn = cos(theta)*(xpn1+1) + sin(theta)*zp - 1;
  ypn = yp;
  zpn = -sin(theta)*(xpn1+1) + cos(theta)*zp;

  clf
  surf(xp,yp,zp,zeros(size(xp)))  % upper hemisphere
  hold on
  surf(xpn,ypn,zpn,ones(size(xp)))  % lower hemisphere
  hold off

  colormap([1 1 1; 1 .8 .8])
  axis([-3 1 -1 1 -2 1])
  daspect([1 1 1])
  axis off
  view([30 30])
  % makeframejpg % make jpg file for this frame (uses claw/matlab/ m-files)
  pause
  end

% If a jpg file was made for each frame, execute    
%      convert -resize 600x600 frame*.jpg spheremovie.gif
% or
%      convert -delay 10 -adjoin frame*.jpg spheremovie.mpeg
% or similar command in unix to create the movie