%
% read S-HIS FBF radiance files
%
% revised: 15 Oct 2004  R. Knuteson, UW-SSEC
%
% 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
%
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('aircraftRollAngle.real4',1,numrec);
altitude = fbfread('aircraftAltitude.real4',1,numrec);
alat = fbfread('aircraftLatitude.real4',1,numrec);
alon = fbfread('aircraftLongitude.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'));
%
scene_ind = find(sta==1 & sceneMirrorAngle > angles(1) & sceneMirrorAngle < angles(2) & ...
                          roll > rollminmax(1) & roll < rollminmax(2));
% 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);
%
% Plot Lat/Lon/Altitude and selected fovs
%
figure
subplot(211)
plot(alon,alat,'b.',selected_scene_lon,selected_scene_lat,'r.')
axis([-160 -145 60 75])
grid
xlabel('longitude')
ylabel('latitude')
title(['MPACE 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([19 24 0 20000])
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);
tmp = apod_meanlwrad;
tmp(find(apod_meanlwrad < 0)) = NaN;
apod_meanlwbt = rad2bt(lwwn,tmp);
%
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);
tmp = apod_meanmwrad;
tmp(find(apod_meanmwrad < 0)) = NaN;
apod_meanmwbt = rad2bt(mwwn,tmp);
%
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);
tmp = apod_meanswrad;
tmp(find(apod_meanswrad < 0)) = NaN;
apod_meanswbt = rad2bt(swwn,apod_meanswrad);
%
% 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(['MPACE 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 280])
xlabel('wavenumber (cm^{-1})')
ylabel('S-HIS Obs. B.T. (K)')
legend('LW','MW','SW') 
title(['MPACE Campaign: ' date_string '  ' num2str(timerange(1)) '-' num2str(timerange(2)) 'UTC  ' ' Scan Angles ' num2str(angles(1)) '-' num2str(angles(2))])
%
% 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,'YLim',[0 2])
xlabel('wavenumber (cm^{-1})')
ylabel('LW Scene STD (RU)')
title(['MPACE 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 1])
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.25])
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
axis([500 3000 0 10])
xlabel('wavenumber (cm^{-1})')
ylabel('S-HIS Scene Std. Dev. (1-sigma) (K)')
legend('LW','MW','SW') 
title(['MPACE 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')
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
xlabel('wavenumber (cm^{-1})')
ylabel('MW')
subplot(313)
plot(swwn, scene_swrad(:,selected_scene_ind), 'r')
grid
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,'YLim',[200 300])
xlabel('Time (UTC)')
ylabel('S-HIS Observed B.T. (K)')
legend('LW 670-690 cm^{-1}','LW 980-985 cm^{-1}',0)
title(['MPACE 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