% Rideout_IndicatorDilution is a compartmental cardiovascular indicator dilution
% model ported from Rideout (ACSL programs IND-DIL, Physiome JSim model 334)
function Rideout_IndicatorDilution()
  close all
  
  % Global Parameters
  Mode = 1;   % see allowed values in function ModeCB below
  F0 = 100;   % [ml/sec] Right Ventricular flow (reference flow)
  Fshunt = 0; % [ml/sec] Shunt flow (diseased heart)
  Tp = 5;     % [sec] Pulmonary Delay
  Ts = 8;     % [sec] Systemic Delay

  % Independent Variable
  t = 0:0.2:100;

  % Solver Parameters
  options = ddeset('RelTol', 1e-7, 'AbsTol', 1e-7);
  
  hGUI = figure('Visible', 'off', 'Position', [360 500 180 200]);
  movegui(hGUI, 'northwest')

  hModeText = uicontrol('Style', 'text', 'String', 'Mode:', ...
                        'BackgroundColor', get(hGUI, 'Color'), ...
                        'HorizontalAlignment', 'Left', ...
                        'Position', [0 188 150 12]);
  hMode = uicontrol('Style', 'popupmenu', 'String', ...
                    {'Healthy', 'Ventricular Septal Defect', ...
                     'Atrial Septal Defect', 'Patent Ductus Arteriosus', ...
                     'Open Loop'}, ...
                    'Position',[0 163 150 25], 'Callback',{@ModeCB});

  hRef = uicontrol('Style', 'text', 'String', ...
                   'Reference Flow: F0 = 100 [ml/sec]', ...
                   'BackgroundColor', get(hGUI, 'Color'), ...
                   'HorizontalAlignment', 'Left', ...
                   'Position', [0 148 180 12]);
  hShuntText = uicontrol('Style', 'text', 'String', ...
                         'Shunt Flow: Fshunt = 0 [ml/sec]', ...
                         'BackgroundColor', get(hGUI, 'Color'), ...
                         'HorizontalAlignment', 'Left', ...
                         'Position', [0 130 150 12]);
  hShunt = uicontrol('Style', 'slider', 'Min', 0, 'Max', F0 * 0.99, ...
                     'Value', 0, ...
                     'Position', [0 105 150 25], 'Callback', {@ShuntCB});

  hRun = uicontrol('Style', 'pushbutton', 'String', 'Run', ...
                   'Position', [0 60 70 25], 'Callback', {@RunCB});

  hLoopText = uicontrol('Style', 'text', 'String', 'VSD Mode only:', ...
                        'BackgroundColor', get(hGUI, 'Color'), ...
                        'HorizontalAlignment', 'Left', ...
                        'Position', [0 35 150 12]);
  hLoop = uicontrol('Style', 'pushbutton', 'String', 'Loop', ...
                    'Position', [0 10 70 25], 'Callback', {@LoopCB});

  set([hGUI hModeText hMode hRef hShuntText hShunt hRun hLoopText hLoop], ...
      'Units', 'normalized')
  set(hGUI, 'NumberTitle', 'off', 'Name', 'Rideout_IndicatorDilution')
  set(hGUI, 'MenuBar','none', 'HandleVisibility', 'off', 'Visible', 'on')

  function ModeCB(hObj, event) 
    str = get(hObj, 'String');
    val = get(hObj, 'Value');
    switch str{val};
    case 'Healthy'
      Mode = 1;
    case 'Ventricular Septal Defect'
      Mode = 2;
    case 'Atrial Septal Defect'
      Mode = 3;
    case 'Patent Ductus Arteriosus'
      Mode = 4;
    case 'Open Loop'
      Mode = 5;
    end
  end % ModeCB

  function ShuntCB(hObj, event)
    Fshunt = get(hObj, 'Value');
    set(hShuntText, 'String', ...
        ['Shunt Flow: Fshunt = ' int2str(Fshunt), ' [ml/sec]'])
  end % ShuntCB

  function RunCB(hObj, event) 
    close all

    if Mode ~= 5
      sol =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);

      figure('name','CONC', 'NumberTitle', 'off');
      subplot(2, 1, 1)
      plot(sol.x, sol.y(2, :), '-k', sol.x, sol.y(4, :), '-r', ...
           sol.x, sol.y(6, :), '-g');
      title('Indicator Concentrations C1, C3, C5')
      xlabel('Time [sec]')
      ylabel('Concentration [mM]')
      legend('C1', 'C3', 'C5')

      subplot(2, 1, 2)
      plot(sol.x, sol.y(6, :), '-k', sol.x, sol.y(7, :), '-r', ...
           sol.x, sol.y(1, :), '-g');
      title('Indicator Concentrations C5, C6, C0')
      xlabel('Time [sec]')
      ylabel('Concentration [mM]')
      legend('C5', 'C6', 'C0')
    else
      sol =  dde23(@dCdt, [Tp Ts], zeros(12, 1), t, options);

      figure('name','OPEN_LOOP', 'NumberTitle', 'off');
      subplot(2, 1, 1)
      plot(sol.x, sol.y(5, :), '-k', sol.x, sol.y(1, :), '-r');
      title('Indicator Concentrations C4 and C0')
      xlabel('Time [sec]')
      ylabel('Concentration [mM]')
      legend('C4', 'C0')

      subplot(2, 1, 2)
      plot(sol.x, sol.y(9, :), '-k', sol.x, sol.y(10, :), '-r');
      title('Time Integrals of C4(t) and C0(t)')
      xlabel('Time [sec]')
      ylabel('Concentration Integral [mM*sec]')
      legend('INT4', 'INT0')

      figure('name','TRANSPORT', 'NumberTitle', 'off');
      subplot(2, 1, 1)
      plot(sol.x, sol.y(5, :) / sol.y(9, end), '-k', ...
           sol.x, sol.y(1, :) / sol.y(9, end), '-r');
      title('Transport Functions H4 and H0')
      xlabel('Time [sec]')
      ylabel('[1/sec]')
      legend('H4', 'H0')

      subplot(2, 1, 2)
      plot(sol.x, sol.y(11, :) / sol.y(9, end), '-k', ...
           sol.x, sol.y(12, :) / sol.y(9, end), '-r');
      title('Average Transit Time TA4 and TA0')
      xlabel('Time [sec]')
      ylabel('[sec]')
      legend('TA4', 'TA0')
    end % if Mode ~= 5
  end % RunCB

  function LoopCB(hObj, event) 
    if Mode == 2
      close all
      figure('name', 'VSD_SWEEP', 'NumberTitle', 'off');
      subplot(2, 1, 1)
      hold(gca, 'on')
      subplot(2, 1, 2)
      hold(gca, 'on')
      Fshunt = 0;
      sol =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
      Fshunt = 10;
      sol10 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
      Fshunt = 20;
      sol20 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
      Fshunt = 30;
      sol30 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
      Fshunt = get(hShunt, 'Value');
      
      subplot(2, 1, 1);
      plot(sol.x, sol.y(1, :), '-k', sol10.x, sol10.y(1, :), '--k', ...
           sol20.x, sol20.y(1, :), '-k', sol30.x, sol30.y(1, :), '--k');
      plot(sol.x, sol.y(8, :), '-r', sol10.x, sol10.y(8, :), '--r', ...
           sol20.x, sol20.y(8, :), '-r', sol30.x, sol30.y(8, :), '--r');
      title('Indicator Concentrations C0, C7 with VSD')
      xlabel('Time [sec]')
      ylabel('Concentration [mM]')
      legend('C0 (Fshunt = 0)', 'C0 (Fshunt = 10)', ...
             'C0 (Fshunt = 20)', 'C0 (Fshunt = 30)', ...
             'C7 (Fshunt = 0)', 'C7 (Fshunt = 10)', ...
             'C7 (Fshunt = 20)', 'C7 (Fshunt = 30)');

      subplot(2, 1, 2);
      plot(sol.x, sol.y(4, :), '-k', sol10.x, sol10.y(4, :), '--k', ...
           sol20.x, sol20.y(4, :), '-k', sol30.x, sol30.y(4, :), '--k');
      title('Indicator Concentration C3 with VSD')
      xlabel('Time [sec]')
      ylabel('Concentration [mM]')
      legend('C3 (Fshunt = 0)', 'C3 (Fshunt = 10)', ...
             'C3 (Fshunt = 20)', 'C3 (Fshunt = 30)');
      subplot(2, 1, 1)
      hold(gca, 'off')
      subplot(2, 1, 2)
      hold(gca, 'off')
    end % if Mode == 2
  end % LoopCB

  function Ret = dCdt(t, C, Cdelay)
    % Parameters
    V0 =  125; % [ml] Right Ventricle
    V1 =  250; % [ml] Pulmonary Artery
    V2 =  125; % [ml] Left Atrium
    V3 =  125; % [ml] Left Ventricle
    V4 =  750; % [ml] Aorta, Large Arteries
    V5 =  200; % [ml] Small Arteries
    V6 = 1000; % [ml] systemic Veins
    V7 =  125; % [ml] Right Atrium
    Fdye = (t <= 0.5) * 0.5; % [umol/sec] Injected dye flow

    if Fshunt < 0.99 * F0 % [ml/sec] Shunt Flow limited to 99% Fshunt
      Fshuntx = Fshunt;
    else
      Fshuntx = 0.99 * F0;
    end
    Fvsdx = Fshuntx * (Mode == 2); % [ml/sec] VSD Flow limited to 99% Fshunt
    Fasdx = Fshuntx * (Mode == 3); % [ml/sec] ASD Flow limited to 99% Fshunt
    Fpda = Fshunt * (Mode == 4); % [ml/sec] Patent Ductus Arteriosus (PDA) Flow

    Fp = F0 + Fpda;   % [ml/sec] Pulmonary Circulation Flow
    F2 = Fp - Fasdx;  % [ml/sec] Left Atrium Flow
    F3 = F2 - Fvsdx;  % [ml/sec] Left Ventricle Flow
    Fs = F3 - Fpda;   % [ml/sec] Systemic Circulation Flow
    F7 = Fs + Fasdx;  % [ml/sec] Right Atrium Flow

    Cp = Cdelay(2,1); % [mM] Pulmonary Circulation Concentration
    Cs = Cdelay(6,2); % [mM] Systemic Circulation Concentration

    % Differential Equations
    if Mode ~= 5
      Ret = [
        (F7 * C(8) + Fvsdx * C(4) - F0 * C(1)) / V0
        (Fdye - Fp * C(2) + Fpda * C(5) + F0 * C(1)) / V1
        (Fp * (Cp - C(3))) / V2
        (F2 * (C(3) - C(4))) / V3
        (F3 * (C(4) - C(5))) / V4
        (Fs * (C(5) - C(6))) / V5
        (Fs * (Cs - C(7))) / V6
        (Fs * C(7) + Fasdx * C(3) - F7 * C(8)) / V7
      ];
    else
      Ret = [
        (F7 * C(8) + Fvsdx * C(4) - F0 * C(1)) / V0
        (Fdye - Fp * C(2) + Fpda * C(5)) / V1
        (Fp * (Cp - C(3))) / V2
        (F2 * (C(3) - C(4))) / V3
        (F3 * (C(4) - C(5))) / V4
        (Fs * (C(5) - C(6))) / V5
        (Fs * (Cs - C(7))) / V6
        (Fs * C(7) + Fasdx * C(3) - F7 * C(8)) / V7
        C(5)
        C(1)
        C(5) * t
        C(1) * t
      ];
    end % if Mode ~= 5
  end % dCdt
end % Rideout_IndicatorDilution
