%
% read S-HIS FBF radiance files
%
% created by: R. Knuteson, UW-SSEC
% revised: 28 Oct 2004  K. Vinson
%
% Read in S-HIS FBF variables and create a mean spectrum (and plot)
% This program should work with any day. Add a path to this script
% and to the "mfiles" directory (e.g. ~/matlab/mfiles).
%
% Prompts for the following inputs:
%   Date String (for plotting)
%   Scene Mirror Angle Range (e.g. 179 183 for nadir; -5 5 for zenith)
%   Min & Max Roll (e.g. -5 5 to remove turns)
%   Time Range in decimal UTC hours (e.g. 21.1 21.5)
%
% set path to this script and the mfiles directory
%
addpath /home/shis/matlab
addpath /home/shis/matlab/mfiles
  FontSize12 = 12;
%
disp('### Read S-HIS FBF Radiance Files ###')
%
lwwn = fbfread('WavenumberLW_Resampled.real8.1287',1,-1);
lwrad = fbfread('SpectrumLW_CalRad_Resampled.real4.1287',1,-1);
numrec = size(lwrad,2);
%
mwwn = fbfread('WavenumberMW_Resampled.real8.1702',1,-1);
mwrad = fbfread('SpectrumMW_CalRad_Resampled.real4.1702',1,-1);
%
swwn = fbfread('WavenumberSW_Resampled.real8.2593',1,-1);
swrad = fbfread('SpectrumSW_CalRad_Resampled.real4.2593',1,-1);
%
sta = fbfread('zfliLW.sta',1,numrec);
instrumentTime = fbfread('instrumentTime.real4',1,numrec);
aircraftTime = fbfread('aircraftTime.real4',1,numrec);
timeHHMMSS = fbfread('timeHHMMSS.real4',1,numrec);
%
sceneMirrorAngle = fbfread('sceneMirrorAngle.real4',1,numrec);
%
roll = fbfread('Roll.real4',1,numrec);
altitude = fbfread('Altitude.real4',1,numrec);
alat = fbfread('Latitude.real4',1,numrec);
alon = fbfread('Longitude.real4',1,numrec);
%
% Prompt for Input
%
date_string = input('Enter Date String for Label: ','s');
angles = str2num(input('Enter the Scene Mirror Angle Range (181 = nadir):','s'));
rollminmax = str2num(input('Enter the Min & Max Roll (e.g. -5 5):','s'));
apodized = input('View apodized spectra? (y or n): ','s')
%
scene_ind = find(sta==1 & sceneMirrorAngle > angles(1) & sceneMirrorAngle < angles(2) & ...
                          roll > rollminmax(1) & roll < rollminmax(2)& ...
                          lwrad(min(find(lwwn>670)),:) ~= 0);
