% analyze_wind_sector_met_city.m
%{
    Reads wind files, computes 1-min flag variable for 'good' in-sector
    time periods at the met city tower.

    Run met_city_pic_10min_avg.m before running this script.
%}

%% initialize
clear
close all
warning('off','MATLAB:MKDIR:DirectoryExists');

%% adjust these parameters as desired, then run

% enter date range, year, leg#
% leg 1: jd 290:347, 2019
% leg 2: jd 348:365, 2019 & jd 001:054, 2020
% leg 3: jd 055:110, 2020
% leg 4: jd 176:211, 2020
% leg 5: jd 240:253, 2020

jd_start = 240;
jd_end = 259;
yyyy = 2020; YYYY = sprintf('%4i',yyyy); YY = YYYY(3:4);
leg = '5';

% root path to data files on this system
root_data_path = '/Users/fluxuser/Documents/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
son2_path = fullfile(root_data_path,'sonic_anemometer_3D_ice_station_tower_10m');
plot_path = fullfile(root_data_path,'_scientific_analysis','flux_plots','wind_sector_plots_met_city');
mkdir(plot_path)
save_path = fullfile(root_data_path,'_scientific_analysis','raw_1min');
mkdir(save_path)
flux_path = fullfile(root_data_path,'_scientific_analysis','flux_10min');
mkdir(flux_path)
path_tower_met = fullfile(root_data_path,'cr1000_ice_station_tower');

plotit = 1; % set to 1 to save plot
N = jd_end - jd_start+ 1; % number of days to process

%% set matlab script path
restoredefaultpath
cd(path_prog);
addpath(genpath(path_prog));
rehash toolboxcache;

%% CO2/CH4 and ship met data

% 10-min co2/ch4 - from ship_picarro_10min_avg.m
dfl = fullfile(flux_path,[YYYY,'_Leg',leg,'_met_city_pic_10min.mat']);
load(dfl,'met_city_pic_10min');
%{
    Output array format: N x 10 arrays
        1   jd     decimal day-of-year timestamp
        2   sec    Picarro timestamp, seconds
        3   Pc     cavity pressure, torr
        4   Tc     cavity temperature, °C
        5   V1     solenoid valve flag
        6   ch4w   methane conc, ppm, wet
        7   ch4d   methane conc, ppm, dry
        8   co2w   carbon dioxide conc, ppm, wet
        9   co2d   carbon dioxide conc, ppm, dry
        10  h20    water vapor conc, %
%}

% 10-min met city logger data - from 2YYYY_LegX_cr1000_10min.mat
dfl = fullfile(flux_path,[YYYY,'_Leg',leg,'_cr1000_10min.mat']);
load(dfl,'cr1000_10min');
%{
    CR1000 output columns:
    1    jd                     day of year
    2    yyyy                   year
    3    mm                     month
    4    dd                     day
    5    hr                     hour
    6    mn                     minute
    7    sec                    second
    8    RECORD                 CR record number
    9    gps_lat_deg            lat deg
    10   gps_lat_min            lat min
    11   gps_lon_deg            lon deg
    12   gps_lon_min            lon min
    13   gps_hdg                heading
    14   gps_alt                altitude, m
    15   gps_qc                 gps QC value
    16   gps_nsat               gps number of sats
    17   gps_hdop               gps hdop
    18   PTemp                  CR1000 temperature, C
    19   batt_volt              battery voltage
    20   call_time_mainscan     iridium call time
    21   apogee_body_T          IR temperature, body, C
    22   apogee_targ_T          IR temperature, target, C
    23   sr50_dist              distance to surface, m
    24   vaisala_RH_2m          2m RH
    25   vaisala_T_2m           2m temperature, C
    26   vaisala_Td_2m          2m dew point, C
    27   vaisala_P_2m           2m pressure, mb
    28   vaisala_RH_6m          6m RH
    29   vaisala_T_6m           6m temperature, C
    30   vaisala_Td_6m          6m dew point, C
    31   vaisala_RH_10m         10m RH
    32   vaisala_T_10m          10m temperature, C
    33   vaisala_Td_10m         10m dew point, C
    34   fp_A_mV                flux plate A, mV
    35   fp_A_Wm2               flux plate A, W/m2
    36   fp_B_mV                flux plate B, mV
    37   fp_B_Wm2               flux plate B, W/m2
    38   mast_T                 Leeds mast, 20m temperature, C
    39   mast_RH                Leeds mast, 20m RH
    40   mast_P                 Leeds mast, 20m pressure, mb
%}

