function test_grad
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Purpose: To test gradient of cost function
%
%  (c) 2004  Data Assimilation Research Centre
%
%  Written  by Amos Lawless
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=1.0d1;       % sigma
b=8.0d0/3.0d0; % rho
r=2.8d1;       % beta
h=0.01d0;        % time step
tstep=10;     % Number of steps
freq=1;
len=12;
x=zeros(tstep+1,1);
y=zeros(tstep+1,1);
z=zeros(tstep+1,1);
alpha_vec=zeros(len,1);
fn_vec=zeros(len,1);
resid_vec=zeros(len,1);
x(1)=1.0d0;
y(1)=2.0d0;
z(1)=1.5d0;
tstep=tstep+1;
%
% Set up data and matrix D
%
[x,y,z]=modeuler(tstep,h,0,z,y,x,a,r,b);
D=zeros(tstep,1);
n_obs = (tstep-1) / freq + 1;
datx=zeros(n_obs,1);
daty=zeros(n_obs,1);
datz=zeros(n_obs,1);
for i=1:freq:tstep
  datx(i)=x(i);
  daty(i)=y(i);
  datz(i)=z(i);
  D(i) = 1.;
end   
%
      x_vec = rand(1,3);
% Make x_vec order one
      scalar1 = 0.0D0;
      for i=1:3
        scalar1 = scalar1 + x_vec(i)*x_vec(i);
      end
      scalar1 = 1.0D0 / sqrt(scalar1);
      for i=1:3
        x_vec(i) = x_vec(i) * scalar1;
      end
%
      x(1) = x_vec(1);
      y(1) = x_vec(2);
      z(1) = x_vec(3);
      [x,y,z]=modeuler(tstep,h,0,z,y,x,a,r,b);
      [f1,grad1] = calcfg(tstep,h,x,y,z,datx,daty,datz,a,r,b,D,freq);
      scalar1 = 0.0D0;
      for i=1:3
        scalar1 = scalar1 + grad1(i)*grad1(i);
      end
      scalar1 = 1.0D0 / sqrt(scalar1);
      for i=1:3
        h_vec(i) = grad1(i) * scalar1;
      end

      alpha = 1.0D0;

      for i_count=1:len
%
      x(1) = x_vec(1) + alpha*h_vec(1);
      y(1) = x_vec(2) + alpha*h_vec(2);
      z(1) = x_vec(3) + alpha*h_vec(3);
      [x,y,z]=modeuler(tstep,h,0,z,y,x,a,r,b);
      [f2,grad2] = calcfg(tstep,h,x,y,z,datx,daty,datz,a,r,b,D,freq);
%
      scalar1 = 0.0D0;
      for i=1:3
        scalar1 = scalar1 + h_vec(i) * grad1(i);
      end
      scalar1 = (f2 - f1)/ (alpha*scalar1);
      scalar2 = abs(scalar1 - 1.0D0);
%
      alpha_vec(i_count) = alpha;
      fn_vec(i_count) = scalar1;
      resid_vec(i_count) = scalar2;
      alpha = alpha * 1.0D-1;

      end % i_count

% Produce plots
h=semilogx(alpha_vec,fn_vec);
title('Verification of gradient calculation')
xlabel('\alpha')
ylabel('\Phi(\alpha)')
%
figure(2)
loglog(alpha_vec,resid_vec)
title('Variation of residual')
xlabel('\alpha')
ylabel('\Phi(\alpha)-1')