%Stats_Ellipse.m
%
%Program to calculate the current statistics for the New York
%Bight.  Requires the input of a station.
%Requires singleJD.m, value2Index.m, gmean.m, gmin.m, gmax.m, gstd.m,
%
%Modified from mid_sc_stat.m by Ben Guitterez
%
%Soupy Alexander, 2/1/2002

station = input('Enter the station of interest: ');


%Select the correct files.
if station == 'A';
    curr_file_h = ['5951adc-a1h.nc'];
    curr_file_lp = ['5951adc-alp.nc'];
    bass_file_h = ['5952v-a1h_d1.nc'];
    bass_file_lp = ['5952v-alp_d1.nc'];
    bass_file2_h = ['5952v-a1h_d2.nc'];
    bass_file2_lp = ['5952v-alp_d2.nc'];
elseif station == 'B';
    curr_file_h = ['5971adc-a1h.nc'];
    curr_file_lp = ['5971adc-alp.nc'];
    bass_file_h = ['5972v-a1h_d1.nc'];
    bass_file_lp = ['5972v-alp_d1.nc'];
    bass_file2_h = ['5972v-a1h_d2.nc'];
    bass_file2_lp = ['5972v-alp_d2.nc'];
elseif station == 'C';
    curr_file_h = ['5991adc-a1h.nc'];
    curr_file_lp = ['5991adc-alp.nc'];
    bass_file_h = ['5992v-a1h.nc'];
    bass_file_lp = ['5992v-alp.nc'];
elseif station == 'D';
    curr_file_h = ['6011adc-a1h.nc'];
    curr_file_lp = ['6011adc-alp.nc'];    
    bass_file_h = ['6012v-a1h_d1.nc'];
    bass_file_lp = ['6012v-alp_d1.nc'];
    bass_file2_h = ['6012v-a1h_d2.nc'];
    bass_file2_lp = ['6012v-alp_d2.nc'];    
elseif station == 'E';
    curr_file_h = ['6031adc-a1h.nc'];
    curr_file_lp = ['6031adc-alp.nc'];
    bass_file_h = ['6032v-a1h.nc'];
    bass_file_lp = ['6032v-alp.nc'];
elseif station == 'F';
    curr_file_h = ['6041adc-a1h.nc'];
    curr_file_lp = ['6041adc-alp.nc']; 
else
    error('Incorrect station ID.  Input in capital letters with single quotes');
end


%Pull out the current data from the adcp
curr_file_h = netcdf(curr_file_h,'nowrite');
curr_file_lp = netcdf(curr_file_lp,'nowrite');

curr_t_h = curr_file_h{'time'}(:);
curr_t2_h = curr_file_h{'time2'}(:);
curr_time_h = singleJD(curr_t_h,curr_t2_h);

curr_t_lp = curr_file_lp{'time'}(:);
curr_t2_lp = curr_file_lp{'time2'}(:);
curr_time_lp = singleJD(curr_t_lp,curr_t2_lp);

u_data_h = curr_file_h{'u_1205'}(:);
v_data_h = curr_file_h{'v_1206'}(:);

u_data_lp = curr_file_lp{'u_1205'}(:);
v_data_lp = curr_file_lp{'v_1206'}(:);

curr_depth_h = curr_file_h{'depth'}(:);
curr_depth_lp = curr_file_lp{'depth'}(:);

%Select out the current data from the proper depth
depth_ind_h(1) = value2Index(curr_depth_h,5);
depth_ind_h(2) = value2Index(curr_depth_h,15);
depth_ind_h(3) = value2Index(curr_depth_h,30);
depth_ind_h(4) = value2Index(curr_depth_h,45);

depth_ind_lp(1) = value2Index(curr_depth_lp,5);
depth_ind_lp(2) = value2Index(curr_depth_lp,15);
depth_ind_lp(3) = value2Index(curr_depth_lp,30);
depth_ind_lp(4) = value2Index(curr_depth_lp,45);

