% Grid, in SI units:
x = [0:5:1000]' * 1000;
y = [0:5:500]' * 1000;

[X,Y] = meshgrid(x,y);

% SI units
% 1m max SSH eddy in midlatitudes
psi0 = 0.30;
g0= 9.8;
f0 = 1e-4;


% Parameters

sig = 30e3;	% width
y0 = max(y)/2;	% center latitude in km

% initialize M
% add to this array
M = zeros(size(X));
count = zeros(size(X));
unity = ones(size(X));


for step = 1:length(x)

% eddy centre in m
   x0 = squeeze(x(end+1-step));

% gaussian anticyclone streamfunction
   psi = psi0*exp(-((X-x0).^2 + (Y-y0).^2)/(2*sig^2));

% velocities
   u =  (Y-y0)/sig^2 .* psi*g0/f0;
   v = -(X-x0)/sig^2 .* psi*g0/f0;

   ke = (u.^2+v.^2);
   num = (u.^2-v.^2);
% only average where the eddy has influence above some threshold
   g = find(ke>0.1^2);
   count(g) = count(g)+unity(g);
   M(g) = M(g) + num(g)./ke(g);

end

% mean of M where there was data
g = find(count>0);
M(g) = M(g) ./ count(g);
b = find(count==0);
M(b) = NaN;


% plot resulting M field
figure(1)
clf
contourf(x/1000,y/1000,M)
set(gca,'XLim',[200 800])
set(gca,'YLim',[150 350])
xlabel('zonal distance in km')
ylabel('meridional distance in km')
title('Westerward propagating Gaussian anticyclone: sigma = 30km')
caxis([-0.8 0.8])