function [c,dcdz]=veldep(z)
%This function gives the velocity and its depth
%derivative as function of depth. The underlying
%velocity model is specified with the global
%variables Z, C and dCdZ which are set in the
%function profile.m. Values outside the nodepoints range
%get the velocity of the nearest nodepoint, except
%for depths greater than the deepest nodepoint where
%the last gradient is positive, then the gradient
%is kept constant.
%Written by Keli Karason, karason@mit.edu

global C Z dCdZ

c=interp1(Z,C,z);
c(z<Z(1))=C(1);
j=find(z>Z(end));
c(j)=C(end)+max(0,(z(j)-Z(end))*dCdZ(end));


nn=length(Z)-1;

if nn>1,
  z(z>Z(end) & dCdZ(end)<0)=-1;
  z(z>Z(end-1))=Z(end-1);
  n=floor(interp1(Z(1:end-1),(1:nn)',z));
  
  j=find(isnan(n));
  n(j)=1;
  dcdz=dCdZ(n);
  dcdz(j)=0;
else
  dcdz=repmat(dCdZ,size(z));
  dcdz(z<Z(1))=0;
  dcdz(z>Z(end) & dcdz<0)=0;
end




 