for index = 1:4;
    if isnan(depth_ind_h(index)) == 0;
u_at_depth_h(:,index) = u_data_h(:,depth_ind_h(index));
v_at_depth_h(:,index) = v_data_h(:,depth_ind_h(index));	
act_depth_h(:,index) = curr_depth_h(depth_ind_h(index));

u_at_depth_lp(:,index) = u_data_lp(:,depth_ind_lp(index));
v_at_depth_lp(:,index) = v_data_lp(:,depth_ind_lp(index));	
act_depth_lp(:,index) = curr_depth_lp(depth_ind_lp(index));

[u_at_depth_h(:,index)] = ridfill(u_at_depth_h(:,index));
[v_at_depth_h(:,index)] = ridfill(v_at_depth_h(:,index));

[u_at_depth_lp(:,index)] = ridfill(u_at_depth_lp(:,index));
[v_at_depth_lp(:,index)] = ridfill(v_at_depth_lp(:,index));

%Calculate principal axes information
[majorax_h(:,index),majoraz_h(:,index),minorax_h(:,index),...
    minoraz_h(:,index),ellip_h(:,index)] = ...
    pcaben(u_at_depth_h(:,index),v_at_depth_h(:,index));
[majorax_lp(:,index),majoraz_lp(:,index),minorax_lp(:,index),...
    minoraz_lp(:,index),ellip_lp(:,index)] = ...
    pcaben(u_at_depth_lp(:,index),v_at_depth_lp(:,index));
end
end
%Assemble titles for currents
for index = 1:4;
    if isnan(depth_ind_h(index)) == 0;
fprintf(['ADCP, Depth = ' num2str(act_depth_h(:,index)) '\r']);
fprintf(['Major (h) = ' num2str(majorax_h(:,index)) '\r']);
fprintf(['Minor (h) = ' num2str(minorax_h(:,index)) '\r']);
fprintf(['Major Orient (h) = ' num2str(majoraz_h(:,index)) '\r']);
fprintf(['Ellip. (h) = ' num2str(ellip_h(:,index)) '\r']);
fprintf(['Major (lp) = ' num2str(majorax_lp(:,index)) '\r']);
fprintf(['Minor (lp) = ' num2str(minorax_lp(:,index)) '\r']);
fprintf(['Major Orient (lp) = ' num2str(majoraz_lp(:,index)) '\r']);
fprintf(['Ellip. (lp)= ' num2str(ellip_lp(:,index)) '\n\n']);
    end
end

if station ~= 'F';
%Pull out the current data from the bass
bass_file_h = netcdf(bass_file_h,'nowrite');
bass_file_lp = netcdf(bass_file_lp,'nowrite');

bass_t_h = bass_file_h{'time'}(:);
bass_t2_h = bass_file_h{'time2'}(:);
bass_time_h = singleJD(bass_t_h,bass_t2_h);

bass_t_lp = bass_file_lp{'time'}(:);
bass_t2_lp = bass_file_lp{'time2'}(:);
bass_time_lp = singleJD(bass_t_lp,bass_t2_lp);

u_data_h = bass_file_h{'u_1205'}(:);
v_data_h = bass_file_h{'v_1206'}(:);

u_data_lp = bass_file_lp{'u_1205'}(:);
v_data_lp = bass_file_lp{'v_1206'}(:);

bass_depth_h = bass_file_h{'depth'}(:);
bass_depth_lp = bass_file_lp{'depth'}(:);

u_at_depth_h = u_data_h;
v_at_depth_h = v_data_h;	
act_depth_h = bass_depth_h;

u_at_depth_lp = u_data_lp;
v_at_depth_lp = v_data_lp;	
act_depth_lp = bass_depth_lp;

[u_at_depth_h] = ridfill(u_at_depth_h);
[v_at_depth_h] = ridfill(v_at_depth_h);

[u_at_depth_lp] = ridfill(u_at_depth_lp);
[v_at_depth_lp] = ridfill(v_at_depth_lp);

