function eddy02
%
%   This subroutines calculates the eddy diffusivities on
% every grid point that the mean velocities were computed
% by the subroutine mapstat, based on the method of 
% integral time.
%
%-------------------------------------------------------------------
%   Written by Vassilis Zervakis, in March 1996, while at
%              AINCO - Interocean, Madrid, Spain.
%
%              Modified by Mike Sparrow in Jan.'97
%                                          Oct.'97
%
%-------------------------------------------------------------------
%
global mapfigure tjump
global minlat minlon maxlat maxlon
global sizelon sizelat
global colors isobaths
global X Y Umean Vmean
global DATAPATH FLOATS
global T11 T22 K11 K22 L11 L22
global dT11 dT22 dL11 dL22 dK11 dK22
%
%  Convert latitudes and longitudes in minutes -------
Ymax = fix(maxlat)*60. + (maxlat-round(maxlat))*0.6;
Ymin = fix(minlat)*60. + (minlat-round(minlat))*0.6;
Xmax = fix(maxlon)*60. + (maxlon-round(maxlon))*0.6;
Xmin = fix(minlon)*60. + (minlon-round(minlon))*0.6;
[rows, cols] = size(Umean) ;
boxes = rows * cols;
%
%  Check if this is the first time T, K and L have been calculated 
%
if T11 ==[]
   T11 = ones(rows,cols) ;
   dT11 = zeros(rows,cols) ;
   dT22 = zeros(rows,cols) ;
end

if K11 ==[]
   K11 = zeros(rows,cols) ;
   K22 = zeros(rows,cols) ;
   dK11 = zeros(rows,cols);
   dK22 = zeros(rows,cols);
end

if L11 ==[]
   L11 = zeros(rows,cols) ;
   L22 = zeros(rows,cols) ;
   dL11 = zeros(rows,cols) ;
   dL22 = zeros(rows,cols) ;
   T22 = ones(rows,cols) ;
end

% ---show which box to compute statistics at: ---------
mbox('OK','Click the cursor where you want','eddy statistics computed');
figure (mapfigure);
hold on;
[x1,y1] = ginput(1);
% ----find which box we have clicked
      xminutes = (x1-minlon)*60.;
      yminutes = (maxlat-y1)*60.;
      ipos = fix(yminutes / sizelat);
      jpos = fix(xminutes / sizelon);
      xmin = minlon + jpos * sizelon/60.;
      xmax = minlon + (jpos+1) * sizelon/60.;
      ymax = maxlat - ipos * sizelat/60.;
      ymin = maxlat - (ipos+1) * sizelat/60.;
      xfill = [xmin xmin xmax xmax ];
      yfill = [ymin ymax ymax ymin ];
NFL = length(FLOATS(:,1));
u1 = zeros(NFL,500);
v1 = zeros(NFL,500);
ur = zeros(NFL,500);
vr = zeros(NFL,500);
Nmax = zeros(NFL,1);
%
%   read a lot of floats
%
ipos = ipos + 1;
jpos = jpos + 1;
size(FLOATS)
% ------Load floats---------------------------
for floatno = 1:NFL,
   filename = [ FLOATS(floatno,:) ] ;
   [float,number,timer,temp,pres,lat,lon,u,v,www,head,start,finish] = rfcload1(DATAPATH,filename);
   NN = length(lat);
   clear number;
   clear temp;
   clear pres;
   clear wwww;
   kl = 1;
   for kk = 1:NN,
    if (lat(kk)<=ymax)&(lat(kk)>=ymin)&(lon(kk)>=xmin)&(lon(kk)<=xmax),
