% ABI_SCALEDRADS    Loads ABI simulation dataset scaled rads
%
% These ABI simulation datasets contain MODIS radiances scaled to
% integers for compression studies.  This m-file script will load the data
% and convert it to radiances in MODIS units (Watts per square meter per
% steradian per micron) and then convert radiances into brightness
% temperatures (K).
%
% Band mapping from MODIS to ABI
%   Note: this is for compression studies only!  Not accurate spectral
%   representation!
% ABI	Center  MODIS   Center  ABI     ABI
% Band  Waveln  Band    Waveln  Res(KM) Bitdepth
% 01    0.47    3       0.47    1       10
% 02    0.64    1       0.659   0.5     12
% 03    0.865   2       0.865   1       10
% 04    1.378   5       1.240   2       10
% 05    1.61    6       1.640   1       10
% 06    2.25    7       2.130   2       10
% 07    3.9     22      3.96    2       14
% 08*   6.19    NA       NA     2       11
% 09    6.95    27      6.78    2       11
% 10    7.34    28      7.34    2       12
% 11    8.5     29      8.55    2       12
% 12    9.61    30      9.72    2       11
% 13*   10.35   NA       NA     2       12
% 14    11.2    31      11.0    2       12
% 15    12.3    32      12.0    2       12
% 16    13.3    33      13.4    2       11
%   % - These bands don't appear in the data files because they're repeated MODIS bands
%
% See also abi_sampledataset, modis_bright

% 03Nov2005 - CIMSS/SSEC/UW-Madison Experimental Software by Mat Gunshor

% Note that variables read in that pertain to 2km data have a 2 in them
% such as abiscaledrads2band##

% Parameters for the user to change:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Mountain Wave Case
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   basefilename='MOD021KM.A2004355.0515.004';
%   satmodis='terra';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CA Fire Case
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     basefilename='MYD021KM.A2003300.2055.003';
%     satmodis='aqua';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Convective Cloud Case
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   basefilename='MOD021KM.A2004329.1640.004';
%   satmodis='terra';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Daytime Clouds Case
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   basefilename='MOD021KM.A2002219.1550.004';
%   satmodis='terra';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Mexican Dust Case
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   basefilename='MYD021KM.A2004064.1955.004';
   satmodis='aqua';
    
    
% Bit depth
bitdepth=[10,12,10,10,10,10,14,11,11,12,12,11,12,12,12,11];
    
% MODIS band number mapping
modisbands(1:6)=[3,1,2,5,6,7];
modisbands(7:16)=[22,27,27,28,29,30,31,31,32,33];

%%%%%%%%%%%%%%%
% Load the data
%%%%%%%%%%%%%%%
eval(['load ',basefilename,'_to_ABIbitdepth_all.mat']);

% Check for visible channel
x=whos('abiscaledradsband01');
s=size(x);
if s(1)>0;
    vis='yes';
    firstband=1;
else
    vis='no';
    firstband=7;
end
clear x s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert scaled radiances to MODIS-style Radiances and Brightness Temps
%   Use the lookup table in scale_band## variables
%   Radiances are Watts per square meter per steradian per micron
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch vis
    case 'yes'
        for loop=[1:3,5]
            eval(['X',num2str(loop,'%02g'),'=[0:(2.^bitdepth(',num2str(loop),'))-1];'])
            eval(['abiradsband',num2str(loop,'%02g'),'=interp1(X',num2str(loop,'%02g'),',scale_band',num2str(loop,'%02g'),',abiscaledradsband',num2str(loop,'%02g'),');'])
        end
end
            
for loop=firstband:16
    if loop==8 || loop==13; continue; end                       % Skip bands 8 and 13 because they're not covered spectrally by MODIS
    if loop==1 || loop==2 || loop==3 || loop==5; continue; end  % Skip bands 1,2,3,5 because they're covered in loop above if visible exists
    eval(['X',num2str(loop,'%02g'),'=[0:(2.^bitdepth(',num2str(loop),'))-1];'])
    eval(['abirads2band',num2str(loop,'%02g'),'=interp1(X',num2str(loop,'%02g'),',scale_band',num2str(loop,'%02g'),',abiscaledrads2band',num2str(loop,'%02g'),');'])
    % Use the function modis_bright to convert to brightness temperatures for non-visible bands
    if loop==4 || loop==6; continue; end                        % Skip converting to temp for visible bands 4 and 6
    eval(['abitemps2band',num2str(loop,'%02g'),'=modis_bright(abirads2band',num2str(loop,'%02g'),',modisbands(',num2str(loop),'),1,satmodis);'])
end
% The interp1 command acts essentially as a lookup table with X## and 
% scale_band## for a given channel number.  (X## is X and scale_band## is Y
% and abiscaledrads2band## is XI).  modis_bright just converts MODIS
% radiances to brightness temperatures.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot a band of interest for this case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxlat=max(latitude2(:))+1.0;
minlat=min(latitude2(:))-1.0;
maxlon=max(longitude2(:))+1.0;
minlon=min(longitude2(:))-1.0;

