function Vq = interp3(varargin)
%MATLAB Code Generation Library Function

%   Limitations and Notes:
%   1. Xq, Yq, and Zq must be the same size. Use MESHGRID to evaluate on on
%      a grid.
%   2. For best performance supply X, Y, and Z as vectors.
%   3. The 'cubic' method requires a uniformly-spaced grid, otherwise an
%      error is issued. Use the 'spline' method explicitly for an
%      unevenly-spaced grid.
%   4. The 'spline' method is optimized for Xq, Yq, and Zq being outputs of
%      MESHGRID or for a small number of interpolation points relative to
%      the dimensions of V. Using the 'spline' method to interpolate over a
%      large set of scattered points may be inefficient.

%   Copyright 1984-2017 The MathWorks, Inc.
%#codegen

narginchk(1,9);
coder.internal.prefer_const(varargin{:});
if nargin == 9
    % interp3(X,Y,Z,V,Xq,Yq,Zq,method,extrap)
    X = varargin{1};
    Y = varargin{2};
    Z = varargin{3};
    V = varargin{4};
    Xq = varargin{5};
    Yq = varargin{6};
    Zq = varargin{7};
    METHOD = StringToMethodID(varargin{8});
    EXTRAP = false;
    interp3_validate(V,Xq,Yq,Zq,METHOD,varargin{9},X,Y,Z);
    extrapval = cast(varargin{9},'like',V);
    Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval, ...
        unmeshgrid(X,2),unmeshgrid(Y,1),unmeshgrid(Z,3));
elseif nargin == 8
    % interp3(X,Y,Z,V,Xq,Yq,Zq,method)
    X = varargin{1};
    Y = varargin{2};
    Z = varargin{3};
    V = varargin{4};
    Xq = varargin{5};
    Yq = varargin{6};
    Zq = varargin{7};
    METHOD = StringToMethodID(varargin{8});
    EXTRAP = default_extrap(METHOD);
    extrapval = cast(coder.internal.nan,'like',V);
    interp3_validate(V,Xq,Yq,Zq,METHOD,extrapval,X,Y,Z);
    Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval, ...
        unmeshgrid(X,2),unmeshgrid(Y,1),unmeshgrid(Z,3));
elseif nargin == 7
    % interp3(X,Y,Z,V,Xq,Yq,Zq)
    X = varargin{1};
    Y = varargin{2};
    Z = varargin{3};
    V = varargin{4};
    Xq = varargin{5};
    Yq = varargin{6};
    Zq = varargin{7};
    METHOD = LINEAR;
    EXTRAP = default_extrap(METHOD);
    extrapval = cast(coder.internal.nan,'like',V);
    interp3_validate(V,Xq,Yq,Zq,METHOD,extrapval,X,Y,Z);
    Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval, ...
        unmeshgrid(X,2),unmeshgrid(Y,1),unmeshgrid(Z,3));
elseif nargin == 6
    % interp3(V,Xq,Yq,Zq,method,extrap)
    V = varargin{1};
    Xq = varargin{2};
    Yq = varargin{3};
    Zq = varargin{4};
    METHOD = StringToMethodID(varargin{5});
    EXTRAP = false;
    interp3_validate(V,Xq,Yq,Zq,METHOD,varargin{6});
    extrapval = cast(varargin{6},'like',V);
    Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval);
elseif nargin == 5
    % interp3(V,Xq,Yq,Zq,method)
    V = varargin{1};
    Xq = varargin{2};
    Yq = varargin{3};
    Zq = varargin{4};
    METHOD = StringToMethodID(varargin{5});
    EXTRAP = default_extrap(METHOD);
    extrapval = cast(coder.internal.nan,'like',V);
    interp3_validate(V,Xq,Yq,Zq,METHOD,extrapval);
    Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval);
elseif nargin == 4
    if coder.internal.isCharOrScalarString(varargin{3})
        % interp3(V,N,method,extrap)
        V = varargin{1};
        N = varargin{2};
        METHOD = StringToMethodID(varargin{3});
        Vq = interp3_ntimes(V,N,METHOD);
    else
        % interp3(V,Xq,Yq,Zq)
        V = varargin{1};
        Xq = varargin{2};
        Yq = varargin{3};
        Zq = varargin{4};
        METHOD = LINEAR;
        EXTRAP = default_extrap(METHOD);
        extrapval = cast(coder.internal.nan,'like',V);
        interp3_validate(V,Xq,Yq,Zq,METHOD,extrapval);
        Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval);
    end
