function lagrspec(Ndays)
%
%   The function lagrspec(Ndays) reads a number of Lagrangian floats,
%   breaks the velocity in records of length Ndays, performs an FFT
%   on each of these, and averages the results.
%
%-------------------------------------------------------------------
%   Written by Vassilis Zervakis, in February 1996, while at
%              AINCO - Interocean, Madrid, Spain.
%-------------------------------------------------------------------
%
global DATAPATH FLOATS
global Ssp1 Ssp2 fr
%
%   read a lot of floats
%
%
% ------Load floats---------------------------
SUM1 = zeros(Ndays/2,1);
SUM2 = zeros(Ndays/2,1);
ans = 1;
numfloats = 0;
nrecords = 0;
floatsread = [];
while ans == 1,
   [float, number, timer, temp, pres, lat, lon, u, v, www] = rfcread1;
%  Make sure this float has not already been read -------------
   flread = 0;
   k = 1;
   while (flread == 0) & (k <= numfloats),
      flread = strcmp( float, floatsread(k,:));
      k = k + 1;        
   end
%  if the float has not been previously read, -----------------
   if flread == 0,
     numfloats = numfloats + 1;
     if floatsread == [],
        floatsread = str2mat(float);
      else
        floatsread = str2mat(floatsread, float);
     end
     NN = length(lat);
%    if the recordlength of the float exceeds the recordlength used ------
%    in the Fourier transform, -------------------------------------------
     if NN > Ndays,
       for nread = 0:Ndays:NN,
         while nread < NN - Ndays,
%    ------break into records of length Ndays-----------------------------
           nrecords = nrecords + 1;
           u1(1:Ndays) = u(nread+1:nread+Ndays);
           v1(1:Ndays) = v(nread+1:nread+Ndays);
           nread = nread + Ndays;
%    ------Apply Fourier transform ---------------------------------------
           S1 = fft(u1);
           S2 = fft(v1);
%    ------Calculate spectrum --------------------------------------------
           SS1 = S1 .* conj(S1) / Ndays;
           SS2 = S2 .* conj(S2) / Ndays;
%    ------Add Spectrum to the sum to calculate average-------------------
           SUM1 = SUM1 + SS1(1:Ndays/2)';
           SUM2 = SUM2 + SS2(1:Ndays/2)';
         end
       end
     end
     clear number;
     clear temp;
     clear pres;
     clear wwww;
%    floats(numfloats) = float;
     charans = mbox('YESNO','Include another float in the statistics? ' ) ;
     ans = strcmp(charans, 'YES');
   else
     charans = mbox('YESNO','THIS FLOAT HAS ALREADY BEEN LOADED.',...
                            'Load another one ? ') ;
     ans = strcmp(charans, 'YES');
   end
end
%
%--Normalize the average spectrum----------------------------------------
SUM1 = SUM1 / nrecords;
SUM2 = SUM2 / nrecords;
%--Normalize the spectrum to cm^2/s^2/cph^2-------------------------------
Ssp1 = SUM1 * (24*2*pi);
Ssp2 = SUM2 * (24*2*pi);
clear SUM1 ;
clear SUM2 ;
%--Calculate frequency vector in cph-------------------------------------
%
Nyquist = Ndays/2 - 1;
fr = [0:Nyquist]'/(2*pi*24.*Ndays);
%--Calculate confidence intervals ---------------------------------------
[upper, lower] = specconf(nrecords);
up = [upper upper];
dn = [lower lower];
fr1 = [ 1 Nyquist] / (2*pi*24*Ndays);
figure('Units','cent',...
       'Position',[2 2 11 14],...
       'Name','K.E. Spectrum',...
       'MenuBar','None');
axes('Units','cent',...
     'Position',[3 2 6 9]);
loglog(fr,Ssp1,fr,Ssp2,'m--');
ylabel('cm^2/s^2/cph');
xlabel('cph');
legend('zonal','meridional');
hold on
loglog(fr1, up, 'w:', fr1, dn, 'w:');
% xx = (fr1(1) + fr1(2)) / 2.;
% yy = (0.7*upper + 1.3*lower) / 2.;
% text(xx,yy,'95%')
uicontrol('Units','cent',...
          'Position',[2 12 3 1],...
          'Style','Pushbutton',...
          'String','Energy Conserv.',...
          'Callback',['global Ssp1 Ssp2 fr,',...
                      'Econserv(Ssp1,Ssp2,fr)']);
