% modified by CAS for OC402 cruise
% draw current position  of drifters (m x m) and track of the 
% central drifter; input file from forecast run. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              DRIFTER_OC391.M
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%run /home/ruoying/Software/opnml_matlab5/startup
% load stw2
    nn=loadgrid('../MESH/stw2');
    nl=nn;
    [nl.y,nl.x]=xy2ll(nn.x,nn.y);

%=================================================================
% Data
% load ../DRIFDAT/drift_2003.dat;
 load /whoinet/oc402_archive/data/shore/0525a/drift_rachels.dat;
 load /whoinet/oc402_archive/data/shore/0525a/drift_argos.dat;
 load /whoinet/oc402_archive/data/shore/0525a/drift_gps.dat;

 a=drift_2003;

 %[a(:,2),a(:,3)]=ll2xy(a(:,3),a(:,2)); % convert obs lon lat to x y

 ff=find(a(:,4)==277);
 a1=a(ff,:);

 ff=find(a(:,4)==87);
 a2=a(ff,:);
 
 ff=find(a(:,4)==213);
 a3=a(ff,:);
 
 ff=find(a(:,4)==629);
 a4=a(ff,:);

 ff=find(a(:,4)==630);
 a5=a(ff,:);

 ff=find(a(:,4)==631);
 a6=a(ff,:);

% load header from *.pth file

t_start=[                        
204839.05876 189044.49597  -1.00   2003    151.37500 %1 277
210916.88167 178201.15055 -14.00   2003    151.32153 %2  87
214366.45684 172728.99464 -14.00   2003    151.28056 %3 213
217898.16475 167602.85161 -14.00   2003    151.25000 %4 629
223893.85492 156452.61221 -14.00   2003    151.20486 %5 630
231203.66896 146569.21463 -14.00   2003    151.16250]%6 631


t_start(:,1:4)=[];

% get drifter paths

   gridname='/MESH/stw2';
   fn='qtc23';

   PTHFILE=['./fcast_archives/OC391_vn_' fn '/CRUISE_Q5.3/DRIFTERS']
%   PTHFILE =['/home/ruoying/MERHAB/OC391_archive/OC391_FC22/CRUISE_Q5.3/DRIFTERS'];
   %PTHFILE =['/home/ruoying/MERHAB/OC391_archive/FCAST_2.0_keston/fcast_archives/RH27/CRUISE_Q5.3/DRIFTERS'];

   [gridname,ndrog,ndts,tsec,pth]=read_pth(PTHFILE,4);
   size(pth);
   ndts=ans(1)/ndrog;
   xdr=pth(:,1);
   xdr=reshape(xdr,ndrog,ndts);
   xdr=xdr';
   ydr=pth(:,2);
   ydr=reshape(ydr,ndrog,ndts);
   ydr=ydr';
   zdr=pth(:,3);
   zdr=reshape(zdr,ndrog,ndts);
   zdr=zdr';
   tdr=pth(:,4);
   tdr=reshape(tdr,ndrog,ndts);
   tdr=tdr';
   [y,x]=xy2ll(xdr,ydr);
   clear pth
%=================================================================
% Observation
     i=6  % 631
     ff=find(x(:,i)==x(1,i));
     d6(:,1)=tdr(ff(length(ff)-1):length(tdr),i);
     d6(:,2)=x(ff(length(ff)-1):length(tdr),i);
     d6(:,3)=y(ff(length(ff)-1):length(tdr),i);
     d6(:,1)=a6(1,1)+(d6(:,1)-d6(1,1))/86400;
     i=5  % 630
     ff=find(x(:,i)==x(1,i));
     d5(:,1)=tdr(ff(length(ff)-1):length(tdr),i);
     d5(:,2)=x(ff(length(ff)-1):length(tdr),i);
     d5(:,3)=y(ff(length(ff)-1):length(tdr),i);
     d5(:,1)=a5(1,1)+(d5(:,1)-d5(1,1))/86400;
     i=4  % 629
     ff=find(x(:,i)==x(1,i));
     d4(:,1)=tdr(ff(length(ff)-1):length(tdr),i);
     d4(:,2)=x(ff(length(ff)-1):length(tdr),i);
     d4(:,3)=y(ff(length(ff)-1):length(tdr),i);
     d4(:,1)=a4(1,1)+(d4(:,1)-d4(1,1))/86400;
     i=3  % 213
     ff=find(x(:,i)==x(1,i));
     d3(:,1)=tdr(ff(length(ff)-1):length(tdr),i);
     d3(:,2)=x(ff(length(ff)-1):length(tdr),i);
     d3(:,3)=y(ff(length(ff)-1):length(tdr),i);
     d3(:,1)=a3(1,1)+(d3(:,1)-d3(1,1))/86400;
     i=2 % 87
     ff=find(x(:,i)==x(1,i));
     d2(:,1)=tdr(ff(length(ff)-1):length(tdr),i);
     d2(:,2)=x(ff(length(ff)-1):length(tdr),i);
     d2(:,3)=y(ff(length(ff)-1):length(tdr),i);
     d2(:,1)=a2(1,1)+(d2(:,1)-d2(1,1))/86400;
     i=1  % 277
     ff=find(x(:,i)==x(1,i));
     d1(:,1)=tdr(ff(length(ff)-1):length(tdr),i);
     d1(:,2)=x(ff(length(ff)-1):length(tdr),i);
     d1(:,3)=y(ff(length(ff)-1):length(tdr),i);
     d1(:,1)=a1(1,1)+(d1(:,1)-d1(1,1))/86400;