elseif nargin == 3
    if coder.internal.isCharOrScalarString(varargin{2})
        % interp3(V,method,extrap)
        V = varargin{1};
        N = coder.internal.indexInt(1);
        METHOD = StringToMethodID(varargin{2});
        Vq = interp3_ntimes(V,N,METHOD);
    else
        % interp3(V,N,method)
        V = varargin{1};
        N = varargin{2};
        METHOD = StringToMethodID(varargin{3});
        Vq = interp3_ntimes(V,N,METHOD);
    end
elseif nargin == 2
    if coder.internal.isCharOrScalarString(varargin{2})
        % interp3(V,method)
        V = varargin{1};
        N = coder.internal.indexInt(1);
        METHOD = StringToMethodID(varargin{2});
        Vq = interp3_ntimes(V,N,METHOD);
    else
        % interp3(V,N)
        V = varargin{1};
        N = varargin{2};
        METHOD = LINEAR;
        Vq = interp3_ntimes(V,N,METHOD);
    end
else
    % interp3(V)
    V = varargin{1};
    N = coder.internal.indexInt(1);
    METHOD = LINEAR;
    Vq = interp3_ntimes(V,N,METHOD);
end

%--------------------------------------------------------------------------

function interp3_validate(V,Xq,Yq,Zq,METHOD,extrapval,X,Y,Z)
AUTOGRID = nargin <= 7;
eml_invariant(ndims(V) == 3, ...
    'Coder:toolbox:interp3_VMustBe3D');
eml_invariant(size(V,1) >= 2 && size(V,2) >= 2 && size(V,3) >= 2, ...
    'Coder:toolbox:NotEnoughPoints');
eml_invariant(isfloat(Xq) && isfloat(Yq) && isfloat(Zq), ...
    'Coder:toolbox:interp3_invalidXYZClass');
eml_invariant(isreal(Xq) && isreal(Yq) && isreal(Zq), ...
    'Coder:toolbox:interp3_nonrealXYZ');
eml_invariant(isa(Yq,class(Xq)) && isa(Zq,class(Xq)), ...
    'Coder:toolbox:interp3_classXIYIandZIMustMatch');
eml_invariant(isequal(size(Xq),size(Yq),size(Zq)), ...
    'Coder:toolbox:interp3_sizeXIYIandZIMustMatch');
if ~AUTOGRID
    eml_invariant(isfloat(X) && isfloat(Y) && isfloat(Z), ...
        'Coder:toolbox:interp3_invalidXYZClass');
    eml_invariant(isreal(X) && isreal(Y) && isreal(Z), ...
        'Coder:toolbox:interp3_nonrealXYZ');
    eml_invariant(isa(Y,class(X)) && isa(Z,class(X)), ...
        'Coder:toolbox:interp3_mixedTypeGrid');
    if eml_is_const(isvector(X)) && isvector(X) && ...
            eml_is_const(isvector(Y)) && isvector(Y) && ...
            eml_is_const(isvector(Z)) && isvector(Z)
        nx = numel(X);
        ny = numel(Y);
        nz = numel(Z);
    else
        eml_invariant(isplaid('meshgrid',X,Y,Z), ...
            'Coder:toolbox:GridMustBeVectorsOrPlaid');
        nx = size(X,2);
        ny = size(Y,1);
        nz = size(Z,3);
    end
    eml_invariant(eml_ndims(V) == 3 && ...
        size(V,1) == ny && ...
        size(V,2) == nx && ...
        size(V,3) == nz, ...
        'Coder:toolbox:interp3_gridToValuesSizeMismatch');
end
if METHOD ~= NEAREST
    eml_invariant(isfloat(V), ...
        'Coder:toolbox:interp3_VMustBeFloat');
    eml_invariant(coder.internal.isBuiltInNumeric(extrapval), ...
        'MATLAB:interp2:extrapvalNotNumeric');
end
eml_invariant(isreal(extrapval) || ~isreal(V), ...
    'Coder:toolbox:interp3_extrapvalNotReal');

%--------------------------------------------------------------------------

