function a_pf = param_Rs_cap_leak_int_t(param_init_vals, id, props) 
  
% param_Rs_cap_leak_int_t - Membrane capacitance and leak integrated over time.
%
% Usage:
%   a_pf = param_Rs_cap_leak_int_t(param_init_vals, id, props)
%
% Parameters:
%   param_init_vals: Array or structure with initial values for 
%     series resistance, Rs [MOhm]; leak conductance, gL [uS]; leak
%     reversal, EL [mV]; cell capacitance, Cm [nF], and a delay [ms]. 
%   id: An identifying string for this function.
%   props: A structure with any optional properties.
% 	   (Rest passed to param_func)
%		
% Returns a structure object with the following fields:
%	param_mult: Holds the inf and tau functions.
%
% Description:
%   Defines a function f(a_pf, struct('v', V [mV], 'dt', dt [ms])) where v is an array of voltage
% values [mV] changing with dt time steps [ms]. The (in)activation is then
% is the result of the integration of dm = (inf-m)/tau for each value of
% v. Initial value is taken from the first voltage value.
%
% See also: param_mult, param_func, param_act, tests_db, plot_abstract
%
% Example:
% >> f_capleak = ...
%    param_Rs_cap_leak_int_t(struct('Rs', 100, 'gL', 3, 'EL', -80, 'Cm', 1e-2), ...
%                        ['Ca chan 3rd instar cap leak']);
%
% $Id: param_Rs_cap_leak_int_t.m 131 2010-06-12 04:02:36Z cengiz $
%
% Author: Cengiz Gunay <cgunay@emory.edu>, 2010/03/02
  
% TODO:
  
  if ~ exist('props', 'var')
    props = struct;
  end

  if ~ exist('id', 'var')
    id = '';
  end

  param_names_ordered = {'Rs', 'gL', 'EL', 'Cm', 'delay'};
  if ~ isstruct(param_init_vals)
    param_init_vals = ...
        cell2struct(num2cell(param_init_vals(:)'), ...
                    param_names_ordered, 2);
  else
    % make sure they're ordered consistently to match the range description
    param_init_vals = ...
        orderfields(param_init_vals, param_names_ordered);
  end
  
  % physiologic parameter ranges
  param_ranges = ...
      [ eps eps -100 eps 0;...
        1e3 1e3 -50 1e3  10];
  
  a_pf = ...
      param_func(...
        {'time [ms]', 'I_{cap+leak} [nA]'}, ...
        param_init_vals, [], ...
        @cap_leak_int, id, ...
        mergeStructs(props, struct('paramRanges', param_ranges)));
  
  function Im = cap_leak_int(p, x)
    Vc = x.v;
    dt = x.dt;
    
    % do the delay as float and interpolate Vc so that the fitting
    % algorithm can move it around
    delay_dt = p.delay/dt;
    delay_dt_int = floor(delay_dt);
    delay_dt_frac = delay_dt - delay_dt_int;

    % prefix some data to reach steady-state
    fixed_delay = round(2/dt);

    % make a new vector for delayed voltage
    Vc_delay = ...
        [ repmat(Vc(1, :), fixed_delay + delay_dt_int + 1, 1); ...
          Vc(2:(end-delay_dt_int), :) - ...
          delay_dt_frac * ...
          diff(Vc(1:(end-delay_dt_int), :)) ];

    [t_tmp, Vm] = ...
        ode15s(@(t, Vm) ...
               (- (Vm - p.EL) * p.gL + ...
                (Vc_delay(round(t/dt) + 1, :)' - Vm) / p.Rs ) / p.Cm, ...
        (0:(size(Vc_delay, 1) - 1))*dt, Vc(1, :));
    
    % after solving Vm, return total input current (except manual
    % correction for leak at -70 mV)
    Im = (Vc_delay - Vm) / p.Rs - ...
         ones(size(Vc_delay, 1), 1) * (Vc_delay(fixed_delay, :) - Vm(fixed_delay, :)) / p.Rs;
    %ones(size(Vc, 1), 1) * (Vm(fixed_delay, :) - p.EL) * p.gL;
    
    % crop the fixed_delay part
    Im = Im((fixed_delay + 1):end, :);
  end

end

