function varargout=forsyth(Te,lbd,f2,r,rc,drho,T,g,xver)
% [G2b,k,l,Zb,Zf]=FORSYTH(Te,lbd,f2,r,rc,drho,T,g,xver)
%
% Calculates the predicted coherence-square and admittance between
% Bouguer/free air anomalies and topography over a range of wavenumbers for a
% lithosphere loaded at the surface and at one other depth which is also
% the depth of compensation. For the coherence it doesn't matter where
% this depth is; it is implicitly included in the loading factor.
% See Forsyth (1985).
%
% Works for a vector of Te OR a vector of f2 OR a vector of r (only one vector).
%
% INPUT:
%
% Te      elastic thickness of plate [m]
% lbd     the wavelengths considered [m] [defaulted to reasonable values]
% f2      spectral ratio of bottom-to-top applied loads [dimensionless, squared]
%         0 for surface loading only (produces unity coherence)
%         1 for equal loads
%         infinity for Moho loading only (produces unity coherence)
% r       initial-load correlation coefficient
% rc      density contrast across topography interface [kg/m^3]
%         rho_crust for loading of continents
% drho    density contrast at compensation interface
%         (usually = rho_mantle - rho_crust) [kg/m^3]
% T       depth to the density contrast/compensation interface [m]
% g       Gravitational acceleration (in m/s^2) [defaulted]
%
% OUTPUT:
%
% G2b       Bouguer coherence-square (between 0 and 1)
% k         wavenumber [rad/m]
% l         wavelength [m]
% Zb        Bouguer addmittance, in mgal/m
% Zf        Free-air addmittance, in mgal/m
%
% See also MCKENZIE, ADMITTANCE, TRANSL
%
% EXAMPLE:
%
% forsyth('demo')
%
% Last modified by fjsimons-at-alum.mit.edu, 03/07/2012

defval('Te','demo')

if ~isstr(Te)
  defval('Te',[20 80]*1e3);
  defval('lbd',linspace(10e3,2000e3,100))
  defval('f2',1);
  defval('rc',2670);
  defval('drho',630);
  defval('T',35e3)
  % Young's modulus
  defval('E',1.4e11);
  % Poisson's ratio
  defval('v',0.25);
  % Gravity, unless specified
  defval('g',fralmanac('GravAcc'));
  G=fralmanac('GravCst');
  defval('xver',0)
  if xver==1
    disp(sprintf('E= %5.3g; v= %5.3f',E,v))
  end

  % Wavenumbers are in radians per meter
  k=2*pi./lbd;

  % Turcotte and Schubert (3-72)
  % Flexural Rigidity [Pa*m^3 or N*m]
  D=(E*Te.^3)/(12*(1-v^2)); 

  % Create grid on which to calculate coherence-square
  if length(Te)>=1 && length(f2)==1 && length(r)==1
    [LL,DD]=meshgrid(lbd,D);
    FF2=f2; RR=r;
  elseif length(f2)>=1 && length(Te)==1 && length(r)==1
    [LL,FF2]=meshgrid(lbd,f2);
    DD=D; RR=r;
  elseif length(r)>=1 && length(Te)==1 && length(f2)==1
    [LL,RR]=meshgrid(lbd,r);
    DD=D; FF2=f2;
  else
    error('Only one vector allowed')
  end
  KK4=(2*pi./LL).^4;
  KK=(2*pi./LL);

  % Forsyth Eqs. (3) and (6)
  xai=1+DD.*KK4/drho/g;
  phi=1+DD.*KK4/rc/g;
  beta=rc./xai/drho;

  % Should you ever need to AVERAGE you would need to average the
  % elements in the numerator and the denominator separately  
  % Olhede and Simons, Eq. (64)
  Ctop=(xai+FF2.*rc^2./drho^2.*phi-RR.*sqrt(FF2)*rc./drho.*[phi.*xai+1]).^2;
  Cbot1=xai.^2+FF2.*rc^2./drho^2-2*RR.*sqrt(FF2)*rc/drho.*xai;
  Cbot2=1+FF2.*rc^2./drho^2.*phi.^2-2*RR.*sqrt(FF2)*rc/drho.*phi;
  % See Forsyth Eq. (25)
  G2b=Ctop./Cbot1./Cbot2;

  if xver==1
    % Simply put another way, Eq. (65)
    Ctop=(xai.*drho^2+FF2.*rc^2.*phi-RR.*sqrt(FF2)*rc.*drho.*[phi.*xai+1]).^2;
    Cbot1=xai.^2.*drho^2+FF2.*rc^2-2*RR.*sqrt(FF2)*rc*drho.*xai;
    Cbot2=drho^2+FF2.*rc^2.*phi.^2-2*RR.*sqrt(FF2)*rc*drho.*phi;
    difer(G2b-Ctop./Cbot1./Cbot2,6,[],NaN);
  end

  % Check the transition wavelength in TRANSL

  % Bouguer admittance in accord with Forsyth's Eqs (11)-(12)
  Zb=-2*pi*G*rc*exp(-KK*T).*...
     (1./xai+phi.*FF2.*beta.^2-RR.*sqrt(FF2)*rc./drho.*(phi.*xai+1)./xai.^2)./...
		   (1+FF2.*beta.^2-2*RR*sqrt(FF2)*rc./drho./xai);

  % Convert to mgal/m
  Zb=Zb/1e-5;

  % Free air admittance also in mgal/m
  Zf=Zb+2*pi*G*rc/1e-5;

  l=lbd;

  % Output
  varns={G2b,k,l,Zb,Zf};
  varargout=varns(1:nargout);
