%look at seatac temeprature records
clear all; close all;
load('C:\Users\thedh\Desktop\research\temperature\stationa_data\seatac.mat', 'date', 'precip','snow_depth', 'snow', 'sun', 'tsun', 'T_AVG','T_MAX', 'T_MIN');

year = rounddown(date/10000);
month=rounddown((date-year*10000)/100);
day = date-year*10000- month*100;
tday = monthday(month)'+day; %day past the start of the year


precip(find(precip == -9999)) = NaN;
snow(find(snow == -9999)) = NaN;
snow_depth(find(snow_depth == -9999)) = NaN;
sun(find(sun == -9999)) = NaN;
tsun(find(tsun == -9999)) = NaN;
T_AVG(find(T_AVG == -9999)) = NaN;
T_MAX(find(T_MAX == -9999)) = NaN;
T_MIN(find(T_MIN == -9999)) = NaN;


%make climatology
for dy=1:365;
    do=find(tday == dy);
    CT_AVG(dy) = nanmean(T_AVG(do,:));
    CT_AVG_length(dy) = length(find(isnan( T_AVG(do,:)) == 0));
    STD_TAVG(dy) = nanstd(T_AVG(do,:));
    
    
    CT_MIN(dy) = nanmean(T_MIN(do,:));
    STD_TMIN(dy) = nanstd(T_MIN(do,:));
      
    CT_MAX(dy) = nanmean(T_MAX(do,:));
    STD_TMAX(dy) = nanstd(T_MAX(do,:));

    
end

LW=2; LW17=1.4; LWA=3;
MS=1.2;
alph=0.2;
FSz= 24; FN = 'Arial';
MS=6; LWM=1.4;


TDAY=1:365;
cut=50; %cutoff period (in days) for the low pass filter
ext=30; % extend the domain by this many days (make periodic) -- for visual purposes only
LOW_ct_ext = low_filter([CT_AVG CT_AVG CT_AVG], 1, cut); %three seasonal cycles to make periodic
LOW_STD_ext = low_filter([STD_TAVG STD_TAVG STD_TAVG], 1, cut); %as above biut for the standard deviation

LOW_CT  = LOW_ct_ext(365+1:365*2); %take middle year
LOW_STD  = LOW_STD_ext(365+1:365*2); 

figure(1); clf;
hold on; plot([tday-365 tday tday+365], [T_AVG T_AVG T_AVG], 'k.','MarkerSize', MS, 'LineWidth', LWM);

hold on; plot([TDAY-365 TDAY TDAY+365], [CT_AVG CT_AVG CT_AVG], 'k--', 'LineWidth', LW); %unsmoothed average at each day
hold on; plot([TDAY-365 TDAY TDAY+365], [LOW_CT LOW_CT LOW_CT], 'k-', 'LineWidth', LW); 

%add in the window of 3 standard deviations around the mean -- compose
%forward and backward upper and lower boundaries of num_std window around
%the smooth average
num_std = 3;
alph = 0.2; %transparency of the patch
X = [TDAY TDAY(end:-1:1)];
Y = [LOW_CT+num_std*LOW_STD LOW_CT(end:-1:1)-num_std*LOW_STD(end:-1:1)];
hold on; patch(X,Y, 'k-', 'EdgeAlpha', alph, 'FaceAlpha', alph)


 YMAX=90; 
 YMIN=20;
win=monthday(12)+21;
sum=monthday(6)+21;
axis([-ext 365+ext YMIN YMAX]); grid on; box on;
tday_min=find(LOW_CT == min(LOW_CT)) %where the smoothed seasonal cycle is a maximum
tday_max=find(LOW_CT == max(LOW_CT)) %where the smoothes seasonal cycle is a minnimum
hold on; plot([0 0], [YMIN YMAX], 'k-','LineWidth', LWA); %mark the new year
hold on; plot([365 365], [YMIN YMAX], 'k-','LineWidth', LWA); % new year
hold on; plot(tday_min*[1 1], [YMIN YMAX], 'b-','LineWidth', LWA);  % plot the time of minnimum
hold on; plot(tday_max*[1 1], [YMIN YMAX], 'r-','LineWidth', LWA); % and maximum
hold on; plot(win*[1 1], [YMIN YMAX], 'b--','LineWidth', LWA); % winter solstice (December 21)
hold on; plot(sum*[1 1], [YMIN YMAX], 'r--','LineWidth', LWA); % summer solstice (June 21)
hold on; plot((win-365)*[1 1], [YMIN YMAX], 'b--','LineWidth', LWA);
hold on; plot((tday_min-365)*[1 1], [YMIN YMAX], 'b-','LineWidth', LWA);

[fund_day] = fund_phase(LOW_CT); % find the phase of the maximum of the first (annual) harmonic fit to the annual cycle (the "fundamental")
[fund_amp] = fund(LOW_CT); % amplitude of fundamental

fund_avg = mean(LOW_CT)+fund_amp*cos(2*pi*(TDAY-fund_day)/365); %recompose the best fit to the annual harmonic

%plot the harmonic along with the time of maximum and minnimum (which
hold on; plot([TDAY-365 TDAY TDAY+365],[fund_avg fund_avg fund_avg], 'm-', 'LineWidth', LW);
hold on; plot(fund_day*[1 1], [YMIN YMAX], 'm--', 'LineWidth', LWA);
hold on; plot((fund_day-365/2)*[1 1], [YMIN YMAX], 'm--', 'LineWidth', LWA);

set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA);
xlabel('Day', 'FontSize', FSz, 'FontName', FN);
ylabel('Temperature (F)', 'FontSize', FSz, 'FontName', FN);