function Vq = interp3_ntimes(V,N,METHOD)
coder.internal.prefer_const(N,METHOD);
eml_invariant(isscalar(N) && isnumeric(N), ...
    'Coder:toolbox:NMustBeANumericScalar');
eml_invariant(ndims(V) == 3, ...
    'Coder:toolbox:interp3_VMustBe3D');
eml_invariant( ...
    size(V,2) >= 2 && size(V,1) >= 2 && size(V,3) >= 2, ...
    'Coder:toolbox:NotEnoughPoints');
eml_invariant(METHOD == NEAREST || isfloat(V), ...
    'Coder:toolbox:interp3_VMustBeFloat');
ny = size(V,1);
nx = size(V,2);
nz = size(V,3);
if isa(V,'single')
    d = 1/pow2(single(floor(N)));
else
    d = 1/pow2(double(floor(N)));
end
[Xq,Yq,Zq] = meshgrid(1:d:nx,1:d:ny,1:d:nz);
Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,false,cast(0,'like',V));

%--------------------------------------------------------------------------

function p = is_unit_colon(x)
% Determine if x is a *constant* colon expression of the form 1:numel(x).
coder.internal.prefer_const(x);
if eml_is_const(x)
    p = coder.const(feval('isequal',x,feval('colon',1,numel(x))));
else
    p = false;
end

%--------------------------------------------------------------------------

function Vq = interp3_dispatch(V,Xq,Yq,Zq,METHOD,EXTRAP,extrapval,X,Y,Z)
coder.internal.prefer_const(METHOD,EXTRAP);
if nargin <= 7
    AUTOGRID = true;
else
    coder.internal.prefer_const(X,Y,Z);
    if coder.const( ...
            is_unit_colon(X) && is_unit_colon(Y) && is_unit_colon(Z))
        AUTOGRID = true;
    else
        AUTOGRID = false;
    end
end
if METHOD == SPLINE
    % Since SPLINE does not really support an automatic grid, allow the
    % X, Y, and Z inputs to pass through, even if they are of the form 1:n.
    if AUTOGRID && nargin <= 7
        Vq = interp3_spline(V,Xq,Yq,Zq,EXTRAP,extrapval);
    else
        Vq = interp3_spline(V,Xq,Yq,Zq,EXTRAP,extrapval,X,Y,Z);
    end
elseif METHOD == CUBIC
    if AUTOGRID
        Vq = interp3_cubic(V,Xq,Yq,Zq,EXTRAP,extrapval);
    else
        Vq = interp3_cubic(V,Xq,Yq,Zq,EXTRAP,extrapval,X,Y,Z);
    end
else
    if AUTOGRID
        Vq = interp3_local(METHOD,V,Xq,Yq,Zq,EXTRAP,extrapval);
    else
        Vq = interp3_local(METHOD,V,Xq,Yq,Zq,EXTRAP,extrapval,X,Y,Z);
    end
end

%--------------------------------------------------------------------------

function Vq = interp3_local(method,V,Xq,Yq,Zq,EXTRAP,extrapval,X,Y,Z)
% The local interpolation methods:  linear and nearest.
% This function assumes numel(X)>=2, numel(Y)>=2, and numel(Z)>=2 and that
% the sizes are consistent with size(V).
coder.internal.prefer_const(method,EXTRAP);
AUTOGRID = nargin <= 7;
if method == LINEAR
    interpf = @scalar_trilinear_interp;
else
    interpf = @scalar_nearest_interp;
end
Vq = coder.nullcopy(zeros(size(Xq),'like',V));
n = coder.internal.indexInt(numel(Xq));
if coder.internal.useParforConst('interp3_orders_0_1',n)
    if AUTOGRID
        ixmax = coder.internal.indexInt(size(V,2) - 1);
        iymax = coder.internal.indexInt(size(V,1) - 1);
        izmax = coder.internal.indexInt(size(V,3) - 1);
        parfor k = 1:n
            Vq(k) = local_autogrid_loop_body(interpf, ...
                Xq(k),Yq(k),Zq(k),V,ixmax,iymax,izmax,EXTRAP,extrapval);
        end
    else
        xmin = X(1);
        ymin = Y(1);
        zmin = Z(1);
        xmax = X(end);
        ymax = Y(end);
        zmax = Z(end);
        parfor k = 1:n
            Vq(k) = local_loop_body(interpf,Xq(k),Yq(k),Zq(k),X,Y,Z,V, ...
                xmin,xmax,ymin,ymax,zmin,zmax,EXTRAP,extrapval);
        end
    end