% 10-min ship met data - from YYYY_LegX_metnav_10min.m
dfl = fullfile(flux_path,[YYYY,'_Leg',leg,'_metnav_10min.mat']);
load(dfl,'metnav_10min');
%{
	metnav is 1440 x 32 array of 1-min met/nav data

    1   jd_min       1-min decimal day of year
    2   wspd0        ship 20m metek true wspd, m/s
    3   wdir0        ship 20m metek true wdir, deg
    4   tson0        ship 20m metek Tsonic, °C
    5   wspd2        met city 10m metek true wspd, m/s
    6   wdir2        met city 10m metek true wdir, deg
    7   tson2        met city 10m metek Tsonic, °C
    8   rwspd1       ship rm young rwspd, m/s
    9   rwdir1       ship rm young rwdir, deg
    10  wspd1        ship rm young true wspd, m/s
    11  wdir1        ship rm young true wdir, deg
    12  yr           year
    13  lat          latitude, decimal deg
    14  lon          longitude, decimal deg
    15  sog          gps speed over ground, m/s
    16  cog          gps course over ground, deg
    17  hed          heading, deg
    18  pitch        pitch angle, deg
    19  roll         roll angle, deg
    20  tair         air temperature, °C
    21  rh           relative humidity, %
    22  dpt          dewpoint, °C
    23  p_mb         air pressure, mb
    24  dir_rad      direct radiation, W/m2
    25  glb_rad      global radiation, W/m2
    26  rwspd        relative wind speed, m/s
    27  rwdir        relative wind direction, deg
    28  wspd         true wind speed, m/s
    29  wdir         true wind direction, deg
    30  vis          visibility, m
    31  tsea         sea surface temperature, °C
    32  ceil         ceiling, m
%}

figure; % vs ship bow tower sonic
subplot(2,1,1); plot(metnav_10min(:,1),metnav_10min(:,2),'.',...
        metnav_10min(:,1),metnav_10min(:,5),'-','markersize',8); grid minor;
    ylabel('wspd m/s'); xlabel(['DOY ',YYYY]); ylim([0,12]);
    legend('ship','m city');
    title('Comparison of Met City 10m & Bow 20m Wind');
subplot(2,1,2); plot(metnav_10min(:,1),metnav_10min(:,3),'.',...
        metnav_10min(:,1),metnav_10min(:,6),'-','markersize',8); grid minor;
    ylabel('wdir deg'); xlabel(['DOY ',YYYY]);
    ylim([0,360]); yticks([0,90,180,270,360]);
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_met_city_wind_v_ship.fig']);
    savefig(ppath);
end

    
% check CO2 vs wind direction
% jj is index array for metnav points that match picarro points
[~,kk,jj] = intersect(met_city_pic_10min(:,1),metnav_10min(:,1));

% matching sets of wind and picarro data
figure; plot(metnav_10min(jj,6),met_city_pic_10min(kk,9),'.','markersize',10); grid;
xlabel('wdir, deg'); xlim([0,360]); xticks(0:30:360);
ylabel('CO_2 ppm'); ylim([390,440]);
title('Met City CO_2 vs. Wind Direction');
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_met_city_co2_v_wdir.fig']);
    savefig(ppath);
end

% ship and tower heading
figure; plot(metnav_10min(:,1),metnav_10min(:,17),'.',...
        cr1000_10min(:,1),cr1000_10min(:,13),'.','markersize',5);
ylim([0,360]); yticks(0:45:360); grid minor;
ylabel('heading, deg'); xlabel(['DOY, ',YYYY]); legend('ship','tower');
title(['MOSAiC, Leg ',leg,', Ship and Tower Heading']);
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_mctower_ship_heading.fig']);
    savefig(ppath);
end

%% compute bearing of ship from tower to assess 'dirty' sector
[~,ia,ib] = intersect(cr1000_10min(:,1),metnav_10min(:,1));
jd_bearing = cr1000_10min(ia,1);
tower_lat = cr1000_10min(ia,9)+cr1000_10min(ia,10)/60;
tower_lon = cr1000_10min(ia,11)+cr1000_10min(ia,12)/60;
ship_lat = metnav_10min(ib,13);
ship_lon = metnav_10min(ib,14);
bearing = NaN(length(ia),1); % compass bearing to ship, deg
dist = NaN(length(ia),1);    % distance to the ship, naut miles
for ii = 1:length(ia)
    lat = [tower_lat(ii);ship_lat(ii)];
    lon = [tower_lon(ii);ship_lon(ii)];
    [bearing(ii),dist(ii)] = legs(lat,lon);
