%=============================================================*
%                                                             *
% MATLAB SCRIPT: main.m                                       *
% Version: 1.0 (created on 12 Aug. 2018)                      *
% AUTHORS: M.U. Qureshi, M.J. Colebank, M.S. Olufsen          *
%          Department of Mathematics                          *
%          North Carolina State University, Raleigh, USA      *
% DATE UPDATED: 12 Aug. 2018.                                 *
%                                                             *
% DESCRIPTION: This script is creates the inteface with the   *
% C++ code for the pulmonary arterial network in mouse. The   *
% script creats and runs the executable by passing selected   *
% model parameters and plots the hemodynamic waveforms        *
%                                                             *
%=============================================================*
%%
clear variables;
close all;
clc;

% These commands need to only be run ONCE in order to run the files. This
% can be run each iteration, but will slow down the run time. It is
% encouraged to comment this out after the first run of this code.
! make clean -f Makefile 
! make -f Makefile
! chmod +x sor06

%% Load in the actual waveform data
P1_dat = load('p_v1.dat');
P2_dat = load('p_v2.dat');
P3_dat = load('p_v3.dat');
Pdata = [P1_dat P2_dat P3_dat];

Q1_dat = load('q_v1.dat');
Q2_dat = load('q_v2.dat');
Q3_dat = load('q_v3.dat');
Qdata = [Q1_dat Q2_dat Q3_dat];

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHANGE k3
k3 = 52021;    % Constant stiffness; valid on [10^4, 10^5]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do not change: Set at nominal values
r1 = 0.20193;  % Scaling factor for proximal resistances
r2 = 0.8889;   % Scaling factor for distal resistances
c1 = 1.4665;   % Scaling factor for preripheral compliances
HB = 9;        % Number of heartbeats for 
cycles = 1;    % Number of cycles to be plotted
id = 1;        % file identfier
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

status = unix(sprintf('./sor06 %f %f %f %f %d %d %d',k3, r1, r2, c1, HB, cycles, id));

if status == 0 %%% Will be zero if the PDE is solved successfully

indx = num2str(id);

%% Loop over all 3 of the vessels and plot the output
for j=1:3
fname = strcat('pu',num2str(j),'_',indx,'.2d');
data = load (fname);

[tt,x,p,q,A,C] = gnuplot(data);

% The subscripts p,m, and d correspond to predictions that are located at
% the proximal, midpoint, and distal ends of the vessels. We only use the
% midpoint data here.

P_p = p(:,1);  Q_p = q(:,1);   A_p = A(:,1); C_p = C(:,1);
P_m = p(:,floor((end)/2));  Q_m = q(:,floor((end)/2));   A_m = A(:,floor((end)/2)); C_m = C(:,floor((end)/2));
P_d = p(:,end); Q_d = q(:,end); A_d = A(:,end); C_d = C(:,end);



time = tt(:,1);

Period = 0.11;

t = time - (HB-cycles)*Period;


%% Plotting Results


% Pressure
figure(j*10);clf;
plot(t,Pdata(:,j),'-k',t,P_m,'--r','LineWidth',2);
grid on;
xlabel('t [s]');
ylabel('p [mmHg]');

figure(j*10 + 1);clf;
plot(t,Qdata(:,j),'-k',t,Q_m,'--r','LineWidth',2);
grid on;
xlabel('t [s]');
ylabel('q [ml/s]');

% Other variables that can be output. Not needed for event purposes.

% figure(3);clf;
% plot(t,C_p,'LineWidth',2);
% grid on;
% xlabel('t [s]');
% ylabel('PWV [cm/s]');
% 
% figure(4);clf;
% plot(Qdat,Pdat,'-k',Q_p,P_p,'-r','LineWidth',2);
% grid on;
% xlabel 'q [ml/s]';
% ylabel 'p [mmHg]';
% 
% figure(5);clf;
% plot(P_p,A_p,'LineWidth',2);
% grid on;
% xlabel 'p [mmHg]';
% ylabel 'A [cm^2]';
end % loop

else
	"No output returned"
end % if status == 0