function varargout=transl(f,Te,r,dr,E,v,g)
% [k12,l12,E,v]=TRANSL(f,Te,r,dr,E,v,g)
%
% Calculates the coherence-square halfway point in function 
% of Forsyth's subsurface-to-surface loading ratio and the elastic
% thickness in km. Returns the wavenumber at which C2=0.5, in rad/km.
%
% INPUT:
%
% f         Amplitude loading ratio [dimensionless, not squared]
% Te        Elastic thickness [km], can be vector
% r         Density [kg/m3]
% dr        Density contrast [kg/m3]
% E         Young's modulus [Pa]
% v         Poisson's ratio [dimensionless]
% g         Gravitational acceleration [m/s2]
%
% OUTPUT:
%
% k12       Wavenumber at which C2=0.5, in rad/km
% l12       Wavelength at which C2=0.5, in km
% E         Young's modulus
% v         Poisson's ratio
%
% SEE ALSO: FORSYTH, COHERETE, COHEREF2
% 
% Last modified by fjsimons-at-alum.mit.edu, 03/05/2012

defval('f',0.5)
defval('Te',17)
defval('E',1.4e11);
defval('v',0.25);
defval('r',2670);
defval('dr',630);
defval('g',fralmanac('GravAcc'));

% Convert to standard units [m]
Te=Te*1000;

disp(sprintf('E= %5.3g; v= %5.3f',E,v))

% How to get to this point
% coh=[ '((1+D*k^4/dr/g)*dr^2+f^2*r^2*(1+D*k^4/r/g))^2/...
%        ((1+D*k^4/dr/g)^2*dr^2+f^2*r^2)/...
%        (dr^2+f^2*r^2*(1+D*k^4/r/g)^2)'];
% ss=[ 'solve(0.5=' coh ',k);'];
% maplesol=parse(maple(ss),',');
% k12=eval(maplesol(1,:))

D=(E*Te.^3)./(12*(1-v^2));

[FF,DD]=meshgrid(f,D);
k12=1./2.*2.^(3./4)./DD./FF.*...
    ((-FF.*dr-FF.*r+FF.^2.*r+dr+...
      (FF.^2.*dr.^2+4.*FF.^2.*dr.*r-2.*FF.^3.*dr.*r+...
       2.*FF.*dr.^2+FF.^2.*r.^2+2.*FF.^3.*r.^2-...
       2.*FF.*r.*dr+FF.^4.*r.^2+dr.^2).^(1./2))...
     .*g.*DD.^3.*FF.^3).^(1./4);

% Now do it exactly as it appears in the EPSL paper:
beta=(FF.^2.*dr.^2+4.*FF.^2.*dr.*r-2.*FF.^3.*dr.*r+...
       2.*FF.*dr.^2+FF.^2.*r.^2+2.*FF.^3.*r.^2-...
       2.*FF.*r.*dr+FF.^4.*r.^2+dr.^2);

% Next line has the typo
%k12=(2*g./DD./FF.*(-FF.*dr-FF.*r+FF.^2.*r+dr+beta.^(1/2))).^(1./4);

% And noting the typo, this should be corrected to
k=(1/2*g./DD./FF.*(-FF.*dr-FF.*r+FF.^2.*r+dr+beta.^(1/2))).^(1./4);

% Check the results
difer(k-k12,[],[],NaN)

% Also check by plugging it in
difer(((1+DD.*k.^4/dr/g)*dr^2+FF.^2*r^2.*(1+DD.*k.^4/r/g)).^2./...
       ((1+DD.*k.^4/dr/g).^2*dr^2+FF.^2*r^2)./...
        (dr^2+FF.^2*r^2.*(1+DD.*k.^4/r/g).^2)-1/2,[],[],NaN)

% Back into the right units
k12=k12*1000;
l12=2*pi./k12;

% Produce output
varns={k12,l12,E,v};
varargout=varns(1:nargout);