clear;
% Example use of the FDPA_Park_Attributes.m code.

% STEP 1 - COLLECT THE DATA
%   Access seismic data to work with. For this example we will use 4 hours
%   of uncorrected data, though one should use acausal (phase-preserving)
%   correction techniques.
%
%   Data used in this example is from seismometer BK ORV, and data is
%   accessed through the International Federation of Digital Seismograph
%   Networks (FDSN). The URL is built through the Northern California
%   Earthquake Data Center (NCEDC) 
%   http://service.ncedc.org/fdsnws/dataselect/1/#urlbuilder
% 
%   Net = BK
%   sta = ORV
%   cha = BHZ/BHN/BHE (iteratively)
%   loc = 00
%   start = 2017-02-15T00:00:00
%   end = 2017-02-15T04:00:00
%
%   At the selected time (February 15th, 2017 at 4pm Pacific Time), 2825
%   cubic meters per second (99791 cubic feet per second) was released
%   through the damaged outflow spillway. Discharge data accessed from the
%   California Department of Water Resources California Data Exchange
%   Center at cdec.water.ca.gov. 
% 
% STEP 2 - READ THE DATA
%   Read in Data using rdmseed (François Beauducel, 2015; accessed at 
%   https://www.mathworks.com/matlabcentral/fileexchange/28803-
%   rdmseed-and-mkmseed--read-and-write-miniseed-files);
%
%   The first column of each channel is in Matlab DateTime (UTC) and the
%   second column is the instrument output in counts.

X = rdmseed('BK_ORV_BHZ.mseed');
ORV_BHZ(:,1) = cat(1,X.t);
ORV_BHZ(:,2) = cat(1,X.d);
clearvars X

X = rdmseed('BK_ORV_BHE.mseed');
ORV_BHE(:,1) = cat(1,X.t);
ORV_BHE(:,2) = cat(1,X.d);

clearvars X

X = rdmseed('BK_ORV_BHN.mseed');
ORV_BHN(:,1) = cat(1,X.t);
ORV_BHN(:,2) = cat(1,X.d);

clearvars X

startdatetime= datestr(ORV_BHZ(1,1)); 
enddatetime = datestr(ORV_BHZ(end,1));

% The Sampling frequency (Hz) for the data is computed from timeseries
Fs = 1/mean(diff(ORV_BHZ(:,1)))/86400; 

% STEP 3 - POLARIZATION ANALYSIS 

%   Call the polarization code- refer to code for details

Z = ORV_BHZ(:,2); 
E = ORV_BHE(:,2);
N = ORV_BHN(:,2);
 
Nwin = 19; %  Select 19 as the number of subwindows
%             (following Koper and Hawley, 2010).

[Evec, Eval, freqs, UD_power, EW_power, NS_power, First_eigenv_power, deg_of_pol, ThetaH, ThetaV, phiHH, phiVH] = FDPA_Park_Attributes(Z,E,N,Fs,Nwin);

% STEP 4
%   Plot the results for a single hour as a scatter plot. Compare to 
%   figure provided with code and results in Goodling et al. (2017).

figure; 
subplot(5,1,1); plot(freqs,real(deg_of_pol),'k.'); ylabel('Degree of Polarization');xlabel('Frequency  (Hz)');
set(gca,'YTick',[0 0.5 1]);box on; grid on;ylim([0 1]);

subplot(5,1,2); plot(freqs,ThetaH,'k.'); ylabel('\Theta_H');xlabel('Frequency (Hz)');
set(gca,'YTick',[0 90 180 270 360]);box on; grid on;ylim([0 360]);

subplot(5,1,3); plot(freqs,ThetaV,'k.'); ylabel('\Theta_V');xlabel('Frequency (Hz)');
set(gca,'YTick',[0 45 90]);box on; grid on;ylim([0 90]);

subplot(5,1,4); plot(freqs,phiVH,'k.'); ylabel('\phi_V_H');xlabel('Frequency (Hz)');
set(gca,'YTick',[-90 -45 0 45 90]);box on; grid on;ylim([-90 90]);

subplot(5,1,5); plot(freqs,phiHH,'k.'); ylabel('\phi_H_H');xlabel('Frequency (Hz)');
set(gca,'YTick',[-180 -90 0 90 180]);box on; grid on;ylim([-180 180]);


suptitle({'Polarization Results, BK ORV', [startdatetime ' to ' enddatetime]});

figure;

plot(freqs,UD_power,'b');hold on
plot(freqs,NS_power,'r');
plot(freqs,EW_power,'g');
plot(freqs,First_eigenv_power,'k');
legend('Up-Down Component','North-South Component',...
       'East-West Component', 'First Eigenvector','location','SouthWest');
xlabel('Frequency (Hz)');ylabel('Power');

title('Power Spectral Density of Components and First Eigenvector');