function [deter,G,GBP] = sismo_Stoneley_det(k,sigma,rhow,rhos,alphaw,alphas,h,betas)
% det is the determinant of the 3 by 3 matrix that gives the amplitudes of the 
% water and crust potentials 
% for a one layer of water over an infinite half-space (crust)
% and G in the transfer function from surface pressure to crust potential amplitude 
% and GP is the TF from surface to bottom pressure
%
% See Ardhuin and Herbers, JFM 2013

l=sqrt(sigma.^2./alphaw.^2-k.^2);
m=sqrt(k.^2-sigma.^2./alphas.^2);
n=sqrt(k.^2-sigma.^2./betas.^2);
fac1=(sigma^2-2*k.^2.*betas^2);
%fac2=4*betas^4.*k.^2.*m.*n-(sigma^2-2*k^2*betas^2)^2;
fac2=rhos.*(4.*betas.^4.*k.^2.*(m.*n)./fac1-fac1);
left=2.*(l.*cos(l.*h).*fac2);
right=2.*rhow.*m.*sin(l.*h).*sigma.^4./fac1; 
deter=(left-right);  % = DET/rhow  with DET given by AH2013 eq. (4.18)
G=2.*l.*m.*(k.^2-n.^2)./((k.^2+n.^2).*deter);

r=-1i*fac2./sigma;
q=m*sigma^2/fac1;


GC=-(l.*r+q.*sigma.*rhow)./(deter*rhow);
GD=-(l.*r-q.*sigma.*rhow)./(deter*rhow);
GBP=(GC+GD)*1i*rhow*sigma;   % p(z=-h)= i rhow * sigma *(C+D)
end

% delta to P = 
% Matrix form:  rhs = [B C A]
%
%r=-1i*fac2./sigma;
%q=m*sigma^2/fac1;
%M3 = [exp(1i*l*h)       exp(-1i*l*h)      0 ;
%      -1i*l             1i*l              q ;
%      1i*rhow*sigma     1i*rhow*sigma    -r;];
%rhs = [1i/rhow/sigma;0;0];
%res3=M3\rhs;