function ncomp = calcNCompMat(depth,freq,model)
%CALCNCOMPMAT Calculate normalized compliance of a layered earth model.
%	Uses M-file for compliance calculation (slower, but universal)
%	ncomp = normalized vertical compliance
%	ncomp = calcNCompMat(wdepth,freqs,model);
%
%	model is a 1-D earth model, with:
%				layer thicknesses (meters) in column 1
%				layer densities (g/cc) in column 2
%				layer compress vels (km/s) in column 3
%				layer shear vels (km/s) in column 4
%	wdepth is the water depth (meters) at the site
%	freqs are the frequencies (Hz) to calculate compliance at

thick=model(:,1);
rho=model(:,2)*1000;
vpsq=1000000*model(:,3).*model(:,3);
vssq=1000000*model(:,4).*model(:,4);
omega = 2*pi*freq;
k = gravd(omega,depth);
p = k ./ omega;
for i = 1:length(p),
	[v u sigzz, sigzx] = raydep(p(i),omega(i),thick,rho,vpsq,vssq);
	ncomp(i) = -k(i)*v(1)/(omega(i)*sigzz(1));
end
