function D = mlm(r1,r2,alpha1,alpha2,dirs)

%MLM Maximum Likelihood Method calculation for directional distribution function
%   D = mlm(R1,R2,ALPHA1,ALPHA2,DIRS) calculates the directional distribution
%   function for a 2D spectrum for NDBC parameters using the Maximum
%   Likelihood Method as outlined by Oltman-Shay & Guza (1984 - JPO). The 
%   inputs are defined on http://www.ndbc.noaa.gov/measdes.shtml - R1 and 
%   R2 are the two spreading parameters (NOTE: 0<r1,r2<1), ALPHA1 and 
%   ALPHA2 are the mean and principal wave directions, and DIRS is the 
%   directional space of the distribution.
%
%   D = mlm(r1,r2,alpha1,alpha2,dirs)
%
% Dave Thompson (dthompson@usgs.gov)

%%
% Set k and c11 to unity. 
k = 1;
c11 = 1;

% Transform directions from Nautical (CW from N) to Cartesian (CCW from E).
dirs = ((270-dirs)<0)*360 + (270-dirs);

% Calculate co and quad-spectrum (from Earle, 1999)
c22 = (1/2)*k^2*c11*(1-r2*cosd(2*alpha2));
c33 = (1/2)*k^2*c11*(1+r2*cosd(2*alpha2));
c23 = (1/2)*k^2*c11*r2*sind(2*alpha2);
q12 = -k*c11*r1*sind(alpha1);
q13 = -k*c11*r1*cosd(alpha1);
c32 = conj(c23);
q21 = -q12;
q31 = -q13;

% cross-spectral data matrix
M = [c11 i*q12 i*q13;
   i*q21 c22 c23;
   i*q31 c32 c33];
M = inv(M); % for Eq. 6 in O-S & G.

% G
for ii=1:length(dirs)
   G(ii,1) = 1;
   G(ii,2) = i*cosd(dirs(ii));
   G(ii,3) = i*sind(dirs(ii));
end

% Calculate direction distribution -> Eq. 6 in O-S & G.
D = zeros(1,length(dirs));
for ii=1:length(dirs)
   for n=1:3
      for m=1:3
         D(ii) = D(ii) + M(n,m) * G(ii,n) * conj(G(ii,m));
      end
   end
   D(ii) = 1/D(ii);
end
D = D/sum(D);