clc;
clear all;
close all;

% This script reads in swaths of HIRAD TB (or excess TB) data and applies
% a cross-track filter.  Technique: The along-track mean is calculated for
% each cross-track position.  That cross-track profile is filtered by first
% performing a least-squares fit to a 4th-order polynomial, and then
% replacing the original along-track mean with that 4th-order polynomial.
% The result is that the artificial along-track stripes in the data are
% removed.  These stripes are evidence that the Level 1 calibration of the
% HIRAD data are still incomplete.  They are non-physical and hence can be
% removed for display and analysis purposes.  It is expected that future,
% improved Level-1 datasets will not require this filtering.

% It is presumed that the user has the Matlab mapping toolbox.  If not, the
% commands axesm, pcolorm, etc., can be replaced by their corresponding
% simpler commands.

% Script developer: Tim Miller, NASA/MSFC, 256-961-7882,
% tim.miller@nasa.gov.  This script is provided for HIRAD data users to see
% exactly how the filtering was performed.

load ../../for_release/Karl_16/data/KARL16_leg2.mat;
%load ../../for_release/Karl_16/data/KARL16_leg10.mat;
%load leg10/TbDATA_Karl_leg10_legV6b.mat;
cname='KARL16-leg2';
cname1='grip_hirad_20100916_';  % Change this date for Earl

tb5=leg.delTB5xx;

lat=leg.swath_lat;
lon=leg.swath_lon;

nalong=size(lat,1);
ncross=size(lat,2);
nskip1=0;   % These parameters are used to process partial legs.
nskip2=0;
na1=nskip1+1;
na2=nalong-nskip2;
width=65;  % Specified parameter.  Use 70 for Earl, 65 for Karl.
ncskip=75-width;
nc1=ncskip+1;
nc2=ncross-ncskip;

tb5=tb5(na1:na2,nc1:nc2);
lat=lat(na1:na2,nc1:nc2);
lon=lon(na1:na2,nc1:nc2);
nalong=size(lat,1);
ncross=size(lat,2);


latmin=min(min(lat))-.2;
lonmin=min(min(lon))-.2;
latmax=max(max(lat))+.2;
lonmax=max(max(lon))+.2;


% S=gshhs('gshhs_h.b',[latmin latmax],[lonmin lonmax]);  % enables plot of coast

scrsz = get(0,'ScreenSize');
x1=1;
y1=.2*scrsz(4);
xwidth=.5*scrsz(3);
ywidth=.7*scrsz(4);
x2=x1+xwidth;  y2=y1+ywidth;
figure('Position',[x1 y1 xwidth ywidth])

axesm('MapProjection','eqdcylin','maplatlimit',[latmin latmax],'maplonlimit',[lonmin lonmax],...
       'grid','on','mlinelocation',.5,'plinelocation',.5,'meridianlabel','on','mlabelround',-1,...
    'labelformat','none','parallellabel','on','plabelround',-1,'flinewidth',1,'fontsize',11,...
    'fontweight','normal');
mlabel('south');
pcolorm (lat, lon, tb5);
set(gca,'clim',[0 100]);
shading flat
colorbar
tightmap;
mlabel('south');

% geoshow(S','facecolor','none');  % "none" gives a clear land mass.  Use 'g' for green, etc.


stime=datestr(leg.tm(na1:na2));
ssec=(stime(:,19:20));
smin=(stime(:,16:17));
shr=(stime(:,13:14));
sec=str2num(ssec);
min=str2num(smin);
hr=str2num(shr);

htime=hr+min/60+sec/3600;

t1=[shr(1,:),':',smin(1,:),':',ssec(1,:)];
t2=[shr(nalong,:),':',smin(nalong,:),':',ssec(nalong,:)];

tightmap;

%title(['HIRAD 5 GHz excess T_B Leg 3 of Karl 16 Sept; ',t1,'-',t2])
title(['HIRAD 5 GHz excess T_B ' cname ' '],'fontsize',12)
print('-djpeg',[cname1 cname,'Before.jpg'],'-r300');


figure
angle=(-width:1:width);

plot(angle,mean(tb5))
hold on

xx=mean(tb5);
xmean=mean(xx);

nfit=5;  % 5 = 4th order polynomial cross-track fit of along-track mean

xbar=ones(nfit,ncross);

n=2;
    xbar(n,:)=angle;
n=3;
    xbar(n,:)=angle.*angle;
n=4;
    xbar(n,:)=angle.*angle.*angle;
n=5;
    xbar(n,:)=angle.*angle.*angle.*angle;

%  Add more terms if needed

fitcoef=xx/xbar;

xfit=fitcoef*xbar;

plot(angle,xfit,'r');
plot(angle,xx-xfit,'g');
grid on;
title([cname ' mean excess T_B, 4th-order fit ']);
print('-djpeg',[cname1 cname 'Mean_4thOrderFit.jpg'],'-r300');

xdif=xx-xfit;
xdif=repmat(xdif,nalong,1);

figure('Position',[x1 y1 xwidth ywidth])

axesm('MapProjection','eqdcylin','maplatlimit',[latmin latmax],'maplonlimit',[lonmin lonmax],...
       'grid','on','mlinelocation',.5,'plinelocation',.5,'meridianlabel','on','mlabelround',-1,...
    'labelformat','none','parallellabel','on','plabelround',-1,'flinewidth',1,'fontsize',11,...
    'fontweight','normal');
mlabel('south');
delTB5filtrd=tb5-xdif;
pcolorm (lat, lon, tb5-xdif);
set(gca,'clim',[0 100]);
shading flat
colorbar
tightmap;
mlabel('south');
title(['HIRAD 5 GHz excess T_B ' cname ' after filtering  '],'fontsize',12)
% geoshow(S','facecolor','none');  % "none" gives a clear land mass.  Use 'g' for green, etc.
print('-djpeg',[cname1 cname 'After.jpg'],'-r300');
filenm=[cname1 cname 'delTB5fltrd'];
xlat=lat; xlon=lon; xtime=leg.tm(na1:na2);
save (filenm, 'delTB5filtrd','xlat','xlon','xtime');

%% Write netcdf file

ncid = netcdf.create([filenm,'.nc'],'NC_NOCLOBBER');

time_dimid = netcdf.defDim(ncid,'TimeElapsed',length(delTB5filtrd(:,1)));
incang_dimid = netcdf.defDim(ncid,'IncAng',length(delTB5filtrd(1,:)));
single_dimid = netcdf.defDim(ncid,'single',1);

delTB5filtrdID = netcdf.defVar(ncid,'delTB5filtrd','double',[time_dimid,incang_dimid]);
swath_latID = netcdf.defVar(ncid,'swath_lat','double',[time_dimid,incang_dimid]);
swath_lonID = netcdf.defVar(ncid,'swath_lon','double',[time_dimid,incang_dimid]);

tmID = netcdf.defVar(ncid,'tm','double',[time_dimid,single_dimid]);

netcdf.endDef(ncid);

netcdf.putVar(ncid,delTB5filtrdID,delTB5filtrd);
netcdf.putVar(ncid,swath_latID,xlat);
netcdf.putVar(ncid,swath_lonID,xlon);
netcdf.putVar(ncid,tmID,xtime);

netcdf.close(ncid);