% CDIP SPECTRUM PLOT - MATLAB

% Plot Wave Energy Density (Ed) and Direction (Dmean) using OPeNDAP service from CDIP THREDDS Server.

%-----------------------------------------------------------------------
%%% PRE-CODE DOWNLOAD INSTRUCTIONS %%%
% You must follow these steps in order to access data directly from NetCDF files
%	- Download the nctoolbox from: https://github.com/nctoolbox
%		1. Download and extract into a folder
%		2. Open matlab and at the matlab prompt cd to the extracted folder where setup_nctoolbox lives. (Use pwd to see where you currently are).
%		3. Run setup_nctoolbox
%		4. Optionally edit startup.m and add lines in as described in the install instructions.
%-----------------------------------------------------------------------

% CHANGE DIRECTORY TO 'NCTOOLBOX' DIRECTORY
cd nctoolbox

% OPEN NCTOOLBOX IN MATLAB
setup_nctoolbox



% USER ENTERS STATION NUMBER AND START DATE FOR PLOT
stn = '029';
startdate = '01/15/2013 00:00';


% CONNECT TO THREDDS SERVER AND OPEN NETCDF FILE
urlbase = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/';  % Set 'base' THREDDS URL and pieces to concatenate with user-defined station number (set above)
p1 = 'p1/';
urlend = 'p1_historic.nc';
dsurl = strcat(urlbase,stn,p1,stn,urlend);
ds = ncdataset(dsurl);

% PRINT LIST OF VARIABLES IN NETCDF FILE
% ds.variables

% GET BUOY NAME AND TRANSPOSE TO HORIZONTAL STRING
buoyname = ds.data('metaStationName');
buoytitle = transpose(buoyname);

% EXTRACT MONTH/YEAR FROM 'STARTDATE' 
monthtitle = datestr(datenum(startdate,'mm/dd/yyyy'),'mmm yyyy'); % Create a string object of Month and Year from startdate


% CONVERT START DATE TO MATLAB SERIAL NUMBER
startser = datenum(startdate, 'mm/dd/yyyy HH:MM');


% CALL TIME VARIABLE 
timevar = ds.data('waveTime');
timeconvert = ds.time('waveTime',timevar); % Convert UNIX timestamps to Matlab serial units

% FIND INDEX NUMBERS OF CLOSEST TIME ARRAY OBJECT TO START/END DATES
diffstart = abs(timeconvert-startser); % Compute the difference between the 'startdate' serial number and every serial number in 'timevar', to determine which difference is smallest (which number most closely matches 'startser')
[startidx startidx] = min(diffstart);  % index of closest object to 'startser', based on the minimum value in 'diffstart' differences array
% closeststart = timeconvert(startidx)  % value of closest object to 'startser' (optional - not used in code)

timestep = 1/48;  % calculate one 30-minute period as a fraction of one Matlab serial day

endidx = startidx + timestep;  % index of enddate object (as 30-minute step from startdate object)
% closestend = timeconvert(endidx);  % value of enddate object (as 30-minute step from startdate object) (optional - not used in code)

startidx2 = [startidx 1]; % Set new start/end index based on 2D [time frequency] variable setup of Ed and Dmean variables
endidx2 = [endidx 64];


% ASSIGN NAMES TO VARIABLES TO BE USED IN SPECTRUM PLOT, AND SPECIFY TIMERANGE
Ed = ds.data('waveEnergyDensity',startidx2,endidx2);
Dmean = ds.data('waveMeanDirection',startidx2,endidx2);
Fq = ds.data('waveFrequency'); % Don't need to specify start/end index values, because you want to call whole array (64 frequency bands)



%% SET UP FIGURE FOR PLOTTING
   
% MAKE SUBPLOTS AND ADD DATA

xvals1 = Fq; % define x-axis variable for Subplot 1
yvals1 = Ed; % define y-variable for Subplot 1
xvals2 = Fq; % define x-axis variable for Subplot 2
yvals2 = Dmean; % define y-variable for Subplot 2

figure % Create new figure

% FIRST SUBPLOT
subplot(2,1,1)
plot(xvals1,yvals1)
ylim([0 10])
grid on
grid minor
title(monthtitle)
xlabel('Frequency (Hz)')
ylabel('Energy Density (m^2/Hz)')

% SECOND SUBPLOT
subplot(2,1,2)
plot(xvals2,yvals2)
ylim([0 360])
grid on
grid minor
xlabel('Period (s)')
ylabel('Direction (deg)')


suptitle(buoytitle) % Set main plot title