function lagrspec(Ndays, rotary)
%
%   The function lagrspec(Ndays,rotary) 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.
%   if the variable rotary is 1, then the spectrum will be rotary
%   if not, it will be in regular zonal and meridional components.
%
%-------------------------------------------------------------------
%   Written by Vassilis Zervakis, in February 1996, while at
%              AINCO - Interocean, Madrid, Spain.
%   Modified by Mike Sparrow, March 1998.
%-------------------------------------------------------------------
%
global DATAPATH FLOATS tjump
global Ssp1 Ssp2 fr rot
%
%   read a lot of floats
%
%
% ------Load floats---------------------------
%
modays = Ndays*tjump
%
% modays is the increased avaraging interval due to tjump.
% tjump is the number of position fixe per day of the float.
%
SUM1 = zeros(modays/2,1);
SUM2 = zeros(modays/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 > modays,
       for nread = 0:modays:NN,
         while nread < NN - modays,
%    ------break into records of length Ndays (modays)---------------------
           nrecords = nrecords + 1;
           u1(1:modays) = u(nread+1:nread+modays);
           v1(1:modays) = v(nread+1:nread+modays);
           nread = nread + modays;
%    ------Apply Fourier transform ---------------------------------------
           S1 = fft(u1);
           S2 = fft(v1);
%    ------Calculate spectrum --------------------------------------------
           if rotary == 1,
%    ----------in rotary cordinates---------------------------------------
               Ucount = ((real(S1) + imag(S2)) + i*(real(S2) - imag(S1))) / 2;
               Uclock = ((real(S1) - imag(S2)) + i*(real(S2) + imag(S1))) / 2; 
               SS1 = Ucount .* conj(Ucount) / modays;
               SS2 = Uclock .* conj(Uclock) / modays;
             else
%    ----------in cartesian coordinates-----------------------------------
               SS1 = S1 .* conj(S1) / modays;
               SS2 = S2 .* conj(S2) / modays;
           end;
%    ------Add Spectrum to the sum to calculate average-------------------
           SUM1 = SUM1 + SS1(1:modays/2)';
           SUM2 = SUM2 + SS2(1:modays/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-------------------------------
Ssp1 = SUM1 * (24*2*pi)/tjump;
Ssp2 = SUM2 * (24*2*pi)/tjump;
clear SUM1 ;
clear SUM2 ;
%--Calculate frequency vector in cph-------------------------------------
%
Nyquist = modays/2 - 1;
fr = [0:Nyquist]'/(24.*Ndays);
%--Calculate confidence intervals ---------------------------------------
[upper, lower] = specconf(nrecords);
up = [upper upper];
dn = [lower lower];
fr1 = [ 1 Nyquist] / (24*Ndays);
figure('Units','cent',...
       'Position',[2 2 11 14],...
       'Name','K.E. Spectrum');
axes('Units','cent',...
     'Position',[3 2 6 9]) ;
loglog(fr,Ssp1,'y',fr,Ssp2,'c--') ;
ylabel('cm^2/s^2/cph') ;
xlabel('cph') ;
if rotary == 1, 
  legend('counterclockwise','clockwise');
else,
  legend('zonal','meridional') ;
end
hold on
loglog(fr1, up, 'w:', fr1, dn, 'w:') ;
rot = rotary;
uicontrol('Units','cent',...
          'Position',[2 12 3 1],...
          'Style','Pushbutton',...
          'String','Energy Conserv.',...
          'Callback',['global Ssp1 Ssp2 fr rot,',...
                      'Econserv(Ssp1,Ssp2,fr, rot)']);
