function cometall(x, y, p, trailcolor)
% COMETALL(x, y, p, trailcolor)  Comet plot of lines whose x 
% coordinates are stored in the rows of the matrix x and the
% y coordinates in the rows of the matrix y.  Every column of 
% these matrices corresponds to a line. p is the length of the 
% trail, in rows of matrix. trailcolor is the color of the
% integrated path of the worm.
% Subroutine derived from the Matlab's subroutine COMET (R),
% by The Mathworks, Inc.
% -------------------------------------------------------------
% Written by Vassilis Zervakis, in February 1996, while at
%            AINCO - Interocean, Madrid, Spain
% -------------------------------------------------------------
%
global firstday firstmonth firstyear
global h1 timehandle
global dd mm yy
global pauseswitch
pauseswitch = 'off';
[M,N] = size(x);
starter = dmy2rd(firstday,firstmonth, firstyear);
% Make a Nx2 Matrix that contains the index of beginning -----------------
% and the index of end of good data (less than 999.)----------------------
ind = zeros(2,N);
ind(1,:) = ind(1,:) + 1;
ind(2,:) = ind(2,:) + M;
for l = 1:N
  for ll = 2: M,
% ---index of start of good data -----------------------------------------
     if (x(ll,l) < 999.) & (x(ll-1,l) == 999.),
        ind(1,l) = ll;
     end
% ---index of end of good data -------------------------------------------
     if (x(ll,l) == 999.) & (x(ll-1,l) < 999.),
        ind(2,l) = ll;
     end
  end
end
%ind
% If the sea is black, make the floats white--------------------------
if trailcolor == 'k',
  bodycolor = 'w';
elseif trailcolor == 'c',
  bodycolor = 'r';
else,
  bodycolor = 'r';
end
% Bodycolor dot head, and line body, and trailcolor line tail.--------------
for l = 1:N,
  head(l) = line('color',bodycolor,'linestyle','.','erase','xor','xdata',x(ind(1,l),l),'ydata',y(ind(1,l),l));
  body(l) = line('color',bodycolor,'linestyle','-','erase','none','xdata',[],'ydata',[]);
  tail(l) = line('color',trailcolor,'linestyle','-','erase','none','xdata',[],'ydata',[]);
end

m = length(x(:,1));
k = p;
size(x);
size(y);
% Grow the body
for i = 2:k+1
  for l = 1:N,
     j = i-1:i;
     if (i < ind(2,l)) & (i >= ind(1,l))
        set(head(l),'xdata',x(i,l),'ydata',y(i,l));
     end
     if (i < ind(2,l)) & (i-1 >= ind(1,l))
        set(body(l),'xdata',x(j,l),'ydata',y(j,l));
     end
  end
% Calculate time --------------------------------------
  tim = starter + i - 1;
  timer = rd2dmy(tim);
  dd = timer(1);
  mm = timer(2);
  yy = timer(3);
  clear timer;
  set(timehandle,'String',[num2str(dd) '/' num2str(mm) '/' num2str(yy)]);
  set(h1,'Enable','off');
% DRAW -----------------------------------------------
  drawnow
end
% Primary loop
for i = k+2:m
   for l = 1:N,
      j = i-1:i;
      if (i < ind(2,l)) & (i >= ind(1,l))
         set(head(l),'xdata',x(i,l),'ydata',y(i,l)) ;
      end
      if (i < ind(2,l)) & (i-1 >= ind(1,l))
         set(body(l),'xdata',x(j,l),'ydata',y(j,l)) ;
      end
      if (i <= ind(2,l)) & (i-k-1 >= ind(1,l))
         set(tail(l),'xdata',x(j-k,l),'ydata',y(j-k,l));
      elseif (i >= ind(2,l)) & (i <= ind(2,l) + k - 1 ),
         j = i-1:i;
         set(tail(l),'xdata',x(j-k,l),'ydata',y(j-k,l));
      end
   end
%  Calculate time --------------------------------------
   tim = starter + i - 1;
   timer = rd2dmy(tim);
   dd = timer(1);
   mm = timer(2);
   yy = timer(3);
   clear timer;
   set(timehandle,'String',[num2str(dd) '/' num2str(mm) '/' num2str(yy)]);
   set(h1,'Enable','off');
%  DRAW -----------------------------------------------
   drawnow
end
set(h1,'Enable','on');

