% The function cpreflex.m implements the cardiopulmonary baroreflex as a
% feedback and control system.
%
% Function arguments:
% S - 3x1 vector containing the desired granularity, duration
% of autonomic impulse responses, and averaging ratio of
% length of averaging window to Sgran
% - [Sgran; Slength; Sratio]
% bPratr - a Slengthx1 column vector whose elements are the Slength most
% recent Pratr samples (Sgran sampling period), where each sample
% corresponds to Pratr data averaged over Sgran*Sratio length intervals
% th - current parameter values
% thsp - 7x1 vector containing parameter setpoint values
% - [Pasp; Fsp; Rasp; Clssp; Crssp; Qvosp; Pratrsp]
%
% Function outputs:
% thc - mandated cardiopulmonary baroreflex adjustment to the controllable parameters
% - [dF; dRa; dCls; dCrs; dQvo]
%
function [thc] = cpreflex(S,bPratr,th,thsp)
% Assigning variables.
Sgran = S(1);
Slength = S(2);
Sratio = S(3);
% Pre-allocating memory for autonomic impulse responses.
sir = zeros(Slength,1);
vir = zeros(Slength,1);
% Creating unit-area beta-sympathetic impulse response.
st1 = ((2/Sgran)-((Sratio*Sgran/2)/Sgran))+1;
en1 = length(Sgran:Sgran:3.0)+st1-1;
st2 = en1+1;
en2 = length(Sgran:Sgran:25.0)+st2-1;
sir(st1:en1) = ((1.0/42.0)*(Sgran:Sgran:3.0))';
sir(st2:en2) = ((-1.0/350.0)*(Sgran:Sgran:25.0)+(1.0/14.0))';
sir = Sgran*sir;
% Creating unit-area parasympathetic impulse response.
numz = (0.5-(Sratio*Sgran/2))/Sgran;
st1 = 1+numz;
en1 = length(Sgran:Sgran:1)+st1-1;
st2 = en1+1;
en2 = length(Sgran:Sgran:1)+st2-1;
vir(st1:en1) = (1)*(Sgran:Sgran:1)';
vir(st2:en2) = (1)*(-(Sgran:Sgran:1)+1)';
vir = Sgran*vir;
% Creating unit-area alpha-sympathetic impulse response for manipulating Ra.
air1 = zeros(Slength,1);
st1 = ((th(50)/Sgran)-((Sratio*Sgran/2)/Sgran))+1;
en1 = length(Sgran:Sgran:3.0)+st1-1;
st2 = en1+1;
en2 = length(Sgran:Sgran:25.0)+st2-1;
air1(st1:en1) = ((1.0/42.0)*(Sgran:Sgran:3.0))';
air1(st2:en2) = ((-(3/42)/25)*(Sgran:Sgran:25.0)+(1.0/14.0))';
air1 = Sgran*air1;
% Creating unit-area alpha-sympathetic impulse response for manipulating Qvo.
air2 = zeros(Slength,1);
st1 = ((th(51)/Sgran)-((Sratio*Sgran/2)/Sgran))+1;
en1 = length(Sgran:Sgran:3.0)+st1-1;
st2 = en1+1;
en2 = length(Sgran:Sgran:25.0)+st2-1;
air2(st1:en1) = ((1.0/42.0)*(Sgran:Sgran:3.0))';
air2(st2:en2) = ((-(3/42)/25)*(Sgran:Sgran:25.0)+(1.0/14.0))';
air2 = Sgran*air2;
% Forming error signal.
mbrealscalar(th(60) == 9999);
if (th(60) == 9999)
% No static saturation.
epd = bPratr-thsp(7);
else
% Static saturation.
epd = th(60)*atan((bPratr-thsp(7))/th(60));
end
% Convolving error signal with each of the impulse responses.
bres = th(38)*0.001*sum(epd.*sir);
pres = th(39)*0.001*sum(epd.*vir);
ares1 = th(92)*th(55)*0.001*sum(epd.*air1);
ares2 = th(92)*th(56)*0.001*sum(epd.*air2);
% Computing mandated change to F.
a = -thsp(2)^2*(bres*th(96)+pres*th(58));
% Computing mandated change to Ra.
b = -(65/9)*ares1;
% Computing mandated change to Cls and Crs.
clsp = ((th(26)-th(9))/th(31))*thsp(4);
crsp = ((th(27)-th(12))/th(32))*thsp(5);
cl = 2*(.2/.5088)*bres*th(94);
cr = (crsp/clsp)*cl;
cl = (th(31)/(th(26)-th(9)))*cl;
cr = (th(32)/(th(27)-th(12)))*cr;
% Computing mandated change to Qvo.
d = (90000/9)*ares2;
% Assigning mandated changes to function output.
thc = [a; b; cl; cr; d];