function eddy01
%
%   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
%                                         Jan. '98
%-------------------------------------------------------------------
%
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 HANDLERR
global L11 L22 K11 K22
global dT11 dT22 dL11 dL22 dK11 dK22
global lag R1 R2 Tu Tv
%  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) ;
   T22 = 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) ;
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 ];
      patch(xfill,yfill,'b')
%drawnow
NFL = length(FLOATS(:,1));
u1 = zeros(NFL,500);
v1 = 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,:)
   DATAPATH
   [float,number,timer,temp,pres,lat,lon,u,v,www] = rfcload1(DATAPATH,filename);
   NN = length(lat);
   clear number;
   clear temp;
   clear pres;
   clear wwww;
   clear start;
   clear finish;
   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----------------------
        u1(floatno,kl) = u(kk) - Umean(ipos,jpos);
        v1(floatno,kl) = v(kk) - Vmean(ipos,jpos);
        Nmax(floatno) = kl;
        kl = kl + 1;
    end
   end
   Nmax
   clear u;
   clear v;
   clear lon;
   clear lat;
end
% -----------Actually calculate correlations -------------------
[mm,nn] = size(u1)
R11 = zeros(mm,nn);
R12 = zeros(mm,nn);
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...)----
    ncount = ncount + Nmax(ifl);
    S1 = xcorr(u1(ifl,1:Nmax(ifl)),'coeff');
    S2 = xcorr(v1(ifl,1:Nmax(ifl)),'coeff');
    n = length(S1);
    m = (n+1)/2;
    size(R11(ifl,:))
    R11(ifl,1:n-m+1) = S1(m:n)';
    R12(ifl,1:n-m+1) = S2(m:n)';
% --Calculate variance of primed velocities --------------
    uvar = uvar + sum((u1(ifl,:) .^ 2));
    vvar = vvar + sum((v1(ifl,:) .^ 2));
  end
end
uvar = uvar / ncount;
vvar = vvar / ncount;
clear S1;
clear S2;
for kk = 1:max(Nmax),
% calculate which and how many floats have data at this lag step
% and fill it in the variable fl -------------------
% --Compute the average correlation for all floats ---
    R1(kk) = 0;
    R2(kk) = 0;
    nf = 0; 
    for kl = 1:NFL,
      if kk <= Nmax(kl),
        R1(kk) = R1(kk) + R11(kl,kk);
        R2(kk) = R2(kk) + R12(kl,kk);
        nf = nf + 1;
      end
    end
    lag(kk) = kk-1;
    R1(kk) = R1(kk)/nf;
    R2(kk) = R2(kk)/nf;
    fl(kk) = nf;
  end
end
clear R11;
clear R12;
clear nhave;
fig1 = figure('Units','cent',...
       'Position',[2 3 12 8]);
axes('Units','cent','Position',[1.5 3.5 10 4],...
     'YLim',[-0.45 1.0],...
     'XTickLabelMode','manual','XTickLabels',[' ']);
line('Xdata',lag(1:max(Nmax)),'Ydata',R1(1:max(Nmax)),'Color','y');
line('Xdata',lag(1:max(Nmax)),'Ydata',R2(1:max(Nmax)),'Color','m');
zoom on;
hold on;
grid on;
ylabel('Autocorrelation');
hold off;
axes('Units','cent','Position',[1.5 1.5 10 2]);
line(lag(1:max(Nmax)), fl(1:max(Nmax)));
t = [ 'Lag (data cycles) [1 day = ',num2str(tjump),' data cycle(s)]'];
ylabel('# of floats');
xlabel(t);
grid on;
drawnow
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 4 12 8]);
YYYmax = max(max(Tu),max(Tv)) + 2;
axes('Units','cent','Position',[1.5 3.5 10 3],...
     'YLim',[0 YYYmax],...
     'XTickLabelMode','manual','XTickLabels',[' ']);
line('Xdata',lag,'Ydata',Tu,'Color','y');
line('Xdata',lag,'Ydata',Tv,'Color','m');
zoom on;
hold on;
grid on;
ylabel('Int. T. Scale (Days)');
hold off;
axes('Units','cent','Position',[1.5 1.5 10 1]);
line(lag(1:max(Nmax)), fl(1:max(Nmax)));
ylabel('# of floats');
t = [ 'Lag (data cycles) [1 day = ',num2str(tjump),' data cycle(s)]'];
xlabel(t);
grid on;
drawnow
%
hold on;
grid on;
clear R1;
clear R2;
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;
%
% ---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);
%
% ---Calculate standard deviation of integral time scales ---
%
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 diffusivities----------------------
%
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)) ];
%text5 = [ text1 '±' num2str(Tsig1) ' days'];
%text6 = [ text2 '±' num2str(Tsig2) ' days'];
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);

set(HANDLERR,'Enable','on') ;
mbox('OK','click on "Errors" in order to obtain','error estimates for K,L etc') ;

