%Stats_Curr_Temp.m
%
%Program to calculate the temperature and 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, 12/17/2001

monthvector = [12 1 2 3 4];

station = input('Enter the station of interest: ');
data_type = input('Hourly Averaged (h) or Low-Pass Filtered (l)?: ');
%%NEED TO SELECT DEPTHS!!!!!!!!!!
hit_depth_type = ...
    input('Enter depth of interest (surface = s, mid(1) = m1, mid(2) = m2, d = deep): ');

if data_type == 'h';
    data_tag = '-a1h.nc';
elseif data_type == 'l';
    data_tag = '-alp.nc';
else
    error('Incorrect data type.  Input "h" or "l" in single quotes');
end


%Select the correct files.
if station == 'A';
    curr_file = ['5951adc' data_tag];
    if hit_depth_type == 's';
        temp_file = ['5941mc' data_tag];
    elseif hit_depth_type == 'd';
        temp_file = ['5952tcp' data_tag];
    elseif hit_depth_type == 'm1'| hit_depth_type == 'm2';
        error('No Mid-Depth Temperature for this station!!!');
    end
elseif station == 'B';
    curr_file = ['5971adc' data_tag];
    if hit_depth_type == 's';
        temp_file = ['5961mc' data_tag];
    elseif hit_depth_type == 'd';
        temp_file = ['5972tc' data_tag];
    elseif hit_depth_type == 'm1';
        temp_file = ['5962sc' data_tag];
    elseif hit_depth_type == 'm2';
        temp_file = ['5963sc' data_tag];
    end
elseif station == 'C';
    curr_file = ['5991adc' data_tag];
    if hit_depth_type == 's';
        temp_file = ['5981mc' data_tag];
    elseif hit_depth_type == 'd';
        temp_file = ['5992tcp' data_tag];
    elseif hit_depth_type == 'm1';
        temp_file = ['5982sc' data_tag];
    elseif hit_depth_type == 'm2';
        temp_file = ['5983sc' data_tag];
    end
elseif station == 'D';
    curr_file = ['6011adc' data_tag];
    if hit_depth_type == 's';
        temp_file = ['6001mc' data_tag];
    elseif hit_depth_type == 'd';
        temp_file = ['6012tcp' data_tag];
    elseif hit_depth_type == 'm1'| hit_depth_type == 'm2';
        error('No Mid-Depth Temperature for this station!!!');
    end
elseif station == 'E';
    curr_file = ['6031adc' data_tag];
    if hit_depth_type == 's';
        temp_file = ['6021mc' data_tag];
    elseif hit_depth_type == 'd';
        temp_file = ['6032tcp' data_tag];
    elseif hit_depth_type == 'm1'| hit_depth_type == 'm2';
        error('No Mid-Depth Temperature for this station!!!');
    end
elseif station == 'F';
    curr_file = ['6041adc' data_tag];
    if hit_depth_type == 's';
        error('No Surface Temperature for this station!!!');
    elseif hit_depth_type == 'd';
        temp_file = ['6043mc' data_tag];
    elseif hit_depth_type == 'm1'| hit_depth_type == 'm2';
        error('No Mid-Depth Temperature for this station!!!');
    end    
else
    error('Incorrect station ID.  Input in capital letters with single quotes');
end


%Pull out the current data
curr_file = netcdf(curr_file,'nowrite');

curr_t = curr_file{'time'}(:);
curr_t2 = curr_file{'time2'}(:);
curr_time = singleJD(curr_t,curr_t2);

u_data = curr_file{'u_1205'}(:);
v_data = curr_file{'v_1206'}(:);

curr_depth = curr_file{'depth'}(:);

%Select out the current data from the proper depth
if hit_depth_type == 's';
    depth_ind = value2Index(curr_depth,1);
elseif hit_depth_type == 'd';
    depth_ind = find(curr_depth == max(curr_depth));
elseif hit_depth_type == 'm1';
    if station == 'B' | station == 'C';
        depth_ind = value2Index(curr_depth,30);
    else
        error('No Mid-Depth Temperature Data from This Station!!!');
    end
elseif hit_depth_type == 'm2';
    if station == 'B';
        depth_ind = value2Index(curr_depth,45);
    elseif station == 'C';
        depth_ind = value2Index(curr_depth,60);
    else
        error('No Mid-Depth Temperature Data from This Station!!!');
    end
