function display_airs
%%%%Program that plots desired AIRS channel with highlighted areas defined by
%%%%lat/lon points in spectra files saved from collocation softare
%%%%Plots are saved to be used in posters, etc.
% Read EOS fields for AIRS-MODIS collocation study
%
% Input
%    filename         : eos_hdf file name; string
%
% Output
%    data : cell list containing field values in double precision
%       radiances: [2378x90x135 double]
%   spectral_freq: [2378x1 double]
%        Latitude: [90x135 double]
%            Time: [90x135 double]
%       Longitude: [90x135 double]
%         scanang: [90x135 double]
%        landFrac: [90x135 double]
%             NeN: [2378x1 double]
%       satheight: [135x1 double]
%          satzen: [90x135 double]
%          solzen: [90x135 double]
%           topog: [90x135 double]
%
%%%%%%%%Variables Used Later - Edit as Desired%%%%%%%%%
%%AIRS filename
filename = ...
'/home/matg/matlab/collocation/20040303_g143/AIRS.2004.03.03.143.L1B.AIRS_Rad.v3.0.10.0.G04065052258.hdf';
%filename = ...
%'/data3/justins/AIRS_DATA/20050828/AIRS.2005.08.28.075.L1B.AIRS_Rad.v4.0.9.0.G05241143908.hdf';
%%Color bar min and max
cbar_min = 220;
cbar_max = 320;
%%Number of colors in colortable
numColors=264;
%%Map boundary colors on plot m is for magenta, k is for black
mapColor='k';
%%AIRS WN to display--notice it will get nearest AIRS channel, watch for coverage gaps!
userwn=890.0;
%%%%%%%%End Variables that need to be editted!%%%%%%%%

%[fn,dn] = uigetfile('AIRS*');
%filename = fullfile(dn,fn);

swathname  = 'L1B_AIRS_Science';

start  = '[]';
stride = '[]';
edge   = '[]';
		
clear flist

% Field name list
flist = {'radiances','spectral_freq','Latitude','Time',...
         'Longitude','scanang','landFrac','NeN',...
	 'satheight','satzen','solzen','topog'};

% Loop through fieldnames
%%%%%%Loading Data
disp('  Load AIRS ...')
for i=1:length(flist)
  s = ['[data.' flist{i} ']=read_eos_field(''' filename ''',''' swathname ''',''' flist{i} ''', start,stride,edge);'];
  eval(s)
end

airs.spectral_freq=data.spectral_freq;
nearest_airs=findnear(airs.spectral_freq,userwn);
airs.channel=nearest_airs(1,2);
airs.wavenumber = data.spectral_freq(airs.channel);
[nchannel, airs.pixels_across, airs.scanlines] = size(data.radiances);
airs.data = zeros(airs.pixels_across, airs.scanlines);
airs.data = rad2bt(airs.wavenumber, ...
                    squeeze(data.radiances(airs.channel, :, :)));
airs.lat = data.Latitude;
airs.lon = data.Longitude;
clear data; clear tmp;

%%%%%Plotting data
%% create latLonLimits
disp('  Plotting data...');
minLon = min(airs.lon(:));
minLat = min(airs.lat(:));
latLonLimits = [minLat max(airs.lat(:)) minLon max(airs.lon(:))];
displayAxis = axesm('mercator', 'maplatlimit', latLonLimits(1:2), ...
            'maplonlimit', latLonLimits(3:4), 'grid', 'off', 'gcolor', mapColor);
surfm(airs.lat,airs.lon,airs.data);
tightmap;
imageTitle = sprintf('AIRS Channel %4d (%7.2f cm^{-1}) Brightness Temp',...
             airs.channel, airs.wavenumber);
title(imageTitle);
load coast; plotm(lat, long, 1, mapColor);
load usalo; plotm(statelat,statelon,mapColor); plotm(uslat,uslon,mapColor);
plotm(gtlakelat,gtlakelon, 1, mapColor); 

tightmap;
%cmap = colormap(bone(numColors));
%colormap(flipud(cmap));
colormap(jet(numColors));
cbar_limits = [cbar_min cbar_max];
caxis(cbar_limits);
cb=colorbar('horizontal');
%%%Plot circles on regions of interest
file1='/home/tims/matlab/colocate/testspectra_dust_land.txt';
file2='/home/tims/matlab/colocate/testspectra_cleardesert3.txt';
file3='/home/tims/matlab/colocate/testspectra_water_clear.txt';
file4='/home/tims/matlab/colocate/testspectra_water_cloudy_water.txt';
file5='/home/tims/matlab/colocate/testspectra_water_cloudy_ice.txt';
plotROI(file1);
plotROI(file2);
plotROI(file3);
plotROI(file4);
plotROI(file5);
%legend

% Print Figure
imagename=['sampleROI_plot' '.png'];
eval(['print -dpng -r100 ',imagename]);
disp('  Complete!...');

function plotROI(filename)
%%%Function that plots the circles on the ROI specified in text file
plotLatLon = textread(filename,'%s',4,'headerlines',6);
lat = plotLatLon(3);
lon = plotLatLon(4);
lat = char(lat); lat = str2double(lat);
lon = char(lon); lon = str2double(lon);
h = plotm(lat,lon,'ko','Linewidth',2');
hold on;
set(h,'markersize',20);
