
% 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 = 'terra';% 'aqua';%
   plttype = "hist";
   year  = 2015;
   yyyy = num2str(year, '%04d');
   inpath0 = ['/data/srb1/wchen/WORK/MODIS/SRB_SW/case_5km_cp_run/input/'];
   dataPath1 = ['/data/srb1/wchen/WORK/MODIS/SRB_SW/case_5km_cp_run/results/'];
   filepath1 = [dataPath1,'/',yyyy,'/',modis,'/'];

   
   dataPath2 = ['/data/srb9/5kmMODIS/Results/'];
   filepath2 = [dataPath2,'/',yyyy,'/'];
   cnt = 0
   for doi=182:212
       doy   = num2str(doi, '%03d');
       datev = datevec(doy2date(doi,year))
       fptn  = [filepath1,'modis.',modis,'.',yyyy,doy,'.*.dat'];
       fList = dir(fptn);
       nFl   = length(fList);
   
      if nFl==0 
         disp('---read_MODIS_swath_new(): Empty directory!')       
      else
         for i=1:nFl
            cnt = cnt+1;
            hhmm  = fList(i).name(end-7:end-4); 
            hrmin = str2num(hhmm);
            [ds2] = read_MODIS_swath_v40(...
                     modis, filepath1, fList(i).name);
            data2 = ds2.varM;
           
            [ds1] = read_MODIS_swath_v40(...
                     modis, filepath2, fList(i).name);
            data1 = ds1.varM;
            [cldLat,cldLon,...
                     ctp, wfrc,wcod,wcre, ifrc,icod,icre, tfrc ] =...
                    get_cloud(...
                              year, datev(2), datev(3), hrmin, doi, modis, inpath0); %0:clear; >0: cloud
            indxcld = find(tfrc~=0);
            indxclr = find(tfrc==0);
            if i==1 && doi==182
               % diff = data2-data1;
               % diffcld = diff(tfrc~=0);
               % diffcld = diffcld(:)
               % diffclr = diff(tfrc==0);
               % diffclr = diffclr(:)
               data2cld = data2(indxcld);
               data2cld = data2cld(:);
               data2clr = data2(indxclr);
               data2clr = data2clr(:);

               data1cld = data1(indxcld);
               data1cld = data1cld(:);
               data1clr = data1(indxclr);
               data1clr = data1clr(:);
            else
               % diff = data2-data1;
               % tmp = diff(tfrc~=0);
               % diffcld = [diffcld;tmp(:)];
               % clear tmp
               % tmp = diff(tfrc==0);
               % diffclr = [diffclr;tmp(:)];
               % clear tmp
               tmp = data2(indxcld);
               data2cld = [data2cld;tmp(:)];
               clear tmp
               tmp = data2(indxclr);
               data2clr = [data2clr;tmp(:)];
               clear tmp

               tmp = data1(indxcld);
               data1cld = [data1cld;tmp(:)];
               clear tmp
               tmp = data1(indxclr);
               data1clr = [data1clr;tmp(:)];
               clear tmp;
            end
            % clear diff;

            % figure()
            % imagesc(ds2.varM)
            % figure()
            % imagesc(ds1.varM)
            % figure()
            % imagesc(tfrc)
            % keyboard
         end%nFL
      end %if nFL==0
	end % doi
   % disp([min(diffcld(:)),max(diffclr(:))])
   % x=-100:100;
   % pd=makedist('Normal','mu',nanmean(diffcld),'sigma',nancov(diffcld));
   % pdf=pdf(pd,x);
   % plot(x,pdf)
   disp("cloud")
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(data1cld,data2cld,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")
   [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(data1clr,data2clr,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(['cnt = ',num2str(cnt)])
   diffcld = data2cld-data1cld;
   diffclr = data2clr-data1clr;
   figure('DefaultTextFontName', 'Times', 'DefaultAxesFontName', 'Times')
   colormap(jet)

   if strcmpi(plttype,'pdf')
      x=-150:5:150;
      pd=makedist('Normal','mu',nanmean(diffcld),'sigma',nancov(diffcld));
      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(diffclr,'BinLimits', [-35 35],'BinEdges' , [-35:1:35],'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',['HIST_diff_ssda_clr_MOD_C07-C061_',modis])
      % h.FaceAlpha = 1 ;
   end
    %keyboard