%
% function [x,y,mag]=provecxy(u,v,yd,[tint2],[lat1],[lon1],[col]);
%
%  to caculate x & y displacements given u & v velocity 
%  and time in yeardays
%
%  If "yd's" length is one, then assumes that is a regular
%  spaced time interval in seconds
%
%  u & v are expected in cm/sec
%  yd is the yearday of observations
%    (or, if its length is one, regular space time in seconds)
%  tint2 is interval to annotate (also in seconds)
%    (if tint2==30*24*3600, annotate day1 of each mth)
%    (if this is not given no plot will be made)
%  x & y are output in Kms unless lat1 and lon1 are included in input
%       in which case x & y are lon, lat
%  lat1 & lon1 (optional) are the starting point in geographic space
%     if they are not specified the starting point is (0,0) in kms
%     (if just "lat1" is specified then this is the y-offset)
%  col is color string

function [x,y,mag]=provecxy(u,v,yd,tint2,lat1,lon1,col);

% initialize some variables
  x(1)=0.0;
  y(1)=0.0;
  sumx=0.0;
  sumy=0.0;

% take only the non-NaN's
  idnn=find(isnan(u)~=1);
  u=u(idnn);
  v=v(idnn);
  yd=yd(idnn);

% figure out the time interval for each observation
  if length(yd)==1;
    tintv=yd*ones(length(u)-1);
  else
% for cases of unequal time use the "yd"
    tintv=diff(yd)*24*3600;
    tint=mean(tintv);
  end;


% for each time step calculate distance covered,
% divide by 100000 to convert from cm to km,
%  and then sum distance
% note that we are using the average velocity between obs.

   for i=2:length(u);
      if tintv(i-1)>10*tint;
        disp('WARNING: time interval >10 the mean');
        yd(i)
      end;
      tx(i)=tintv(i-1)*(u(i-1)+u(i))/2/100000;
      ty(i)=tintv(i-1)*(v(i-1)+v(i))/2/100000;
      x(i)=sumnnan(tx(1:i));
      y(i)=sumnnan(ty(1:i));
      tmag(i)=sqrt(tx(i)^2+ty(i)^2);
      mag(i)=sum(tmag(1:i));
   end;

%************************* PLOT IN XY SPACE *******************

% offset y if requested
  if nargin==5;
%    cc=input('color?','s');
    cc='w';
    y=y+lat1;
  else;
    cc='w';
  end;

% now plot
  if nargin<6;
    figure(1);
    plot(x,y,cc);
    hold on;
  end;

% annotate plot only when the "tint2" interval is specified
if nargin>=4;

%  calculate a displacement of yearday label
  dyd=(max(y)-min(y))/20;

  if tint2~=30*24*3600;

%   add regularly timed labels
    for i=1:floor(tint2/tint):length(x);
      plot(x(i),y(i),'.','color','r','markersize',24)
      date=calday(95,yd(i));
      h2=text(x(i),y(i)+dyd,[num2str(date(2)),'/',num2str(date(3))]);
%     case where hourly lables needed
      if tint2<3600*24;
        h3=text(x(i),y(i)+dyd*2,[num2str(date(4)*100+date(5)),'GMT']);
        set(h3,'fontsize',10)
      end;
    end;

  else;
%   add monthly labels
    cmth=['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];

    for i=1:length(x);
      date=calday(95,yd(i));
      if date(3)==1&date(4)==0;
        plot(x(i),y(i),'.','color','w','markersize',24)
        h2=text(x(i),y(i),cmth(date(2)),...
         'horizontalalignment','center',...
         'verticalalignment','middle',...
         'fontsize',6,'color','k');
      end;
    end;      
  end;
      
      set(h2,'horizontalalignment','center');
else;
     plot(x(1),y(1),'.','color','r','markersize',24)
    h2=text(x(1),y(1),'start');
end;

% fix axis
  axis('equal');
  ylabel('Kilometers (N/S)');
  xlabel('Kilometers (E/W)');

%***********************   CALCULATE LAT/LON *****************

% now calulcate lat lon results if asked for
if nargin>=6;
  [latd]=dist(lat1,lon1,lat1-0.5,lon1); 
  [lond]=dist(lat1,lon1,lat1,lon1+0.5); 
  llratio=latd/lond

% figure lat step to next "adjusted" station
  latstep=y/1.81/60;

% figure lon step by using the lon to lat distance ratio 
  lonstep=(x/1.81/60)*llratio;

% figure new lat lon time series
  for k=1:length(x);
    nlat(k)=lat1+latstep(k);
    nlon(k)=lon1+lonstep(k);
  end;

% now plot lat lon in another window if posted "tint2" is given

%************************  PLOT LAT/LON  *********************

  if nargin>=6;
%    figure(2);
    if nargin~=7;col='g';end;

    disp(['col=',col])
    plot(nlon,nlat,'-','color',col,'linewidth',2);

 
%   calculate a displacement of yearday label
    hold on;
    dyd=(max(nlat)-min(nlat))/20;
    for i=1:floor(tint2/tint):length(nlon);
      plot(nlon(i),nlat(i),'.','color','r','markersize',24)
%     h1=text(nlon(i),nlat(i),'.','color','r','fontsize',24)
%     set(h1,'verticalalignment','middle','horizontalalignment','center');
     date=calday(95,yd(i));
     h2=text(nlon(i),nlat(i)+dyd,[num2str(date(2)),'/',num2str(date(3))]);
%    case where hourly lables needed
     if tint2<3600*24;
       h3=text(nlon(i)+dyd,nlat(i)+dyd,[num2str(date(4)*100+date(5)),'GMT']);
       set(h3,'fontsize',10)
     end;
      
     set(h2,'horizontalalignment','center');
    end;


%   fix axis
    set(gca,'aspectratio',[1 .78]);
    ylabel('Latitude (N/S)');
    xlabel('Longitude (E/W)');
  end;
  x=nlon;
  y=nlon;
end;

 