%% pwint2m.m is modified version of pwint2.m. It is designed to 
%% interpolated within one or more two dimensional arrays of the
%% same size at the same points of interest. The user inputs the 
%% required two dimensional arrays in the form of a single three
%% dimensional array which has dimensions of length(yv) by 
%% length(xv) by n, were n is the number of 2-D arrays. The output
%% is in the form of a single 3-D array with dimensions of size(xi)
%% by n.

function zres=pwint2m(xv,yv,zmat,xi,yi)

% xv = x-coordinates vector
% yv = y-coordinates vector
% zmat = matrix of z values (size(yv,xv) by n were n equals number of arrays to lookup)
% xi = x-coordinates matrix points of interest
% yi = y-coordinates matrix points of interest
% zres = calculated z values (size(xi) by n)

[szy,szx,szz]=size(zmat);                       %size of zmat
sxi=size(xi);                                   %size of xi matrix

xv=xv(:);                                       %make xv is a column vector
yv=yv(:);                                       %make yv is a column vector
xi=xi(:)';                                      %make xi is a row vector
yi=yi(:)';                                      %make yi is a row vector

xil=length(xi);                                 %length of xi 

xvm=xv*ones(1,xil);                             %required xv matrix
yvm=yv*ones(1,xil);                             %required yv matrix
xim=ones(szx,1)*xi;                             %required xi matrix
yim=ones(szy,1)*yi;                             %required yi matrix

xaic=sum(xvm<=xim);                             %row index of required x values below or equal to xi
xbic=((szx+1)-sum(xvm>=xim));                   %row index of required x values above or equal to xi
yaic=sum(yvm<=yim);                             %row index of required y values below or equal to xi
ybic=((szy+1)-sum(yvm>=yim));                   %row index of required y values above or equal to xi

xai=xaic+szx*(0:(xil-1));                       %index of required x values below or equal to xi
xbi=xbic+szx*(0:(xil-1));                       %index of required x values above or equal to xi
yai=yaic+szy*(0:(xil-1));                       %index of required y values below or equal to xi
ybi=ybic+szy*(0:(xil-1));                       %index of required y values above or equal to xi

xa=xvm(xai);                                    %required x values below or equal to xi
xb=xvm(xbi);                                    %required x values above or equal to xi
ya=yvm(yai);                                    %required y values below or equal to yi
yb=yvm(ybi);                                    %required y values above or equal to yi

t=zeros(size(xi));                              %t parameter
u=zeros(size(xi));                              %u parameter

t(xb~=xa)=(xi(xb~=xa)-xa(xb~=xa))./(xb(xb~=xa)-xa(xb~=xa));     %correct t parameter values 
u(yb~=ya)=(yi(yb~=ya)-ya(yb~=ya))./(yb(yb~=ya)-ya(yb~=ya));     %correct u parameter values

%Calculate the zmat values
zres=zeros(sxi(1),sxi(2),szz);                                  %results matrix
for i=1:szz
 zmati=zmat(:,:,i);                                             %i'th zmat 2-D array
 zaa=zmati(yaic+(xaic-1)*szy);                                  %1st required z values
 zab=zmati(ybic+(xaic-1)*szy);                                  %2nd required z values
 zba=zmati(yaic+(xbic-1)*szy);                                  %3rd required z values
 zbb=zmati(ybic+(xbic-1)*szy);                                  %4th required z values
 zi=(1-t).*(1-u).*zaa+t.*(1-u).*zba+(1-t).*u.*zab+t.*u.*zbb;    %calculate z using bilinear interpolation
 zi=reshape(zi,sxi);                                            %reshape zi to same size as input xi
 zres(:,:,i)=zi;                                                %record i'th array results
end 