elseif strcmp(Te,'demo')
  % Illustrates the functions FORSYTH, MCKENZIE, ADMITTANCE, AND TRANSL
  % For EQUAL loading at both interfaces 
  % For UNCORRELATED loading only 
  % Elastic thickness [km]
  Te=[20 80];
  % Density contrasts in kg/m^3
  DEL=[2670,630];
  % Depth to the second loading interface [km]
  T=40;
  % Calculates coherence and admittance a la Forsyth
  [G2bF,k,l,ZbF,ZfF]=forsyth(Te*1e3,[],1,0,DEL(1),DEL(2),T*1e3);
  % Figures out the transitional wavelength
  [k12,l12]=transl(1,Te,DEL(1),DEL(2));
  % Calculates coherence and admittance a la McKenzie
  [l,ZbTe1,G2bTe1,ZfTe1]=mckenzie([0 DEL(1) DEL(1)+DEL(2)],[0 T],[1 1],Te(1));
  [l,ZbTe2,G2bTe2,ZfTe2]=mckenzie([0 DEL(1) DEL(1)+DEL(2)],[0 T],[1 1],Te(2));
  % Calculates Bouguer admittance a la old-fashioned version
  [Q1,kQ]=admittance(Te(1)*1e3,1,0,T*1e3,DEL(1),DEL(2));
  [Q2,kQ]=admittance(Te(2)*1e3,1,0,T*1e3,DEL(1),DEL(2));
  
  % Make the plot
  clf
  [ah,ha,H]=krijetem(subnum(3,1));
  axes(ah(1))
  pF=semilogx(k*1000,G2bF,'Color','g','LineW',2);
  hold on
  pMTe1=semilogx(2*pi./l*1000,G2bTe1,'b');
  pMTe2=semilogx(2*pi./l*1000,G2bTe2,'r');
  set([pMTe1 pMTe2],'MarkerS',4); 
  grid on; openup(gca,6); 
  xl(1)=xlabel('Wavenumber (rad/km)');
  yl(1)=ylabel('Bouguer Coherence \gamma^2');
  pl(1)=plot(k12(1),0.5,'x');
  pl(2)=plot(k12(2),0.5,'o');
  hold off

  axes(ah(2))
  pZf=semilogx(k*1000,ZfF*1000); hold on
  pZfTe1=semilogx(2*pi./l*1000,ZfTe1*1000,'b');
  pZfTe2=semilogx(2*pi./l*1000,ZfTe2*1000,'r');
  grid on; openup(gca,6); hold off
  xl(3)=xlabel('Wavenumber (rad/km)');
  yl(2)=ylabel('Free-air Admittance Z_f (mgal/km)');

  axes(ah(3))
  pZb=semilogx(k*1000,ZbF*1000); hold on
  pZbTe1=semilogx(2*pi./l*1000,ZbTe1*1000,'b');
  pZbTe2=semilogx(2*pi./l*1000,ZbTe2*1000,'r');
  pQ1=semilogx(kQ*1000,Q1*1000,'g');
  pQ2=semilogx(kQ*1000,Q2*1000,'y');
  grid on; openup(gca,6); hold off
  xl(4)=xlabel('Wavenumber (rad/km)');
  yl(3)=ylabel('Bouguer Admittance Z_b (mgal/km)');
  
  % Cosmetics
  set(ah,'xlim',[3e-3 1e-1])
  set(ah(1),'ytick',[0:0.25:1])
  xx(1)=xtraxis1d(ah(1)); 
  xl(2)=xlabel('Wavelength (km)');
  xx(2)=xtraxis1d(ah(2)); 
  xx(3)=xtraxis1d(ah(3)); 
  longticks([ah xx],2)
  fig2print(gcf,'tall')
  figdisp
end