% Enter the Time Range
timerange = str2num(input('Enter the Time Range (decimal UTC):','s'));
%
scene_time = instrumentTime(scene_ind);
selected_scene_ind = find(scene_time > timerange(1) & scene_time < timerange(2));
if (length(selected_scene_ind) <= 0) disp('ERROR: No Match Found'); end
selected_scene_time = scene_time(selected_scene_ind);
%
scene_lat = alat(scene_ind);
scene_lon = alon(scene_ind);
selected_scene_lat = scene_lat(selected_scene_ind);
selected_scene_lon = scene_lon(selected_scene_ind);
scene_altitude = altitude(scene_ind);
selected_scene_altitude = scene_altitude(selected_scene_ind);
%
mean_selected_scene_lat = mean(selected_scene_lat)
mean_selected_scene_lon = mean(selected_scene_lon)
mean_selected_scene_altitude = mean(selected_scene_altitude)
%
% Plot Lat/Lon/Altitude and selected fovs
%
figure
subplot(211)
plot(alon,alat,'b.',selected_scene_lon,selected_scene_lat,'r.')
%axis([-102 -92 28 43])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
grid
xlabel('Longitude')
ylabel('Latitude')
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])
subplot(212)
plot(instrumentTime,altitude,'b.',selected_scene_time,selected_scene_altitude,'r.')
%axis([16 23 0 20000])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
grid
xlabel('Time (UTC)')
ylabel('Altitude (m)')
%
% compute mean
%
scene_lwrad = lwrad(:,scene_ind);
tmp = scene_lwrad(:,selected_scene_ind);
unap_meanlwrad = mean(tmp,2);
tmp = unap_meanlwrad;
tmp(find(unap_meanlwrad < 0)) = NaN;
unap_meanlwbt = rad2bt(lwwn,tmp);
apod_meanlwrad = apodizer_rok(unap_meanlwrad,2);
lw_rolloff = [610 1100]';
[lwwn_zeropad apod_meanlwrad_zeropad] = zeropad4(lwwn,apod_meanlwrad,15,lw_rolloff);
[lwwn_zeropad unap_meanlwrad_zeropad] = zeropad4(lwwn,unap_meanlwrad,15,lw_rolloff);
tmp = apod_meanlwrad;
tmp(find(apod_meanlwrad < 0)) = NaN;
apod_meanlwbt = rad2bt(lwwn,tmp);
tmp2 = apod_meanlwrad_zeropad;
tmp2(find(apod_meanlwrad_zeropad < 0)) = NaN;
apod_meanlwbt_zeropad = rad2bt(lwwn_zeropad,tmp2);
unap_meanlwbt_zeropad = rad2bt(lwwn_zeropad,unap_meanlwrad_zeropad);
%
scene_mwrad = mwrad(:,scene_ind);
unap_meanmwrad = mean(scene_mwrad(:,selected_scene_ind),2);
tmp = unap_meanmwrad;
tmp(find(unap_meanmwrad < 0)) = NaN;
unap_meanmwbt = rad2bt(mwwn,tmp);
apod_meanmwrad = apodizer_rok(unap_meanmwrad,2);
mw_rolloff = [1042 1812];
[mwwn_zeropad apod_meanmwrad_zeropad] = zeropad4(mwwn,apod_meanmwrad,15,mw_rolloff);
[mwwn_zeropad unap_meanmwrad_zeropad] = zeropad4(mwwn,unap_meanmwrad,15,mw_rolloff);
tmp = apod_meanmwrad;
tmp(find(apod_meanmwrad < 0)) = NaN;
apod_meanmwbt = rad2bt(mwwn,tmp);
tmp2 = apod_meanmwrad_zeropad;
tmp2(find(apod_meanmwrad_zeropad < 0)) = NaN;
apod_meanmwbt_zeropad = rad2bt(mwwn_zeropad,tmp2);
unap_meanmwbt_zeropad = rad2bt(mwwn_zeropad,unap_meanmwrad_zeropad);
%
scene_swrad = swrad(:,scene_ind);
unap_meanswrad = mean(scene_swrad(:,selected_scene_ind),2);
tmp = unap_meanswrad;
tmp(find(unap_meanswrad < 0)) = NaN;
unap_meanswbt = rad2bt(swwn,unap_meanswrad);
apod_meanswrad = apodizer_rok(unap_meanswrad,2);
sw_rolloff = [1812 3000];
[swwn_zeropad apod_meanswrad_zeropad] = zeropad4(swwn,apod_meanswrad,15,sw_rolloff);
[swwn_zeropad unap_meanswrad_zeropad] = zeropad4(swwn,unap_meanswrad,15,sw_rolloff);
tmp = apod_meanswrad;
tmp(find(apod_meanswrad < 0)) = NaN;
apod_meanswbt = rad2bt(swwn,apod_meanswrad);
apod_meanswbt_zeropad = rad2bt(swwn_zeropad,apod_meanswrad_zeropad);
unap_meanswbt_zeropad = rad2bt(swwn_zeropad,unap_meanswrad_zeropad);
%
%
% find cut-off indices
index.lwwn_cutoff = find(lwwn_zeropad > 610 & lwwn_zeropad < 1100);
index.mwwn_cutoff = find(mwwn_zeropad > 1042 & mwwn_zeropad < 1812);
index.swwn_cutoff = find(swwn_zeropad > 1812 & swwn_zeropad < 3000);
%
%
%
if apodized == 'n'
% Unapodized Radiance Plot with zeropadding
figure
plot(lwwn_zeropad(index.lwwn_cutoff), unap_meanlwrad_zeropad(index.lwwn_cutoff), 'b', mwwn_zeropad(index.mwwn_cutoff), unap_meanmwrad_zeropad(index.mwwn_cutoff), 'g', swwn_zeropad(index.swwn_cutoff), unap_meanswrad_zeropad(index.swwn_cutoff), 'r')
grid
axis([600 3000 0 130])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
text(1600,120,'Unapodized','FontSize',FontSize12,'Fontweight','bold');
xlabel('Wavenumber (cm^{-1})')
ylabel('S-HIS Obs. Radiance (mW/(m^2 str cm^{-1})')
legend('LW','MW','SW') 
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])

