%Ellipse_Calc.m
%
%Program which calculated the ellipse parameters from a given site at a given depth.
%Prints out a result which can be copied into notepad to save as a .txt file; user
%also chooses if low-pass filtered or hour-averaged data is used.
%Calls singleJD.m, value2Index.m, ridfill.m, gregorian.m, & pcaben.m
%Modified from a program (bot_elp.m) received from Ben Guitterez
%
%Soupy Alexander
%12/7/2001

ncclear

station = input('Choose a station: ');
hit_depth = input('Choose a depth (B = Bottom): ');
data_type = input('Choose a data type (h = hourly averaged, l = low-pass filtered: ');

%Differentation of hour averaged vs. low-pass filtered
if data_type == 'h';
    dataTag = '-a1h';
elseif data_type == 'l';
    dataTag = '-alp';
else
    error('Non-existent data type.  Be sure to use lower case letter and single quotes')
end

%Selection of data files
if station == 'A' & hit_depth ~= 'B';
    fileName = ['5951adc' dataTag '.nc'];
elseif station == 'A' & hit_depth == 'B';
    fileName = ['5952v' dataTag '_d1.nc'];
elseif station == 'B' & hit_depth ~= 'B';
    fileName = ['5971adc' dataTag '.nc'];
elseif station == 'B' & hit_depth == 'B';
    fileName = ['5972v' dataTag '_d1.nc'];
elseif station == 'C' & hit_depth ~= 'B';
    fileName = ['5991adc' dataTag '.nc'];
elseif station == 'C' & hit_depth == 'B';
    fileName = ['5992v' dataTag '.nc'];
elseif station == 'D' & hit_depth ~= 'B';
    fileName = ['6011adc' dataTag '.nc'];
elseif station == 'D' & hit_depth == 'B';
    fileName = ['6012v' dataTag '_d1.nc'];
elseif station == 'E' & hit_depth ~= 'B';
    fileName = ['6031adc' dataTag '.nc'];
elseif station == 'E' & hit_depth == 'B';
    fileName = ['6032v' dataTag '.nc'];
elseif station == 'F';
    fileName = ['6041adc' dataTag '.nc'];
else
    error('Non-existent station.  Be sure to use capital letter and single quotes')
end    

%Extract useful information
file = netCDF(fileName,'nowrite');
time = file{'time'};
time2 = file{'time2'};
jdTime = singleJD(time,time2);
u_data = file{'u_1205'}(:);
v_data = file{'v_1206'}(:);
depth_vector = file{'depth'}(:);

%Determine the index of the depth closest to the desired depth
if hit_depth ~= 'B';
    [depth_index,dist_to_depth] = value2Index(depth_vector,hit_depth);
        if dist_to_depth > 5;
            error('No depth data available within 5 meters of selected depth')
        end
elseif hit_depth == 'B';
    depth_index = 1;
end    
        
u_data_hit = u_data(:,depth_index);
v_data_hit = v_data(:,depth_index);
act_depth = depth_vector(depth_index);

u_data_hit = ridfill(u_data_hit);
v_data_hit = ridfill(v_data_hit);

%Calculate the current speed and direction
cspd = sqrt(u_data_hit.^2 + v_data_hit.^2);
cdir = atan2(v_data_hit,u_data_hit)*(180/pi);

%Find the gregorian date vectors
greg_time = gregorian(jdTime);

Y = greg_time(:,1);
M = greg_time(:,2);
D = greg_time(:,3);
H = greg_time(:,4);
MI = greg_time(:,5);
S = greg_time(:,6);

%define start and stop dates for output.
date1 = datestr(datenum(Y(1),M(1),D(1),H(1),MI(1),S(1)),0);
date2 = datestr(datenum(Y(end),M(end),D(end),H(end),MI(end),S(end)),0);

%Start printing outputs
fprintf('\nEllipse Parameters\n\n\n')
if data_type == 'h',
   
   fprintf(' Location: New York Bight\n')
   fprintf(' Data: Hour-averaged\n')
   fprintf([' Depth: ' num2str(act_depth,'%.2f') ' m\n\n'])
   
elseif data_type == 'l',
   
   fprintf(' Location: New York Bight\n')
   fprintf(' Data: Low-pass filtered\n')
   fprintf([' Depth: ' num2str(act_depth,'%.2f') ' m\n\n'])
end

fprintf('Start date: %15s\nStop date: %15s\n\n',date1,date2)
fprintf('Month\t\tPoints\t\tMajor\t\tMinor\t\tMajor\t\tEllipticity\n')
fprintf('\t\t\t\tAxis\t\tAxis\t\tOrientation\n')

if data_type == 'h';
fprintf('\t\t(Hours)\t\t(cm/s)\t\t(cm/s)\t\t(degrees)\n\n')
elseif data_type == 'l';
fprintf('\t\t(Points)\t\t(cm/s)\t\t(cm/s)\t\t(degrees)\n\n')
end

for j = [12 1 2 3 4]
   thismo = find(M==j);
   if isempty(thismo)| thismo ==1, 
   fprintf('%3.0f\t%6s\n',j,' No Data')
   else
            dlen = length(u_data_hit(thismo));
   [A,B,C,D3,E] = pcaben(u_data_hit(thismo),v_data_hit(thismo));
   fprintf('%3.0f\t\t%3.0f\t\t%5.2f\t\t%4.2f\t\t%3.0f\t\t%5.2f\n',j,dlen,A,C,B,E)
   end
end
fprintf('\n')

dlen = length(u_data_hit);

[A,B,C,D3,E] = pcaben(u_data_hit,v_data_hit);

fprintf('All\t\t%3.0f\t\t%5.2f\t\t%4.2f\t\t%3.0f\t\t%5.2f\n',dlen,A,C,B,E)