end
dist = dist*1852;  % nautical miles to meters

figure;
subplot(2,1,1);
    plot(jd_bearing,bearing,'.','markersize',5);
    ylim([0,360]); yticks(0:45:360); grid minor;
    ylabel('bearing to ship, deg'); xlabel(['DOY, ',YYYY]);
    title(['MOSAiC, Leg ',leg,', Ship Bearing and Distance from Met City']);
subplot(2,1,2);
    plot(jd_bearing,dist,'.','markersize',5); grid minor;
    ylabel('distance to ship, m'); xlabel(['DOY, ',YYYY]);
    ylim([0,800]);
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_MCtower_to_ship.fig']);
    savefig(ppath);
end

% copy bearing to length of metnav_10min time base
theta = NaN(length(metnav_10min),1);
theta(ib) = bearing;

%% 10-min sector flag for flux calculations
uL = theta + 20; lL = theta - 20;  % range of bad wind sector
wdir = metnav_10min(:,6);
sFlag = zeros(length(metnav_10min),1);
for ii = 1:length(sFlag)
    if lL(ii)>0 && uL(ii)<=360
        sFlag(ii) = wdir(ii)<lL(ii) | wdir(ii)>uL(ii);
    elseif uL(ii)>360   % catch case where upper limit exceeds 360
        sFlag(ii) = wdir(ii)<lL(ii) & wdir(ii)>(uL(ii)-360);
    elseif lL(ii)<0     % catch case where lower limit is negative
        sFlag(ii) = wdir(ii)>uL(ii) & wdir(ii)<(lL(ii)+360);
    else
        % nothing - sFlag is zero
    end
end

% flux flag covering entire wind data set
fluxFlag_met_city_10min = [metnav_10min(:,1),sFlag,metnav_10min(:,5:6)];
sfl = fullfile(flux_path,[YYYY,'_Leg',leg,'_fluxFlag_met_city_10min.mat']);
save(sfl,'fluxFlag_met_city_10min'); % save 10-min flux flag

%% plot wind

xmin = nanmin1(metnav_10min(:,1)); xmax = nanmax1(metnav_10min(:,1));
f1 = figure('position',[1000,600,700,900]);
left_color = [0 0 0]; right_color = [0 0 0];
set(f1,'defaultAxesColorOrder',[left_color; right_color]);
subplot(3,1,1)
    plot(metnav_10min(:,1)+5/1440,metnav_10min(:,6),'ro','markersize',2); hold on;
    plot(metnav_10min(sFlag==1,1)+5/1440,sFlag(sFlag==1),'go','markersize',2);
    xlim([xmin,xmax]); ylim([0,360]); yticks([0,90,180,270,360]);grid;
    ylabel('Wind Dir'); xlabel('day-of-year');
    title(['Met City Wind Sector Analysis, 10-min Average, Leg',leg,'  (In-Sector = Green)'],...
           'fontweight','bold','fontsize',14);
subplot(3,1,2)
    h4 = plot(metnav_10min(:,1)+5/1440,metnav_10min(:,5),'r.',metnav_10min(:,1),...
              metnav_10min(:,2),'b.','markersize',6); grid; hold on;
    h5 = plot(metnav_10min(sFlag==1,1)+5/1440,sFlag(sFlag==1)*0,'go','markersize',2);
    xlim([xmin,xmax]); ylim([0,20]);
    ylabel('Wind Spd, m/s'); xlabel('day-of-year');
    legend(h4,{'met city','bow-3D'},'location','best');
subplot(3,1,3)
    plot(met_city_pic_10min(:,1),met_city_pic_10min(:,9),'b.','markersize',6); hold on;
    plot(metnav_10min(sFlag==1,1)+5/1440,sFlag(sFlag==1)*410,'go','markersize',2);
    ylabel('CO_2, ppm'); xlim([xmin,xmax]); ylim([395,415]); grid;
    xlabel('day-of-year');
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_met_city_wind_sector.fig']);
    savefig(f1, ppath);
end

