% The function resp_act.m computes respiratory-related waveforms
% over the entire integration period.
%
% Function arguments:
%	th - current parameter values
%	at - column vector containing the arrival times of each
%            respiratory cycle
%       deltaT - half of the integration step size (s)
%	stime - start time
%	etime - end time
%
% Function outputs:
%	thbreathe - 4xtime vector containing respiratory-related data
%	          - [Qlu; Pth; dPth; Ppac], where time is length of 
%                   the vector stime:deltaT:etime
%

function thbreathe = resp_act(th,at,deltaT,stime,etime)

% Establishing times in which the respiratory-related waveforms
% are to be computed.
time = stime:deltaT:etime;

% Pre-allocating memory for the function outputs.
thbreathe = zeros(4,length(time));

% No breathing.  
if (length(at) == 1)

	% Setting respiratory-related waveforms to their functional
	% reserve values.
	thbreathe(1,:) = th(98)*ones(1,length(time));
	thbreathe(2,:) = th(23)*ones(1,length(time));
	thbreathe(3,:) = zeros(1,length(time));
	thbreathe(4,:) = th(25)*ones(1,length(time));

% Step change in breathing.
elseif (length(at) == 2)

	  % Parametrization for step Qlu function  -- hyperbolic tangent.

	     % Width of transition band of step with small values
	     % indicating a small transition band.
	     tau = 0.25;

	     % Height of step.
	     sigmailv = at(2);

	  time1 = stime:deltaT:th(103);
	  time2 = time1(end)+deltaT:deltaT:etime;

	  Q = (sigmailv./(1+exp(-(time-th(103))/tau)))+th(98);	

	  dQ = ((sigmailv/tau)*exp(-(abs(time-th(103)))/tau))./((1+exp(-(abs(time-th(103)))/tau)).^2);

	  d2Q1 = -((sigmailv/tau^2)*exp(-(abs(time1-th(103)))/tau).*(exp(-(abs(time1-th(103)))/tau)-1))./((1+exp(-(abs(time1-th(103)))/tau)).^3);
	
	  d2Q2 = ((sigmailv/tau^2)*exp(-(abs(time2-th(103)))/tau).*(exp(-(abs(time2-th(103)))/tau)-1))./((1+exp(-(abs(time2-th(103)))/tau)).^3);

	  d2Q = [d2Q1 d2Q2];
	  thbreathe(1,:) = Q;
	  thbreathe(2,:) = th(91)-th(100)*dQ-(1/th(101))*(Q-th(102));
	  thbreathe(3,:) = -th(100)*d2Q-(1/th(101))*dQ;
	  thbreathe(4,:) = th(91)-th(100)*dQ;

% Impulse change in breathing.
elseif (length(at) == 3)

          % Parametrization for impulse Qlu function -- derivative of
          % hyperbolic tangent.

	     % Width of transition band of step with small values
	     % indicating a small transition band.
	     tau = 0.25;

	     % Area of impulse.
	     sigmailv = at(2);

	  time1 = stime:deltaT:th(103);
	  time2 = time1(end)+deltaT:deltaT:etime;

	  Q = ((sigmailv/tau)*exp(-(abs(time-th(103)))/tau))./((1+exp(-(abs(time-th(103)))/tau)).^2)+th(98);
	  dQ1 = -((sigmailv/tau^2)*exp(-(abs(time1-th(103)))/tau).*(exp(-(abs(time1-th(103)))/tau)-1))./((1+exp(-(abs(time1-th(103)))/tau)).^3);
	  dQ2 = ((sigmailv/tau^2)*exp(-(abs(time2-th(103)))/tau).*(exp(-(abs(time2-th(103)))/tau)-1))./((1+exp(-(abs(time2-th(103)))/tau)).^3);
	  dQ = [dQ1 dQ2];
	  d2Q = (sigmailv/tau^3)*(exp(-3*(abs(time-th(103)))/tau)-4*exp(-2*(abs(time-th(103)))/tau)+exp(-(abs(time-th(103)))/tau))./((1+exp(-(abs(time-th(103)))/tau)).^4);
	
	  thbreathe(1,:) = Q;
	  thbreathe(2,:) = th(91)-th(100)*dQ-(1/th(101))*(Q-th(102));
	  thbreathe(3,:) = -th(100)*d2Q-(1/th(101))*dQ;
	  thbreathe(4,:) = th(91)-th(100)*dQ;