else
    error('Improper depth type selected.');
end

u_at_depth = u_data(:,depth_ind);
v_at_depth = v_data(:,depth_ind);	
act_depth = curr_depth(depth_ind);

[u_at_depth] = ridfill_nan(u_at_depth);
[v_at_depth] = ridfill_nan(v_at_depth);

%Convert to speed and direction
c_speed = sqrt(u_at_depth.^2 + v_at_depth.^2);
c_dir = atan2(v_at_depth,u_at_depth)*(180/pi);
    


%Assemble titles for wind and currents
fprintf('Monthly and Yearly Current and Temperature Statistics\n\n\n')
fprintf('  Location: New York Bight.\n')
if data_type == 'h';
    fprintf('  Data: Hour-averaged\n')
elseif data_type == 'l';
    fprintf('  Data: Low-pass Filtered\n')
end
fprintf('     Depth: %2.1f',act_depth)
fprintf(' m\n\n')
   

%Deal with the time parameters for the input data
greg_vect = gregorian(curr_time);
Y = greg_vect(:,1);
M = greg_vect(:,2);
D = greg_vect(:,3);
H = greg_vect(:,4);
MI = greg_vect(:,5);
S = greg_vect(:,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);

curr_data = [u_at_depth v_at_depth c_speed];

%Calculates and output Monthly statistics to the screen

fprintf('Start date: %15s\nStop date: %15s\n\n',date1,date2)

%predefine output files

		st = zeros(1,5);
      smean = zeros(1,5);
      sstd = zeros(1,5);
      smax = zeros(1,5);
      smin = zeros(1,5);

for i = 1:3,
   if i == 1,
      	ivarname = 'u_1205 Eastward Velocity';
    elseif i == 2,
      	ivarname = 'v_1206 Northward Velocity';
    elseif i == 3,
      	ivarname = 'CS_300 Current Speed';
   	end
   
    fprintf('Variable: %12s\n\n',ivarname)
    fprintf('Month\t\tPoints\t\t   Mean\t\tStandard\tMaximum\tMinimum\n')
    fprintf('\t\t\t\t\t\tDeviation\n')
    if data_type == 'h';
        fprintf('\t\t(Hours)\t\t  (cm/s)\t\t(cm/s)\t\t(cm/s)\t\t(cm/s)\n\n')
    elseif data_type == 'l';
        fprintf('\t\t(Points)\t\t  (cm/s)\t\t(cm/s)\t\t(cm/s)\t\t(cm/s)\n\n')
    end
    
    data = curr_data(:,i);
   		for j = 1:5;
            thismo = find(M==monthvector(j));
            if isempty(thismo) | length(thismo)==1,
               fprintf('%4.0f\t %8s\n',monthvector(j),' No Data')
               dmean = NaN;
               dstd = NaN;
               dmax = NaN;
               dmin = NaN;
               jt = NaN;
            else
               %Calculate and prep statistics for output
      		dlen = length(find(isnan(data(thismo))==0));   
      		dmax = gmax(data(thismo));   
      		dmin = gmin(data(thismo));
            dmean = gmean(data(thismo));			%ignores nans
            dstd = gstd(data(thismo));				%ignores nans
            dmeanday = gmean(D(thismo));
            jt = monthvector(j);
            %Print out results
   			fprintf('%4.0f\t\t%6.0f\t\t%9.1f\t\t%5.1f\t\t%6.1f\t\t%6.1f\n',...
               jt,dlen,dmean,dstd,dmax,dmin)
          	end
         %Patch results into new variables
            st(j) = jt;
            smean(j) = dmean;
            sstd(j) = dstd;
            smax(j) = dmax;
            smin(j) = dmin;

   		end
         fprintf('\n')
         %Reassign the output variables according to the value of i.
         if i == 1,
            ut = st;
            umean = smean;
            ustd = sstd;
            umax = smax;
            umin = smin;
        elseif i == 2,
            vt = st;
            vmean = smean;
            vstd = sstd;
            vmax = smax;
            vmin = smin;
        elseif i == 3,
            cst = st;
            csmean = smean;
            csstd = sstd;
            csmax = smax;
            csmin = smin;
         end
         
