function gammak=gammakos(k,th,params,Hk)
% gammak=GAMMAKOS(k,th,params,Hk)
%
% Computes the first partial derivatives of the likelihood function for
% the isotropic Forsyth model in Olhede & Simons. The numbers between
% brackets indicate the order in which the derivatives are being
% presented. 
%
% INPUT:
%
% k        Wavenumber(s) at which this is to be evaluated [1/m]
% D        (1) Isotropic flexural rigidity [Nm]
% f2       (2) The sub-surface to surface initial loading ratio
% DEL      Two density contrasts, surface and subsurface [kg/m^3] 
% g        Gravitational acceleration [m/s^2]
% s2       (3) The first Matern parameter, aka sigma^2
% nu       (4) The second Matern parameter
% rho      (5) The third Matern parameter
% Hk       A complex matrix of Fourier-domain observations
%
% OUTPUT:
%
% gammak   A 5-column vector with the wavenumbers unwrapped, containing
%          the five partials of the likelihood function as columns.
%
% Last modified by fjsimons-at-alum.mit.edu, 2/10/2011

% Extract the needed parameters from the input
D=th(1);
s2=th(3);
nu=th(4);
rho=th(5);
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 get the special matrices etc.
[m,A]=mAos(k,th,params,phi,xi,pxm);

% Then get the power spectrum
S11=maternos(k,s2,nu,rho);

% Now compute the score properly speaking
gammak=nan(length(k(:)),5);
for j=1:5
  gammak(:,j)=-2*m{j}-hformos(S11,A{j},Hk);
end