%Pull principal axes info.
[majorax_h,majoraz_h,minorax_h,minoraz_h,ellip_h] = ...
    pcaben(u_at_depth_h,v_at_depth_h);
[majorax_lp,majoraz_lp,minorax_lp,minoraz_lp,ellip_lp] = ...
    pcaben(u_at_depth_lp,v_at_depth_lp);

%Assemble titles for wind and currents
fprintf(['BASS, Depth = ' num2str(act_depth_h) '\r']);
fprintf(['Major (h) = ' num2str(majorax_h) '\r']);
fprintf(['Minor (h) = ' num2str(minorax_h) '\r']);
fprintf(['Major Orient (h) = ' num2str(majoraz_h) '\r']);
fprintf(['Ellip. (h) = ' num2str(ellip_h) '\r']);
fprintf(['Major (lp) = ' num2str(majorax_lp) '\r']);
fprintf(['Minor (lp) = ' num2str(minorax_lp) '\r']);
fprintf(['Major Orient (lp) = ' num2str(majoraz_lp) '\r']);
fprintf(['Ellip. (lp)= ' num2str(ellip_lp) '\n\n']);
end

if station == 'A' | station == 'B' | station == 'D';
%Pull out the current data from the bass
bass_file2_h = netcdf(bass_file2_h,'nowrite');
bass_file2_lp = netcdf(bass_file2_lp,'nowrite');

bass_t_h = bass_file2_h{'time'}(:);
bass_t2_h = bass_file2_h{'time2'}(:);
bass_time_h = singleJD(bass_t_h,bass_t2_h);

bass_t_lp = bass_file2_lp{'time'}(:);
bass_t2_lp = bass_file2_lp{'time2'}(:);
bass_time_lp = singleJD(bass_t_lp,bass_t2_lp);

u_data_h = bass_file2_h{'u_1205'}(:);
v_data_h = bass_file2_h{'v_1206'}(:);

u_data_lp = bass_file2_lp{'u_1205'}(:);
v_data_lp = bass_file2_lp{'v_1206'}(:);

bass_depth_h = bass_file2_h{'depth'}(:);
bass_depth_lp = bass_file2_lp{'depth'}(:);

u_at_depth_h = u_data_h;
v_at_depth_h = v_data_h;	
act_depth_h = bass_depth_h;

u_at_depth_lp = u_data_lp;
v_at_depth_lp = v_data_lp;	
act_depth_lp = bass_depth_lp;

[u_at_depth_h] = ridfill(u_at_depth_h);
[v_at_depth_h] = ridfill(v_at_depth_h);

[u_at_depth_lp] = ridfill(u_at_depth_lp);
[v_at_depth_lp] = ridfill(v_at_depth_lp);

%Pull principal axes info.
[majorax_h,majoraz_h,minorax_h,minoraz_h,ellip_h] = ...
    pcaben(u_at_depth_h,v_at_depth_h);
[majorax_lp,majoraz_lp,minorax_lp,minoraz_lp,ellip_lp] = ...
    pcaben(u_at_depth_lp,v_at_depth_lp);

%Assemble titles for wind and currents
fprintf(['BASS (2), Depth = ' num2str(act_depth_h) '\r']);
fprintf(['Major (h) = ' num2str(majorax_h) '\r']);
fprintf(['Minor (h) = ' num2str(minorax_h) '\r']);
fprintf(['Major Orient (h) = ' num2str(majoraz_h) '\r']);
fprintf(['Ellip. (h) = ' num2str(ellip_h) '\r']);
fprintf(['Major (lp) = ' num2str(majorax_lp) '\r']);
fprintf(['Minor (lp) = ' num2str(minorax_lp) '\r']);
fprintf(['Major Orient (lp) = ' num2str(majoraz_lp) '\r']);
fprintf(['Ellip. (lp)= ' num2str(ellip_lp) '\n\n']);
end