% 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;
 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='qtc24'

for N=1:3
   PTHFILE=['./fcast_archives/OC391_vn_' fn '/CRUISE_Q5.' num2str(N) '/DRIFTERS']
   %PTHFILE =['/home/ruoying/MERHAB/OC391_archive/OC391_FC22/CRUISE_Q5.' num2str(N) '/DRIFTERS'];
 %PTHFILE =['/home/ruoying/MERHAB/OC391_archive/FCAST_2.0_keston/fcast_archives/RH27/CRUISE_Q5.' num2str(N) '/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';
   x=xdr; y=ydr;
   clear pth

% 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',[.4 .4 .4],'linewidth',0.2)
%    xlabel('Longitude','FontWeight','bold','FontSize',14)
%    ylabel('Latitude','FontWeight','bold','FontSize',14)
%    title('Observed and simulated (OC391-FC22) drifter paths ',...
%    'FontWeight','bold','FontSize',14,'color','k');
%    set(gcf,'paperposition',[-0.5 1.5 9.5 8]);

%    hold on
%    plot(x(1,6:11),y(1,6:11),'o','markersize',8,'markerfacecolor',...
%     'c','markeredgecolor','r')

%=================================================================
% 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
 
      nmin=min([d1(1,1) d2(1,1) d3(1,1) d4(1,1) d5(1,1) d6(1,1)]);
      nmax=max([d1(length(d1),1) d2(length(d2),1) d3(length(d3),1) ...
                d4(length(d4),1) d5(length(d5),1) d6(length(d6),1)]);
 
  ii=find(diff(a1(:,1))==0); a1(ii,:)=[];    
  cx=interp1(a1(:,1),a1(:,2),d1(:,1));
  cy=interp1(a1(:,1),a1(:,3),d1(:,1));
  c1=sqrt((cx-d1(:,2)).^2+(cy-d1(:,3)).^2);

  ii=find(diff(a2(:,1))==0); a2(ii,:)=[];
  cx=interp1(a2(:,1),a2(:,2),d2(:,1));
  cy=interp1(a2(:,1),a2(:,3),d2(:,1));
  c2=sqrt((cx-d2(:,2)).^2+(cy-d2(:,3)).^2);

  ii=find(diff(a3(:,1))==0); a3(ii,:)=[];
  cx=interp1(a3(:,1),a3(:,2),d3(:,1));
  cy=interp1(a3(:,1),a3(:,3),d3(:,1));
  c3=sqrt((cx-d3(:,2)).^2+(cy-d3(:,3)).^2);

  ii=find(diff(a4(:,1))==0); a4(ii,:)=[];
  cx=interp1(a4(:,1),a4(:,2),d4(:,1));
  cy=interp1(a4(:,1),a4(:,3),d4(:,1));
  c4=sqrt((cx-d4(:,2)).^2+(cy-d4(:,3)).^2);

  ii=find(diff(a5(:,1))==0); a5(ii,:)=[];
  cx=interp1(a5(:,1),a5(:,2),d5(:,1));
  cy=interp1(a5(:,1),a5(:,3),d5(:,1));
  c5=sqrt((cx-d5(:,2)).^2+(cy-d5(:,3)).^2);

  ii=find(diff(a6(:,1))==0); a6(ii,:)=[];
  cx=interp1(a6(:,1),a6(:,2),d6(:,1));
  cy=interp1(a6(:,1),a6(:,3),d6(:,1));
  c6=sqrt((cx-d6(:,2)).^2+(cy-d6(:,3)).^2);

  figure(N)
  set(gcf,'paperposition',[0.75 3 7 5]);
  set(gcf,'color','w');


  plot(d1(:,1),c1*10e-4,'b-.','linewidth',1);
  axis([nmin nmax 0 45])
  grid
  gtext('#1')
  hold on
  plot(d2(:,1),c2*10e-4,'k-','linewidth',1);
  gtext('#2')
  plot(d3(:,1),c3*10e-4,'k:','linewidth',2);
  gtext('#3')
  plot(d4(:,1),c4*10e-4,'g-','linewidth',1);
  gtext('#4')
  plot(d5(:,1),c5*10e-4,'m-','linewidth',1);
  gtext('#5')
  plot(d6(:,1),c6*10e-4,'r--','linewidth',1);
  gtext('#6')
  hold off
  xlabel('Time [days]');
  ylabel('Distance [km]');
  title('Distance between observed and simulated drifters at the end of the cruise');

  ds(1)=c1(length(c1))*10e-4/(max(d1(:,1))-min(d1(:,1)));
  ds(2)=c2(length(c2))*10e-4/(max(d2(:,1))-min(d2(:,1)));
  ds(3)=c3(length(c3))*10e-4/(max(d3(:,1))-min(d3(:,1)));
  ds(4)=c4(length(c4))*10e-4/(max(d4(:,1))-min(d4(:,1)));
  ds(5)=c5(length(c5))*10e-4/(max(d5(:,1))-min(d5(:,1)));
  ds(6)=c6(length(c6))*10e-4/(max(d6(:,1))-min(d6(:,1)));

  dis=[c1(length(c1))*10e-4 c2(length(c2))*10e-4 c3(length(c3))*10e-4 ...
       c4(length(c4))*10e-4 c5(length(c5))*10e-4 c6(length(c6))*10e-4];
  
  day=[max(d1(:,1))-min(d1(:,1)) max(d2(:,1))-min(d2(:,1)) ...
       max(d3(:,1))-min(d3(:,1)) max(d4(:,1))-min(d4(:,1)) ...
       max(d5(:,1))-min(d5(:,1)) max(d6(:,1))-min(d6(:,1))];

  ds=mean(dis./day);

  for k=1:6
   text(0.05,0.95-(k-1)*0.04,['Drouge:' num2str(k) ' Dist:' num2str(dis(k)) ' Days:' num2str(day(k))],'unit','norm')
  end
  text(0.05,0.95-6*0.04,['Average Divergence :' num2str(ds) '(km/day)'],'unit','norm','fontweight','bold')

 eval(['print -djpeg ' fn '_' num2str(N) '.jpg']); close
end