% Unapodized Brightness Temperature with zeropadding
figure
plot(lwwn_zeropad(index.lwwn_cutoff), unap_meanlwbt_zeropad(index.lwwn_cutoff), 'b', mwwn_zeropad(index.mwwn_cutoff), unap_meanmwbt_zeropad(index.mwwn_cutoff), 'g', swwn_zeropad(index.swwn_cutoff), unap_meanswbt_zeropad(index.swwn_cutoff), 'r')
grid
axis([500 3000 190 305])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
text(1500,300,'Unapodized','FontSize',FontSize12,'Fontweight','bold');
xlabel('Wavenumber (cm^{-1})')
ylabel('S-HIS Obs. B.T. (K)')
legend('LW','MW','SW') 
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])

else

% Apodized Radiance Plot
%figure
%plot(lwwn(1:end-5), apod_meanlwrad(1:end-5), 'b', mwwn(1:end-5), apod_meanmwrad(1:end-5), 'g', swwn(1:end-5), apod_meanswrad(1:end-5), 'r')
%grid
%xlabel('Wavenumber (cm^{-1})')
%ylabel('S-HIS Obs. Radiance (mW/(m^2 str cm^{-1})')
%legend('LW','MW','SW') 
%title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])

% Apodized Radiance Plot with zeropadding
figure
plot(lwwn_zeropad(index.lwwn_cutoff), apod_meanlwrad_zeropad(index.lwwn_cutoff), 'b', mwwn_zeropad(index.mwwn_cutoff), apod_meanmwrad_zeropad(index.mwwn_cutoff), 'g', swwn_zeropad(index.swwn_cutoff), apod_meanswrad_zeropad(index.swwn_cutoff), 'r')
grid
axis([600 3000 0 130])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
text(1600,120,'Apodized','FontSize',FontSize12,'Fontweight','bold');
xlabel('Wavenumber (cm^{-1})')
ylabel('S-HIS Obs. Radiance (mW/(m^2 str cm^{-1})')
legend('LW','MW','SW') 
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])


% Apodized Brightness Temperature
%figure
%plot(lwwn(1:end-5), apod_meanlwbt(1:end-5), 'b', mwwn(5:end-5), apod_meanmwbt(5:end-5), 'g', swwn(1:end-5), apod_meanswbt(1:end-5), 'r')
%grid
%axis([500 3000 210 310])
%xlabel('Wavenumber (cm^{-1})')
%ylabel('S-HIS Obs. B.T. (K)')
%legend('LW','MW','SW') 
%title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])
%
% Apodized Brightness Temperature with zeropadding
figure
plot(lwwn_zeropad(index.lwwn_cutoff), apod_meanlwbt_zeropad(index.lwwn_cutoff), 'b', mwwn_zeropad(index.mwwn_cutoff), apod_meanmwbt_zeropad(index.mwwn_cutoff), 'g', swwn_zeropad(index.swwn_cutoff), apod_meanswbt_zeropad(index.swwn_cutoff), 'r')
grid
axis([500 3000 210 310])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
text(1500,300,'Apodized','FontSize',FontSize12,'Fontweight','bold');
xlabel('Wavenumber (cm^{-1})')
ylabel('S-HIS Obs. B.T. (K)')
legend('LW','MW','SW') 
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])
%
end

% compute standard deviation
%
%scene_lwrad = lwrad(:,scene_ind);
unap_stdlwrad = std(scene_lwrad(:,selected_scene_ind),1,2);
tmp = unap_stdlwrad;
tmp(find(unap_stdlwrad < 0)) = NaN;
unap_stdlwbt = rad2bt(lwwn,tmp);
apod_stdlwrad = apodizer_rok(unap_stdlwrad,2);
tmp = apod_stdlwrad;
tmp(find(apod_stdlwrad < 0)) = NaN;
apod_stdlwbt = rad2bt(lwwn,apod_meanlwrad+tmp)-apod_meanlwbt;
%
%scene_mwrad = mwrad(:,scene_ind);
unap_stdmwrad = std(scene_mwrad(:,selected_scene_ind),1,2);
tmp = unap_stdmwrad;
tmp(find(unap_stdmwrad < 0)) = NaN;
unap_stdmwbt = rad2bt(mwwn,tmp);
apod_stdmwrad = apodizer_rok(unap_stdmwrad,2);
tmp = apod_stdmwrad;
tmp(find(apod_stdmwrad < 0)) = NaN;
apod_stdmwbt = rad2bt(mwwn,apod_meanmwrad+tmp)-apod_meanmwbt;
%
%scene_swrad = swrad(:,scene_ind);
unap_stdswrad = std(scene_swrad(:,selected_scene_ind),1,2);
tmp = unap_stdswrad;
tmp(find(unap_stdswrad < 0)) = NaN;
unap_stdswbt = rad2bt(swwn,unap_stdswrad);
apod_stdswrad = apodizer_rok(unap_stdswrad,2);
tmp = apod_stdswrad;
tmp(find(apod_stdswrad < 0)) = NaN;
apod_stdswbt = rad2bt(swwn,apod_meanswrad+tmp)-apod_meanswbt;