else
    if AUTOGRID
        ixmax = coder.internal.indexInt(size(V,2) - 1);
        iymax = coder.internal.indexInt(size(V,1) - 1);
        izmax = coder.internal.indexInt(size(V,3) - 1);
        for k = 1:n
            Vq(k) = local_autogrid_loop_body(interpf, ...
                Xq(k),Yq(k),Zq(k),V,ixmax,iymax,izmax,EXTRAP,extrapval);
        end
    else
        xmin = X(1);
        ymin = Y(1);
        zmin = Z(1);
        xmax = X(end);
        ymax = Y(end);
        zmax = Z(end);
        for k = 1:n
            Vq(k) = local_loop_body(interpf,Xq(k),Yq(k),Zq(k),X,Y,Z,V, ...
                xmin,xmax,ymin,ymax,zmin,zmax,EXTRAP,extrapval);
        end
    end
end

%--------------------------------------------------------------------------

function Vqk = local_loop_body(interpf,Xqk,Yqk,Zqk,X,Y,Z,V, ...
    xmin,xmax,ymin,ymax,zmin,zmax,EXTRAP,extrapval)
coder.inline('always');
coder.internal.prefer_const(interpf,EXTRAP);
if EXTRAP || ...
        Xqk >= xmin && Xqk <= xmax && ...
        Yqk >= ymin && Yqk <= ymax && ...
        Zqk >= zmin && Zqk <= zmax
    ix = coder.internal.bsearch(X,Xqk);
    iy = coder.internal.bsearch(Y,Yqk);
    iz = coder.internal.bsearch(Z,Zqk);
    x1 = X(ix);
    x2 = X(ix + 1);
    y1 = Y(iy);
    y2 = Y(iy + 1);
    z1 = Z(iz);
    z2 = Z(iz + 1);
    v111 = V(iy,ix,iz);
    v211 = V(iy,ix + 1,iz);
    v121 = V(iy + 1,ix,iz);
    v221 = V(iy + 1,ix + 1,iz);
    v112 = V(iy,ix,iz + 1);
    v212 = V(iy,ix + 1,iz + 1);
    v122 = V(iy + 1,ix,iz + 1);
    v222 = V(iy + 1,ix + 1,iz + 1);
    Vqk = interpf(x1,x2,y1,y2,z1,z2, ...
        v111,v211,v121,v221,v112,v212,v122,v222, ...
        Xqk,Yqk,Zqk);
else
    Vqk = extrapval;
end

%--------------------------------------------------------------------------

function Vqk = local_autogrid_loop_body(interpf, ...
    Xqk,Yqk,Zqk,V,ixmax,iymax,izmax,EXTRAP,extrapval)
coder.inline('always');
coder.internal.prefer_const(interpf,ixmax,iymax,izmax,EXTRAP,extrapval);
if EXTRAP || ...
        Xqk >= 1 && Xqk <= size(V,2) && ...
        Yqk >= 1 && Yqk <= size(V,1) && ...
        Zqk >= 1 && Zqk <= size(V,3)
    ix = autogrid_lookup(Xqk,ixmax);
    iy = autogrid_lookup(Yqk,iymax);
    iz = autogrid_lookup(Zqk,izmax);
    x1 = cast(ix,'like',Xqk);
    x2 = x1 + 1;
    y1 = cast(iy,'like',Yqk);
    y2 = y1 + 1;
    z1 = cast(iz,'like',Zqk);
    z2 = z1 + 1;
    v111 = V(iy,ix,iz);
    v211 = V(iy,ix + 1,iz);
    v121 = V(iy + 1,ix,iz);
    v221 = V(iy + 1,ix + 1,iz);
    v112 = V(iy,ix,iz + 1);
    v212 = V(iy,ix + 1,iz + 1);
    v122 = V(iy + 1,ix,iz + 1);
    v222 = V(iy + 1,ix + 1,iz + 1);
    Vqk = interpf(x1,x2,y1,y2,z1,z2, ...
        v111,v211,v121,v221,v112,v212,v122,v222, ...
        Xqk,Yqk,Zqk);
