%% pwint2.m is my second attempt to produce
%% my own version of MATLAB's interp2 command
%% (see Philint2.m for my first go a year and a half ago).
%% The reason for the new code is that I have learnt
%% a few new matrix index tricks which should speed
%% up the code. Also I have read the relativent section
%% in Numerical Recipes (3.6 Interpolation in Two or
%% More Dimensions) which gives the simple formulas
%% for bilinear interpolation. 
%% The output from pwint2.m is the same size as the input xi
%% matrix.
%%

function zi=pwint2(xv,yv,zmat,xi,yi)

% xv = x-coordinates vector
% yv = y-coordinates vector
% zmat = matrix of z values (size(yv,xv))
% xi = x-coordinates matrix points of interest
% yi = y-coordinates matrix points of interest
% zi = calculated z values 

[szy,szx]=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

zaa=zmat(yaic+(xaic-1)*szy);                    %1st required z values
zab=zmat(ybic+(xaic-1)*szy);                    %2nd required z values
zba=zmat(yaic+(xbic-1)*szy);                    %3rd required z values
zbb=zmat(ybic+(xbic-1)*szy);                    %4th required z values

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 z using bilinear interpolation
zi=(1-t).*(1-u).*zaa+t.*(1-u).*zba+(1-t).*u.*zab+t.*u.*zbb;

zi=reshape(zi,sxi);                             %reshape zi to same size as input xi

