function projectile_normal ( )

%*****************************************************************************80
%
%% projectile_normal uses normal equations on the projectile function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2019
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'projectile_normal:\n' );
  fprintf ( 1, '  Use normal equations to seek minimizer \n' );
  fprintf ( 1, '  of sum-of-squares error for projectile data.\n' );
%
%  Get the data.
%
  [ t, h ] = projectile_data ( );
  n = length ( t );
%
%  Set up the matrix A.
%
  A = zeros ( n, 3 );
  A(:,1) = ones(n,1);
  A(:,2) = t;
  A(:,3) = t.^2;
%
%  Set the system matrix AtA and right hand side Aty:
%
  AtA = A' * A;
  Aty = A' * h;
%
%  Solve 3x3 system for c;
%
  c = AtA \ Aty;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  c = (%20.16g,%20.16g,%20.16g)\n', c )
  fprintf ( 1, '  f(c) = %20.16g\n', projectile ( c ) )
  fprintf ( 1, '  f''(c) = (%g,%g,%g)\n', projectile_df ( c ) )
%
%  Compare with Professor Layton's values:
%
  c2 = [ 1097.714286, 4.573428573, -17.14285714 ];
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compare values from Professor Layton''s notes:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  c2 = (%20.16g,%20.16g,%20.16g)\n', c2 )
  fprintf ( 1, '  f(c2) = %20.16g\n', projectile ( c2 ) )
  fprintf ( 1, '  f''(c2) = (%g,%g,%g)\n', projectile_df ( c2 ) )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'projectile_normal:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  return
end
function e = projectile ( c )
  [ t, h ] = projectile_data ( );
  e = sum ( ( h - c(1) - c(2) * t - c(3) * t.^2 ).^2 );
  return
end
function dedc = projectile_df ( c )
  [ t, h ] = projectile_data ( );
  dedc1 = sum ( - 2.0 * ( h - c(1) - c(2) * t - c(3) * t.^2 ) );
  dedc2 = sum ( - 2.0 * ( h - c(1) - c(2) * t - c(3) * t.^2 ) .* t );
  dedc3 = sum ( - 4.0 * ( h - c(1) - c(2) * t - c(3) * t.^2 ) .* t.^2 );
  dedc = [ dedc1, dedc2, dedc3 ];
  return
end
function [ t, h ] = projectile_data ( )
  data = load ( 'projectile_data.txt' );
  t = data(:,1);
  h = data(:,2);
  return
end


