%WIND_TRACK  Script to calculate ADCIRC tides and low-passed tide gauge (wind effect) 

echo on

track=load('95007_bathy.dat');
depth=track(:,9);
%
% Run Median Filter over depth data - 15 filter window
%
depth_filt=medfilt(depth,15);

% find the tidal correction from ADCIRC tidal model
% ztide is the tidal correction

[ztide,lon,lat,jd_track]=tide_track(track);  


% Below are the 'routines' for determining the wind compenent of 
% the two tidal stations

% Load the observed tide gauge data

a=load('95007_SandyHook.dat'); 
b=load('95007_Montauk.dat'); 

% form the 6 column Gregorian date vector
g=a(:,1:5);
g(:,6)=g(:,5)*0;

% convert to decimal Julian Day 
jd_tide=julian(g);

% in this case, the tide gauge data is
% relative to MLLW, and we want MSL, so instead
% of downloading the data again, we simply figured
% out what the offset was between MSL and MLLW
% 
%  By default - will use MLLW as datum for bathy
%  corresponds to all of USGS older bathy
%  Can change this at a later date if necessary
%  by applying a constant offset to the data
%
% aoffset=0.780;  % MSL-MLLW at Sandy Hook
% boffset=0.387;  % MSL-MLLW at Montauk
a_lat=40+28/60;
a_lon=-(74+.6/60);
b_lat=41+2.9/60;
b_lon=-(71+57.6/60);

% low-pass the tide gauge data: low-pass filter these data
% at a greater time period than the tides; effectively 'filters out'
% the tides - leaving the wind component at the tide gauge

alp=pl33tn(a(:,6),6/60,33); 
blp=pl33tn(b(:,6),6/60,33);

% convert to MSL

% alp=alp-aoffset; 
% blp=blp-boffset; 

% interpolate tide gauge data onto 
% ship track time base (ship track = 10 sec; tide gauge = 6 min)

atrack=interp1(jd_tide,alp,jd_track);
btrack=interp1(jd_tide,blp,jd_track);

% calculate relative position of the track points 
% between Sandy Hook and Montauk
% so that we can weight the wind component from
% the nearest tide gauge more.

% simple distance calc

dx=cos(a_lat*pi/180)*(lon-a_lon);
dy=(lat-a_lat);
dist1=sqrt(dx.^2+dy.^2);

dx=cos(b_lat*pi/180)*(lon-b_lon);
dy=(lat-b_lat);
dist2=sqrt(dx.^2+dy.^2);

% wind component

z_wind=atrack+(btrack-atrack).*dist1./(dist2+dist1);


% The final corrections will be depth_filt - tide_corr - wind_component
% depth_corr = ((depth_filt - ztide) - z_wind)

depth_corr = ((depth_filt - ztide) - z_wind);

%  Now, add 1 meter to the corrected depth values to account
% for distance of transducer below water line (draft correction)
%

depth_draft = (depth_corr + 1.0);

% The bathy data need to be filtered - working with Rich - looks
% like mean filter (boxcar 10, e.g.) works well - noise is 
% removed, but filtered data remain true to original
%
% Chuck has written an M file that will look at the cross points
% (where track lines cross) in order to evaluate the Z values
% at these overlapping areas.  This will enable us to compare
% the depth_corr with the depth_filt to see how effective the
% tidal correction was.
%

echo off