function out = running_I0poly(UTC_time, field_val1, field_val3, bin_width_min, init_ps)

% This function was written to calculate the running value(s) of various
% sets of data with 2 dependent variables, such as fits
% Inputs: 1) UTC_time should be a vector, 2) field_val1 (x), field_val2
% (y) - is taken as UTC_time, and field_val3 (z)
% should be vectors 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 (i.e. fittypes equations). This input can be a single string, or a 
% cell. Example calculation types include:
% ft1 = % fittype('c-x*(p0+p1*y)','options',fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[-4 0 0]));
% ft2 = % fittype('c-x*(p0+p1*y+p2*y^2)','options',fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[-4 0 0 0]));
% see matlab 'fit' page and related pages for more info

% 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=curr_day.UTC_time; 
% field_val1=curr_day.AM;
% field_val3=curr_day.ln_avg_intens_i1;
% bin_width_min=set_raw.deltime(2);
% init_ps = [-4.2, 0.2]; %initial p's (ln(I0) and tau)

%% fixed inputs
% for these fits, parameter c is the ln(I0), parameter p0 is tau_tot,
% parameter p1 is the delta(tau_tot)/min, and p2 is the change per minute
% (per minute)...etc.
% ft1 = fittype('c-x*(p0+p1*y)','independent',{'x','y'},'dependent',{'z'},'options', ...
%   fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[init_ps(1) init_ps(2) 0], 'Lower', [-6.5, 0, -0.1], 'Upper', [-0.5, 5, 0.1]));
ft2 = fittype('c-x*(p0+p1*y+p2*y^2)','independent',{'x','y'},'dependent',{'z'},'options', ...
  fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[init_ps(1) init_ps(2) 0 0], 'Lower', [-6.5, 0, -0.1, -0.01], 'Upper', [-0.5, 5, 0.1, 0.01]));
ft3 = fittype('c-x*(p0+p1*y+p2*y^2+p3*y^3)','independent',{'x','y'},'dependent',{'z'},'options', ...
  fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[init_ps(1) init_ps(2) 0 0 0], 'Lower', [-6.5, 0, -0.1, -0.01, -0.005], 'Upper', [-0.5, 5, 0.1, 0.01, 0.005]));
ft4 = fittype('c-x*(p0+p1*y+p2*y^2+p3*y^3+p4*y^4)','independent',{'x','y'},'dependent',{'z'},'options', ...
  fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[init_ps(1) init_ps(2) 0 0 0 0], 'Lower', ...
  [-6.5, 0, -0.1, -0.01, -0.005, -0.002], 'Upper', [-0.5, 5, 0.1, 0.01, 0.005, 0.002]));
ft5 = fittype('c-x*(p0+p1*y+p2*y^2+p3*y^3+p4*y^4+p5*y^5)','independent',{'x','y'},'dependent',{'z'},'options', ...
  fitoptions('Method','NonlinearLeastSquares', 'StartPoint',[init_ps(1) init_ps(2) 0 0 0 0 0], 'Lower', ...
  [-6.5, 0, -0.1, -0.01, -0.005, -0.001, -5E-4], 'Upper', [-0.5, 5, 0.1, 0.01, 0.005, 0.001, 5E-4]));
fit_types = {ft2, ft3, ft4, ft5}; %, 'lineb', 'r'};

field_val2=UTC_time;

%% data check
if min(size(field_val1))~=1 || min(size(field_val2))~=1 || min(size(field_val3))~=1 || ...
     length(field_val1)~=length(UTC_time) || length(field_val2)~=length(UTC_time) ...
     || length(field_val3)~=length(UTC_time);
  error('Variables field_vals(1,2,3) should be vectors the same length as UTC_time')
end
n_tpoints = length(UTC_time);

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

for i1=1:length(fit_types);
    if ~isequal(class(fit_types{i1}), 'fittype');
        warning(['type number ', num2str(i1), ' is not a recognized calculation type...quitting.'])
        fit_types{i1} = [];
    end
end
if any(cellfun(@isempty,fit_types));
  error('Use fittype, see matlab page on fittype (and fitoptions) for more info')
end
  

%% sometimes more outputs are calculated than what is called for. This subsection adds these in


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


%% pre-allocate some values
pre_out.n = NaN*zeros(n_tpoints, 1); %pre-allocate, number of points in calculation
cof_names = cell(length(fit_types),1);
for i2=1:length(fit_types);
   cof_names{i2} = coeffnames(fit_types{i2}); %save coefficient names
   pre_out.(['calc_' num2str(i2)]) = NaN*zeros(length(field_val1),length(cof_names{i2})); 
   pre_out.(['gof_' num2str(i2)]) = NaN*zeros(length(field_val1), 1); %goodness of fit indicators to save
end
 

%%
% tic
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
  
  pre_out.n(i1) = sum(i1_indx_indx); %number of points in running bin
  
  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
  i1_vals3 = field_val3(i1_indx, :); %pulls values close to current time point
  i1_vals2 = (i1_vals2-nanmean(i1_vals2))*1440; %convert to minutes and center the time
  
  if pre_out.n(i1)>10; %requiring at least 10 points for a fit
    for i2=1:length(fit_types);
      [coeffs_i2, fitind_i2] = fit([i1_vals1,i1_vals2],i1_vals3,fit_types{i2}); 

      for i3=1:length(cof_names{i2})
        pre_out.(['calc_' num2str(i2)])(i1,i3) = coeffs_i2.(cof_names{i2}{i3}); %save the fit coefficients
      end

      pre_out.(['gof_' num2str(i2)])(i1) = fitind_i2.adjrsquare; %save the adjusted R2

    end
  end
  
%   if mod(i1,100)==0; %display the time
%     toc
%   end
end

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



