function [s,a,chk] = mem(a1,a2,b1,b2,ang)
% function [s,a,chk] = mem(a1,a2,b1,b2,ang)
%
% Maximum entropy method (MEM) for the estimation of a directional
% distribution from pitch-and-roll data.
%
% (Lygre and Krogstad, JPO v16 1986: NOTE - there is a typo in the
% paper...BOTH exponetials in eq. 13 should have negative signs.
% This is fixed in Krogstad, 1991- THH has the book)
%
% edited by T.Hilmer,UH.  I think J.Aucan wrote the original code. added
% support for vectors and variable angular resolution and averaging. 
%
% Since the MEM is a fcn(coefficients), it doesn't make sense to have ab12
% as matrices.  Vectorize them, run them thru 'mem', then reshape back to
% matrix. WARNING: since this script breaks the calculation down to
% sub-degree angular resolution, then sums over angle bins, it can easily
% overload your available memory.
%
% a1,b1,b1,b2   fourier coef. of directional distribution
%               (normalized ala Long [1980]) - TRIG COORDINATES
%  
% ang   optional. ang(1) is the # of angle bins returned.   
%                 ang(2) is the # of angles within each bin
%       so the resolution is ang(1)*ang(2) and the output is an average
%       over ang(2) samples.
%       if you want the frequency bins to be truly centered, ang(1)*ang(2)
%       needs to be divisible by 4
%       TO MATCH CDIP & Aucan, use ang=[72 50] 
%
% s     MEM estimate of normalized directional spectrum, 1 deg resolution
%       in COMPASS COORDINATES (direction waves are coming from)
%
% a     radian angles of bin centers
%
% chk   check factor: MEM likes to make narrow spikes for directional
%       spectra if it can (it fits the a1,a2, b1 abd b2's exactly). If the
%       directional peak is extremely narrow then the angular resolution
%       may be insufficient (rare), and chk will be significantly .lt. 1.
%       TH: I assume the chk factor is essentially a sum of energy or
%       power.

if nargin < 5
    ang = 360;
    anr = 10;
else
    anr = ang(2);
    ang = ang(1);
end

dim = [length(a1) ang*anr];

%... switch to Lygre & Krogstad notation c
d1 = a1;
d2 = b1;
d3 = a2;
d4 = b2;

c1 = d1+i*d2;
c2 = d3+i*d4;

p1 = (c1-c2.*conj(c1))./(1-abs(c1).^2);
p2 = c2-c1.*p1;

x = 1-p1.*conj(c1)-p2.*conj(c2);

a = (0:ang*anr-1)+ceil(anr/2);% ceil or floor, will be left or right of center if not divisible by 2
a = 2*pi*a/(ang*anr);

e1 = cos(a)-i*sin(a);
e2 = cos(2*a)-i*sin(2*a);

y = abs(1-p1*e1-p2*e2).^2;% These are the memory hogs
x = repmat(x,1,dim(2));

s = abs(x./y)/dim(2);

s = reshape(s,[dim(1) anr ang]);
s = squeeze(sum(s,2));% sum across degree bins

a = reshape(a,[anr ang]);
a = squeeze(mean(a));% get angle centers also:

chk = sum(s,2);
% TH: I got 1E-15 difference between the 'chk' (power sum I assume) values
% for the loop vs vectorized version.  Strangely enough, all values for 's'
% were EXACTLY the same.
if any(abs(chk-1) > 0.01), warning('check MEM accuracy'), end

% switch from whatever the hell coordinate system we're in to cartesian:
a = 3/2*pi - a;

% %... switch from trig to compass coordinates
% ndir = 270 - ndir;
% if ndir > 360
%     ndir = ndir-360;
% end
% if ndir < 1
%     ndir = ndir+360;
% end