function analemma ( ecc, lon_deg, obliq_deg )

%*****************************************************************************80
%
%% analemma() computes the analemma.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 November 2018
%
%  Author:
%
%    Original C version by Brian Tung.
%    Octave version by John Burkardt.
%
%  Input:
%
%    real ECC, the orbital eccentricity.
%
%    real LON_DEG, the longitude of the perihelion in degrees.
%
%    real OBLIQ_DEG, the obliquity in degrees.
%
  days = 365.242;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'analemma():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Compute and plot the analemma, equation of time, and declination.\n' );
  fprintf ( 1, '  This program is based on a C program by Brian Tung.\n' );
%
%  Set defaults.
%  Internally, longitude and obliquity are in radians.
%
  if ( nargin < 1 )
    ecc = 0.01671;
  elseif ( ischar ( ecc ) )
    ecc = str2num ( ecc );
  end

  if ( nargin < 2 )
    lon = 1.347;
    log_deg = lon * 180.0 / pi;
  elseif ( ischar ( lon_deg ) )
    lon = str2num ( lon_deg ) * pi / 180.0;
  else
    lon = lon_deg * pi / 180.0;
  end

  if ( nargin < 3 )
    obliq = 0.4091;
    obliq_deg = obliq * 180.0 / pi;
  elseif ( ischar ( obliq_deg ) )
    obliq = str2num ( obliq_deg ) * pi / 180.0;
  else
    obliq = obliq_deg * pi / 180.0;
  end
%
%  Report data.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %g, the orbital eccentricity\n', ecc );
  fprintf ( 1, '  %g, the longitude of the perihelion in degrees\n', lon_deg );
  fprintf ( 1, '  %g, the obliquity in degrees.\n', obliq_deg );
%
%  Compute the data using vector operations.
%
  n = 10001;

  t01 = zeros ( 1, n );
  eot = zeros ( 1, n );
  dec = zeros ( 1, n );

  i = 1 : n;
  f = ( i - 1 ) / ( n - 1 );
  tau = 2.0 * pi * f;
%
%  Set theta to the current longitude. 
%
  theta = atan2 ( sqrt ( 1.0 - ecc * ecc ) * sin ( tau ), cos ( tau ) - ecc );
%
%  Rotate clockwise in XY plane by theta, corrected by lon.
%
  x1 = cos ( theta - ( lon - pi / 2.0 ) );
  y1 = sin ( theta - ( lon - pi / 2.0 ) );
  z1 = 0.0;
% 
%  Rotate counter-clockwise in XZ plane by obliq.
%
  x2 = cos ( obliq ) * x1 + sin ( obliq ) * z1;
  y2 = y1;
  z2 = - sin ( obliq ) * x1 + cos ( obliq ) * z1;
% 
%  Set t equal to real time from tau and
%  rotate counter-clockwise by t, corrected by lon 
%
  t = tau - ecc * sin ( tau );
  x3 =   cos ( t - ( lon - pi / 2.0 ) ) .* x2 + sin ( t - ( lon - pi / 2.0 ) ) .* y2;
  y3 = - sin ( t - ( lon - pi / 2.0 ) ) .* x2 + cos ( t - ( lon - pi / 2.0 ) ) .* y2;
  z3 = z2;

  t01 = t / 2.0 / pi;
  eot = - atan2 ( y3, x3 ) * 4.0 * 180 / pi * days / ( days + 1.0 );
  dec = asin ( z3 ) * 180.0 / pi;
%
%  Plot the equation of time.
%
  figure ( 1 );
  plot ( t01, eot, 'b-' );  
  grid ( 'on' );
  xlabel ( '<---Normalized Date--->' );
  ylabel ( '<---Minutes Early/Late--->' );
  title ( 'Equation of Time' );
  print ( '-dpng', 'eot.png' );
%
%  Plot the declination.
%
  figure ( 2 );
  plot ( t01, dec, 'b-' );
  grid ( 'on' );
  xlabel ( '<---Normalized Date--->' );
  ylabel ( '<---Degrees North/South--->' );
  title ( 'Declination' );
  print ( '-dpng', 'declination.png' );
%
%  Plot the analemma.
%
  figure ( 3 );
  plot ( eot, dec, 'b-' ); 
  grid ( 'on' );
  xlabel ( '<---Minutes Early/Late--->' );
  ylabel ( '<---Degrees North/South--->' );
  title ( 'Analemma' );
  print ( '-dpng', 'analemma.png' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'analemma()\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return;
end

