% Rideout_PressureFlow0 is an uncontrolled CV loop pressure-flow model ported
% from Rideout (ACSL program PF-0, Physiome JSim model 338)
function Rideout_PressureFlow0()
  close all
  
  % Parameters
  Mode = 1;       % see allowed values in function ModeCB below
  R40 = 250;      % [gr/cm^4/sec] VSD Resistance
  R0  =   10;     % [gr/cm^4/sec] Right Ventricle Resistance
  R1  =   10;     % [gr/cm^4/sec] Pulmonary Artery Resistance
  R2  =   26.6;   % [gr/cm^4/sec] Pulmonary Artery 2 Resistance
  R3  =  149.663; % [gr/cm^4/sec] Lung Capillaries Resistance
  R4  =   10;     % [gr/cm^4/sec] Left Ventricle Resistance
  R5  =   13.32;  % [gr/cm^4/sec] Aortic Resistance
  R6  =    6.66;  % [gr/cm^4/sec] Descending Aorta Resistance
  R7  =   26.64;  % [gr/cm^4/sec] Abdominal Aorta Resistance
  R8  = 5328;     % [gr/cm^4/sec] Leg Arteries Resistance
  R9  =   13.32;  % [gr/cm^4/sec] Veins Resistance
  R10 =   13.32;  % [gr/cm^4/sec] Vena Cava Resistance
  R11 = 1332;     % [gr/cm^4/sec] Upper Body Resistance
  R12 = 1332;     % [gr/cm^4/sec]
  R13 = 5328;     % [gr/cm^4/sec]

  L2 = 1.332;     % [gr/cm^4] Pulmonary Artery 2 Inertance
  L6 = 1.332;     % [gr/cm^4] Descending Aorta Inertance
  L7 = 1.332;     % [gr/cm^4] Abdominal Aorta Inertance

  C1  = .001001;  % [cm^4*sec^2/gr] Pulmonary Artery Compliance
  C2  = .003003;  % [cm^4*sec^2/gr] Pulmonary Artery 2 Compliance
  C3  = .005005;  % [cm^4*sec^2/gr] Lung Capillaries Compliance
  C5  = .000501;  % [cm^4*sec^2/gr] Aortic Compliance
  C6  = .001001;  % [cm^4*sec^2/gr] Descending Aorta Compliance
  C7  = .001336;  % [cm^4*sec^2/gr] Abdominal Aorta Compliance
  C8  = .008314;  % [cm^4*sec^2/gr] Leg Arteries Compliance
  C9  = .01001 ;  % [cm^4*sec^2/gr] Veins Compliance
  C10 = .005005;  % [cm^4*sec^2/gr] Vena Cava Compliance
  C11 = .002002;  % [cm^4*sec^2/gr] Upper Body Compliance

  PLIM = 2000 * 1332; % [gr/cm/sec^2] Valve Pressure limit
  
  LD =   43.956;  % [gr/cm^4/sec^2] Diastolic Left Ventricle Stiffness
  LS = 1998;      % [gr/cm^4/sec^2] Systolic Left Ventricle Stiffness
  RD =   58.608;  % [gr/cm^4/sec^2] Diastolic Right Ventricle Stiffness
  RS =  399.6;    % [gr/cm^4/sec^2] Systolic Right Ventricle Stiffness

  TH = 0.8;       % [sec] Heart period (75 beats per minute) 

  TCO   = 4;      % [sec] Time Constant for Cardiac Output
  TFP   = 4;      % [sec] Time Constant for Pulmonary Flow
  TFVSD = 4;      % [sec] Time Constant for VSD Flow
  
  % Independent Variable
  dt = 0.02;
  tspan = (0:dt:3.2)';

  % Solver Parameters
  options = ddeset('RelTol', 1e-7, 'AbsTol', 1e-7);
  
  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', {'Healthy', 'VSD'}, ...
                    '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_PressureFlow0')
  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 'Healthy'
      Mode = 1;
    case 'VSD'
      Mode = 2;
    end
  end % ModeCB

  % Run button callback function
  function RunCB(hObj, event) 
    close all

    % Dependent Variables
    PINIT = [14 14 6 104 104 104 9 8 7 60] * 1332 * (Mode == 1) + ...
            [19 18 6  79  79  79 9 9 8 47] * 1332 * (Mode ~= 1);
    FINIT = [19 17 8] * (Mode == 1) + [28 1 12] * (Mode ~= 1);
    QINIT = [154 158] * (Mode == 1) + [170 162] * (Mode ~= 1);
    COINIT = 103 * (Mode == 1) + 76 * (Mode ~= 1);
    FPINIT = 102 * (Mode == 1) + 152 * (Mode ~= 1);
    FVSDINIT = 70 * (Mode ~= 1);
    GD = (1 / R40) * (Mode ~= 1);

    SSW = sin(2 * pi * tspan / TH);
    ACT = min(max(SSW, 0), 1);
    S0 = RS * ACT + RD;
    S4 = LS * ACT + LD;

    [tout res] = ode23(@myode, tspan, ...
                       [PINIT FINIT QINIT COINIT FPINIT FVSDINIT], options);
    Q4 = res(:, 15);
    Q0 = res(:, 14);
    [P0 P4 F4 F5 F40] = myode2(res(:, 3), res(:, 4), Q0, Q4);
    
    figure('name', 'VOLUME_FLOW', 'NumberTitle', 'off');
    subplot(2, 2, 1)
    plot(tout, Q4, '-k');
    title('Left Ventricular Volume')
    xlabel('Time [sec]')
    ylabel('Volume [ml]')
    axis([0 3.5 30 180])
    
    subplot(2, 2, 2)
    plot(tout, F4, '-b', tout, F5, '-r');
    title('Left Ventricular and Aortic Flows')
    xlabel('Time [sec]')
    ylabel('Flow [ml/sec]')
    legend('F4', 'F5')
    axis([0 3.5 0 1600])

    subplot(2, 2, 3)
    plot(tout, Q0, '-r');
    title('Right Ventricular Volume')
    xlabel('Time [sec]')
    ylabel('Volume [ml]')

    subplot(2, 2, 4)
    plot(tout, F40, '-g');
    title('VSD Flow')
    xlabel('Time [sec]')
    ylabel('Flow [ml/sec]')

    figure('name', 'PRESSURE_LOCUS', 'NumberTitle', 'off');
    subplot(2, 2, 1)
    plot(tout, P0 / 1332, '-k', tout, res(:, 1) / 1332, '-r');
    title('Right Ventricular and Pulmonary Artery Pressures')
    xlabel('Time [sec]')
    ylabel('Pressure [mmHg]')
    legend('P0', 'P1')

    subplot(2, 2, 2)
    plot(tout, P4 / 1332, '-r', tout, res(:, 4) / 1332, '-b');
    title('Left Ventricular and Aortic Pressures')
    xlabel('Time [sec]')
    ylabel('Pressure [mmHg]')
    legend('P4', 'P5')
    axis([0 3.5 0 160])

    subplot(2, 2, 3)
    plot(Q0, P0 / 1332, '-g');
    title('Right Ventricle Locus')
    xlabel('Volume [ml]')
    ylabel('Pressure [mmHg]')

    subplot(2, 2, 4)
    plot(Q4, P4 / 1332, '-r');
    title('Left Ventricle Locus')
    xlabel('Volume [ml]')
    ylabel('Pressure [mmHg]')
    axis([30 165 0 150])

    % Solve ODEs
    function myodeout = myode(t, myodein)
      P    = [myodein(14) * S0(round(t / dt) + 1); myodein(1:3); ...
             myodein(15) * S4(round(t / dt) + 1); myodein(4:10)];
      F2   = myodein(11);
      F6   = myodein(12);
      F7   = myodein(13);
      CO   = myodein(16);
      FP   = myodein(17);
      FVSD = myodein(18);

      F0 = min(max(P(11) - P(1), 0), PLIM) / R0;
      F1 = min(max(P(1) - P(2), 0), PLIM) / R1;
      F3 = (P(3) - P(4)) / R3;
      F4 = min(max(P(4) - P(5), 0), PLIM) / R4;
      F5 = min(max(P(5) - P(6), 0), PLIM) / R5;
      F8 = (P(8) - P(9)) / R8;
      F9 = (P(9) - P(10)) / R9;
      F10 = (P(10) - P(11)) / R10;
      F11 = (P(7) - P(12)) / R11;
      F12 = (P(12) - P(10)) / R12;
      F13 = (P(8) - P(10)) / R13;
      F40 = GD * (P(5) - P(1));
      
      % Differential Equations
      myodeout = [
        (F1 - F2) / C1               % P1:t
        (F2 - F3) / C2               % P2:t
        (F3 - F4) / C3               % P3:t
        (F5 - F6) / C5               % P5:t
        (F6 - F11 - F7) / C6         % P6:t
        (F7 - F8 - F13) / C7         % P7:t
        (F8 - F9) / C8               % P8:t
        (F9 - F10 + F12 + F13) / C9  % P9:t
        (F10 - F0) / C10             % P10:t
        (F11 - F12) / C11            % P11:t
        (P(2) - P(3) - F2 * R2) / L2 % F2:t
        (P(6) - P(7) - F6 * R6) / L6 % F6:t
        (P(7) - P(8) - F7 * R7) / L7 % F7:t
        F0 - F1 + F40                % Q0:t
        F4 - F5 - F40                % Q4:t
        (F12 + F13 + F8 - CO) / TCO  % CO:t
        (F3 - FP) / TFP              % FP:t
        (F40 - FVSD) / TFVSD         % FVSD:t
      ];
    end % myode
    
    % Calculate rest of dependent variables after solving ODEs
    function [P0 P4 F4 F5 F40] = myode2(P3, P5, Q0, Q4)
      P0 = S0 .* Q0;
      P4 = S4 .* Q4;
      F4 = min(max(P3 - P4, 0), PLIM) / R4;
      F5 = min(max(P4 - P5, 0), PLIM) / R5;
      F40 = GD * (P4 - P0);
    end % myode2
  end % RunCB
end % Rideout_PressureFlow0