Fisherkos.m

function mcF=Fisherkos(k,th,params,xver)
% mcF=Fisherkos(k,th,params,xver)
%
% Calculates the entries in the Fisher matrix for Olhede & Simons (2013)
% for the UNCORRELATED initial-loading model prior to wavenumber averaging. 
%
% INPUT:
%
% k        Wavenumber(s) at which this is to be evaluated [1/m]
% th       The five-parameter vector with elements:
%          th(1)=D    Isotropic flexural rigidity 
%          th(2)=f2   The sub-surface to surface initial loading ratio 
%          th(3)=s2   The first Matern parameter, aka sigma^2 
%          th(4)=nu   The second Matern parameter 
%          th(5)=rho  The third Matern parameter 
% params   A structure with AT LEAST these constants that are known:
%          DEL   surface and subsurface density contrast [kg/m^3]
%          g     gravitational acceleration [m/s^2]
%
% OUTPUT:
%
% mcF      The 15-column Fisher matrix, in this order:
%          [1]  FDD    [2]  Ff2f2
%          [3]  Fs2s2  [4]  Fnunu [5] Frhorho
%          [6]  FDf2   [7]  FDs2  [8]  FDnu  [9] FDrho
%          [10] Ff2s2  [11] Ff2nu [12] Ff2rho
%          [13] Fs2nu  [14] Fs2rho
%          [15] Fnurho
%
% EXAMPLE:
% 
% [~,~,th0,p,k]=simulos([],[],[],1);
% mcF=Fisherkos(k,th0,p,1);
%
% Last modified by fjsimons-at-alum.mit.edu, 10/25/2014

defval('xver',1)

% Extract the parameters from the input
D=th(1);
f2=th(2);
DEL=params.DEL;
g=params.g;

% The number of parameters to solve for
np=length(th);
% The number of unique entries in an np*np symmetric matrix
npp=np*(np+1)/2;

% First the auxiliary quantities
phi=phios(k,D,DEL,g);
xi = xios(k,D,DEL,g);
% Note that this has a zero at zero wavenumber
pxm=(phi.*xi-1);

% First compute the "means" parameters
m=mAos(k,th,params,phi,xi,pxm,xver);

% Forcefully set f2 to a positive number even if it means a throw back
f2=abs(f2);

% Then compute the various entries in the order of the paper
lk=length(k(:));
% Some of them depend on the wave vectors, some don't
mcF=cellnan([npp 1],...
	    [lk 1 1 repmat(lk,1,6) 1  repmat(lk,1,5)],...
	    repmat(1,1,npp));

% FDD
warning off MATLAB:divideByZero
fax=2*k(:).^8/g^2*dpos(DEL,-2,-2)./pxm.^2.*xi.^(-2);
warning on MATLAB:divideByZero
mcF{1}=fax.*(2*dpos(DEL,0,2)*xi.^4+2*dpos(DEL,2,0)*phi.^2.*xi.^2+...
	     4*dpos(DEL,2,0)-4*dpos(DEL,1,1)*phi.*xi.^3+6*dpos(DEL,1,1)*xi.^2-... 
	     4*dpos(DEL,2,0)*phi.*xi+f2*dpos(DEL,2,0)*xi.^2+...
	     dpos(DEL,0,2)*xi.^2/f2);
% Ff2f2
mcF{2}=1/f2^2;

% Fthsths
for j=3:np
  mcF{j}=2*m{j}.^2;
end

% FDf2 might write it in explicitly, compare to FISHERKROS
mcF{np+1}=ddTos(k,th,params,phi,xi,pxm);

% All further combinations
jcombo=combntns(1:np,2);
for j=2:length(jcombo)
  mcF{np+j}=2*m{jcombo(j,1)}.*m{jcombo(j,2)};
end

% Verification step
if xver==1
  % The "linear" tracecheck has already happed in MAOS
  % This should be twice the negative sum of the eigenvalues of
  % the product of this...
  [M,N]=dTos(k,th,params,phi,xi,pxm);
  % and this...
  [invT,detT,L]=Tos(k,th,params,phi,xi,pxm);
  % ... which we'll be again checking at a random wavenumber
  tracecheck(L,{M N},{m{1:2}},9,1)

  % ... and then these also for ueber-completeness
  for ind=3:np
    tracecheck(L,{-repmat(m{ind},[size(invT,1)/size(m{ind},1)],3).*invT},...
	       {m{ind}},9,1)
  end

  % And then we check the sum of the squares of the eigenvalues which
  % hasn't happened so far
  tracecheck(L,{M N},{mcF{1:2}},8,2)
end