% The function abreflex.m implements the arterial baroreflex as a negative
% 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] 
%	bP - a Slengthx1 column vector whose elements are the Slength most 
%	     recent Pa samples (Sgran sampling period), where each sample 
%            corresponds to Pa averaged over Sgran*Sratio length intervals
%	th - current parameter values
%	thsp - 7x1 vector containing parameter setpoint values
%	     - [Pasp; Fsp; Rasp; Clssp; Crssp; Qvosp; Pratrsp]
%	deltaPasp - a Slengthx1 vector containing the values of fluctuations in Pasp
%
% Function outputs:
%	thc - mandated arterial baroreflex adjustment to the controllable parameters
%	    - [dF; dRa; dCls; dCrs; dQvo]
%

function [thc] = abreflex(S,bP,th,thsp,deltaPsp)

% 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(45)/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(49)/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(46) == 9999);
if (th(46) == 9999)

        % No static saturation.
	epd = bP-(thsp(1)+deltaPsp);

else

	% Static saturation.
	epd = th(46)*atan((bP-(thsp(1)+deltaPsp))/th(46));

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(37)*0.001*sum(epd.*air1);
ares2 = th(92)*th(54)*0.001*sum(epd.*air2);

% Computing mandated change to F.
a = -thsp(2)^2*(bres*th(95)+pres*th(97));

% Computing mandated change to Ra.
b = -(2/th(3))*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(93);
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 = 2*(750/.5088)*ares2;

% Assigning mandated changes to function output.
thc = [a; b; cl; cr; d];