function [z,dz]=scoord_rps(s,Cs,hc,zeta,h); 
% SCOORD_RPS returns 2D z values for 1d arrays of h and zeta
% Usage: z=scoord_rps(s,Cs,hc,zeta,h); 
% take 1D vectors of zeta and h and returns z and dz values
% arrange so we can matrix multiply
Cs=Cs(:).';
s=s(:).';
onez=ones(size(h(:)));
z = zeta(:)*(1.0+s) + onez*hc*s + (h(:)-hc)*Cs;
dz=diff(z,1,2);  %1st order difference along dimension 2