function mcF=Fisherkos(k,th,params,xver)
% mcF=Fisherkos(k,th,params,xver)
%
% Calculates the entries in the Fisher matrix for Olhede & Simons
% the UNCORRELATED initial-loading model prior to wavenumber averaging. 
%
% INPUT:
%
% k        Wavenumber(s) at which this is to be evaluated [1/m]
% th       The parameter vector with elements:
%          D   Isotropic flexural rigidity [Nm]
%          f2  The sub-surface to surface initial loading ratio
%          s2  The first Matern parameter, aka sigma^2
%          nu  The second Matern parameter
%          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
%
% Last modified by fjsimons-at-alum.mit.edu, 02/10/2011

defval('xver',0)

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

% 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);

% Then compute the various entries in the order of the paper
% WATCH OUT TO MAKE SURE THIS INITIALIZATION IS OK
%mcF=cellnan(10,length(k(:)),10);
lk=length(k(:));
mcF=cellnan([15 1],[lk 1 1 repmat(lk,1,6) 1  repmat(lk,1,5)],repmat(1,1,15));

% FDD
warning off MATLAB:divideByZero
fax=2*k(:).^8/g^2*dpos(DEL,-2,-2).*xi.^(-2)./pxm.^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:5
  mcF{j}=2*m{j}.^2;
end

% FDf2
mcF{6}=ddTos(k,th,params,phi,xi,pxm);

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

% Verification step, use TRACECHECK HERE
if xver==1
  % This should be the sum of the squared eigenvalues of 
  [M,N]=dTos(k,th,params,phi,xi,pxm);
  [invT,detT,L]=Tos(k,th,params,phi,xi,pxm);
  randi=round(rand*length(invT));
  Mrand=[M(randi,1) M(randi,2) ; M(randi,2) M(randi,3)];
  Nrand=[N(randi,1) N(randi,2) ; N(randi,2) N(randi,3)];  
  Lrand=[L(randi,1)    0       ; L(randi,2) L(randi,3)];
  % Check RELATIVE precision
  % Already done elsewhere but do again to be sure
  difer((-2*m{1}(randi)-sum(eig(Lrand'*Mrand*Lrand)))/m{1}(randi))
  difer((-2*m{2}-sum(eig(Lrand'*Nrand*Lrand)))/m{2})
  difer(mcF{1}(randi)-sum([eig(Lrand'*Mrand*Lrand)].^2))
  difer((mcF{2}-sum([eig(Lrand'*Nrand*Lrand)].^2))/mcF{2})
end