else
    Vqk = extrapval;
end

%--------------------------------------------------------------------------

function idx = autogrid_lookup(x,maxidx)
% Substitute for bsearch when the grid is automatic.
coder.inline('always')
if x <= 1
    idx = coder.internal.indexInt(1);
elseif x <= maxidx
    idx = coder.internal.indexInt(floor(x));
else
    idx = coder.internal.indexInt(maxidx);
end

%--------------------------------------------------------------------------

function Vq = interp3_cubic(V,Xq,Yq,Zq,EXTRAP,extrapval,X,Y,Z)
coder.inline('always');
T = zeros('like',V);
RT = zeros('like',real(V));
AUTOGRID = nargin <= 6;
if AUTOGRID
    nx = coder.internal.indexInt(size(V,2));
    ny = coder.internal.indexInt(size(V,1));
    nz = coder.internal.indexInt(size(V,3));
    ixmax = nx - 1;
    iymax = ny - 1;
    izmax = nz - 1;
else
    dx = cast(EqualSpacing(X),'like',RT);
    dy = cast(EqualSpacing(Y),'like',RT);
    dz = cast(EqualSpacing(Z),'like',RT);
    eml_invariant(dx > 0 && dy > 0 && dz > 0, ...
        'Coder:toolbox:CubicGridMustBeUniform');
    nx = coder.internal.indexInt(numel(X));
    ny = coder.internal.indexInt(numel(Y));
    nz = coder.internal.indexInt(numel(Z));
    xmin = X(1);
    ymin = Y(1);
    zmin = Z(1);
    xmax = X(end);
    ymax = Y(end);
    zmax = Z(end);
end
eml_invariant(nx >= 3 && ny >= 3 && nz >= 3, ...
    'MATLAB:interp3:cubic:Vsize')
Vq = coder.nullcopy(zeros(size(Xq),'like',T));
% Expand V so interpolation is valid at the boundaries.
VV = zeros(size(V)+2,'like',T);
for iz = 1:nz
    for ix = 1:nx
        for iy = 1:ny
            VV(iy+1,ix+1,iz+1) = V(iy,ix,iz);
        end
    end
end
for iz = 1:nz+2
    for ix = 1:nx+2
        VV(1,ix,iz) = 3*VV(2,ix,iz) - 3*VV(3,ix,iz) + VV(4,ix,iz); % Y edges
        VV(ny+2,ix,iz) = 3*VV(ny+1,ix,iz) - 3*VV(ny,ix,iz) + VV(ny-1,ix,iz);
    end
end
for iz = 1:nz+2
    for iy = 1:ny+2
        VV(iy,1,iz) = 3*VV(iy,2,iz) - 3*VV(iy,3,iz) + VV(iy,4,iz); % X edges
        VV(iy,nx+2,iz) = 3*VV(iy,nx+1,iz) - 3*VV(iy,nx,iz) + VV(iy,nx-1,iz);
    end
end
for ix = 1:nx+2
    for iy = 1:ny+2
        VV(iy,ix,1) = 3*VV(iy,ix,2) - 3*VV(iy,ix,3) + VV(iy,ix,4); % Z edges
        VV(iy,ix,nz+2) = 3*VV(iy,ix,nz+1) - 3*VV(iy,ix,nz) + VV(iy,ix,nz-1);
    end
end
nout = coder.internal.indexInt(numel(Xq));
if coder.internal.useParforConst('interp3_orders_2_3',nout)
    if AUTOGRID
        parfor k = 1:nout
            Vq(k) = cubic_autogrid_loop_body(Xq(k),Yq(k),Zq(k),VV, ...
                ixmax,iymax,izmax,EXTRAP,extrapval);
        end
    else
        parfor k = 1:nout
            Vq(k) = cubic_loop_body(Xq(k),Yq(k),Zq(k),X,Y,Z,VV, ...
                xmin,xmax,ymin,ymax,zmin,zmax,dx,dy,dz,EXTRAP,extrapval);
        end
    end
