function errstat
%
%   This subroutines calculates the errorbars for the estimates
% of the mean velocities that were computed
% by the subroutine mapstat 
%
%-------------------------------------------------------------------
%   Written by Vassilis Zervakis, in June  1996, while at
%              AINCO - Interocean, Madrid, Spain.
%-------------------------------------------------------------------
%
%              Modified by Mike Sparrow Dec. 1996
%
%-------------------------------------------------------------------
%
global mapfigure
global minlat minlon maxlat maxlon
global sizelon sizelat sizevect gridsize
global colors isobaths skale
global X Y Umean Vmean
global DATAPATH FLOATS
global T11 T22 K11 K22
global L11 L22 HANDLERR deltau deltav
global dL11 dL22 dK11 dK22 dT11 dT22
%
%  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;
NFL = length(FLOATS(:,1));
u1 = zeros(NFL,500);
v1 = zeros(NFL,500);
dudv = zeros(NFL,500);
dudvav = zeros(NFL);
Nmax = zeros(NFL,1);
header0 = '               Results of statistical analysis               ';
header0b= '               -------------------------------               ';
header1 = ' Lon   Lat   Umean Vmean  dU    dV    U*    V*   dU*   dV*   ';
header2 = ' dec.  deg.  cm/s  cm/s  cm/s  cm/s  cm/s  cm/s  cm/s  cm/s  ';
headers = '                                                             ';
header3 = 'T11   dT11  T22   dT22  L11   dL11   L22  dL22               ' ;
header4 = '---------days---------  ------kilometers------               ' ;
header5 = '    K11        dK11        K22        dK22 R.Str. Npoints ' ;
header6 = '    --------------cm^2 s^-1--------------- (/den)         ' ;
fid1 = fopen('errorbar.asc','wt');
figure(mapfigure);
hold on;
for ipos = 1:rows,
  for jpos = 1:cols,
    Ntotalbox = 0;
    if (Umean(ipos,jpos) ~= 0.) & (Umean(ipos,jpos)< 990.),
      xmin = minlon + (jpos-1) * sizelon/60.;
      xmax = minlon + (jpos) * sizelon/60.;
      ymax = maxlat - (ipos-1) * sizelat/60.;
      ymin = maxlat - (ipos) * sizelat/60.;
      xfill = [xmin xmin xmax xmax ];
      yfill = [ymin ymax ymax ymin ];
%
%     read a lot of floats
%
% --------Load floats---------------------------
      Nmax = zeros(NFL,1);
      uvar = 0.;
      vvar = 0.;
      ncount = 0;
      for floatno = 1:NFL,
        filename = [ FLOATS(floatno,:) ] ;
        [float, number, timer, temp, pres, lat, lon, u, v, www] = rfcload1(DATAPATH,filename);
        NN = length(lat);
        clear number;
        clear temp;
        clear pres;
        clear wwww;
        kl = 1;
        datainbox = 'no ';
        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);
             dudv(floatno,kl) = u1(floatno,kl)*v1(floatno,kl);
             Nmax(floatno) = kl;
             kl = kl + 1;
             datainbox = 'yes';
             dudvav(floatno) = sum(dudv(floatno,:));
          end
        end
        clear u;
        clear v;
        clear lon;
        clear lat;
        ifl = floatno;
        if (Nmax(ifl) > 0) & (datainbox == 'yes'),
% ------( if there are actually points of this float in this box...)----
% ------Calculate variance of primed velocities --------------
          ncount = ncount + Nmax(ifl);
          uvar = uvar + sum((u1(ifl,:) .^ 2));
          vvar = vvar + sum((v1(ifl,:) .^ 2));
          Ntotalbox = Ntotalbox + Nmax(ifl);
        end
      end
% ----Having read all the floats for this box, estimate error bars....
%
%  Check if this is the first time T11 and T22 have been calculated 
%
if deltau ==[]
   deltau = zeros(rows,cols) ;
   deltav = zeros(rows,cols) ;
