% Rideout_IndicatorDilution_Modified is a compartmental cardiovascular indicator
% dilution model based on Rideout_IndicatorDilution with systemic flow used as
% constant reference. (Physiome JSim model 335)
function Rideout_IndicatorDilution_Modified()
  close all
  
  % Global Parameters
  Fs = 100;   % [ml/sec] Right Ventricular flow (reference flow)
  Fvsd = 0;   % [ml/sec] Ventricular Septal Defect (VSD) Flow
  Fasd = 0;   % [ml/sec] Atrial Septal Defect (ASD) Flow
  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')

  hRef = uicontrol('Style', 'text', 'String', ...
                   'Reference Flow: Fs = 100 [ml/sec]', ...
                   'BackgroundColor', get(hGUI, 'Color'), ...
                   'HorizontalAlignment', 'Left', ...
                   'Position', [0 188 180 12]);

  hFvsdText = uicontrol('Style', 'text', 'String', ...
                        'VSD Flow: Fvsd = 0 [ml/sec]', ...
                        'BackgroundColor', get(hGUI, 'Color'), ...
                        'HorizontalAlignment', 'Left', ...
                        'Position', [0 168 150 12]);
  hFvsd = uicontrol('Style', 'slider', 'Min', 0, 'Max', Fs * 2, 'Value', 0, ...
                    'Position', [0 143 150 25], 'Callback', {@FvsdCB});

  hFasdText = uicontrol('Style', 'text', 'String', ...
                        'ASD Flow: Fasd = 0 [ml/sec]', ...
                        'BackgroundColor', get(hGUI, 'Color'), ...
                        'HorizontalAlignment', 'Left', ...
                        'Position', [0 125 150 12]);
  hFasd = uicontrol('Style', 'slider', 'Min', 0, 'Max', Fs * 2, 'Value', 0, ...
                    'Position', [0 100 150 25], 'Callback', {@FasdCB});

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

  hLoopVsd = uicontrol('Style', 'pushbutton', 'String', 'Loop VSD', ...
                       'Position', [0 40 70 25], 'Callback', {@LoopVsdCB});

  hLoopAsd = uicontrol('Style', 'pushbutton', 'String', 'Loop ASD', ...
                       'Position', [0 10 70 25], 'Callback', {@LoopAsdCB});

  set([hGUI hRef hFvsdText hFvsd hFasdText hFasd hRun hLoopVsd hLoopAsd], ...
      'Units', 'normalized')
  set(hGUI, 'NumberTitle', 'off', 'Name', 'Rideout_IndicatorDilution_Modified')
  set(hGUI, 'MenuBar','none', 'HandleVisibility', 'off', 'Visible', 'on')

  function FvsdCB(hObj, event)
    Fvsd = get(hObj, 'Value');
    set(hFvsdText, 'String', ...
        ['VSD Flow: Fvsd = ' int2str(Fvsd), ' [ml/sec]'])
  end % FvsdCB

  function FasdCB(hObj, event)
    Fasd = get(hObj, 'Value');
    set(hFasdText, 'String', ...
        ['ASD Flow: Fasd = ' int2str(Fasd), ' [ml/sec]'])
  end % FasdCB

  function RunCB(hObj, event) 
    close all

    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')
  end % RunCB

  function LoopVsdCB(hObj, event) 
    close all
    figure('name', 'SWEEP', 'NumberTitle', 'off');
    subplot(2, 1, 1)
    hold(gca, 'on')
    subplot(2, 1, 2)
    hold(gca, 'on')
    Fvsd = 0;
    sol =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fvsd = 10;
    sol10 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fvsd = 20;
    sol20 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fvsd = 30;
    sol30 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fvsd = get(hFvsd, '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 Shunt Aperture Sweep')
    xlabel('Time [sec]')
    ylabel('Concentration [mM]')
    legend('C0 (Fvsd = 0)', 'C0 (Fvsd = 10)', ...
           'C0 (Fvsd = 20)', 'C0 (Fvsd = 30)', ...
           'C7 (Fvsd = 0)', 'C7 (Fvsd = 10)', ...
           'C7 (Fvsd = 20)', 'C7 (Fvsd = 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 Shunt Aperture Sweep')
    xlabel('Time [sec]')
    ylabel('Concentration [mM]')
    legend('C3 (Fvsd = 0)', 'C3 (Fvsd = 10)', ...
           'C3 (Fvsd = 20)', 'C3 (Fvsd = 30)');
    subplot(2, 1, 1)
    hold(gca, 'off')
    subplot(2, 1, 2)
    hold(gca, 'off')
  end % LoopVsdCB

  function LoopAsdCB(hObj, event) 
    close all
    figure('name', 'SWEEP', 'NumberTitle', 'off');
    subplot(2, 1, 1)
    hold(gca, 'on')
    subplot(2, 1, 2)
    hold(gca, 'on')
    Fasd = 0;
    sol =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fasd = 10;
    sol10 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fasd = 20;
    sol20 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fasd = 30;
    sol30 =  dde23(@dCdt, [Tp Ts], zeros(8, 1), t, options);
    Fasd = get(hFasd, '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 Shunt Aperture Sweep')
    xlabel('Time [sec]')
    ylabel('Concentration [mM]')
    legend('C0 (Fasd = 0)', 'C0 (Fasd = 10)', ...
           'C0 (Fasd = 20)', 'C0 (Fasd = 30)', ...
           'C7 (Fasd = 0)', 'C7 (Fasd = 10)', ...
           'C7 (Fasd = 20)', 'C7 (Fasd = 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 Shunt Aperture Sweep')
    xlabel('Time [sec]')
    ylabel('Concentration [mM]')
    legend('C3 (Fasd = 0)', 'C3 (Fasd = 10)', ...
           'C3 (Fasd = 20)', 'C3 (Fasd = 30)');
    subplot(2, 1, 1)
    hold(gca, 'off')
    subplot(2, 1, 2)
    hold(gca, 'off')
  end % LoopAsdCB

  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

    F7 = Fs + Fasd;   % [ml/sec] Right Atrium Flow
    F0 = F7 + Fvsd;   % [ml/sec] Right Ventricle Flow
    F2 = Fs - Fasd;   % [ml/sec] Left Atrium Flow
    F3 = F2 - Fvsd;   % [ml/sec] Left Ventricle Flow

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

    % Differential Equations
    Ret = [
      (F7 * C(8) + Fvsd * C(4) - F0 * C(1)) / V0
      (Fdye - F0 * C(2) + F0 * C(1)) / V1
      (F0 * (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) + Fasd * C(3) - F7 * C(8)) / V7
    ];
  end % dCdt
end % Rideout_IndicatorDilution_Modified