else
    if AUTOGRID
        for k = 1:nout
            Vq(k) = cubic_autogrid_loop_body(Xq(k),Yq(k),Zq(k),VV, ...
                ixmax,iymax,izmax,EXTRAP,extrapval);
        end
    else
        for k = 1:nout
            Vq(k) = cubic_loop_body(Xq(k),Yq(k),Zq(k),X,Y,Z,VV, ...
                xmin,xmax,ymin,ymax,zmin,zmax,dx,dy,dz,EXTRAP,extrapval);
        end
    end
end

function X = localevaluate(x,iter)
coder.inline('always');
coder.internal.prefer_const(iter);
if iter == 0
    X = ((2-x).*x-1).*x;
elseif iter == 1
    X = (3*x-5).*x.*x+2;
elseif iter == 2
    X = ((4-3*x).*x+1).*x;
else
    X = (x-1).*x.*x;
end

%--------------------------------------------------------------------------

function Vqk = cubic_loop_body(Xqk,Yqk,Zqk,X,Y,Z,VV, ...
    xmin,xmax,ymin,ymax,zmin,zmax,dx,dy,dz,EXTRAP,extrapval)
coder.inline('always');
coder.internal.prefer_const(EXTRAP);
realtype = zeros('like',real(VV));
if EXTRAP || ( ...
        Xqk >= xmin && Xqk <= xmax && ...
        Yqk >= ymin && Yqk <= ymax && ...
        Zqk >= zmin && Zqk <= zmax)
    ix = coder.internal.bsearch(X,Xqk);
    iy = coder.internal.bsearch(Y,Yqk);
    iz = coder.internal.bsearch(Z,Zqk);
    s = (cast(Xqk,'like',realtype) - cast(X(ix),'like',realtype))/dx;
    t = (cast(Yqk,'like',realtype) - cast(Y(iy),'like',realtype))/dy;
    w = (cast(Zqk,'like',realtype) - cast(Z(iz),'like',realtype))/dz;
    vik = zeros('like',VV);
    for iw = coder.unroll(0:3)
        ww = localevaluate(w,iw);
        for is = coder.unroll(0:3)
            ss = localevaluate(s,is);
            for it = coder.unroll(0:3)
                tt = localevaluate(t,it);
                vik = vik + VV(iy+it,ix+is,iz+iw)*ss*tt*ww;
            end
        end
    end
    Vqk = vik/8;
else
    Vqk = extrapval;
end

%--------------------------------------------------------------------------

function Vqk = cubic_autogrid_loop_body(Xqk,Yqk,Zqk,VV, ...
    ixmax,iymax,izmax,EXTRAP,extrapval)
coder.inline('always');
coder.internal.prefer_const(ixmax,iymax,izmax,EXTRAP);
realtype = zeros('like',real(VV));
if EXTRAP || ( ...
        Xqk >= 1 && Xqk <= (ixmax + 1) && ...
        Yqk >= 1 && Yqk <= (iymax + 1) && ...
        Zqk >= 1 && Zqk <= (izmax + 1))
    ix = autogrid_lookup(Xqk,ixmax);
    iy = autogrid_lookup(Yqk,iymax);
    iz = autogrid_lookup(Zqk,izmax);
    s = cast(Xqk,'like',realtype) - cast(ix,'like',realtype);
    t = cast(Yqk,'like',realtype) - cast(iy,'like',realtype);
    w = cast(Zqk,'like',realtype) - cast(iz,'like',realtype);
    vik = zeros('like',VV);
    for iw = coder.unroll(0:3)
        ww = localevaluate(w,iw);
        for is = coder.unroll(0:3)
            ss = localevaluate(s,is);
            for it = coder.unroll(0:3)
                tt = localevaluate(t,it);
                vik = vik + VV(iy+it,ix+is,iz+iw)*ss*tt*ww;
            end
        end
    end
    Vqk = vik/8;
else
    Vqk = extrapval;
end

%--------------------------------------------------------------------------

function vi = scalar_trilinear_interp(x1,x2,y1,y2,z1,z2, ...
    v111,v211,v121,v221,v112,v212,v122,v222, ...
    xi,yi,zi)
coder.inline('always');
T = zeros('like',v111);
if xi == x1
    v11 = cast(v111,'like',T);
    v21 = cast(v121,'like',T);
    v12 = cast(v112,'like',T);
    v22 = cast(v122,'like',T);
elseif xi == x2
    v11 = cast(v211,'like',T);
    v21 = cast(v221,'like',T);
    v12 = cast(v212,'like',T);
    v22 = cast(v222,'like',T);
