%
% matlab script that illustrates extraction and plotting of the
% BYU-created ASCAT/TRMM PR/ECMWF collocated dataset
%
% written 13 Jan 2022 by DGL at BYU 
%

clf;

% input directory
input_dir='ATE_col/';

% read names in input director
D=dir([input_dir,'ATE*.dat']);

% read each file and plot information
nD=length(D);
for n=1:nD
  
  %file name
  fname=[D(n).folder,'/',D(n).name];

  %load files
  vals=load(fname);

  % definition of columns in vals
  %  time difference in seconds
  % year/month/day/hour/m/s
  % lat/lon values for TRMM PR rain measurements (deg)
  % lat/lon values for SIR-interpolated ASCAT sigma0 measurements (deg)
  % forward, middle, aft values of ASCAT sigma0 (dB)
  % forward, middle, aft values of ASCAT azimuth observation angles (deg)
  % forward, middle, aft values of ASCAT incidence obseration angles (deg)
  % TRMM PR reported rain (mm/hr) 
  % u and v components of the from bilineary interpolated from ECMWF

  % select data within +/-'thres' collocation time
  thres=15*60;  % time in seconds
  ind=find(abs(vals(:,1)) < thres);
  
  % extract data within time interval for plotting
  tdif=vals(ind,1);
  tlat=vals(ind,8);
  tlon=vals(ind,9);
  alat=vals(ind,10);
  alon=vals(ind,11);
  sigf=vals(ind,12);
  sigm=vals(ind,13);
  siga=vals(ind,14);
  train=vals(ind,21);
  espd=sqrt(vals(ind,22).^2+vals(ind,23).^2); % ECMWF wind speed
  edir=atan2d(-vals(ind,22),-vals(ind,23)); % ECMWF meteorological wind direction 

  % data start date/time (first data entry)
  datestr=sprintf('%4d-%2.2d-%2.2d:%2.2d:%2.2d:%2.2d',vals(1,2:7));
  
  % prepare to plot
  t20=18+zeros(size(tdif));

  sigr=double([min([min(sigf(:)),min(sigm(:)),min(siga(:))]),max([max(sigf(:)),max(sigm(:)),max(siga(:))])]);

  % plot arrays
  
  myfigure(1)
  scatter3(tlon,tlat,tdif,t20,tdif,'filled');
  view(0,90);
  colorbar;
  xlabel('longitude (deg)')
  ylabel('latitude (deg)')
  title(sprintf('%s (%d) time difference (s)',datestr,thres));

  if 0 % no need to plot lat and lon arrays
    myfigure(10)
    scatter3(tlon,tlat,tlat,t20,tlat,'filled');
    view(0,90);
    colorbar;
    xlabel('longitude (deg)')
    ylabel('latitude (deg)')
    title(sprintf('%s TRMM latitude (deg)',datestr));

    myfigure(11)
    scatter3(tlon,tlat,tlon,t20,tlon,'filled');
    view(0,90);
    colorbar;
    xlabel('longitude (deg)')
    ylabel('latitude (deg)')
    title(sprintf('%s TRMM longitude (deg)',datestr));
  end
  
  myfigure(2)
  scatter3(tlon,tlat,espd,t20,espd,'filled');
  view(0,90);
  colorbar;
  xlabel('longitude (deg)')
  ylabel('latitude (deg)')
  title(sprintf('%s ECMWF surface wind speed (m/s)',datestr));
  
  myfigure(3)
  scatter3(alon,alat,sigf,t20,sigf,'filled');
  caxis(sigr);
  view(0,90);
  colorbar;
  xlabel('longitude (deg)')
  ylabel('latitude (deg)')
  title(sprintf('%s ASCAT fore sigma-0 (dB)',datestr));

  myfigure(4)
  scatter3(alon,alat,sigm,t20,sigm,'filled');
  caxis(sigr);
  view(0,90);
  colorbar;
  xlabel('longitude (deg)')
  ylabel('latitude (deg)')
  title(sprintf('%s ASCAT mid sigma-0 (dB)',datestr));

  myfigure(5)
  scatter3(alon,alat,siga,t20,siga,'filled');
  caxis(sigr);
  view(0,90);
  colorbar;
  xlabel('longitude (deg)')
  ylabel('latitude (deg)')
  title(sprintf('%s ASCAT aft sigma-0 (dB)',datestr));
  
  myfigure(6)
  scatter3(tlon,tlat,train,t20,train,'filled');
  view(0,90);
  colorbar;
  xlabel('longitude (deg)')
  ylabel('latitude (deg)')
  title(sprintf('%s TRMM PR rain (mm/hr)',datestr));

  fprintf('%s (return to continue)',datestr);
  pause;
  fprintf('\n');
  
end