% compare ship main mast T/RH with bow tower T/RH
% for evidence of interference from the ship's exhaust stacks

%% initialize
clear
close all

% start date for good PIR data
stjd = 85;

% matlab script path
path_prog = fullfile('c:\Users\Sergio\Documents\MATLAB\2016_CAPRICORN\check_ship_data');
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

fig_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\stack_plume_analysis';
mkdir(fig_path);

%% concat 1-min ship & psd met files
path_proc_data = fullfile('D:','Data','CAPRICORN_2016','Investigator','flux','Processed_Raw');
filecat(fullfile(path_proc_data,'CAPRICORN_2016_proc_met_1min_*.txt'),fullfile(path_proc_data,'CAPRICORN_2016_proc_met_1min_all.txt'),'Overwrite');
filecat(fullfile(path_proc_data,'CAPRICORN_2016_proc_scs_1min_*.txt'),fullfile(path_proc_data,'CAPRICORN_2016_proc_scs_1min_all.txt'),'Overwrite');
filecat(fullfile(path_proc_data,'CAPRICORN_2016_proc_rad_1min_*.txt'),fullfile(path_proc_data,'CAPRICORN_2016_proc_rad_1min_all.txt'),'Overwrite');

%% Load ship and psd met data & define variables - same timebase for all
psd = load(fullfile(path_proc_data,'CAPRICORN_2016_proc_met_1min_all.txt'));
rad = load(fullfile(path_proc_data,'CAPRICORN_2016_proc_rad_1min_all.txt'));
scs = load(fullfile(path_proc_data,'CAPRICORN_2016_proc_scs_1min_all.txt'));

jd = psd(:,1);
ta = psd(:,7);                % noaa air temp
rh  = psd(:,8);               % noaa rh
pir = (rad(:,2)+rad(:,3))/2;  % mean noaa pir

ta_p = scs(:,19);             % ship port/stbd Ta, RH, PIR
ta_s = scs(:,20);
rh_p = scs(:,21);
rh_s = scs(:,22);
pir_p = rad(:,6);
pir_s = rad(:,7);

% use mean of ships anemometers for relative wind direction
rwdir = (unwrap(scs(:,10)*pi/180) + unwrap(scs(:,14)*pi/180))/2;
rwdir = mod(rwdir+2*pi, 2*pi)*180/pi;
rwdir(rwdir>180) = rwdir(rwdir>180)-360;  %+/- 180 format

%% compute differences & plot
Dta_p = ta_p - ta;
Dta_s = ta_s - ta;
Drh_p = rh_p - rh;
Drh_s = rh_s - rh;
Dpir_p = pir_p - pir;
Dpir_s = pir_s - pir;

figure; plot(rwdir,Dta_p,'r.',rwdir,Dta_s,'b.');
hold on; plot([0,0],[0,5],'k-',[-180,180],[0,0],'k-');
grid; xlabel('Relative Wind Direction'); ylabel('Ship Air Temp - Bow Tower Air Temp');
axis([-180,180,-1,5]); legend('port','stbd','location','north');
title('Influence of ship exhaust on main mast air temperature');
saveas(gcf,fullfile(fig_path,['Temp.fig']));

figure; plot(rwdir,Drh_p,'r.',rwdir,Drh_s,'b.');
hold on; plot([0,0],[0,100],'k-',[-180,180],[0,0],'k-');
grid; xlabel('Relative Wind Direction'); ylabel('Ship RH - Bow Tower RH');
axis([-180,180,-20,40]); legend('port','stbd','location','north');
title('Influence of ship exhaust on main mast humidity');
saveas(gcf,fullfile(fig_path,['RH.fig']));

ii = jd > stjd;  % select PIR data after fix for bias problem
figure; plot(rwdir(ii),Dpir_p(ii),'r.',rwdir(ii),Dpir_s(ii),'b.');
hold on; plot([0,0],[-50,50],'k-',[-180,180],[0,0],'k-');
grid; xlabel('Relative Wind Direction'); ylabel('Ship PIR - NOAA PIR');
axis([-180,180,-50,50]); legend('port','stbd','location','north');
title('Influence of ship exhaust on PIR radiation');
saveas(gcf,fullfile(fig_path,['PIR.fig']));