else
    r = (xi - x1)/(x2 - x1);
    onemr = 1 - r;
    if v111 == v211
        v11 = cast(v111,'like',T);
    else
        v11 = cast(onemr*v111 + r*v211,'like',T);
    end
    if v121 == v221
        v21 = cast(v121,'like',T);
    else
        v21 = cast(onemr*v121 + r*v221,'like',T);
    end
    if v112 == v212
        v12 = cast(v112,'like',T);
    else
        v12 = cast(onemr*v112 + r*v212,'like',T);
    end
    if v122 == v222
        v22 = cast(v122,'like',T);
    else
        v22 = cast(onemr*v122 + r*v222,'like',T);
    end
end
if yi == y1
    v1 = v11;
    v2 = v12;
elseif yi == y2
    v1 = v21;
    v2 = v22;
else
    r = (yi - y1)/(y2 - y1);
    onemr = 1 - r;
    if v11 == v21
        v1 = v11;
    else
        v1 = cast(onemr*v11 + r*v21,'like',T);
    end
    if v12 == v22
        v2 = v12;
    else
        v2 = cast(onemr*v12 + r*v22,'like',T);
    end
end
if zi == z1 || v1 == v2
    vi = v1;
elseif zi == z2
    vi = v2;
else
    r = (zi - z1)/(z2 - z1);
    vi = cast((1 - r)*v1 + r*v2,'like',T);
end

%--------------------------------------------------------------------------

function zi = scalar_nearest_interp(x1,x2,y1,y2,z1,z2, ...
    v111,v211,v121,v221,v112,v212,v122,v222, ...
    xi,yi,zi)
coder.inline('always');
dx1 = xi - x1;
dx2 = x2 - xi;
dy1 = yi - y1;
dy2 = y2 - yi;
dz1 = zi - z1;
dz2 = z2 - zi;
if dx1 < dx2
    if dy1 < dy2
        if dz1 < dz2
            zi = v111;
        else
            zi = v112;
        end
    else
        if dz1 < dz2
            zi = v121;
        else
            zi = v122;
        end
    end
else
    if dy1 < dy2
        if dz1 < dz2
            zi = v211;
        else
            zi = v212;
        end
    else
        if dz1 < dz2
            zi = v221;
        else
            zi = v222;
        end
    end
end

%--------------------------------------------------------------------------

function Vq = interp3_spline(V,Xq,Yq,Zq,EXTRAP,extrapval,X,Y,Z)
%3-D spline interpolation
if nargin < 7
    X = ones('like',real(V)):size(V,2);
    Y = ones('like',real(V)):size(V,1);
    Z = ones('like',real(V)):size(V,3);
end
Vq = TensorInterp23(@spline,Y,X,Z,V,Yq,Xq,Zq);
if ~EXTRAP
    Vq = MaskOutOfRangeValues(extrapval,Vq,X,Y,Z,Xq,Yq,Zq);
end

%--------------------------------------------------------------------------

function m = LINEAR
coder.inline('always');
m = uint8(0);

function m = NEAREST
coder.inline('always');
m = uint8(1);

function m = CUBIC
coder.inline('always');
m = uint8(3);

function m = SPLINE
coder.inline('always');
m = uint8(4);

%--------------------------------------------------------------------------

function m = StringToMethodID(method)
coder.internal.prefer_const(method);
eml_invariant(eml_is_const(method), ...
    'Coder:toolbox:MethodMustBeConstant');
n = coder.const(eml_max(1,strlength(method)));
islinear = coder.const(strncmpi(method,'linear',n));
isnearest = coder.const(strncmpi(method,'nearest',n));
iscubic = coder.const(strncmpi(method,'cubic',n));
isspline = coder.const(strncmpi(method,'spline',n));
eml_invariant(islinear || isnearest || iscubic || isspline, ...
    'Coder:toolbox:interp3_unsupportedMethod');
if isnearest
    m = NEAREST;
elseif iscubic
    m = CUBIC;
elseif isspline
    m = SPLINE;
else
    m = LINEAR;
end

%--------------------------------------------------------------------------

function p = default_extrap(METHOD)
if METHOD == SPLINE
    p = true;
else
    p = false;
end

%--------------------------------------------------------------------------

