% compute statistics on model/data time series pairs
% step 1.  >> hindcast_gom3_uv (save output in mat file)
% step 2.  >> hindcast_stats  (load mat file from Step 1)
% step 3.  >> hindcast_plot
%load fvcom_gom3_uv.mat
load month.mat
isub=3;
for i=1:length(lon_mod),
  [tidestruc_obs{i}]=t_tide(wobs{i},'start_time',jd{i}(1),'error','wboot');
  [tidestruc_mod{i}]=t_tide(wmod{i},'start_time',jd{i}(1),'error','wboot');
  
  wobslp{i}=pl33tn(wobs{i},1,33); % low-pass filter to 33 hours (remove tides)
  wmodlp{i}=pl33tn(wmod{i},1,33);
  wobslp{i}=wobslp{i}(1:isub:end);
  wmodlp{i}=wmodlp{i}(1:isub:end);
  jdlp{i}=jd{i}(1:isub:end);
  %[amp,theta,trans]=comcor(wilp,wmodlp); %complex correlation of low-pass data
  %[th,maj,min,wr]=princax(wmodlp);
  [mean1(i), std1(i), mean2(i), std2(i), corr(i), transfn(i)]=stats_v(wobslp{i},wmodlp{i},1);
  
  % plot up data & model time series
  if(0)
    figure;
    set(gcf,'color','white','pos',[10 100 1200 500]);
    ii=1:200;
    plot(jdi(ii),real([wobs{i}(ii) wmod{i}(ii)]));datetick
    legend('DATA','MODEL')
    toc
    shg
  end
  
  if(0)
    subplot(121);
    plot(jdi,umod,'-k+', jdi,vmod, '-r+');
    hold on
    plot(jdi, ui, '-b^', jdi, vi, '-g^');
    hold off
    datetick;
    grid
    legend('Model U Velocity','Model V Velocity', 'Data U Velocity', 'Data V Velocity');
    %     iext=strfind(nci.location,'dodsC/'); %edited by acrosby
    ext=nci.location((strfind(nci.location,'dodsC/')+6):end); %edited by acrosby
    title(sprintf('Temp comparison at %6.2f,%6.2f, Depth=%5.1f',...
      loni,lati,depi),'interpreter','none');
    
    subplot(122);
    %     xlim([-74 -65]);
    %     ylim([38 45]);
    hold on
    pcolor(double(subs.grid.lon),double(subs.grid.lat),double(subs.data));% changed acrosby
    
    shading flat
    scatter(loni,lati,'MarkerEdgeColor','black','Marker','o','LineWidth',1.6,'SizeData',65);
    %     dasp(42);
    caxis([-300 0])
    hold off
    axis xy;
    title(ext,'interpreter','none');
    set(gca,'tickdir','out');
    drawnow
    figure(gcf); %same as shg?
  end
end
%%  t_tide output for vector input
%    [fmaj,emaj,fmin,emin,finc,einc,pha,epha]
%        1    2    3    4    5    6   7    8
for i=1:length(lon_mod)
  ind=strmatch('M2  ',tidestruc_mod{i}.name)
  m2maj_mod(i)=tidestruc_mod{i}.tidecon(ind,1);
  m2pha_mod(i)=tidestruc_mod{i}.tidecon(ind,7);
  m2inc_mod(i)=tidestruc_mod{i}.tidecon(ind,5);
  m2inc_obs(i)=tidestruc_obs{i}.tidecon(ind,5);
  m2maj_obs(i)=tidestruc_obs{i}.tidecon(ind,1);
  m2pha_obs(i)=tidestruc_obs{i}.tidecon(ind,7);
end
ii=find(m2pha_mod>180);
m2pha_mod(ii)=m2pha_mod(ii)-360;
ii=find(m2pha_obs>180);
m2pha_obs(ii)=m2pha_obs(ii)-360;