%{
 process_pir.m  -  reads Investigator raw NMEA meteo files and
    reprocesses PIR data to W/m2.  PIR data is approx 1 Hz.
%}

% initializde
clear
close all

% 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;

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

% move to raw met directory
cd ('D:\Data\CAPRICORN_2016\Investigator\flux\SCS\met_nmea\');
files = dir('*.METEO');

ppir = []; spir = [];

for ii = 1:length(files)
    dfl = files(ii).name;
    disp(['processing ',dfl]);
    
    met = importdata(dfl);  % load as cell array
    
    % scan for port pir data
    raw = {}; praw = []; temp = [];
    temp = met(cellfun(@(x) ~isempty(strfind(x,'$PPPIR')),met));
    for j=1:length(temp)
        raw =textscan(temp{j,:},'%4f-%2f-%2f%*c%2f:%2f:%6.3f %*6c %f %f %f %f %*[^\n]','delimiter',',*');
        try
            praw(j,:)=[raw{1},raw{2},raw{3},raw{4},raw{5},raw{6},raw{7},raw{8},raw{9},raw{10}];
        catch
            praw(j,:)=NaN(1,10);
        end;
    end;
    praw(isnan(praw(:,1)),:) = [];
    jd_p = datenum(praw(:,1:6))-datenum(2015,12,31);
    praw = [jd_p,praw(:,7:10)];
    
    % repeat for stbd pir
    raw = {}; sraw = []; temp = [];
    temp = met(cellfun(@(x) ~isempty(strfind(x,'$PSPIR')),met));
    for j=1:length(temp)
        raw =textscan(temp{j,:},'%4f-%2f-%2f%*c%2f:%2f:%6.3f %*6c %f %f %f %f %*[^\n]','delimiter',',*');
        try
            sraw(j,:)=[raw{1},raw{2},raw{3},raw{4},raw{5},raw{6},raw{7},raw{8},raw{9},raw{10}];
        catch
            sraw(j,:)=NaN(1,10);
        end;
    end;
    sraw(isnan(sraw(:,1)),:) = [];
    jd_s = datenum(sraw(:,1:6))-datenum(2015,12,31);
    sraw = [jd_s,sraw(:,7:10)];
    
    ppir = [ppir;praw];
    spir = [spir;sraw];
end;

clear praw sraw temp raw met jd_p jd_s

% reference 5 sec timestamp
stjd = floor(ppir(1,1));
endjd = ceil(ppir(end,1));
jdref = stjd:5/86400:endjd; jdref = jdref(:); % 5-sec bin edges

% compute 5 sec mean of raw pir data and save to single array
% cols: jd, W/m2_port, uV_port, Tc_port, Td_port, W/m2_stbd, uV_stbd, Tc_stbd, Td_stbd, W/m2_port_recalc, W/m2_stbd_recalc
%        1,     2,        3,       4,       5,        6,        7,       8,       9,           10,               11
pir5 = NaN(length(jdref)-1,11);
pir5(:,1) = jdref(1:end-1);
avg = interval_avg(ppir(:,1),ppir(:,2:5),jdref);
pir5(:,2:5) = avg(:,2:5);
avg = interval_avg(spir(:,1),spir(:,2:5),jdref);
pir5(:,6:9) = avg(:,2:5);

clear ppir spir jdref avg

% constants
Sb = 5.67e-8;     % Stephan Boltzman Const.
CC = 3.5;         % PIR case constant
CCpsd = 4;
sens_port = 3.42; % uV per W/m2 for port PIR
sens_stbd = 2.85; % uV per W/m2 for stbd PIR

% W/m2 = uV/sens + Sb*Tc^4 + CC*Sb*(Tc^4 - Tc^4)
pir5(:,10) = pir5(:,3)/sens_port + Sb*pir5(:,4).^4 - CC*Sb*(pir5(:,5).^4 - pir5(:,4).^4);
pir5(:,11) = pir5(:,7)/sens_stbd + Sb*pir5(:,8).^4 - CC*Sb*(pir5(:,9).^4 - pir5(:,8).^4);

pir5(:,12) = pir5(:,3)/sens_port + Sb*pir5(:,4).^4 - CCpsd*Sb*(pir5(:,5).^4 - pir5(:,4).^4);
pir5(:,13) = pir5(:,7)/sens_stbd + Sb*pir5(:,8).^4 - CCpsd*Sb*(pir5(:,9).^4 - pir5(:,8).^4);

jdship = pir5(:,1);

zz = find(isfinite(pir5(:,2)));
figure; plot(jdship,pir5(:,2)-pir5(:,12),'b.');grid;
xlim([jdship(zz(1)),jdship(zz(end))]); xlabel('day-of-year'); ylabel('PIR(CC=3.5) - PIR(cc=4), W/m2');
title('Port PIR Difference due to Case Constant, 3.5 vs 4');

% load psd data for jd 75-76
%{
1        decimal day of year
2        pir1 W/m2
3        psp1 W/m2
4        PIR1 case temperature C
5        PIR1 dome temperature C
6        SST snake C
7        air temperature C
8        RH
9        rain rate, org, mm/hr
10       org carrier voltage
11       T/RH aspirator fan
12       atmospheric pressure mb
%}
psd = load('D:\Data\CAPRICORN_2016\Investigator\flux\Processed_In_Field\CAPRICORN_2016_proc_met_1min_2016_03_24_084.txt');
jdpsd = psd(:,1);

% load ship 1-min data for jd 75-76
%{
1       jd                  decimal day of year
2       lat                 seapath latitude
3       lon                 seapath longitude, deg E
4       sog                 seapath sog, m/s
5       cog                 seapath cog, deg
6       hed                 gyro heading, deg
7       wspd_port           true wspd, port, m/s @ 25m
8       wdir_port           true wdir, port, deg
9       rwspd_port          rel wspd, port, m/s @ 25m
10      rwdir_port          rel wdir, port, 0-360 deg
11      wspd_stbd           true wspd, stbd, m/s @ 25m
12      wdir_stbd           true wdir, stbd, deg
13      rwspd_stbd          rel wspd, stbd, m/s @ 25m
14      rwdir_stbd          rel wdir, stbd, 0-360 deg
15      wspd_sonic          true wspd, sonic, m/s @ 22.5m
16      wdir_sonic          true wdir, sonic, deg
17      rwspd_sonic         rel wspd, sonic, m/s @ 22.5m
18      rwdir_sonic         rel wdir, sonic, 0-360 deg
19      ta_port             air temp, port, C @ 31.4m
20      ta_stbd             air temp, stbd, C 31.4m
21      rh_port             rel humidity, port, % 31.4m
22      rh_stbd             rel humidity, stbd, % 31.4m
23      Pmb                 atm preessure, mb 20m
24      rain                accumulated rain, RM Young, mm
25      rain_rate           rain rate, RM Young, mm/hr
26      rain_optical        accumulated rain, optical disdrometer, mm
27      rain_rate_optical   rain rate, optical disdrometer, mm/hr
28      PSP_port            solar rad flux, port, W/m2
29      PSP_stbd            solar rad flux, stbd, W/m2
30      PIR_port            IR rad flux, port, W/m2
31      PIR_stbd            IR rad flux, stbd, W/m2
32      SST_IR              ISAR infrared sea surface temp, C
33      SST                 TSG water temp, C @ depth of port drop keel
34      Sal                 TSG salinity, psu
35      uship_long          longitudinal speed through water, m/s
36      uship_trans         transverse speed through water, m/s
37      SST_depth           depth of SST on port drop keel, m
38      sog_n               N-S component of ship speed
39      sog_w               E-W component of ship speed
40      hed_n               N-S component of ship heading
41      hed_w               E-W component of ship heading
%}
scs = load('D:\Data\CAPRICORN_2016\Investigator\flux\Processed_In_Field\CAPRICORN_2016_proc_scs_1min_2016_03_24_084.txt');
jdscs = scs(:,1);

% plot temperatures
figure; plot(jdpsd,psd(:,7),'r.',jdscs,scs(:,19),'g.');
hold on; plot(jdpsd,psd(:,4),'r-',jdpsd,psd(:,5),'b-');
plot(jdship,pir5(:,4)-273.15,'g-',jdship,pir5(:,5)-273.15,'k-',jdship,pir5(:,8)-273.15,'m-',jdship,pir5(:,9)-273.15,'m--');
grid; xlim([jdpsd(1),jdpsd(end)]); xlabel('day-of-year'); ylabel(' Temperature, C');
ylim([10,25]); title('Air and PIR Case/Dome Temperatures, C');
legend('NOAA Tair','Ship Tair Port','NOAA PIR1 Tc','NOAA PIR1 Td','Ship Port Tc','Ship Port Td','Ship Stbd Tc','Ship Stbd Td','location','north');

% interp noaa Tc/Td to ship pir timebase
% Tc2 = interp1(jdpsd,psd(:,4),pir5(:,1))+273.15;
% Td2 = interp1(jdpsd,psd(:,5),pir5(:,1))+273.15;
% figure; plot(jdpsd,psd(:,7),'r.',jdscs,scs(:,19),'g.');
% hold on; plot(jdpsd,psd(:,4),'m.',jdpsd,psd(:,5),'b.',jdship,Tc2-273.15,'g-',jdship,Td2-273.15,'y-'); grid;
% xlim([jdpsd(1),jdpsd(end)]); xlabel('day-of-year'); ylabel(' Temperature, C');
% title('Air and PIR Case/Dome Temperatures, C');
% legend('NOAA Tair','Ship Tair Port','NOAA PIR1 Tc','NOAA PIR1 Td','Ship Port Tc2','Ship Port Td2','location','north');

% recompute ship pir with adjusted Tc/Td
pir5(:,14) = pir5(:,3)/sens_port + Sb*Tc2.^4 - CC*Sb*(Td2.^4 - Tc2.^4);
pir5(:,15) = pir5(:,7)/sens_stbd + Sb*Tc2.^4 - CC*Sb*(Td2.^4 - Tc2.^4);

% compare noaa pir1 and stbd pir
ii = find(jdship>84.4 & jdship<84.7); % nighttime
jj = find(jdpsd>84.4 & jdpsd<84.7); % nighttime
dIR = mean(pir5(ii,6)) - mean(psd(jj,2));
figure; plot(jdpsd,psd(:,2),'r-',jdship,pir5(:,6),'b-'); grid;
xlim([jdpsd(1),jdpsd(end)]); ylim([290,400]);
xlabel('day-of-year'); ylabel('Downwelling IR, W/m2');
title(['IR Flux with adj Tc/Td, mean night diff = ',sprintf('%2.2f',dIR),' W/m2']);
legend('NOAA PIR1', 'Ship PIR Stbd','location','north');




