function ShowGN(A0,P0,del,m)
%
%% showgn() illustrates the Gauss-Newton retrieval of information from noise.
%
% Pre:  A0 and P0 define an untilted orbit.
%       del is the relative error in a radius vector sampling.
%       m is the number of noisy samplings of the radius vector.
%
% Post: Illustrates the Gauss-Newton ``retrieval'' of A0 and
%       P0 from the noisy data.
%  
% Usage: ShowGN(5,1,.1,50)     
%      
  fprintf ( 1, '\n' );
  fprintf ( 1, 'showgn()\n' );

   con_size = 20;        % Size of contour plot
   kmax     = 25 ;       % Maximum number of GN steps.
   Lmax     = 8;         % Maximum step length factor
   tol      = .000001;   % Size of gradient for termination.
   
   % Generate a noisy orbit
   
   theta = linspace(0,2*pi,m)';
   c = cos(linspace(0,2*pi,m)');
   r = ((2*A0*P0)./(P0*(1-c) + A0*(1+c))).*(1 + del*randn(m,1));
   plist = [r c];
   
   % Plot the true orbit and the radius vector samples
   
   close all
   Orbit(linspace(0,2*pi,200),A0,P0,0,'-');
   hold on
   x = r.*c;
   y = r.*sin(theta);
   plot(x,y,'o')
   hold off   
   pause(.01)

   % Generate Contour Plot

   range = linspace(.9*min(r),1.5*max(r),con_size);
   Z = zeros(con_size,con_size);
   for i=1:con_size
      for j=1:con_size
	   	f = rho([range(j);range(i)],plist);
	    Z(i,j) = (f'*f)/2;
      end
   end
   figure
   C = contour(range,range,Z,20);
   title('Enter Contour Labels (Optional). To quit, strike <return>.')
   xlabel('A'); ylabel('P'); Clabel(C,'manual')
  
   % Initialize for iteration
  
   title('Enter Starting Point for Gauss-Newton')
   zc = zeros(2,1);
   [zc(1),zc(2)] = ginput(1); 
   Jc = Jrho(zc,plist);
   Fc =  rho(zc,plist);
   wc = (Fc'*Fc)/2;
   k=0;
   gnorm = norm(Jc'*Fc);
   clc
   home
   disp(sprintf('A0 = %5.3f,  P0 = %5.3f,   Realtive Noise = %5.3f,  Samples = %3.0f\n',A0,P0,del,m))
   disp('Iteration     wc         A          P          norm(grad)')
   disp('----------------------------------------------------------')
   disp(sprintf('   %2.0f   %10.4f %10.6f %10.6f       %5.3e ',k,wc,zc(1),zc(2),gnorm))
   figure(3)
   
   % Gauss-Newton with linesearch
   
   while (gnorm > tol ) & (k<kmax)
      if k<=2
         [zc,Fc,wc,Jc] = GNStep('rho','Jrho',zc,Fc,wc,Jc,plist,Lmax,0);
      else 
	     [zc,Fc,wc,Jc] = GNStep('rho','Jrho',zc,Fc,wc,Jc,plist,Lmax,1);
      end
	  k=k+1;
	  gnorm = norm(Jc'*Fc);
	  disp(sprintf('   %2.0f   %10.4f %10.6f %10.6f       %5.3e ',k,wc,zc(1),zc(2),gnorm))
   end
   figure(1)
   hold on
   Orbit(linspace(0,2*pi),zc(1),zc(2),0,'.');
   hold off
  
  return
end
   
	 
  