else

	flag = at(1);
	at = at(2:length(at));

	% Fixed-rate breathing -- sinusoidal Qlu variations plus
	% functional reserve volume.
	if (flag == 0)
   
	   % Calculating the respiratory-related waveforms all at once.
	   thbreathe(1,:) = (th(77)/2)*(1-cos(2*pi*time/th(42)))+th(98);
	   thbreathe(2,:) = th(91)-th(100)*(th(77)*pi/th(42))*sin(2*pi*time/th(42)) - (1/th(101))*(thbreathe(1,:)-th(102));
	   thbreathe(3,:) = -th(100)*(th(77)/2)*(2*pi/th(42))^2*cos(2*pi*time/th(42))-(1/th(101))*(th(77)*pi/th(42))*sin(2*pi*time/th(42));
	   thbreathe(4,:) = th(91)-th(100)*(th(77)*pi/th(42))*sin(2*pi*time/th(42));		  
	% Random-interval breathing -- each Qlu respiratory cycle is one period
	% of a sinusoid plus the functional reserve volume, but the duration of
	% the periods are random.
	else

           % Prior to initiation of first respiratory cycle.  Setting
           % respiratory-related waveforms to functional reserve values.
	   startm = 1;
	   m = 1;
	   while (time(m) < at(1))

		m=m+1;	

	   end
	   endm=m-1;
	   thbreathe(1,startm:endm) = th(98)*ones(1,endm);
	   thbreathe(2,startm:endm) = th(23)*ones(1,endm);
	   thbreathe(3,startm:endm) = zeros(1,endm);
	   thbreathe(4,startm:endm) = th(25)*ones(1,endm);

	   % After initiation of first respiratory cycle, calculating
	   % respiratory-related waveforms one respiratory cycle at
	   % a time.
	   for counting = 2:length(at)

	       startm = endm+1;
	       m = startm;
	       while (m <= length(time) & time(m) < at(counting))
	
         	     m=m+1;

	       end
	       endm=m-1;

	       period = at(counting)-at(counting-1);
			

               % Establishing tidal volume for each respiratory cycle.

	          % Constant instantaneous alveolar ventilation rate
	          % set to mean rate of fixed-rate breathing (ml/s).
	          mbrealscalar(flag == 1);
	          if (flag == 1)

		       meanavr = (th(77)-th(99))/th(42);
		       Qtr = meanavr*period + th(99);

	          % Constant tidal volume resulting in a nearly constant
                  % instantaneous alveolar ventilation rate.
	          else

	               Qtr = 236.36*(period^(1/2)); 

	          end

	       thbreathe(1,startm:endm) = (Qtr/2)*(1-cos(2*pi*(time(startm:endm)-at(counting-1))/period))+th(98);
	       thbreathe(2,startm:endm) = th(91)-th(100)*(Qtr*pi/period)*sin(2*pi*(time(startm:endm)-at(counting-1))/period) - (1/th(101))*(thbreathe(1,startm:endm)-th(102));
	       thbreathe(3,startm:endm) = -th(100)*(Qtr/2)*(2*pi/period)^2*cos(2*pi*(time(startm:endm)-at(counting-1))/period)-(1/th(101))*(Qtr*pi/period)*sin(2*pi*(time(startm:endm)-at(counting-1))/period);
	       thbreathe(4,startm:endm) = th(91)-th(100)*(Qtr*pi/period)*sin(2*pi*(time(startm:endm)-at(counting-1))/period);

	   end

	end

end