function mcs = trmmplots(hdfpath,hdffile)
% function mcs = trmmplots(hdfpath,hdffile)
% version 1.0
% 
% This matlab function and subfunctions plot some of the data from the
% tamuu / tamucc MCS database.  
% 
% The code is mainly a conversion from Chuntao Liu's IDL script. 
% Edit it to suit your needs.
% 
% caveats...
%  the color scaling is not as well done as I would like
% -- 
% EXAMPLE
% 
% hdfpath = '~/tamuu/tamuu.data/data/';
% hdffile = '86182_1.6_-1.4_2.00000_6.HDF'; 
%    note the above file produced from mcs database using on line tool
% plot the results...
% mcs = trmmplots(hdfpath,hdffile)
% 
% ------------------------------------------
% Dr Greg King                             &
% Institute of Marine Sciences (ICM-CSIC)  &
% Passeig Maritim Barceloneta, 37 - 49     &
% 08003 Barcelona, Spain                   &
% Questions to: gking@icm.csic.es          &
% ------------------------------------------

% for convenience, I put the path and filename here
%  for regular use, better to supply these as arguments in the 
%  call to trmmplots

% ----------------------------------------------------------------------
usemon2 = false;
[top1pos, top2pos, top3pos, bot1pos, bot2pos, bot3pos] = getfigtopbotpositions(usemon2);
% ----------------------------------------------------------------------
mcs = read_sds([hdfpath,hdffile]);
% ----------------------------------------------------------------------

id = strfind(hdffile,'_');
mcs.orbitS = hdffile(1:id(1)-1);
mcs.xlatS = hdffile(id(1)+1:id(2)-1);
mcs.xlonS = hdffile(id(2)+1:id(3)-1);
mcs.xlonF = str2num(mcs.xlonS);
mcs.xlatF = str2num(mcs.xlatS);
mcs.lonmin = min(mcs.PR_LON(:));
mcs.lonmax = max(mcs.PR_LON(:));
mcs.latmin = min(mcs.PR_LAT(:));
mcs.latmax = max(mcs.PR_LAT(:));
mcs.axis = [mcs.lonmin mcs.lonmax mcs.latmin mcs.latmax];

% calc mcstime
yearF = mean(mcs.YEAR(1));
monthF = mean(mcs.MONTH(1));
dayF = mean(mcs.DAY(1));
hrF = mean(double(mcs.HOUR) + double(mcs.MINUTE)/60);

mcs.yearS = num2str(yearF);
mcs.monthS = num2str(monthF,'%02d');
mcs.dayS = num2str(dayF,'%02d');
mcs.yymmddS = [mcs.yearS,'-',mcs.monthS,'-',mcs.dayS];
mcs.timeS = num2str(hrF,3);

fig.num = 1;
fig.pos = top1pos;
plot_PCT85(mcs,fig)

fig.num = 2;
fig.pos = bot1pos;
plot_nearsurfz(mcs,fig)

fig.num = 3;
fig.pos = top2pos;
plot_virs_ch4(mcs,fig)

fig.num = 4;
fig.pos = bot2pos;
plot_nearsurfrain(mcs,fig)

fig.num = 5;
fig.pos = top3pos;
plot_stratconv(mcs,fig)

fig.num = 6;
fig.pos = bot3pos;
plot_2a12rainrate(mcs,fig)

end

% ----------------------------------------------------------------------
% subfunctions follow
% ----------------------------------------------------------------------

function fdata = read_sds(hdffile)
% ; GENERIC PROGRAM TO READ SDS ARRAYS FROM AN HDF FILE
% ; ALL OF THE ARRAYS ARE PUT INTO AN OUTPUT STRUCTURE, fdata

fdata = struct([]);

fileinfo = hdfinfo(hdffile);
numsds = numel(fileinfo.SDS);
for i = 1:numsds,   
    sdsinfo = fileinfo.SDS(i);
    field = sdsinfo.Name;
    fdata(1).(field) = hdfread(sdsinfo);
end

end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function [top1pos, top2pos, top3pos, bot1pos, bot2pos, bot3pos] = getfigtopbotpositions(usemon2)
% usemon2 = false;
% [top1pos, top2pos, top3pos, bot1pos, bot2pos, bot3pos] = getfigtopbotpositions(usemon2)