%   ----Subtract the mean velocity of that box----------------------
%   ----(Calculate the prime velocities ----------------------------
        u1(floatno,kl) = u(kk) - Umean(ipos,jpos);
        v1(floatno,kl) = v(kk) - Vmean(ipos,jpos);
%   ----(Save the raw velocities within the box---------------------
        ur(floatno,kl) = u(kk);
        vr(floatno,kl) = v(kk);
        Nmax(floatno) = kl;
        kl = kl + 1;
    end
   end
   Nmax
   clear u;
   clear v;
   clear lon;
   clear lat;
end
% -----------Compute the structure functions -------------------
[mm,nn] = size(u1)
D1 = zeros(nn,1);
D2 = zeros(nn,1);
npoints = zeros(nn,1);
R1 = zeros(nn,1);
R2 = zeros(nn,1);
uvar = 0.;
vvar = 0.;
ncount = 0;
for ifl = 1:NFL,
  if Nmax(ifl) > 1,
% --( if there are actually points of this float in this box...)----
    for t2 = 1: Nmax(ifl)-1,
      for t1 = 1:Nmax(ifl)-t2,
         npoints(t1) = npoints(t1) + 1;
         D1(t1) = D1(t1) + (ur(ifl,t1+t2)-ur(ifl,t2)).^2 ;
         D2(t1) = D2(t1) + (vr(ifl,t1+t2)-vr(ifl,t2)).^2 ;
      end
    end
% --Calculate variance of primed velocities --------------
    ncount = ncount + Nmax(ifl);
    uvar = uvar + sum((u1(ifl,1:Nmax(ifl)) .^ 2));
    vvar = vvar + sum((v1(ifl,1:Nmax(ifl)) .^ 2));
  end
end
uvar = uvar / ncount;
vvar = vvar / ncount;
kk = 1;
while npoints(kk) > 0,
  D1(kk) = D1(kk) ./ npoints(kk);
  D2(kk) = D2(kk) ./ npoints(kk);
  kk = kk + 1;
end
lag = [0:kk-1];
h5 = figure
plot(lag,D1(1:kk));
xlabel('time step');
ylabel('Structure Function D');
hold on
plot(lag,D2(1:kk));
pause
close(h5);
% --Calculate correlation coefficients from structure functions---
R1 = 1 - D1 ./ (2*uvar);
R2 = 1 - D2 ./ (2*vvar);
clear D1;
clear D2;
lag = [1:max(Nmax)];
fig1 = figure('Units','cent',...
       'Position',[2 2 12 10]);
plot(lag(1:max(Nmax)), R1(1:max(Nmax)), lag(1:max(Nmax)), R2(1:max(Nmax)))
hold on;
grid on;
t = [ 'Lag (data cycles) [1 day = ',num2str(tjump),' data cycle(s)]'];
xlabel(t);
ylabel('Autocorrelation');
hold off;
pause
close(fig1);
% ---Compute integral time scales ------
for kk = 1:max(Nmax),
  Tu(kk) = (sum(R1(1:kk-1)))/tjump;
  Tv(kk) = (sum(R2(1:kk-1)))/tjump;
end
fig2 = figure('Units','cent',...
              'Position',[3 3 12 10]);
plot(lag,Tu,lag,Tv);
hold on;
grid on;
clear R1;
clear R2;
t = [ 'Lag (data cycles) [1 day = ',num2str(tjump),' data cycle(s)]'];
xlabel(t);
ylabel('Integral time scale (days) ');
hold off
pause
% -----Compute Eddy diffusivities-------
mbox('OK','Use the mouse to show beginning and end of integration',...
          'for computation of K11 (E-W comp., yellow line) ');
figure(fig2);
A = ginput(2);
tx1 = A(1,1);
tx2 = A(2,1);
if tx2 > length(Tu),
   tx2 = length(Tu);
end
clear A;
mbox('OK','Use the mouse to show beginning and end of integration',...
          'for computation of K22 (N-S comp., magenta line) ');
figure(fig2);
A = ginput(2);
ty1 = A(1,1);
ty2 = A(2,1);
if ty2 > length(Tv),
   ty2 = length(Tv);
end
clear A;
close;
% ---Integrate to compute integral time scales------------
T11(ipos,jpos) = sum(Tu(tx1:tx2))/(tx2 - tx1 + 1);
T22(ipos,jpos) = sum(Tv(ty1:ty2))/(ty2 - ty1 + 1);
%
dT11(ipos,jpos) = sqrt(sum( (Tu(tx1:tx2) - T11(ipos,jpos)) .^2 ) / (tx2-tx1+1));
dT22(ipos,jpos) = sqrt(sum( (Tv(ty1:ty2) - T22(ipos,jpos)) .^2 ) / (ty2-ty1+1));
%
% --- Computation of eddy diffusicities ------------------
K11(ipos,jpos) = (T11(ipos,jpos)*24.*3600.) * uvar ;
K22(ipos,jpos) = (T22(ipos,jpos)*24.*3600.) * vvar ;
L11(ipos,jpos) = T11(ipos,jpos).*sqrt(uvar);
L22(ipos,jpos) = T22(ipos,jpos).*sqrt(vvar);
yst = (maxlat - minlat) / 20;
xst = (maxlon - minlon) / 20;
yy1 = (ymax + ymin)/2 - 2.0*yst;
yy2 = (ymax + ymin)/2 - 1.0*yst;
yy3 = (ymax + ymin)/2;
yy4 = (ymax + ymin)/2 + 1.0*yst;
yy7 = (ymax + ymin)/2 + 2.0*yst;
yy8 = (ymax + ymin)/2 + 3.0*yst;
xx1 = (xmax + xmin)/2;
text1 = [ 'T11 = ' num2str(T11(ipos,jpos)) ' +/- ' num2str(dT11(ipos,jpos))];
text2 = [ 'T22 = ' num2str(T22(ipos,jpos)) ' +/- ' num2str(dT22(ipos,jpos))];
text3 = [ 'K11 = ' num2str(K11(ipos,jpos)) ];
text4 = [ 'K22 = ' num2str(K22(ipos,jpos)) ];
text7 = [ 'L11 = ' num2str(L11(ipos,jpos)) ];
text8 = [ 'L22 = ' num2str(L22(ipos,jpos)) ];
figure(mapfigure);
hold on
if colors == 'yes',
  textcolor = 'k';
 else
  textcolor = 'w';
end
pos = [xx1,yy1];
text('Units','data',...
     'FontSize',9,...
     'HorizontalAlignment','center',...
     'EraseMode','none',...
     'Color',textcolor,...
     'String',text1,...
     'Position',pos);
pos = [xx1,yy2];
text('Units','data',...
     'FontSize',9,...
     'HorizontalAlignment','center',...
     'EraseMode','none',...
     'Color',textcolor,...
     'String',text2,...
     'Position',pos);
pos = [xx1,yy3];
text('Units','data',...
     'FontSize',9,...
     'HorizontalAlignment','center',...
     'EraseMode','none',...
     'Color',textcolor,...
     'String',text3,...
     'Position',pos);
pos = [xx1,yy4];
text('Units','data',...
     'FontSize',9,...
     'HorizontalAlignment','center',...
     'EraseMode','none',...
     'Color',textcolor,...
     'String',text4,...
     'Position',pos);
pos = [xx1,yy7];
text('Units','data',...
     'FontSize',9,...
     'HorizontalAlignment','center',...
     'EraseMode','none',...
     'Color',textcolor,...
     'String',text7,...
     'Position',pos);
pos = [xx1,yy8];
text('Units','data',...
     'FontSize',9,...
     'HorizontalAlignment','center',...
     'EraseMode','none',...
     'Color',textcolor,...
     'String',text8,...
     'Position',pos);

drawnow
