function [N, Mu, dN, dMu, dL, D, dD] = ieshapefun(ietype, xi, P, ptype, a, b)
%IESHAPEFUN Calculation of shape function values for infinite elements
% Inputs:
%   ietype - Type of infinite element, valid values:
%       122 - 2D infinite element with 2 base points (iequad)
%       133 - 3D infinite element with 3 base points (iepenta)
%       134 - 3D infinite element with 4 base points (iehexa)
%   xi - Quadrature points in (s,t,u) coordinate system (dim x Nxi matrix)
%   P  - polynomial order (scalar)
%   ptype - Polynomial type, valid values
%       1, 'l', 'lag', or 'lagrange' - Lagrange polynomials
%       2, 'j', 'jac', or 'jacobi'   - Jacobi polynomials
%       3, 'leg' or 'legendre' - Legendre polynomials
%   a - alpha parameter of Jacobi polynomials (unused in Lagrange case) 
%   b - beta parameter of Jacobi polynomials (unused in Lagrange case)
%
% Outputs:
%   N - pressure shape function values (Nxi x dim P matrix)
%   Mu - phase function values (Nxi x dim P x Nf matrix)
%   dN - derivative of pressure shape functions (dim Nxi x dim P matrix)
%
% See also:
%   ielem_mkc, ieshape22m, ieshape33m, ieshape34m, jacpoly, lagpolys, 
%   model_mk

% Peter Rucz
% 2009-2015

% Check argument number
narginchk(4,6);
nargoutchk(1,7);

% Dimensionality
ndim = floor(mod(ietype, 100) / 10);

% Switch on polynomial type
switch ptype
    % Lagrange polynomials and their derivatives
    case {1, 'l','lag','lagrange'}
        z = linspace(-1, 1-2/P, P);         % Create vector of zeroes
        Lag = lagpolys(z);                  % Create Lagrange polynomials
        dLag = cell(P,1);                   % Preallocate derivatives
        % Scale the shape functions
        switch ndim
            case 2 % Two-dimensional elements (scale by sqrt((1-z)/2))
                for c1 = 1 : P
                    Lag{c1} = Lag{c1}/polyval(Lag{c1},z(c1))/sqrt((1-z(c1))/2);
                    dLag{c1} = Lag{c1}(1:end-1).*(length(Lag{c1})-1:-1:1);
                end
            case {1, 3} % Three-dimensional elements (scale by (1-z)/2)
                for c1 = 1 : P
                    Lag{c1} = Lag{c1}/polyval(Lag{c1},z(c1))/((1-z(c1))/2);
                    dLag{c1} = Lag{c1}(1:end-1).*(length(Lag{c1})-1:-1:1);
                end
            otherwise
                error('ieshapefun:ndim', 'Invalid dimensionality: %d', ndim);
        end
        % Fill arguments to be passed
        arg1 = Lag;
        arg2 = dLag;
    % Jacobi polynomials
    case {2, 'j','jac','jacobi'}
        % Fill arguments to be passed
        arg1 = a;
        arg2 = b;
    case {3, 'leg', 'legendre'}
        x = xi(:, 1);
        Nx = numel(x);
        LegValues = zeros(Nx, P);
        dLegValues = zeros(Nx, P);
        for iP = 1 : P
            [LegValues(:, iP), dLegValues(:, iP)] = legpoly(iP-1, x);
        end
        arg1 = LegValues;
        arg2 = dLegValues;
    otherwise
        error('ieshapefun:invalidArg','Unknown polynomial family.');
end % of switch on ptype

% Switch ietype
switch ietype
    case 111 % 1D infinite element with 1 base point
        [N, Mu, dN, dMu, dL, D, dD] = ieshape11m(xi, P, ptype, arg1, arg2);
    case 122 % 2D infinite element with 2 base points (iequad)
        [N, Mu, dN, dMu, dL, D, dD] = ieshape22m(xi, P, ptype, arg1, arg2);
    case 133 % 3D infinite element with 3 base points (iepenta)
        [N, Mu, dN, dMu, dL, D, dD] = ieshape33m(xi, P, ptype, arg1, arg2);
    case 134 % 3D infinite element with 4 base points (iehexa)
        [N, Mu, dN, dMu, dL, D, dD] = ieshape34m(xi, P, ptype, arg1, arg2);
end % of switch on ietype
        
end % of ieshapefun
