% Rideout_PressureFlowNP is a CV loop nonpulsatile model ported from Rideout
% (ACSL program PF-NP, Physiome JSim model 360)
function Rideout_PressureFlowNP()
  close all

  % Parameters
  Mode =    1;            % see allowed values in function ModeCB below

  G1    = 24 / 1332;      % [cm^4*sec/gr] L. Vent. Preload Conductance
  G2IC  =   .821 / 1332;  % [cm^4*sec/gr] L. Vent. Afterload Init. Conductance
  G2NEW =   .4105 / 1332; % [cm^4*sec/gr] L. Vent. Afterload New Conductance
  G3    = 40 / 1332;      % [cm^4*sec/gr] R. Vent. Preload Conductance
  G4    =  3.889 / 1332;  % [cm^4*sec/gr] R. Vent. Afterload Conductance

  RS =  1.0111 * 1332;    % [gr/cm^4/sec] Systemic Arterial Resistance
  RP =   .12222 * 1332;   % [mmHg*sec/ml] Pulmonary Arterial Resistance

  CS =   2.6316 / 1332;   % [cm^4*sec^2/gr] Systemic Arterial Compliance
  CR = 225 / 1332;        % [cm^4*sec^2/gr] Systemic Venous Compliance
  CP =   6.9444 / 1332;   % [cm^4*sec^2/gr] Pulmonary Arterial Compliance
  CL =  42.857 / 1332;    % [cm^4*sec^2/gr] Pulmonary Venous Compliance

  QSIC = 1000;            % [ml] Systemic Arterial Initial Volume
  QRIC = 5400;            % [ml] Systemic Venous Initial Volume
  QPIC =  500;            % [ml] Pulmonary Arterial Initial Volume
  QLIC = 1800;            % [ml] Pulmonary Venous Initial Volume

  QSU  =  750;            % [ml] Systemic Arterial Unstressed Volume
  QRU  = 4500;            % [ml] Systemic Venous Unstressed Volume
  QPU  =  375;            % [ml] Pulmonary Arterial Unstressed Volume
  QLU  = 1500;            % [ml] Pulmonary Venous Unstressed Volume

  A = 25;                 % [ml/sec] Infusion Flow Amplitude

  TSTT = 1;               % [sec] Infusion flow at t > TSTT
  WID = 2;                % [sec] Infusion flow pulse width
  TCH = 1;                % [sec] Stiffness change at t > TCH
  
  % Solver Parameters
  options = ddeset('RelTol', 1e-7, 'AbsTol', 1e-7);

  % GUI Setup
  hGUI = figure('Visible', 'off', 'Position', [360 500 180 80]);
  movegui(hGUI, 'northwest')

  hModeText = uicontrol('Style', 'text', 'String', 'Mode:', ...
                        'BackgroundColor', get(hGUI, 'Color'), ...
                        'HorizontalAlignment', 'Left', ...
                        'Position', [0 68 150 12]);
  hMode = uicontrol('Style', 'popupmenu', 'String', ...
                    {'Infusion', 'Left Ventricular Stiffness Change'}, ...
                    'Position', [0 43 150 25], 'Callback', {@ModeCB});

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

  set([hGUI hModeText hMode hRun], 'Units', 'normalized')
  set(hGUI, 'NumberTitle', 'off', 'Name', 'Rideout_PressureFlowNP')
  set(hGUI, 'MenuBar','none', 'HandleVisibility', 'off', 'Visible', 'on')

  % Mode popup menu callback function
  function ModeCB(hObj, event) 
    str = get(hObj, 'String');
    val = get(hObj, 'Value');
    switch str{val};
    case 'Infusion'
      Mode = 1;
    case 'Left Ventricular Stiffness Change'
      Mode = 2;
    end
  end % ModeCB

  % Run button callback function
  function RunCB(hObj, event)
    close all
    
    % Independent Variable
    dt = 0.1;
    tmax = 12 * (Mode == 1) + 9 * (Mode == 2);
    tspan = (0:dt:tmax)';

    % Dependent Variables
    FI = A * (tspan >= TSTT & tspan < (TSTT + WID)) .* (Mode == 1);
    
    [tout Q] = ode23(@myode, tspan, [QSIC QRIC QPIC QLIC], options);
    [QT P F] = myode2(tspan, Q);

    if Mode == 1
      figure('name', 'INFUSION', 'NumberTitle', 'off');
      subplot(2, 1, 1)
      plot(tout, FI, '-r');
      title('Infusion Flow Pulse')
      xlabel('Time [sec]')
      ylabel('Flow [ml/sec]')

      subplot(2, 1, 2)
      plot(tout, QT / 1000, '-g');
      title('Total Volume')
      xlabel('Time [sec]')
      ylabel('Volume [liter]')
      axis([0 tmax 8.64 8.77])

      figure('name', 'VOLUME', 'NumberTitle', 'off');
      subplot(2, 2, 1)
      plot(tout, Q(:, 1) / 1000, '-k');
      title('Systemic Arterial Volume')
      xlabel('Time [sec]')
      ylabel('Volume [liter]')
      axis([0 tmax 1 1.012])

      subplot(2, 2, 2)
      plot(tout, Q(:, 2) / 1000, '-r');
      title('Systemic Venous Volume')
      xlabel('Time [sec]')
      ylabel('Volume [liter]')

      subplot(2, 2, 3)
      plot(tout, Q(:, 3) / 1000, '-g');
      title('Pulmonary Arterial Volume')
      xlabel('Time [sec]')
      ylabel('Volume [liter]')
      axis([0 tmax .5 .512])

      subplot(2, 2, 4)
      plot(tout, Q(:, 4) / 1000, '-b');
      title('Pulmonary Venous Volume')
      xlabel('Time [sec]')
      ylabel('Volume [liter]')
      axis([0 tmax 1.8 1.812])
    end % if Mode == 1
    
    figure('name', 'PRESSURE', 'NumberTitle', 'off');
    subplot(2, 2, 1)
    plot(tout, P(:, 1) / 1332, '-k');
    title('Systemic Arterial Pressure')
    xlabel('Time [sec]')
    ylabel('Pressure [mmHg]')
    
    subplot(2, 2, 2)
    plot(tout, P(:, 2) / 1332, '-r');
    title('Systemic Venous Pressure')
    xlabel('Time [sec]')
    ylabel('Pressure [mmHg]')

    subplot(2, 2, 3)
    plot(tout, P(:, 3) / 1332, '-g');
    title('Pulmonary Arterial Pressure')
    xlabel('Time [sec]')
    ylabel('Pressure [mmHg]')

    subplot(2, 2, 4)
    plot(tout, P(:, 4) / 1332, '-b');
    title('Pulmonary Venous Pressure')
    xlabel('Time [sec]')
    ylabel('Pressure [mmHg]')

    figure('name', 'FLOW', 'NumberTitle', 'off');
    plot(tout, F(:, 2), '-k', tout, F(:, 3), '-r', tout, F(:, 4), '-g', tout, F(:, 1), '-b');
    title('Systemic Arterial and Venous, Pulmonary Arterial and Venous Flows')
    xlabel('Time [sec]')
    ylabel('Flow [ml/sec]')
    legend('FR', 'FP', 'FL', 'FS')
    
    % Solve ODEs
    function myodeout = myode(t, Q)
      G2 = G2NEW * (t >= TCH & Mode == 2) + G2IC * (t < TCH | Mode == 1);
      
      PS = (Q(1) - QSU) / CS;
      PR = (Q(2) - QRU) / CR;
      PP = (Q(3) - QPU) / CP;
      PL = (Q(4) - QLU) / CL;

      FS = (PS - PR) / RS;
      FR = G3 * PR - G4 * PP;
      FP = (PP - PL) / RP;
      FL = G1 * PL - G2 * PS;

      % Differential Equations
      myodeout = [
        FL - FS                         % QS:t
        FS - FR + FI(round(t / dt) + 1) % QR:t
        FR - FP                         % QP:t
        FP - FL                         % QL:t
      ];
    end % myode
    
    % Calculate rest of dependent variables after solving ODEs
    function [QT P F] = myode2(t, Q)
      G2 = G2NEW * (t >= TCH & Mode == 2) + G2IC * (t < TCH | Mode == 1);
      QT = Q(:, 1) + Q(:, 2) + Q(:, 3) + Q(:, 4);
      P(:, 1) = (Q(:, 1) - QSU) / CS;
      P(:, 2) = (Q(:, 2) - QRU) / CR;
      P(:, 3) = (Q(:, 3) - QPU) / CP;
      P(:, 4) = (Q(:, 4) - QLU) / CL;
      F(:, 1) = (P(:, 1) - P(:, 2)) / RS;
      F(:, 2) = G3 * P(:, 2) - G4 * P(:, 3);
      F(:, 3) = (P(:, 3) - P(:, 4)) / RP;
      F(:, 4) = G1 * P(:, 4) - G2 .* P(:, 1);
    end % myode2
  end % RunCB
end % Rideout_PressureFlowNP