%Stats_Current.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_nan(u_at_depth_h(:,index));
[v_at_depth_h(:,index)] = ridfill_nan(v_at_depth_h(:,index));

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

%Convert to speed and direction
c_speed_h(:,index) = sqrt(u_at_depth_h(:,index).^2 + v_at_depth_h(:,index).^2);
c_dir_h(:,index) = atan2(v_at_depth_h(:,index),u_at_depth_h(:,index))*(180/pi);
    
c_speed_lp(:,index) = sqrt(u_at_depth_lp(:,index).^2 + v_at_depth_lp(:,index).^2);
c_dir_lp(:,index) = atan2(v_at_depth_lp(:,index),u_at_depth_lp(:,index))*(180/pi);
    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(['Points = ' num2str(length(find(isnan(u_at_depth_h(:,index))==0))) '\r']);
fprintf(['Mean U = ' num2str(nanmean(u_at_depth_h(:,index))) '\r']);
fprintf(['St. Dev. U = ' num2str(nanstd(u_at_depth_h(:,index))) '\r']);
fprintf(['St. Dev. U (lp) = ' num2str(nanstd(u_at_depth_lp(:,index))) '\r']);
fprintf(['Mean V = ' num2str(nanmean(v_at_depth_h(:,index))) '\r']);
fprintf(['St. Dev. V = ' num2str(nanstd(v_at_depth_h(:,index))) '\r']);
fprintf(['St. Dev. V (lp) = ' num2str(nanstd(v_at_depth_lp(:,index))) '\r']);
fprintf(['Mean Speed = ' num2str(nanmean(c_speed_h(:,index))) '\r']);
fprintf(['St. Dev. Speed = ' num2str(nanstd(c_speed_h(:,index))) '\r']);
fprintf(['St. Dev. Speed (lp) = ' num2str(nanstd(c_speed_lp(:,index))) '\r']);
fprintf(['Max Speed = ' num2str(max(c_speed_h(:,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_nan(u_at_depth_h);
[v_at_depth_h] = ridfill_nan(v_at_depth_h);

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

%Convert to speed and direction
c_speed_h = sqrt(u_at_depth_h.^2 + v_at_depth_h.^2);
c_dir_h = atan2(v_at_depth_h,u_at_depth_h)*(180/pi);
    
c_speed_lp = sqrt(u_at_depth_lp.^2 + v_at_depth_lp.^2);
c_dir_lp = atan2(v_at_depth_lp,u_at_depth_lp)*(180/pi);

%Assemble titles for wind and currents
fprintf(['BASS, Depth = ' num2str(act_depth_h) '\r']);
fprintf(['Points = ' num2str(length(find(isnan(u_at_depth_h)==0))) '\r']);
fprintf(['Mean U = ' num2str(nanmean(u_at_depth_h)) '\r']);
fprintf(['St. Dev. U = ' num2str(nanstd(u_at_depth_h)) '\r']);
fprintf(['St. Dev. U (lp) = ' num2str(nanstd(u_at_depth_lp)) '\r']);
fprintf(['Mean V = ' num2str(nanmean(v_at_depth_h)) '\r']);
fprintf(['St. Dev. V = ' num2str(nanstd(v_at_depth_h)) '\r']);
fprintf(['St. Dev. V (lp) = ' num2str(nanstd(v_at_depth_lp)) '\r']);
fprintf(['Mean Speed = ' num2str(nanmean(c_speed_h)) '\r']);
fprintf(['St. Dev. Speed = ' num2str(nanstd(c_speed_h)) '\r']);
fprintf(['St. Dev. Speed (lp) = ' num2str(nanstd(c_speed_lp)) '\r']);
fprintf(['Max Speed = ' num2str(max(c_speed_h)) '\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_nan(u_at_depth_h);
[v_at_depth_h] = ridfill_nan(v_at_depth_h);

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

%Convert to speed and direction
c_speed_h = sqrt(u_at_depth_h.^2 + v_at_depth_h.^2);
c_dir_h = atan2(v_at_depth_h,u_at_depth_h)*(180/pi);
    
c_speed_lp = sqrt(u_at_depth_lp.^2 + v_at_depth_lp.^2);
c_dir_lp = atan2(v_at_depth_lp,u_at_depth_lp)*(180/pi);

%Assemble titles for wind and currents
fprintf(['BASS (2), Depth = ' num2str(act_depth_h) '\r']);
fprintf(['Points = ' num2str(length(find(isnan(u_at_depth_h)==0))) '\r']);
fprintf(['Mean U = ' num2str(nanmean(u_at_depth_h)) '\r']);
fprintf(['St. Dev. U = ' num2str(nanstd(u_at_depth_h)) '\r']);
fprintf(['St. Dev. U (lp) = ' num2str(nanstd(u_at_depth_lp)) '\r']);
fprintf(['Mean V = ' num2str(nanmean(v_at_depth_h)) '\r']);
fprintf(['St. Dev. V = ' num2str(nanstd(v_at_depth_h)) '\r']);
fprintf(['St. Dev. V (lp) = ' num2str(nanstd(v_at_depth_lp)) '\r']);
fprintf(['Mean Speed = ' num2str(nanmean(c_speed_h)) '\r']);
fprintf(['St. Dev. Speed = ' num2str(nanstd(c_speed_h)) '\r']);
fprintf(['St. Dev. Speed (lp) = ' num2str(nanstd(c_speed_lp)) '\r']);
fprintf(['Max Speed = ' num2str(max(c_speed_h)) '\n\n']);
end