end
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
%
%  Check to see if T has been estimated for this box...
%
%T11
%T22
%dT11
%dT22
%K11
%K22
%dK11
%dK22
rstress = -sum(dudvav(:));
urms = sqrt(uvar/ncount);
vrms = sqrt(vvar/ncount);
deltau(ipos,jpos) = (2*urms) / (sqrt(Ntotalbox/T11(ipos,jpos)));
deltav(ipos,jpos) = (2*vrms) / (sqrt(Ntotalbox/T22(ipos,jpos)));
durms  = 1.96*urms/(sqrt(2*Ntotalbox/T11(ipos,jpos)));
dvrms  = 1.96*vrms/(sqrt(2*Ntotalbox/T22(ipos,jpos)));
dL11(ipos,jpos)=L11(ipos,jpos)*sqrt((durms/urms).^2+(dT11(ipos,jpos)/T11(ipos,jpos)).^2) ;
dL22(ipos,jpos)=L22(ipos,jpos)*sqrt((dvrms/vrms).^2+(dT22(ipos,jpos)/T22(ipos,jpos)).^2) ;
dK11(ipos,jpos)=K11(ipos,jpos)*sqrt(2*((durms/urms).^2)+(dT11(ipos,jpos)/T11(ipos,jpos).^2)) ;
dK22(ipos,jpos)=K22(ipos,jpos)*sqrt(2*((dvrms/vrms).^2)+(dT22(ipos,jpos)/T22(ipos,jpos).^2)) ;
xcent = (xmin+xmax) /2.;
ycent = (ymin+ymax) /2.;
times = 100 / (sizevect * gridsize) ;
fprintf(fid1,'%62s\n',headers);
fprintf(fid1,'%62s\n',header0);
fprintf(fid1,'%62s\n',header0b);
fprintf(fid1,'%62s\n',headers);
fprintf(fid1,'%62s\n',header1);
fprintf(fid1,'%62s\n',header2);
printvar = [ xcent; ycent; Umean(ipos,jpos); Vmean(ipos,jpos);...
            deltau(ipos,jpos); deltav(ipos,jpos); urms; vrms; durms; dvrms] ;
fprintf(fid1,'%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f', printvar);
fprintf(fid1,'%12s\n','            ');
fprintf(fid1,'%62s\n',headers);
fprintf(fid1,'%62s\n',header3);
fprintf(fid1,'%62s\n',header4);
printvar2 = [ T11(ipos,jpos); dT11(ipos,jpos); T22(ipos,jpos); dT22(ipos,jpos);...
            L11(ipos,jpos); dL11(ipos,jpos); L22(ipos,jpos); dL22(ipos,jpos)] ;
fprintf(fid1,'%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f', printvar2);
fprintf(fid1,'%27s\n','                           ');
fprintf(fid1,'%62s\n',headers);
fprintf(fid1,'%62s\n',header5);
fprintf(fid1,'%62s\n',header6);
printvar3 = [ K11(ipos,jpos); dK11(ipos,jpos); K22(ipos,jpos); dK22(ipos,jpos); rstress; Ntotalbox] ;
fprintf(fid1,'%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5i\n', printvar3);
fprintf(fid1,'%27s\n','                           ');
% .......and plot them as ellipses around the origin of the velocity vector --
%   if (rows > 1) & (cols > 1),
%
%     ik = 0 ;
%     for theta = 0 : 30 : 360 ;
%       ik = ik + 1 ;
%       radangle = theta * pi / 180 ;
%
%       uuu(ik)=X(ipos,jpos)+skale*deltau*cos(radangle)/(sizevect*gridsize) ;
%       vvv(ik)=Y(ipos,jpos)+skale*deltav*sin(radangle)/(sizevect*gridsize) ;
%
%     end
%     line(uuu,vvv)
%     drawnow
%    end
    end
  end
end
fprintf(fid1,'%12s\n','            ');
fprintf(fid1,'%12s\n','            ');
fprintf(fid1,'%36s\n','U* = root of the mean square of uï [= U-Umean]');
fclose(fid1);
set(HANDLERR,'Enable','off');
refresh(mapfigure);
mbox('OK','The results are stored in the text file','errorbar.asc in the current directory');
