
% read MODIS (5-km) swath outputs
% The file name is designated by the start time of 5-minutes swath. 
% The file is in binary format. 
% The geo-location information and flux for every pixel within a 5-minutes swath 
% are stored in the order of latitude (270*406), longitude (270*406) 
% and radiative flux (270*406). 
% Usually, the output grid size is 270*406 pixels. 
% However, there are swaths that are of 270*408 pixels or 270*410 pixels.
% Those cases are denoted by the additional suffix "_408" or "_410" to
% the file names. 
addpath('/data/srb1/wchen/tools/m_fx//','-end')
clear all;clc; 
   modis = 'aqua';%'terra';% 
   if strcmp(modis,'terra') satStr='MOD';
   elseif strcmp(modis,'aqua')  satStr='MYD';
   end
   plttype = "hist";
   year  = 2015;
   yyyy = num2str(year, '%04d');
   load aqua_fluxes_mod5km_Jul2015_all.mat;
   %all sky flx @ cloudy pixels
   indflx = find(~isnan(flx));%keyboard;
   hist_plt(flxcld(indflx)-c6flx(indflx),'aqua_hist_diffflxcldflxc6',-600,300,10)
   keyboard
   diffflxc6 = flx(indflx)-c6flx(indflx);
   diffflxc7 = flx(indflx)-c7flx(indflx);
   diffc7c6  = c7flx(indflx)-c6flx(indflx);

   %all sky flx @ clear pixels
   indclc = find(isnan(flx));
   diffclcc6 = flxclc(indclc)-c6flx(indclc);
   diffclcc7 = flxclc(indclc)-c7flx(indclc);
   diffclcc7c6 = c7flx(indclc)-c6flx(indclc);

   disp("cloudy pixels")
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(c6flx(indflx),flx(indflx),1,0,'no');  
   statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
   disp(statistics)
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(c7flx(indflx),flx(indflx),1,0,'no');  
   statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
   disp(statistics)
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(c6flx(indflx),c7flx(indflx),1,0,'no');  
   statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
   disp(statistics)

   disp("clear pixels")
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(c6flx(indclc),flxclc(indclc),1,0,'no');  
   statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
   disp(statistics)
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(c7flx(indclc),flxclc(indclc),1,0,'no');  
   statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
   disp(statistics)
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(c6flx(indclc),c7flx(indclc),1,0,'no');  
   statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
   disp(statistics)
   %
   hist_plt(diffflxc6,'aqua_hist_diffflxc6',-300,600,10)
   hist_plt(diffflxc7,'aqua_hist_diffflxc7',-300,600,10)
   hist_plt(diffc7c6,'aqua_hist_diffflxc6c7',-35,35,5)
   %
   hist_plt(diffclcc6,'aqua_hist_diffclcc6',-600,400,10)
   hist_plt(diffclcc7,'aqua_hist_diffclcc7',-600,400,10)
   hist_plt(diffclcc7c6,'aqua_hist_diffclcc6c7',-35,35,5)

   function pdf_plt(x,y,figname)
      figure('DefaultTextFontName', 'Times', 'DefaultAxesFontName', 'Times')
      colormap(jet)
      pd=makedist('Normal','mu',nanmean(y),'sigma',nancov(y));
      ypdf=pdf(pd,x);
      plot(x,ypdf,'k','LineWidth',2.75)
      axis tight
      %set(gca,'ylim',[0 12e-3],'FontSize', 28,'linewidth',2.5);
      ylabel('PDF','FontName', 'Times', 'FontSize', 30)
      set(gcf,'unit','normalized','Position',[.2 .2 .6 .8]);%[x y width height]
      print(gcf,'-djpeg','-r300',figname)
    end
    function hist_plt(y,figname,limitmin,limitmax,limiti)
      figure('DefaultTextFontName', 'Times', 'DefaultAxesFontName', 'Times')
      colormap(jet)
      histogram(y,'BinLimits', [limitmin limitmax],'BinEdges' , [limitmin:limiti:limitmax],'EdgeColor','k','FaceColor','b');
      axis tight
      set(gca,'FontSize', 28,'linewidth',2.5);
      ylabel('Numbers','FontName', 'Times', 'FontSize', 28)
      set(gcf,'unit','normalized','Position',[.2 .2 .6 .8]);%[x y width height]
      print(gcf,'-djpeg','-r300',figname)
    end
   %{
   figure('DefaultTextFontName', 'Times', 'DefaultAxesFontName', 'Times')
   colormap(jet)
   if strcmpi(plttype,'pdf')
      x=-150:5:150;
      pd=makedist('Normal','mu',nanmean(diffflxc6),'sigma',nancov(diffflxc6));
      pdf=pdf(pd,x);
      plot(x,pdf,'k','LineWidth',2.75)
      axis tight
      set(gca,'ylim',[0 12e-3],'FontSize', 28,'linewidth',2.5);
      ylabel('PDF','FontName', 'Times', 'FontSize', 30)
      set(gcf,'unit','normalized','Position',[.2 .2 .6 .8]);%[x y width height]
      print(gcf,'-djpeg','-r300',['PDF_diff_ssda_cld_MOD_C07-C061_',modis])
   end
   if strcmpi(plttype,'hist')
      histogram(diffflxc6,'BinLimits', [-100 100],'BinEdges' , [-100:10:100],'EdgeColor','k','FaceColor','b');
      axis tight
      set(gca,'FontSize', 28,'linewidth',2.5);
      ylabel('Numbers','FontName', 'Times', 'FontSize', 28)
      set(gcf,'unit','normalized','Position',[.2 .2 .6 .8]);%[x y width height]
      print(gcf,'-djpeg','-r300',['aqua_hist_diff_ssda_clr_MOD_C07-C061_',modis])
      % h.FaceAlpha = 1 ;
   end
   %}
    %keyboard