function out = running_corr_vals(UTC_time, field_val1, field_val2, calc_types, bin_width_min)

% This function was written to calculate the running value(s) of various
% pairs of data (e.g. lines & slopes)
% Inputs: 1) UTC_time should be a vector, 2) field_val1 and field_val2 
% should be a vector that are the same size as UTC_time
% 3) calc_types is a list of % the different calculations that can be run 
% for a correlated running value. This input can be a single string, or a 
% cell. Calculation types include:
% linem - calculate the running slope of a linear fit line
% lineb - calculate the running intercept of a linear fit line
% r - linear correlation coefficient
% 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_vals1 = rand(length(UTC_time),1);
% field_vals2 = rand(length(UTC_time),1);
% calc_types = {'linem', 'lineb', 'r'};
% bin_width_min = 10;

% UTC_time=atf_filter.UTC_time; 
% field_val1=atf_filter.AM;
% field_val2=atf_filter.ln_avg_intens(:,1);
% calc_types= {'tsrem', 'tsreb'}; %, 'lineb', 'r'};
% bin_width_min=set_raw.deltime(2);

%% data check
if min(size(field_val1))~=1 || min(size(field_val2))~=1 || ...
     length(field_val1)~=length(UTC_time) || length(field_val2)~=length(UTC_time);
  error('Variables field_vals1 and field_vals2 should be vectors the same length as UTC_time')
end
n_tpoints = length(UTC_time);

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 = {'linem', 'lineb', 'r', 'tsrem', 'tsreb'}; %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('linem', calc_types) || ismember('lineb', calc_types);
    calc_types2(end+1:end+2) = {'linem', 'lineb'};
end
if ismember('tsreb', calc_types) && ~ismember('tsrem', calc_types);
    calc_types2(end+1) = {'tsrem'}; % add in Theil-Sen robust estimator slope
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_val1)); 
end
    
%% prep for TheilSen estimator
if ismember('tsrem', calc_types2) || ismember('tsreb', calc_types2);
    tsre_mat = sparse(length(field_val1),length(field_val1)); %pre-allocate sparse matrix. Only about 2% of size of filled matrix
        
    for i1=1:size(field_val1); %prepare the slopes matrix 
        % accumulate slopes
        timeind_i1 = false([1,length(field_val1)]); %pre-allocate/form index
        timeind_i1(i1:end) = (UTC_time(i1:end)-UTC_time(i1))*24*60<(bin_width_min/2); %index points within time index
        tsre_mat(i1,timeind_i1) = (field_val2(i1)-field_val2(timeind_i1))./(field_val1(i1) - field_val1(timeind_i1));
    end
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_vals1 = field_val1(i1_indx, :); %pulls values close to current time point
  i1_vals2 = field_val2(i1_indx, :); %pulls values close to current time point
  
  if running_npoints(i1)>=2; %require at least 2 points
    for i2=1:length(calc_types2);
      if isequal(calc_types2{i2}, 'r');
          pre_out.(calc_types2{i2})(i1) = nancorrcoef(i1_vals1, i1_vals2); %calculates correlation coefficient
      elseif isequal(calc_types2{i2}, 'linem');

          linefit = nanpolyfit(i1_vals1, i1_vals2,1);

          pre_out.linem(i1) = linefit(1);
          pre_out.lineb(i1) = linefit(2);
      elseif isequal(calc_types2{i2}, 'lineb');
          %do nothing for now
      elseif isequal(calc_types2{i2}, 'tsrem'); % Theil-Sen Robust slope estimator
          pre_out.tsrem(i1) = nanmedian(nonzeros([tsre_mat(i1,i1_indx), tsre_mat(i1_indx,i1)']));

          if ismember('tsreb', calc_types2); %also include intercept
              pre_out.tsreb(i1) = nanmedian(i1_vals2-pre_out.tsrem(i1)*i1_vals1);
          end
      end

    end
  end
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



