function out = running_vals(UTC_time, field_vals, calc_types, bin_width_min)

% This function was written to calculate the running value(s) of various
% fields.
% Inputs: 1) UTC_time should be a vector, 2) field_vals should be a vector
% or a matrix with one side equal to the UTC_time length, these are values
% that the running product is calculated over, 3) calc_types is a list of
% the different calculations that can be run for a running value. This 
% input can be a single string, or a cell. Calculation types include:
% min - calculate the running minimum
% max - calculate the running maximum
% range - calculate the running range
% mean - calculate the running mean
% median - calculate the running median
% std - calculate the running std
% 4) bin_width_min is the width of the running bin in minutes. Note: it is
% assumed that there are no more than 11 values per minute.


%% example inputs
% UTC_time = datenum(2015,1,1,10,0,0):(1/(60*24)):datenum(2015,1,1,15,0,0);
% field_vals = rand(length(UTC_time),2);
% calc_types = {'min', 'max', 'range', 'mean', 'median', 'std'};
% bin_width_min = 10;
% UTC_time = filtered.utc_time;
% field_vals = filtered.airmass;
% calc_types = {'min', 'max', 'range', 'mean', 'median', 'std'};
% bin_width_min = 10;


%% data check
if size(field_vals,1)==length(UTC_time); 
    %do nothing
elseif size(field_vals,2)==length(UTC_time); 
    field_vals = field_vals'; %transpose matrix
else
  error('UTC_time and spec_fields appear to have different lengths')
end
[n_tpoints, ~] = size(field_vals);

if iscell(calc_types);
    % do nothing
elseif ischar(calc_types);
    calc_types = {calc_types};
else
    error('Calc_types needs to be a cell or string of the possible calculation types')
end

internal_ct_list = {'min', 'max', 'range', 'mean', 'median', 'std', 'sqdf'}; %internal list of calculation types
for i1=1:length(calc_types);
    if isempty(calc_types{i1}) || ~ismember(calc_types{i1}, internal_ct_list);
        warning([calc_types{i1}, ' is not a recognized calculation type...skipping.'])
        calc_types{i1} = [];
    end
end
calc_types = calc_types(~cellfun(@isempty,calc_types)); %remove unrecognized calculation types

%% sometimes more outputs are calculated than what is called for. This subsection adds these in
calc_types2 = calc_types; %values calculated internally (includes ancillary types)
if ismember('range', calc_types);
    calc_types2(end+1:end+2) = {'min', 'max'};
end
calc_types2 = unique(calc_types2); %remove duplicates 


%% 
t_spacing = abs(diff(UTC_time))*24*60;
t_spacing = t_spacing(t_spacing>0.0001); 
min_t_space = min(t_spacing); %minimum time spacing in minutes - this is likely a slight underestimate in I2S data

self_index = findclosestdaily_n(UTC_time, UTC_time, floor(bin_width_min/min_t_space)); %index closest neighboring data points
running_npoints = NaN*zeros(n_tpoints, 1); %pre-allocate, number of points in calculation

%% pre-allocate some values
for i1=1:length(calc_types2);
   pre_out.(calc_types2{i1}) = NaN*zeros(size(field_vals)); 
end
    

%%
for i1=1:n_tpoints; %over each timestep
  i1_indx_indx = abs(self_index.closest_indx_diff(i1,:))<(bin_width_min/2)/(60*24); %within specified number of minutes of self
  
  running_npoints(i1) = sum(i1_indx_indx);
  
  i1_indx = self_index.closest_indx(i1, i1_indx_indx); %index from main dataset
  i1_vals = field_vals(i1_indx, :); %pulls values close to current time point
  
  for i2=1:length(calc_types2);
    if isequal(calc_types2{i2}, 'min');
        pre_out.(calc_types2{i2})(i1,:) = nanmin(i1_vals);
    elseif isequal(calc_types2{i2}, 'max');
        pre_out.(calc_types2{i2})(i1,:) = nanmax(i1_vals);
    elseif isequal(calc_types2{i2}, 'mean');
        pre_out.(calc_types2{i2})(i1,:) = nanmean(i1_vals);
    elseif isequal(calc_types2{i2}, 'median');
        pre_out.(calc_types2{i2})(i1,:) = nanmedian(i1_vals);
    elseif isequal(calc_types2{i2}, 'std') && running_npoints(i1)>1;
        pre_out.(calc_types2{i2})(i1,:) = nanstd(i1_vals);
    elseif isequal(calc_types2{i2}, 'sqdf') && running_npoints(i1)>1; %rms of differences from center point
        pre_out.(calc_types2{i2})(i1,:) = sqrt(running_npoints(i1)*sum((i1_vals-i1_vals(1)).^2));    
    elseif isequal(calc_types2{i2}, 'range');
        %do nothing for now
    end
    
  end
end
if ismember('range', calc_types);
    pre_out.range = pre_out.max - pre_out.min;
end

out.n = running_npoints;
out.UTC_time = UTC_time;
for i2=1:length(calc_types2);
    out.(calc_types2{i2}) = pre_out.(calc_types2{i2});
end

%% plot check
% figure(1); clf; hold on;
% plot(out.UTC_time, field_vals(:,1), '.')
% plot(out.UTC_time, out.range(:,1), '.-')



