function measure = qmesh(N,x,rho,k,m)
%QMESH computes the mesh quality measure.
%
% MEASURE = QMESH(N,X,RHO) computes the maximum norm of the
% equidistribution quality measure for mesh X of N points.
% RHO is the adaptation function defined on X.
%
% Usage with optional arguments:
% MEASURE = QMESH(N,X,RHO,K,M)
%
% The optional arguments are K and M. It is required that M <= K.
% When these optional arguments are specified, the overall
% mesh quality in the weighted L^2 norm is computed. 
%
% K is the degree of interpolation polynomials and M has the meaning that
% the interpolation error is measured in the semi-H^M norm.
%
% Example 1: (K,M) = (1,0) for linear interpolation with L^2 norm
% being used for measuring the error.
%
% Example 2: (K,M) = (1,1) for linear interpolation with H^1 seminorm
% being used for measuring the error.
%

%
% 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.
%

% Initialization

if nargin < 3
   error('The first three arguments must be given.')
end
sigma = 0;
for i=2:N
   sigma = sigma + (x(i)-x(i-1))*(rho(i)+rho(i-1));
end
sigma = 0.5*sigma;

% Compute the quality measure in the maximum norm

measure = 0;
if nargin == 3
   for i=2:N
      qadp = (N-1)*(x(i)-x(i-1))*0.5*(rho(i)+rho(i-1))/sigma;
      if qadp > measure
         measure = qadp;
      end
   end
   return;
end

% Compute the quality measure in the L^2 norm

if nargin < 5
   error('All of the arguments must be given.')
end
if k < m
   error('m must be less than or equal to k.')
end

temp = 0;
temp1 = 0;
for i=2:N
   qadp = (N-1)*(x(i)-x(i-1))*0.5*(rho(i)+rho(i-1))/sigma;
   temp = temp + (x(i)-x(i-1))*0.5*(rho(i)+rho(i-1))*qadp^(2*(k-m+1));
   temp1 = temp1 + (x(i)-x(i-1))*0.5*(rho(i)+rho(i-1));
end
measure = sqrt(temp/temp1);
% end of qmesh
