% Illustration of planar wavefronts going through a lens and getting focused
% into a converging spherical wave

function main ()

  % lens index
   n=1.5; 

  % number of points, used for plotting
   N = 100;

  % radii of lens surfaces
   R1 = 0.5; 
   R2 = 1.5;

  % centers of circles (y coord is 0)
   O1 = -2.9;
   O2 = -O1;

  % focal length
   f = (n-1)*(1/R1+1/R2); f = 1/f;
   
   % theta0 determines the width of the lens
   theta0=pi/6;
   Theta = linspace(-theta0, theta0, N);

  % right face of the lens
   L1x = R1*cos(Theta)+O1;
   L1y =R1*sin(Theta); 

   % left size of the lens
   L2x=-R2*cos(Theta)+O2;
   L2y = R2*sin(Theta);

   % flat top part
   Topx = [L1x(N), L2x(N)];
   Topy = [L1y(N), L2y(N)];

   % flat bottom part
   Botx = [L1x(1) L2x(1)];
   Boty = [L1y(1), L2y(1)];

   % the lens
   Lensx = [L1x rv_vec(Topx), rv_vec(L2x), Botx];
   Lensy = [L1y rv_vec(Topy), rv_vec(L2y), Boty];

   % Parameters for graphing
   Lens_color  = [204, 226, 239]/256;
   Lens_border = 0.3*[1, 1, 1];
   lbw = 1.3; % lens border width
   wavefr_color = [1, 0, 0];
   wavefr_bdw   = 2;
   
   % spacing between wavefronts (both plane and spherical ones)
   spacing = 0.25;

   % 2*H is the height of the plane wavefronts
   H = L1y(N); 

   % theta2 = slope of the line going from the upper-right
   % end of the lens to the focus point
   theta2 = atan(L1y(N)/(f-L1x(N)));

   % Shape of the spherical wavefronts.
   Theta = linspace(-theta2, theta2, N);
   X = -cos(Theta);
   Y =  sin(Theta);

   S = -f; % start ploting waves from here to the right

   % number of frames in the movie
   num_frames = 10;
   Shifts = linspace(0, spacing, num_frames+1);

   % start at S+shift, plot the wavefronts
   for frame_no = 1:num_frames

      shift = Shifts(frame_no);
      
      s = S+shift;

      % plotting window
      figure(1); clf; hold on; axis equal; axis off;

      % plot the plane wavefronts
      while s < 0
         plot([s, s], [-H, H], 'color', wavefr_color, 'linewidth', wavefr_bdw);
         s = s + spacing;
      end
      
   
      % plot the spherical wavefronts
      s = s - 10*spacing; % backtrack a bit
      while s < f
         
         rho = f-s;
         
         if rho*Y(N) <= L1y(N)
            plot(rho*X+f, rho*Y, 'color', wavefr_color, 'linewidth', wavefr_bdw);
         end
         
         s = s + spacing;
         
      end

      % plot the lens
      fill(Lensx, Lensy, Lens_color, 'EdgeColor', Lens_border, 'LineWidth', lbw);
%      get(H)
%      return
      
      % Invisible points to force MATLAB to keep the
      % plotting window fixed.
      tiny = 0.15*spacing;
      white = 0.999*[1, 1, 1];
      plot(S-tiny,   H+tiny, 'color', white);
      plot(S-tiny,  -H-tiny, 'color', white);
      plot(f+tiny,   H+tiny, 'color', white);
      plot(f+tiny,  -H-tiny, 'color', white);

      % Rotate by 90 degrees
      set(gca, 'View', [90, 90])

      % save current file
      frame_file = sprintf('Frame%d.eps', 1000+frame_no);
      disp(frame_file);
      saveas(gcf, frame_file, 'psc2');
      pause(0.07)
   end

% The frames were converted to a movie with the command
% convert -antialias -loop 10000  -delay 8 -compress LZW Frame100* Lens_and_wavefronts.gif
   
function W = rv_vec(V)

   K = length(V);

   W = V;
   for i=1:K
      W(i) = V(K-i+1);
   end