switch basefilename
    case 'MOD021KM.A2004355.0515.004'
        %%% Mountain Wave Case %%%
        figure
        axesm('MapProjection','mercator','maplatlimit',[minlat maxlat],'maplonlimit',[minlon maxlon])
        set(gca,'FontWeight','Bold','FontSize',11)
        surfm(latitude2,longitude2,abitemps2band09);
        displaym(worldlo('POline'))
        h=displaym(usalo('stateborder'));
        set(h,'LineStyle','-');
        set(gca,'CLim',[220 260])   % Enhanced so you can see the mountain waves over Colorado
        hc=colorbar;
        set(hc,'FontSize',11,'FontWeight','Bold')
        tightmap
        titletxt=[basefilename,' ABI Band 8 (MODIS Band ',num2str(modisbands(8)),' - 6.78\mum)'];
        title(titletxt)
	cimsslabelvert2
        %eval(['print -dpng -r100 ',basefilename,'_band09.png'])
        %eval(['! convert ',basefilename,'_band09.png ',basefilename,'_band09.jpg'])
    case 'MYD021KM.A2003300.2055.003'
        %%% CA Fire Case %%%
        figure
        axesm('MapProjection','mercator','maplatlimit',[minlat maxlat],'maplonlimit',[minlon maxlon])
        set(gca,'FontWeight','Bold','FontSize',11)
        surfm(latitude2,longitude2,abitemps2band07);
        displaym(worldlo('POline'))
        h=displaym(usalo('stateborder'));
        set(h,'LineStyle','-');
        set(gca,'CLim',[270 340])   % Enhanced so you can see fires in southern CA
        hc=colorbar;
        set(hc,'FontSize',11,'FontWeight','Bold')
        tightmap
        titletxt=[basefilename,' ABI Band 7 (MODIS Band ',num2str(modisbands(7)),' - 3.96\mum)'];
        title(titletxt)
	cimsslabelvert2
        %eval(['print -dpng -r100 ',basefilename,'_band07.png'])
        %eval(['! convert ',basefilename,'_band07.png ',basefilename,'_band07.jpg'])
    case 'MOD021KM.A2004329.1640.004'
        %%% Convective Cloud Case %%%
        figure
        axesm('MapProjection','mercator','maplatlimit',[minlat maxlat],'maplonlimit',[minlon maxlon])
        set(gca,'FontWeight','Bold','FontSize',11)
        surfm(latitude2,longitude2,abitemps2band14);
        displaym(worldlo('POline'))
        h=displaym(usalo('stateborder'));
        set(h,'LineStyle','-');
        set(gca,'CLim',[200 290])   % Enhanced so you can see cold convective clouds
        hc=colorbar;
        set(hc,'FontSize',11,'FontWeight','Bold')
        tightmap
        titletxt=[basefilename,' ABI Band 14 (MODIS Band ',num2str(modisbands(14)),' - 11.2\mum)'];
        title(titletxt)
	cimsslabelvert2
        %eval(['print -dpng -r100 ',basefilename,'_band14.png'])
        %eval(['! convert ',basefilename,'_band14.png ',basefilename,'_band14.jpg'])
    case 'MOD021KM.A2002219.1550.004'
        %%% Daytime Clouds Case %%%
        figure
        axesm('MapProjection','mercator','maplatlimit',[minlat maxlat],'maplonlimit',[minlon maxlon])
        set(gca,'FontWeight','Bold','FontSize',11)
        surfm(latitude2,longitude2,abitemps2band14);
        displaym(worldlo('POline'))
        h=displaym(usalo('stateborder'));
        set(h,'LineStyle','-');
        set(gca,'CLim',[230 310])   % Enhanced so you can see clouds
        hc=colorbar;
        set(hc,'FontSize',11,'FontWeight','Bold')
        tightmap
        titletxt=[basefilename,' ABI Band 14 (MODIS Band ',num2str(modisbands(14)),' - 11.2\mum)'];
        title(titletxt)
	cimsslabelvert2
        eval(['print -dpng -r100 ',basefilename,'_band14.png'])
        eval(['! convert ',basefilename,'_band14.png ',basefilename,'_band14.jpg'])
    case 'MYD021KM.A2004064.1955.004'
        %%% Mexican Dust Case %%%
        figure
        axesm('MapProjection','mercator','maplatlimit',[minlat maxlat],'maplonlimit',[minlon maxlon])
        set(gca,'FontWeight','Bold','FontSize',11)
        surfm(latitude2,longitude2,abitemps2band11);
        displaym(worldlo('POline'))
        h=displaym(usalo('stateborder'));
        set(h,'LineStyle','-');
        set(gca,'CLim',[230 310])   % Enhanced so you can see the dust and ground details
        hc=colorbar;
        set(hc,'FontSize',11,'FontWeight','Bold')
        tightmap
        titletxt=[basefilename,' ABI Band 11 (MODIS Band ',num2str(modisbands(11)),' - 8.5\mum)'];
        title(titletxt)
	cimsslabelvert2
        eval(['print -dpng -r100 ',basefilename,'_band11.png'])
        eval(['! convert ',basefilename,'_band11.png ',basefilename,'_band11.jpg'])
end