%
% Apodized Standard Deviation Plot
figure
subplot(311)
plot(lwwn(1:end-5), apod_stdlwrad(1:end-5), 'b')
grid
set(gca,'FontSize',FontSize12,'FontWeight','bold');
set(gca,'YLim',[0 4.2])
xlabel('wavenumber (cm^{-1})')
ylabel('LW Scene STD (RU)')
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Number of Samples: ' num2str(length( selected_scene_ind))])
subplot(312)
plot(mwwn(1:end-5), apod_stdmwrad(1:end-5), 'g')
grid
set(gca,'YLim',[0 3.5])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
xlabel('wavenumber (cm^{-1})')
ylabel('MW Scene STD (RU)')
subplot(313)
plot(swwn(1:end-5), apod_stdswrad(1:end-5), 'r')
grid
set(gca,'YLim',[0 0.4])
set(gca,'FontSize',FontSize12,'FontWeight','bold');
xlabel('wavenumber (cm^{-1})')
ylabel('SW Scene STD (RU)')
%
% Apodized Brightness Temperature
figure
plot(lwwn(1:end-5), apod_stdlwbt(1:end-5), 'b', mwwn(1:end-5), apod_stdmwbt(1:end-5), 'g', swwn(1:end-5), apod_stdswbt(1:end-5), 'r')
grid
set(gca,'FontSize',FontSize12,'FontWeight','bold');
axis([500 3000 0 10])
xlabel('wavenumber (cm^{-1})')
ylabel('S-HIS Scene Std. Dev. (1-sigma) (K)')
legend('LW','MW','SW') 
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Number of Samples: ' num2str(length( selected_scene_ind))])
%
% Diagnostic Plots
%
% Plot of all the scenes used in the average to QC for bad spectra
%
% Unapodized Radiance Plot
figure
subplot(311)
plot(lwwn, scene_lwrad(:,selected_scene_ind), 'b')
set(gca,'FontSize',FontSize12,'FontWeight','bold');
grid
xlabel('wavenumber (cm^{-1})')
ylabel('LW')
title(['Selected Spectra: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])
subplot(312)
plot(mwwn, scene_mwrad(:,selected_scene_ind), 'g')
grid
set(gca,'FontSize',FontSize12,'FontWeight','bold');
xlabel('wavenumber (cm^{-1})')
ylabel('MW')
subplot(313)
plot(swwn, scene_swrad(:,selected_scene_ind), 'r')
grid
set(gca,'FontSize',FontSize12,'FontWeight','bold');
xlabel('wavenumber (cm^{-1})')
ylabel('SW')
%
% Timeseries of IR channels
%
lw_670_690_ind = find(lwwn > 670 & lwwn < 690);
lw_670_690_meanrad_ts = mean(scene_lwrad(lw_670_690_ind,selected_scene_ind),1);
lw_670_690_meanbt_ts = rad2bt(680, lw_670_690_meanrad_ts);
%
lw_980_985_ind = find(lwwn > 980 & lwwn < 985);
lw_980_985_meanrad_ts = mean(scene_lwrad(lw_980_985_ind,selected_scene_ind),1);
lw_980_985_meanbt_ts = rad2bt(982.5, lw_980_985_meanrad_ts);
%
%
figure
plot(selected_scene_time,lw_670_690_meanbt_ts,'b.',...
     selected_scene_time,lw_980_985_meanbt_ts,'r.')
grid
set(gca,'FontSize',FontSize12,'FontWeight','bold');
set(gca,'YLim',[190 300],'XLim',[timerange(1) timerange(2)]);
xlabel('Time (UTC)')
ylabel('S-HIS Observed B.T. (K)')
legend('LW 670-690 cm^{-1}','LW 980-985 cm^{-1}',0)
title(['AVE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])

% save a matlab "mat" file containing the mean and std spectra
numSamples = length( selected_scene_ind);
save shis_mean_spectra  date_string timerange angles rollminmax numSamples selected_scene_time ...
       lwwn unap_meanlwrad mwwn unap_meanmwrad swwn unap_meanswrad ...
       lwwn unap_stdlwrad mwwn unap_stdmwrad swwn unap_stdswrad