if nargin == 0,
    usemon2 = false;
end

if usemon2,
    mon2 = [1300 0 0 0];
else
    mon2 = [0 0 0 0];
end

%  positions for figure windows on the screen
left = mon2 + [14   100   560   420];

top1pos = left+[0 600 0 0];
top2pos = left+[600 600 0 0];
top3pos = left+[1200 600 0 0];

bot1pos = left;
bot2pos = left+[600 0 0 0];
bot3pos = left+[1200 0 0 0];

end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function plot_PCT85(mcs,fig)

mytitle=[mcs.orbitS,' / 85GHz PCT (K) / ',mcs.yymmddS];
fn = fig.num;
fpos = fig.pos;
levels = 10;
clims = [160,300];

figure(fn)
set(gcf,'Position',fpos);

contour(mcs.TMI_LON,mcs.TMI_LAT,mcs.PCT85,levels)
hold on
plot([mcs.lonmin, mcs.lonmax],[mcs.xlatF mcs.xlatF],'k--')
plot([mcs.xlonF, mcs.xlonF],[mcs.latmin,mcs.latmax],'k--')
plot(mcs.PR_LON(:,1),mcs.PR_LAT(:,1),'k--','LineWidth',1)
plot(mcs.PR_LON(:,end),mcs.PR_LAT(:,end),'k--','LineWidth',1)
axis(mcs.axis);
caxis(clims);
cmap = colormap;
colormap(1-cmap);
colorbar
xlabel('LONGITUDE')
ylabel('LATITUDE')
title(mytitle)

end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function plot_nearsurfz(mcs,fig)

mytitle=[mcs.orbitS,' / 2A25 near surface reflectivity / ',mcs.yymmddS];
fn = fig.num;
fpos = fig.pos;

levels = (0:9)*6;

figure(fn)
set(gcf,'Position',fpos);

contour(mcs.PR_LON,mcs.PR_LAT,mcs.NEARSURFZ,levels);
hold on
plot([mcs.lonmin, mcs.lonmax],[mcs.xlatF mcs.xlatF],'k--')
plot([mcs.xlonF, mcs.xlonF],[mcs.latmin,mcs.latmax],'k--')
plot(mcs.PR_LON(:,1),mcs.PR_LAT(:,1),'k--','LineWidth',1)
plot(mcs.PR_LON(:,end),mcs.PR_LAT(:,end),'k--','LineWidth',1)
axis(mcs.axis);
colorbar
xlabel('LONGITUDE')
ylabel('LATITUDE')
title(mytitle)
end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function plot_virs_ch4(mcs,fig)
% on VIRS grid

mytitle=[mcs.orbitS,' / CH4 T_{B} (K) / ',mcs.yymmddS];
fn = fig.num;
fpos = fig.pos;
clims = [188,304];
levels = 10;

figure(fn)
set(gcf,'Position',fpos);

contour(mcs.VIRS_LON,mcs.VIRS_LAT,mcs.VIRS_SWATH_CH4,levels);
hold on
plot([mcs.lonmin, mcs.lonmax],[mcs.xlatF mcs.xlatF],'k--')
plot([mcs.xlonF, mcs.xlonF],[mcs.latmin,mcs.latmax],'k--')
plot(mcs.PR_LON(:,1),mcs.PR_LAT(:,1),'k--','LineWidth',1)
plot(mcs.PR_LON(:,end),mcs.PR_LAT(:,end),'k--','LineWidth',1)
axis(mcs.axis);
caxis(clims);
cmap = colormap;
colormap(1-cmap);
colorbar
xlabel('LONGITUDE')
ylabel('LATITUDE')
title(mytitle)

end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function plot_nearsurfrain(mcs,fig)

%        colors=[0,indgen(59)/60.0*45+15]
%        levelind=indgen(10)*6

mytitle=[mcs.orbitS,' / 2A25 near surface rainrate (mm/hr) / ',mcs.yymmddS];
fn = fig.num;
fpos = fig.pos;

levels = 0:60;