%Compute yearly statistics for currents
      	dlen = length(find(isnan(data)==0));   
      	dmax = gmax(data);   
      	dmin = gmin(data);   
      	dmean = gmean(data);
      	dstd = gstd(data);   
     	fprintf('  All\t\t%6.0f\t\t%9.1f\t\t%5.1f\t\t%6.1f\t\t%6.1f\n',...
         dlen,dmean,dstd,dmax,dmin)
         fprintf('\n\n\n')
                  
         yt = 2.5;
         ymean = dmean;
         ystd = dstd;
         ymax = dmax;
         ymin = dmin;
         
         %Reassign the output variables according to the data type.
         if i == 1,
            uyt = yt;
            uymean = ymean;
            uystd = ystd;
            uymax = ymax;
            uymin = ymin;
        elseif i == 2,
            vyt = yt;
            vymean = ymean;
            vystd = ystd;
            vymax = ymax;
            vymin = ymin;
        elseif i == 3,
            csyt = yt;
            csymean = ymean;
            csystd = ystd;
            csymax = ymax;
            csymin = ymin;
         end
end
      
%Now calculate monthly statistics for temperature
tempfile = netcdf(temp_file,'nowrite');
temp_t = tempfile{'time'}(:);
temp_t2 = tempfile{'time2'}(:);
temp_depth = tempfile{'depth'};

temp_time = singleJD(temp_t,temp_t2);

greg_vect = gregorian(temp_time);

M2 = greg_vect(:,2);
D2 = greg_vect(:,3);

temp = tempfile{'T_20'}(:);
temp = ridfill_nan(temp);

ivarname2 = 'T_20 Temperature';

tt = zeros(1,5);
tmean = zeros(1,5);
tstd = zeros(1,5);
tmax = zeros(1,5);
tmin = zeros(1,5);

fprintf('Variable: %12s\n\n',ivarname2)
fprintf('Month\t\tPoints\t\t  Mean\t\tStandard\tMaximum\tMinimum\n')
fprintf('\t\t\t\t\t\tDeviation\n')

if data_type == 'h';
    fprintf('\t\t(Hours)\t\t(Deg. C)\t(Deg. C)\t(Deg. C)\t(Deg. C)\n\n')
elseif data_type == 'l';
    fprintf('\t\t(Points)\t\t(Deg. C)\t(Deg. C)\t(Deg. C)\t(Deg. C)\n\n')
end

data2 = temp;

for j = 1:5;
    thismo2 = find(M2==monthvector(j));
    if isempty(thismo2) | length(thismo2)==1,
        fprintf('%4.0f\t\t%8s\n',monthvector(j),' No Data')
        dmean2 = NaN;
        dstd2 = NaN;
        dmax2 = NaN;
        dmin2 = NaN;
        jt2 = NaN;
    else
    %calculate and prep statistics for output
    dlen2 = length(find(isnan(data2(thismo2))==0));   
    dmax2 = gmax(data2(thismo2));   
    dmin2 = gmin(data2(thismo2));
    dmean2 = gmean(data2(thismo2));			
    dstd2 = gstd(data2(thismo2));				
    dmeanday2 = mean(D2(thismo2));
    jt2 = monthvector(j);
    %print out results
   	fprintf('%4.0f\t\t%6.0f\t\t%9.1f\t\t%5.1f\t\t%6.1f\t\t%6.1f',...
    jt2,dlen2,dmean2,dstd2,dmax2,dmin2)
    end
    %patch results into new variables
    tt(j) = jt2;
    tmean(j) = dmean2;
    tstd(j) = dstd2;
    tmax(j) = dmax2;
    tmin(j) = dmin2;
    fprintf('\n')
end     

%Take out data not from the specified 5 months;
%creepy_data = find(M2 > 4 & M2 < 12);
%
%if isempty(creepy_data) == 0;
%data2(creepy_data) = NaN;
%end

%Calculate yearly statistics for temperature.
dlen2 = length(find(isnan(data2)==0));   
dmax2 = gmax(data2);   
dmin2 = gmin(data2);   
dmean2 = gmean(data2);
dstd2 = gstd(data2);   
fprintf('  All\t\t%6.0f\t\t%9.1f\t\t%5.1f\t\t%6.1f\t\t%6.1f\n',...
    dlen2,dmean2,dstd2,dmax2,dmin2)
fprintf('\n\n\n')
                  
tyt = 2.5;
tymean = dmean2;
tystd = dstd2;
tymax = dmax2;
tymin = dmin2;

PlotStatsNY