clear all;
close all;

%weather types
LineWidth = 1.5;
WeatherTypes = {'Clear', ...                                        %00
                'Haze/smoke/dust, vis=>1km', ...                    %04
                'Haze/smoke/dust, vis<1km', ...                     %05
                'Mist', ...                                         %10
                'Fog*', ...                                         %20
                'PRECIPITATION*', ...                               %21
                'Drizzle (not freezing) or snow grains*', ...       %22
                'Rain (not freezing)*', ...                         %23
                'Snow*', ...                                        %24
                'Freezing rain or freezing drizzle*', ...           %25
                'Fog', ...                                          %30
                'Fog or ice fog, in patches', ...                   %31
                'Fog or ice fog, thinning', ...                     %32
                'Fog or ice fog, no change', ...                    %33
                'Fog or ice fog, thickening', ...                   %34
                'PRECIPITATION', ...                                %40
                'Precipitation, slight or moderate', ...            %41
                'Precipitation, heavy', ...                         %42
                'DRIZZLE', ...                                      %50
                'Drizzle, not freezing, slight', ...                %51
                'Drizzle, not freezing, moderate', ...              %52
                'Drizzle, not freezing, heavy', ...                 %53
                'Drizzle, freezing, slight', ...                    %54
                'Drizzle, freezing, moderate', ...                  %55
                'Drizzle, freezing, heavy', ...                     %56
                'RAIN', ...                                         %60
                'Rain, light', ...                                  %61
                'Rain, moderate', ...                               %62
                'Rain, heavy', ...                                  %63
                'Rain, freezing, light', ...                        %64
                'Rain, freezing, moderate', ...                     %65
                'Rain, freezing, heavy', ...                        %66
                'Rain (or drizzle) and snow, light', ...            %67
                'Rain (or drizzle) and snow, moderate or heavy', ...%68
                'SNOW', ...                                         %70
                'Snow, light', ...                                  &71
                'Snow, moderate', ...                               %72
                'Snow, heavy', ...                                  %73
                'Ice pellets, light', ...                           %74
                'Ice pellets, moderate', ...                        %75
                'Ice pellets, heavy', ...                           %76
                'SHOWERS OR INTERMITTENT PRECIPITATION', ...        %80
                'Rain showers, light', ...                          %81
                'Rain showers, moderate', ...                       %82
                'Rain showers, heavy', ...                          %83
                'Rain showers, violent ( >32 mm/h )', ...           %84
                'Snow showers, light', ...                          %85
                'Snow showers, moderate', ...                       %86
                'Snow showers, heavy', ...                          %87
                };
SYNOPCodes = [00, 04, 05, 10, ...
              20, 21, 22, 23, 24, 25, ...
              30, 31, 32, 33, 34, ...
              40, 41, 42, ...
              50, 51, 52, 53, 54, 55, 56, ...
              60, 61, 62, 63, 64, 65, 66, 67, 68, ...
              70, 71, 72, 73, 74, 75, 76, ...
              80, 81, 82, 83, 84, 85, 86, 87];


%parameters
% FilePath = 'G:\Shared drives\ASTraL\pilot2023\data\pwd\';%'G:\Shared drives\ASTraL\pilot2023\codes\readpwd\';
FilePath = '/Users/eliz/DATA/ASTRAL_2023/Revelle/pwd/';
ds = datastore(FilePath,'IncludeSubfolders',1,'Type','tabulartext');
ds.Files = ds.Files(contains(ds.Files,'PWD'));
ppath = FilePath;

%load file
figure
for ifile = 1:length(ds.Files)
    fileID = fopen(ds.Files{ifile});

    %read file
    A = textscan(fileID,'%s%s%s%f%f%s%f%f%f%f%f%s');%, '%s%f%f%f%f%s%f%f%f%f%f%f', 'HeaderLines', 0, 'Delimiter', ' ');

    %close file
    fclose(fileID);

    %get variables from matrix
    t = A{1};
    VisStatus = A{2};
    HWStatus = A{3};
    Vis1Min = A{4};
    Vis10Min = A{5};
    CodeNWS = A{6};
    CodePW = A{7};
    CodePW15Min = A{8};
    CodePW1hr = A{9};
    WaterIntensity = A{10};
    WaterSum = A{11};
    SnowSum = A{12};
    t0 = datenum(t, 'yyyymmddTHHMMSS');
    t01 = datetime(t0,'ConvertFrom','datenum');

    if ifile == 1
        t0big = t0;
        WaterIntensitybig = WaterIntensity;
        WaterSumbig = WaterSum;
        CodePWbig = CodePW;
    else
        t0big = [t0big;t0];
        WaterIntensitybig = [WaterIntensitybig;WaterIntensity];
        WaterSumbig = [WaterSumbig;WaterSum];
        CodePWbig = [CodePWbig;CodePW];
    end

    clear A;

end

%% PLOT the entire time series, or move up into loop to plot each day

t0 = t0big;
CodePW = CodePWbig;
WaterIntensity = WaterIntensitybig;
WateSum = WaterSumbig;

%%%%% visibility
% range = 2750:4000;
figure(1)
subplot(3,1,1);
semilogy(t0, Vis1Min, 'LineWidth', LineWidth)
xlim([min(t0) max(t0)]);
hold on
% semilogy(t0, Vis10Min)
% plot([min(t0),max(t0)), [1000,1000], 'LineWidth', LineWidth)
% xlim([t0(2500) t0(3500)])
hold off
box on
datetick('x','DD','keeplimits')
grid on;
% xlabel('t')
ylabel('visibility (m)')
% yticks([1e3, 1e4])
% yticklabels({'10^3','10^4'})

%%%%% rain intensity
% figure(2)
subplot(3,1,2)
hold on
plot(t0, WaterIntensity, 'LineWidth', LineWidth)
xlim([min(t0) max(t0)]);
% xlim([t0(2500) t0(3500)])
datetick('x','DD','keeplimits')
grid on;
% xlabel('t')
ylabel('Rain intensity')

%%%%% codes
% figure(3)
subplot(3,1,3)
FaceColor = {'r', 'g', 'b', 'y', 'm', 'c', 'k'};
Codes = unique(CodePW);
Codes(Codes==0) = [];
LegendString = cell(size(Codes));
hold on
for i=1:length(Codes)
    area(t0(CodePW==Codes(i)), zeros(sum(CodePW==Codes(i)),1), 'FaceColor', FaceColor{i})
    LegendString{i} = WeatherTypes{SYNOPCodes==Codes(i)};
end
parts = [1; find(diff(CodePW)~=0) + 1];
for i=1:length(parts)-1
    if CodePW(parts(i))~=0
        Color = FaceColor{CodePW(parts(i))==Codes};
        area(t0(parts(i):parts(i+1)), ones(parts(i+1)-parts(i)+1,1), 'FaceColor', Color, 'EdgeColor', 'none')
    end
end
box on
xlim([min(t0) max(t0)]);
ylabel('Weather code')
if ~isempty(Codes)
    legend(LegendString)%, 'Location', 'best')
end
datetick('x','DD','keeplimits')
grid on;
hold off
% set(findall(gcf,'-property','FontSize'),'FontSize',15)

print('-dpng',[ppath 'rain.png']);

%% save data

pwd.t = t0;
pwd.prate = WaterIntensity;
pwd.paccum = WaterSum;
pwd.code = CodePW;
pwd.codelegend = Codes;
pwd.weathercode = {'mist';'rain, light';'rain, moderate'};
save([ppath 'ASTRAL_rain.mat'],'pwd');