% The following shifts the data (where non-zero), presumeably to produce
% the green areas in Chuntau's plots.  We need a different colormap to
% reproduce that.
data = mcs.NEARSURFRAIN;
id = data > 0;
if sum(id(:)) > 0, 
    data(id) = data(id)/60*45 + 15;
end
% ----------------------------------------------------------------------

figure(fn)
set(gcf,'Position',fpos);

contour(mcs.PR_LON,mcs.PR_LAT,data,levels);
hold on
plot([mcs.lonmin, mcs.lonmax],[mcs.xlatF mcs.xlatF],'k--')
plot([mcs.xlonF, mcs.xlonF],[mcs.latmin,mcs.latmax],'k--')
plot(mcs.PR_LON(:,1),mcs.PR_LAT(:,1),'k--','LineWidth',1)
plot(mcs.PR_LON(:,end),mcs.PR_LAT(:,end),'k--','LineWidth',1)
axis(mcs.axis);
colorbar
xlabel('LONGITUDE')
ylabel('LATITUDE')
title(mytitle)
end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function plot_2a12rainrate(mcs,fig)

mytitle=[mcs.orbitS,' / 2A12 rainrate (mm/hr) / ',mcs.yymmddS];
fn = fig.num;
fpos = fig.pos;
clims = [0 60];

levels = 0:60;
% colors=(0:58)/60*45 + 15
%        cid=where(data gt 0,ccc)
%        if(ccc gt 0) then data[where(data gt 0)]=dtmp[where(data gt 0)]/60.0*45.0+15

% The following shifts the data (where non-zero), presumeably to produce
% the green areas in Chuntau's plots.  We need a different colormap to
% reproduce that.
data = mcs.TMIRAIN;
id = data > 0;
if sum(id(:)) > 0, 
    data(id) = data(id)/60*45 + 15;
end

figure(fn)
set(gcf,'Position',fpos);

contour(mcs.TMI_LON,mcs.TMI_LAT,data,levels)
hold on
plot([mcs.lonmin, mcs.lonmax],[mcs.xlatF mcs.xlatF],'k--')
plot([mcs.xlonF, mcs.xlonF],[mcs.latmin,mcs.latmax],'k--')
plot(mcs.PR_LON(:,1),mcs.PR_LAT(:,1),'k--','LineWidth',1)
plot(mcs.PR_LON(:,end),mcs.PR_LAT(:,end),'k--','LineWidth',1)
axis(mcs.axis);
caxis(clims)
colorbar
xlabel('LONGITUDE')
ylabel('LATITUDE')
title(mytitle)

end
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

function plot_stratconv(mcs,fig)
% PR 2A23 Rain type 
% stratiform: 100-153
% convective: 200-293

% IDL...
% colorbar,
% levels,
% colors,
% levelind=[1,2,3]
% tics=['Norain','N/A','Strat','Conv']

mytitle=[mcs.orbitS,' / 2A23 Rain Types / ',mcs.yymmddS];
fn = fig.num;
fpos = fig.pos;
clims = [0 55];

dat = mcs.PR_2A23_RAINTYPE;
data=NaN(size(dat));
% stratiform: 100-153
id = dat >= 100 & dat <= 153;
if(sum(id(:)) >= 0),
    data(id)=30;
end
% convective: 200-293
id = dat >= 200 & dat <= 293;
if(sum(id(:)) >= 0),
    data(id)=45;
end
levels=[0,29,43,51];
cmap = colormap;
%        colors=[0,30,45]
%        levelind=[0]
figure(fn)
set(gcf,'Position',fpos);

contourf(mcs.PR_LON,mcs.PR_LAT,data,levels);
hold on
plot([mcs.lonmin, mcs.lonmax],[mcs.xlatF mcs.xlatF],'k--')
plot([mcs.xlonF, mcs.xlonF],[mcs.latmin,mcs.latmax],'k--')
plot(mcs.PR_LON(:,1),mcs.PR_LAT(:,1),'k--','LineWidth',1)
plot(mcs.PR_LON(:,end),mcs.PR_LAT(:,end),'k--','LineWidth',1)
axis(mcs.axis);
caxis(clims)
colorbar
xlabel('LONGITUDE')
ylabel('LATITUDE')
title(mytitle)

end


