function x = meshgen_mmpde5xi(N,y,rho,xi,dt,tau)
%MESHGEN_MMPDE5xi generates a new mesh of N points using the MMPDE5xi algorithm.
%
% X = MESHGEN_MMPDE5XI(N,Y,RHO) generates a new mesh of N points using
% the MMPDE5 algorithm. Y is the old mesh of N points. RHO is the adaptation
% function defined on the old mesh Y. 
%
% Usage with optional arguments:
% X = MESHGEN_MMPDE5XI(N,Y,RHO,XI)
% X = MESHGEN_MMPDE5XI(N,Y,RHO,XI,DT)
% X = MESHGEN_MMPDE5XI(N,Y,RHO,XI,DT,TAU)
%
% The optional arguments are XI, DT and TAU. Their default values are
% XI(i) = (i-1)/(N-1), i = 1,...,N,
% DT = 0.1, and TAU = 1.
%
% XI specifies the computational mesh of N points.
%
% DT sets the pseudo-time step size.
%
% TAU sets the value of parameter tau.
%

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

% declare arrays

x = zeros(N,1);
xi_new = zeros(N,1);
A = sparse(N,N);

% Initialization

if nargin < 3
   error('The first three arguments must be given.')
end
if nargin < 4
   xi = linspace(0,1,N)';
end
if nargin < 5
   dt = 0.1;
end
if nargin < 6
   tau = 1.0;
end

% Assemble the coefficient matrix and the right-hand side term 

dt1 = dt/tau;
A(1,1) = 1; xi_new(1) = xi(1);
for i = 2:N-1
   dx = 2/(y(i+1)-y(i-1));
   dx1 = 1/(y(i)-y(i-1));
   dx2 = 1/(y(i+1)-y(i));
   A(i,i-1) = -2.0/(rho(i)+rho(i-1))*dx*dx1*dt1;
   A(i,i+1) = -2.0/(rho(i)+rho(i+1))*dx*dx2*dt1;
   A(i,i) = 1.0 - A(i,i-1) - A(i,i+1);
   xi_new(i) = xi(i);
end
A(N,N) = 1; xi_new(N) = xi(N);

% solve the system for the new mesh

xi_new = A\xi_new;
xi_new(1) = xi(1); xi_new(N) = xi(N);

% use linear interpolation to find new physical coordinates

for i=1:N
   if xi(i)<=xi_new(1)
      m = 1;
   elseif xi(i)>=xi_new(N-1)
      m=N-1;
   else
      while m<N
	     if (xi(i)>=xi_new(m)) && (xi(i)<xi_new(m+1))
		    break;
		 end
		 m = m + 1;
      end
   end
   temp = (xi(i)-xi_new(m))/(xi_new(m+1)-xi_new(m));
   x(i) = y(m) + temp*(y(m+1)-y(m));
end

% end of meshgen_mmpde5xi
