% Rideout_PressureFlow1 is an uncontrolled CV loop pressure-flow model ported
% from Rideout (ACSL program PF-1, Physiome JSim model 343)
function Rideout_PressureFlow1()
  close all
  
  % Parameters
  Mode =    1;           % see allowed values in function ModeCB below
  
  RP1  =   10;           % [gr/cm^4/sec] Pulm. Art. 1 Resistance
  RP2  =   40;           % [gr/cm^4/sec] Pulm. Art. 2 Resistance
  RP3  =   80;           % [gr/cm^4/sec] Pulm. Art. 3 Resistance
  RL1  =   30;           % [gr/cm^4/sec] Pulm. Veins 1 Resistance
  RL2  =   10;           % [gr/cm^4/sec] Pulm. Veins 2 Resistance
  RLA  =    5;           % [gr/cm^4/sec] Left Atrium Resistance
  RLV  =    5;           % [gr/cm^4/sec] Left Ventricle Resistance
  RA1  =   10;           % [gr/cm^4/sec] Aortic 1 Resistance
  RA2  =  160;           % [gr/cm^4/sec] Aortic 2 Resistance
  RA3  = 1000;           % [gr/cm^4/sec] Syst. Art. Resistance
  RV1  =   90;           % [gr/cm^4/sec] Syst. Veins 1 Resistance
  RV2  =   10;           % [gr/cm^4/sec] Syst. Veins 2 Resistance
  RRA  =    5;           % [gr/cm^4/sec] Right Atrium Resistance
  RRV  =    5;           % [gr/cm^4/sec] Right Ventricle Resistance
  RPW1 =   10;           % [gr/cm^4/sec]
  RPW2 =   10;           % [gr/cm^4/sec]

  LP1  =    1;           % [gr/cm^4] Pulm. Art. 1 Inertance
  LL2  =    1;           % [gr/cm^4] Pulm. Veins 2 Inertance
  LLA  =    1;           % [gr/cm^4] Left Atrium Inertance
  LLV  =    1;           % [gr/cm^4] Left Ventricle Inertance
  LA1  =    1;           % [gr/cm^4] Aortic 1 Inertance
  LV2  =    1;           % [gr/cm^4] Syst. Veins 2 Inertance
  LRA  =    1;           % [gr/cm^4] Right Atrium Inertance
  LRV  =    1;           % [gr/cm^4] Right Ventricle Inertance

  CP1 = .0001;           % [cm^4*sec^2/gr] Pulm. Art. 1 Compliance
  CP2 = .0003;           % [cm^4*sec^2/gr] Pulm. Art. 2 Compliance
  CP3 = .0027;           % [cm^4*sec^2/gr] Pulm. Art. 3 Compliance
  CL1 = .001;            % [cm^4*sec^2/gr] Pulm. Veins 1 Compliance
  CL2 = .001;            % [cm^4*sec^2/gr] Pulm. Veins 2 Compliance
  CLA = .01176;          % [cm^4*sec^2/gr] Left Atrium Compliance
  CA1 = .00018;          % [cm^4*sec^2/gr] Aortic 1 Compliance
  CA2 = .00023;          % [cm^4*sec^2/gr] Aortic 2 Compliance
  CA3 = .00182;          % [cm^4*sec^2/gr] Syst. Art. Compliance
  CV1 = .021;            % [cm^4*sec^2/gr] Syst. Veins 1 Compliance
  CV2 = .045;            % [cm^4*sec^2/gr] Syst. Veins 2 Compliance
  CRA = .045;            % [cm^4*sec^2/gr] Right Atrium Compliance

  PP1EDM =  7.2 * 1332;  % [gr/cm/sec^2] Pulm. Art. 1 End-Diastolic Pressure
  PP2EDM =  7 * 1332;    % [gr/cm/sec^2] Pulm. Art. 1 End-Diastolic Pressure
  PP3EDM =  6.6 * 1332;  % [gr/cm/sec^2] Pulm. Art. 1 End-Diastolic Pressure
  PL1EDM =  4.45 * 1332; % [gr/cm/sec^2] Pulm. Veins 1 End-Diastolic Pressure
  PL2EDM =  3.62 * 1332; % [gr/cm/sec^2] Pulm. Veins 2 End-Diastolic Pressure
  PLAEDM =  3.45 * 1332; % [gr/cm/sec^2] Left Atrium End-Diastolic Pressure
  PLVEDM =  4 * 1332;    % [gr/cm/sec^2] Left Ventricle End-Diastolic Pressure
  PA1EDM = 64.3 * 1332;  % [gr/cm/sec^2] Aortic 1 End-Diastolic Pressure
  PA2EDM = 64 * 1332;    % [gr/cm/sec^2] Aortic 2 End-Diastolic Pressure
  PA3EDM = 63 * 1332;    % [gr/cm/sec^2] Syst. Art. End-Diastolic Pressure
  PV1EDM = 13.5 * 1332;  % [gr/cm/sec^2] Syst. Veins 1 End-Diastolic Pressure
  PV2EDM =  7.2 * 1332;  % [gr/cm/sec^2] Syst. Veins 2 End-Diastolic Pressure
  PRAEDM =  6.64 * 1332; % [gr/cm/sec^2] Right Atrium End-Diastolic Pressure
  PRVEDM =  7.4 * 1332;  % [gr/cm/sec^2] Right Ventricle End-Diastolic Pressure

  FP1IC =  0;            % [ml/sec] Pulm. Art. 1 Initial Flow
  FL2IC = 33;            % [ml/sec] Pulm. Veins 2 Initial Flow
  FLAIC =  0;            % [ml/sec] Left Atrium Initial Flow
  FLVIC =  0;            % [ml/sec] Left Ventricle Initial Flow
  FA1IC =  4.6;          % [ml/sec] Aortic 1 Initial Flow
  FV2IC = 95;            % [ml/sec] Syst. Veins 2 Initial Flow
  FRAIC =  0;            % [ml/sec] Right Atrium Initial Flow
  FRVIC =  6;            % [ml/sec] Right Ventricle Initial Flow
  FIS   = 10;            % [ml/sec] Sudden change in FLV at t > TIS

  QP1U =    7.8;         % [ml] Pulm. Art. 1 Unstressed Volume
  QP2U =   23.4;         % [ml] Pulm. Art. 2 Unstressed Volume
  QP3U =  210.5;         % [ml] Pulm. Art. 3 Unstressed Volume
  QL1U =   69;           % [ml] Pulm. Veins 1 Unstressed Volume
  QL2U =   69;           % [ml] Pulm. Veins 2 Unstressed Volume
  QLAU =  814.5;         % [ml] Left Atrium Unstressed Volume
  QLVU =   10;           % [ml] Left Ventricle Unstressed Volume
  QA1U =   35.1;         % [ml] Aortic 1 Unstressed Volume
  QA2U =   85;           % [ml] Aortic 2 Unstressed Volume
  QA3U =  710;           % [ml] Syst. Art. Unstressed Volume
  QV1U =  909;           % [ml] Syst. Veins 1 Unstressed Volume
  QV2U = 1948;           % [ml] Syst. Veins 2 Unstressed Volume
  QRAU = 1948;           % [ml] Right Atrium Unstressed Volume
  QRVU =   10;           % [ml] Right Ventricle Unstressed Volume

  LD =    45;            % [gr/cm^4/sec^2] Diastolic L Ventricle Stiffness
  LSI = 2500;            % [gr/cm^4/sec^2] Systolic L Ventricle Init. Stiffness
  RD =    72;            % [gr/cm^4/sec^2] Diastolic R Ventricle Stiffness
  RSI =  350;            % [gr/cm^4/sec^2] Systolic R Ventricle Init. Stiffness
  DLS = -625;            % [gr/cm^4/sec^2] Sudden change in LS at t > THI
  DRS =  -87.5;          % [gr/cm^4/sec^2] Sudden change in RD at t > THI
  
  SV1 = .9;	             % Stiffness first harmonic factor
  SV2 = .25;	           % Stiffness second harmonic factor

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

  THI = 2.4;             % [sec] Time of Ventricular Infarction
  TIS = 3;               % [sec] Time of Ventricular Infusion
  
  % 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', {'Healthy', ...
                    'Myocardial Infarction', 'Infusion or Bleeding'}, ...
                    '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_PressureFlow1')
  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 'Myocardial Infarction'
      Mode = 2;
    case 'Infusion or Bleeding'
      Mode = 3;
    end
  end % ModeCB

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

    % Independent Variable
    dt = 0.02;
    tmax = 3.2 * (Mode == 1) + 9 * (Mode ~= 1);
    tspan = (0:dt:tmax)';

    % Dependent Variables
    FINIT = [FP1IC FL2IC FLAIC FLVIC FA1IC FV2IC FRAIC FRVIC];
    QP1IC = QP1U + PP1EDM * CP1;
    QP2IC = QP2U + PP2EDM * CP2;
    QP3IC = QP3U + PP3EDM * CP3;
    QL1IC = QL1U + PL1EDM * CL1;
    QL2IC = QL2U + PL2EDM * CL2;
    QLAIC = QLAU + PLAEDM * CLA;
    QLVIC = QLVU + PLVEDM / LD;
    QA1IC = QA1U + PA1EDM * CA1;
    QA2IC = QA2U + PA2EDM * CA2;
    QA3IC = QA3U + PA3EDM * CA3;
    QV1IC = QV1U + PV1EDM * CV1;
    QV2IC = QV2U + PV2EDM * CV2;
    QRAIC = QRAU + PRAEDM * CRA;
    QRVIC = QRVU + PRVEDM / RD;
    QINIT = [QP1IC QP2IC QP3IC QL1IC QL2IC QLAIC QLVIC ...
             QA1IC QA2IC QA3IC QV1IC QV2IC QRAIC QRVIC];

    LS = LSI + DLS * (Mode == 2) .* (tspan > THI);
    RS = RSI + DRS * (Mode == 2) .* (tspan > THI);
    X = rem(tspan, TH);
    XX = (X <= TS + dt / 2);
    SSW = SV1 * sin(pi * XX .* X / TS) - SV2 * sin(2 * pi * XX .* X / TS);
    ACT = min(max(SSW, 0), 1);
    SLV = LD * (1 - ACT) + LS .* ACT;
    SRV = RD * (1 - ACT) + RS .* ACT;

    [tout res] = ode23(@myode, tspan, [FINIT QINIT], options);
    [PRV PP1 PLV PA1] = myode2(tspan, res);
    F = res(:, 1:8);
    Q = res(:, 9:22);

    figure('name', 'STIFF_ACT_VOLUME', 'NumberTitle', 'off');
    subplot(2, 2, 1)
    plot(tout, SRV, '-k');
    title('Right Stiffness')
    xlabel('Time [sec]')
    ylabel('Stiffness [gr/cm^4/sec^2]')
    axis([0 tmax 0 400])
    
    subplot(2, 2, 2)
    plot(tout, ACT, '-b');
    title('Activity Function')
    xlabel('Time [sec]')
    axis([0 tmax 0 1.1])

    subplot(2, 2, 3)
    plot(tout, SLV, '-r');
    title('Left Stiffness')
    xlabel('Time [sec]')
    ylabel('Stiffness [gr/cm^4/sec^2]')
    axis([0 tmax 0 2700])

    subplot(2, 2, 4)
    plot(tout, Q(:, 14), '-k', tout, Q(:, 7), '-r');
    title('Ventricular Volumes')
    xlabel('Time [sec]')
    ylabel('Volume [ml]')
    legend('QRV', 'QLV')
    axis([0 tmax 60 165])

    if Mode == 1
      figure('name', 'PRESSURE_FLOW_LOCUS', 'NumberTitle', 'off');
      subplot(2, 2, 1)
      plot(tout, PRV / 1332, '-k', tout, PP1 / 1332, '-r');
      title('Right Ventricular and Pulm. Artery Pressures')
      xlabel('Time [sec]')
      ylabel('Pressure [mmHg]')
      legend('PRV', 'PP1')
      axis([0 tmax 0 30])

      subplot(2, 2, 2)
      plot(tout, F(:, 5), '-k', tout, F(:, 3), '-r');
      title('Aortic and Left Atrial Flows')
      xlabel('Time [sec]')
      ylabel('Flow [ml/sec]')
      legend('FA1', 'FLA')
      axis([0 tmax 0 900])

      subplot(2, 2, 3)
      plot(tout, PLV / 1332, '-k', tout, PA1 / 1332, '-r');
      title('Left Ventricular and Aortic Pressures')
      xlabel('Time [sec]')
      ylabel('Pressure [mmHg]')
      legend('PLV', 'PA1')
      axis([0 tmax 0 150])

      subplot(2, 2, 4)
      plot(Q(:, 7), PLV / 1332, '-r');
      title('Left Ventricle Locus')
      xlabel('Volume [ml]')
      ylabel('Pressure [mmHg]')
      axis([0 150 0 150])
    end % if Mode == 1
    
    % Solve ODEs
    function myodeout = myode(t, myodein)
      F = myodein(1:8);
      Q = myodein(9:22);
      
      PP1  = (Q(1)  - QP1U) / CP1 + RPW1 * (F(8) - F(1));
      PP2  = (Q(2)  - QP2U) / CP2;
      PP3  = (Q(3)  - QP3U) / CP3;
      PL1  = (Q(4)  - QL1U) / CL1;
      PL2  = (Q(5)  - QL2U) / CL2;
      PLA  = (Q(6)  - QLAU) / CLA;
      PLV  = (Q(7)  - QLVU) * SLV(round(t / dt) + 1);
      PA1  = (Q(8)  - QA1U) / CA1 + RPW2 * (F(4) - F(5));
      PA2  = (Q(9)  - QA2U) / CA2;
      PA3  = (Q(10) - QA3U) / CA3;
      PV1  = (Q(11) - QV1U) / CV1;
      PV2  = (Q(12) - QV2U) / CV2;
      PRA  = (Q(13) - QRAU) / CRA;
      PRV  = (Q(14) - QRVU) * SRV(round(t / dt) + 1);

      FP2  = (PP2 - PP3) / RP2;
      FP3  = (PP3 - PL1) / RP3;
      FL1  = (PL1 - PL2) / RL1;
      FLAD = (PLA - PLV - RLA * F(3)) / LLA;
      FLVD = (PLV - PA1 - RLV * F(4)) / LLV;
      FA2  = (PA2 - PA3) / RA2;
      FA3  = (PA3 - PV1) / RA3;
      FV1  = (PV1 - PV2) / RV1;
      FRAD = (PRA - PRV - RRA * F(7)) / LRA;
      FRVD = (PRV - PP1 - RRV * F(8)) / LRV;

      FI   = FIS * (Mode == 3 && t > TIS);

      % Differential Equations
      myodeout = [
        (PP1 - PP2 - RP1 * F(1)) / LP1                            % FP1:t
        (PL2 - PLA - RL2 * F(2)) / LL2                            % FL2:t
        FLAD * (FLAD >= 0 && F(3) < 1e4 || FLAD <= 0 && F(3) > 0) % FLA:t
        FLVD * (FLVD >= 0 && F(4) < 1e5 || FLVD <= 0 && F(4) > 0) % FLV:t
        (PA1 - PA2 - RA1 * F(5)) / LA1                            % FA1:t
        (PV2 - PRA - RV2 * F(6)) / LV2                            % FV2:t
        FRAD * (FRAD >= 0 && F(7) < 1e4 || FRAD <= 0 && F(7) > 0) % FRA:t
        FRVD * (FRVD >= 0 && F(8) < 1e5 || FRVD <= 0 && F(8) > 0) % FRV:t
        F(8) - F(1)                                               % QP1:t
        F(1) - FP2                                                % QP2:t
        FP2 - FP3                                                 % QP3:t
        FP3 - FL1                                                 % QL1:t
        FL1 - F(2)                                                % QL2:t
        F(2) - F(3)                                               % QLA:t
        F(3) - F(4)                                               % QLV:t
        F(4) - F(5) + FI                                          % QA1:t
        F(5) - FA2                                                % QA2:t
        FA2 - FA3                                                 % QA3:t
        FA3 - FV1                                                 % QV1:t
        FV1 - F(6)                                                % QV2:t
        F(6) - F(7)                                               % QRA:t
        F(7) - F(8)                                               % QRV:t
      ];
    end % myode
    
    % Calculate rest of dependent variables after solving ODEs
    function [PRV PP1 PLV PA1] = myode2(t, myodein)
      F = myodein(:, 1:8);
      Q = myodein(:, 9:22);
      
      PP1  = (Q(:, 1)  - QP1U) / CP1 + RPW1 * (F(:, 8) - F(:, 1));
      PLV  = (Q(:, 7)  - QLVU) .* SLV(round(t / dt) + 1);
      PA1  = (Q(:, 8)  - QA1U) / CA1 + RPW2 * (F(:, 4) - F(:, 5));
      PRV  = (Q(:, 14) - QRVU) .* SRV(round(t / dt) + 1);
    end % myode2
  end % RunCB
end % Rideout_PressureFlow1