% try another one version:

      a1=sortrows(a1,1);
      a2=sortrows(a2,1);
      a3=sortrows(a3,1);
      a4=sortrows(a4,1);
      a5=sortrows(a5,1);
      a6=sortrows(a6,1);

% determine ( visually) the earlest and the latest time
 
  ii=find(diff(a1(:,1))==0); a1(ii,:)=[];    
  ax1=interp1(a1(:,1),a1(:,2),d1(:,1));
  ay1=interp1(a1(:,1),a1(:,3),d1(:,1));

  ii=find(diff(a2(:,1))==0); a2(ii,:)=[];
  ax2=interp1(a2(:,1),a2(:,2),d2(:,1));
  ay2=interp1(a2(:,1),a2(:,3),d2(:,1));

  ii=find(diff(a3(:,1))==0); a3(ii,:)=[];
  ax3=interp1(a3(:,1),a3(:,2),d3(:,1));
  ay3=interp1(a3(:,1),a3(:,3),d3(:,1));

  ii=find(diff(a4(:,1))==0); a4(ii,:)=[];
  ax4=interp1(a4(:,1),a4(:,2),d4(:,1));
  ay4=interp1(a4(:,1),a4(:,3),d4(:,1));

  ii=find(diff(a5(:,1))==0); a5(ii,:)=[];
  ax5=interp1(a5(:,1),a5(:,2),d5(:,1));
  ay5=interp1(a5(:,1),a5(:,3),d5(:,1));

  ii=find(diff(a6(:,1))==0); a6(ii,:)=[];
  ax6=interp1(a6(:,1),a6(:,2),d6(:,1));
  ay6=interp1(a6(:,1),a6(:,3),d6(:,1));

% start drawing 

    clf
    set(gcf,'defaultaxesfontsize',12,'color','w')
    whitebg('w')
    set(gcf,'color','w')

    hb=plotbnd(nl);
    set(hb,'Linewidth',2,'color','k')
    axis([-69.2000  -68.1000   43.3000   44.1000])
    hh=lcontour(nl,'z',[30:30:210])
    set(hh,'color',0.9*[1 1 1])%,'linewidth',0.2)
    xlabel('Longitude','FontWeight','bold','FontSize',14)
    ylabel('Latitude','FontWeight','bold','FontSize',14)
    title('Observed(dash) and simulated(solid)drifter paths ',...
    'FontWeight','bold','FontSize',14,'color','k');
    set(gcf,'paperposition',[-0.5 1.5 9.5 8]);

    hold on
    plot(x(1,1:6),y(1,1:6),'o','markersize',12,'markerfacecolor',...
     'r','markeredgecolor','k');

    plot(ax1,ay1,'b--','linewidth',3); plot(d1(:,2),d1(:,3),'b','linewidth',3);
    plot(ax2,ay2,'r--','linewidth',3); plot(d2(:,2),d2(:,3),'r','linewidth',3);
    plot(ax3,ay3,'g--','linewidth',3); plot(d3(:,2),d3(:,3),'g','linewidth',3);
    plot(ax4,ay4,'k--','linewidth',3); plot(d4(:,2),d4(:,3),'k','linewidth',3);
    plot(ax5,ay5,'m--','linewidth',3); plot(d5(:,2),d5(:,3),'m','linewidth',3);
    plot(ax6,ay6,'c--','linewidth',3); plot(d6(:,2),d6(:,3),'c','linewidth',3);

eval(['print -djpeg Q3_' fn '.jpg'])
