function [CO,flow] = est11_mf(abp, onset, MAP, HR, age, gender)
% EST11_mf  CO estimator 11: Wesseling's non-linear, time-varying 3-element model
%
%   In:  ABP       <kx1>     vector  --- abp waveform
%        ONSET     <(n+1)x1> vector  --- beat-to-beat onset time in samples
%        MAP       <nx1>     vector  --- beat-to-beat mean arterial pressure
%        HR        <nx1>     vector  --- beat-to-beat heart rate
%        AGE       <1x1>     scalar  --- age
%        GENDER    <1x1>     scalar  --- male=1, female=2
%
%   Out: CO        <nx1>     vector  --- estimated beat-by-beat CO
%        FLOW      <kx1>     vector  --- instantaneous flow (125 Hz)
%
%   Written by James Sun (xinsun@mit.edu) on Nov 19, 2005.
%   - v2.0 - 3/29/06 - uses state-space implementation, various bug fixes
%   - v2.1 - 5/10/06 - stroke volume is sum of positive flow only

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% constants
rho = 1.05;            % density of blood: 1.05 [g/cc]
l   = 80;              % aortic effective length [cm]
P1  = 57 - 0.44*age;   % pressure for arc-tangent relationship
switch gender
    case 1
        P0   = 76 - 0.89*age;   % male
        Amax = 5.62;
    case 2
        P0   = 72 - 0.89*age;   % female
        Amax = 4.12;
end

% initial conditions
R  = 100/50; % 2 [mmHg*s/cc]
x  = 0;      % initial condition for state variable

%% minor business
flow   = zeros(length(abp),1);
CO     = zeros(length(MAP),1);
abp    = double(abp);   % change to double FP precision
HR     = double(HR);
MAP    = double(MAP);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% main loop (runs Wesseling's model, updating nonlinear parameters for each iteration)
k=1;
for i=1:length(abp)
    
    %% first parameters
    P  = abp(i);        % ABP sample [mmHg]
    m  = (P-P0)/P1;
    A  = Amax*(0.5+atan(m)/pi); % aortic cross-sectional area [cm^2]
    Cp = Amax/(pi*P1*(1+m^2));  % aortic compliance per unit length [cm^2/mmHg]

    %% circuit parameters
    Z  = sqrt(rho/(1333*A*Cp)); % aortic characteristic impedance  [mmHg*s/cc]
    C  = l*Cp;                  % arterial compliance (windkessel) [cc/mmHg]

    %% state-space setup
    AA = -(1/R+1/Z)/C;
    BB = -1/(Z^2*C);
    CC = 1;
    DD = 1/Z;

    %% run simulation
    [y,x]   = SSsolve(AA,BB,CC,DD,P,[0 0.008],x(end));
    flow(i) = y(end);

    %% Update MAP, CO, R
    if k+1 > length(onset)
        break
    elseif i >= onset(k+1)
        pos_flow = flow(onset(k):end);
        pos_flow = pos_flow(pos_flow>0);
        SV    = sum(pos_flow)/125;
        CO(k) = SV*HR(k)/60;  % preserve CO in units of [cc/s]
        R     = MAP(k)/CO(k); % peripheral resistance [mmHg*s/cc]
        
        k = k+1; % point to next index for R,HR,onset,CO,etc.
    end
end

CO = 60/1000 * CO; % convert final CO into L/min