function err = interpolation_error(FnctV,jmax,npde,t,x,v,p,k,m,FnctDV)
%
% err = interpolation_error(FnctV,jmax,npde,t,x,v,p,k,m,FnctDV)
%
% the function computes the W^(m,p) (0 <= m <= k) seminorm of
% constant (k=0) or linear (k=1) interpolation error.
%
% The first 6 arguments must be specified by the user. The other 4
% arguments (p,k,m,FnctDV) are optional. Particularly, FnctDV is
% used only for the case k=1 and m=1.
%
% variables:
%
% FnctV:    the exact solution provided by the user. Its syntax is
%
%                  y = FnctV(t,x)
%
%           Its output should be a column vector of size npde.
% npde:     the number of solution components.
% jmax:     the number of mesh points.
% t:        current time instant.
% x:        a column vector of size jmax containing the mesh points.
% v:        a array of size npde by jmax, v(i,j): i=1:npde, j=1:jmax,
%           containing nodal values of the solution.
% p:        p (0 < p <= inf) indicates L^p is used for computing the error.
%           default value: p = 2.
% k:        k=0 (constant interppolation) and k=1 (linear interpolation).
%           default value: k = 1.
% m:        indicates if the error or its derivative is comouted.
%           m=0 if k=0
%           m=0 or m=1 if k=1.
%           default value: m = 0.
% FnctDV:   its syntax is
%
%                  dy = FnctDV(t,x)
%
%           Its output should be a column vector of size npde. This function is needed
%           only for the case m = 1.
%

%
% Copyright (C) 2010 Weizhang Huang and Robert D. Russell
% all rights reserved.
%
% This program is provided "as is", without warranty of any kind.
% Permission is hereby granted, free of charge, to use this program
% for personal, research, and education purposes. Distribution or use
% of this program for any commercial purpose is permissible
% only by direct arrangement with the copyright owner.
%

   % set the default values for the parameters

   if nargin < 6
       error('The first 6 arguments must be defined.');
   elseif nargin < 7
      p=[];
	  k=[];
	  m=[];
	  FnctDV=[];
   elseif nargin < 8
	  k=[];
	  m=[];
	  FnctDV=[];
   elseif nargin < 9
	  m=[];
	  FnctDV=[];
   elseif nargin < 10
      FnctDV=[];
   elseif nargin > 10
      error('Two many arguments are specified.');
   end    

   if isempty(p) 
      p = 2;
   end
   if isempty(k)
      k = 1;
   end
   if isempty(m)
      m = 0;
   end

   % check parameters
   
   if p<eps
      error('p must be: p > 0.');
   end
   if (k~=0) && (k~=1)
      error('k must be 0 or 1.');
   end
   if (k==0)
      if (m~=0)
         error('m must be: 0 <= m <= k.');
      end
   end
   if (k==1)
      if (m~=0) && (m~=1)
         error('m must be: 0 <= m <= k.');
      end
   end
   if (m==1)
      if isempty(FnctDV)
	     error('The user must specify FnctDV for the case m = 1.');
      end
   end

% define weights for composite simpson's quadrature rule

npts1 = 20;
npts = 2*npts1+1;
w = zeros(1,npts);
for i=1:npts1, w(2*i)=4.0; end
for i=1:(npts1+1), w(2*i-1)=2.0; end
w(1)=1.0; w(npts)=1.0;
w=w/(3*(npts-1));

% declare arrays

vexact = zeros(npde,1);
vcompt = zeros(npde,1);
temp = zeros(npde,1);

% compute the L^p norm of the constant interpolation error

if m == 0 && k == 0
   err = 0;
   for i=2:jmax
      xxx = linspace(x(i-1),x(i),npts);
	  h = x(i)-x(i-1);
      for j=1:npts
         xx = xxx(j);
         vexact = feval(FnctV, t, xx);
         vcompt = (v(:,i)+v(:,i-1))*0.5;
         temp = abs(vexact-vcompt);
         if p == inf
            if max(temp) > err, err = max(temp); end
         else
            err = err + w(j)*h*norm(temp,p)^p;
         end
      end
   end
   if p ~= inf, err = (err/(x(jmax)-x(1)))^(1.0/p); end
end


% compute the L^p norm of the linear interpolation error

if m == 0 && k == 1
   err = 0;
   for i=2:jmax
      xxx = linspace(x(i-1),x(i),npts);
      h = x(i)-x(i-1);
      for j=1:npts
         xx = xxx(j);
         vexact = feval(FnctV, t, xx);
         vcompt = ((xx-x(i-1))*v(:,i)+(x(i)-xx)*v(:,i-1))/h;
         temp = abs(vexact-vcompt);
         if p == inf
            if max(temp) > err, err = max(temp); end
         else
            err = err + w(j)*h*norm(temp,p)^p;
         end
      end
   end
   if p ~= inf, err = (err/(x(jmax)-x(1)))^(1.0/p); end
end
      
% compute the L^p norm of the gradient of the linear interpolation error

if m == 1 && k == 1
   err = 0;
   for i=2:jmax
      xxx = linspace(x(i-1),x(i),npts);
      h = x(i)-x(i-1);
      for j=1:npts
         xx = xxx(j);
         vexact = feval(FnctDV, t, xx);
         vcompt = (v(:,i)-v(:,i-1))/h;
         temp = abs(vexact-vcompt);
         if p == inf
            if max(temp) > err, err = max(temp); end
         else
			err = err + w(j)*h*norm(temp,p)^p;
         end
      end
   end
   if p ~= inf, err = (err/(x(jmax)-x(1)))^(1.0/p); end
end

% end of interpolation_error
