% The function simulate.m implements various computational models of the
% human cardiovascular system.  The models include the following three
% major components:
%
% 1) Lumped parameter, pulsatile heart and circulation
%    (various preparations)
% 2) Short-term regulatory system
%    a) Arterial baroreflex system
%    b) Cardiopulmonary baroreflex system
%    c) Direct neural coupling mechanism between Qlu and F
% 3) Resting physiologic perturbations
%    a) respiratory activity
%    b) eExogenous disturbances to
%         i) Ra
%        ii) F
%	iii) Qvo
%	 iv) Pasp
%
% Function arguments:
%	flag - 8x1 vector indicating the status of the model
%	       (overrides parameter values assigned in th)
%	     - [preparation breathing dncm baro dra dpasp df dqvo]
% 
%	     - preparation - status of pulsatile heart and circulation
%			   - 0 - intact circulation (nominal)
%			   - 1 - heart-lung unit designed for measuring cardiac function curves 
%			   - 2 - systemic circulation designed for measuring venous return curves
%			   - 3 - left ventricle
%			         (automated step changes in Ppv and Pa NOT implemented)
%			   - 4 - intact circulation for measurement of single ventricular
%				 contraction Pa response
%			   - 5 - intact circulation with only linear elements
%			   - 6 - intact circulation with third-order systemic arteries
%			   - 7 - intact circulation with nonlinear systemic arterial compliance
%			   - 8 - intact circulation with third-order systemic arteries
%				 and nonlinear large, elastic arterial compliance
%			   - 9 - intact circulation with an arterial pressure reservoir
%			   - 10 - intact circulation with contracting atria
%	     -	breathing - status of breathing (Qlu)
%			  - 0 - no breathing (Qlu constant)
%			  - 1 - fixed-rate breathing (Qlu periodic)
%			  - 2 - random-interval breathing (Qlu random; 5 s mean period)
%			        instantaneous alveolar ventilation rate constant over time
%			  - 3 - random-interval breathing (Qlu random; 5 s mean period)
%				tidal volume fixed; instantaneous alveolar ventilation rate
%				not constant
%			  - 4 - step increase in breathing (Qlu step)
%			  - 5 - impulse increase in breathing (Qlu impulse)
%	     - dncm - status of direct neural coupling mechanism between Qlu and F
%		    - 0 - dncm off
%		    - 1 - dncm on (not applicable for step increase in Qlu)
%	     - baro - status of baroreflexes
%		    - 0 - arterial and cardiopulmonary baroreflexes off
%		    - 1 - arterial baroreflex on only
%		    - 2 - cardiopulmonary baroreflex on only
%		    - 3 - arterial and cardiopulmonary baroreflexes on
%	     - dra - status of bandlimited disturbance to Ra
%		   - 0 - dra off
%		   - 1-3 - dra on
%		         - 1 - 1/f fluctuations
%			 - 2 - bandlimited, white fluctuations
%			 - 3 - sinusoidal fluctuations
%	     - dpasp - status of sinusoidal disturbance to Pasp
%		  - 0 - dpasp off
%		  - 1 - dpasp on
%	     - df - status of 1/f disturbance to F
%		  - 0 - df off
%		  - 1 - df on
%	     - dqvo - status of bandlimited disturbance to Qvo
%		    - 0 - dqvo off
%		    - 1 - dqvo on
%	th - 113x1 vector containing the parameter values characterizing
%            the simulation
%       outputfile - prefix of the output files in MIT format generated by the model
%	parameterfile - name of a file which contains the parameter values
%                       characterizing the model and its execution
%	signals - the signals to be viewed by WAVE as the data are generated
%
%	NOTE: If this program is executed in the MATLAB environment, it will
%             not be possible to implement on-line viewing and parameter updating
%             Thus, the last three arguments should be omitted.
%             
% Function outputs:
%	P - matrix representing simulated pressures [mmHg]
%	  = [Pl Pe Pv Pr Ppa Ppv Pth Palv Pra Pm] for preparation == 6,8
%	  = [Pl Pa Pv Pr Ppa Ppv Pth Palv Pra Pla] for preparation == 10
%	  = [Pl Pa Pv Pr Ppa Ppv Pth Palv Pra] otherwise
%	Q - matrix representing simulated volumes [ml]
%	  = [Ql Qe Qv Qr Qpa Qpv Qlu Qm] for preparation == 6,8
%	  = [Ql Qa Qv Qr Qpa Qpv Qlu Qra Qla] for preparation == 10
%	  = [Ql Qa Qv Qr Qpa Qpv Qlu] otherwise
%	q - matrix representing simulated flow rates [ml/s]
%	  = [qpv ql qm qv qr qpa qe] for preparation == 6,8
%	  = [qla ql qa qra qr qpa qv qpv] for preparation == 10
%	  = [qpv ql qa qv qr qpa] otherwise
%	ap - matrix representing simulated adjustable parameters
%	   = [Cls Crs Qvo Ra F Tsys]
%	ve - matrix representing the ventricular variable elastances
%            [mmHg/ml]
%	   = [El Er Ela Era] for preparation == 10
%	   = [El Er] otherwise
%	qrs - vector representing times of ventricular contractions [s]
%	t - vector representing the associated time [s]
%	num - matrix representing either cardiac function or venous
%	      return curve numerics
%	    = [mra mpa mco] for preparation == 1
%	    = [mra mco] for preparation == 2
%	    = [0, 0] otherwise
%
% Function usage (in MATLAB):
%
%	[P,Q,q,ap,ve,qrs,t,num] = simulate(flag,th);
%

function [P,Q,q,ap,ve,qrs,t,num] = simulate(flag,th,outputfile,parameterfile,signals)

% Declaring C functions.
%#function wave_remote read_key write_param

% Checking number of arguments.  If there are only two arguments,
% this indicates MATLAB execution and thus on-line viewing is 
% turned off.  (If it were left on, MATLAB would produce an error.)
if (nargin == 2)
   signals = '-1';
end

% Checking dimensionality of the argument flag -- status parameters.
% If flag is of incorrect dimension, the program will exit in error.
if (length(flag) ~= 8)
	error('flag is not of correct dimension');
end

% Assigning status parameters to new variable names.

preparation = flag(1);
breathing = flag(2);
dncm = flag(3);
baro = flag(4);
dra = flag(5);
dpasp = flag(6);
df = flag(7);
dqvo = flag(8);

% Assigning numerical integration variables.

   % Number of integration steps.
   tics = th(106);

   % Integration step size (s).
   deltaT = 1/th(113);

% Assigning integrate and fire variables.

   % Fraction of current cardiac cycle that integration step is at.
   ipfm = 0; 

   % Current number of ventricular contractions.
   qrs_index = 1;

   % Time in current cardiac cycle that integration is at. 
   sumdeltaT = 0; 

% Short-term regulatory systems and resting physiologic 
% perturbation variables.
% Note there are restrictions on altering Sgran and Sratio.
% 1) Sgran*Sratio <= 1.0.
% 2) (1.5/Sgran) must be an integer.
% 3) ((2/Sgran)-((Sratio*Sgran/2)/Sgran)) must be an integer. 
% 4) (0.5-(Sratio*Sgran/2))/Sgran must be an integer. 

  % Granularity (s).
  Sgran = 0.0625;

  % Ratio of length of averaging window to Sgran
  % (must be an integer).
  Sratio = 4;

  % Duration (samples).
  Slength = 40/Sgran;


% Generation of a vector of arrival times for each respiratory cycle.  

   % at(1) = 0 & length(at) ~= 1: Qlu periodic
   if (breathing == 1)

	at = zeros(ceil((deltaT*tics)+1)+1,1);      
	i = 2;
	x = (at(i) <= deltaT*tics+1);
	mbintscalar(x);
	while (x)
		i = i+1;
		at(i) = at(i-1)+th(42);
		x = (at(i) <= deltaT*tics+1);
		mbintscalar(x);
	end	
	at = at(1:i); 

   % at(1) = 1: Qlu random with constant alveolar ventilation rate (2)
   % at(1) = 2: Qlu random with fixed tidal volume (3)
   elseif (breathing == 2 | breathing == 3)

	at = zeros(ceil((deltaT*tics)+1)+1,1);
	if (breathing == 2)
		at(1) = 1;
	else
		at(1) = 2;
	end

	i = 2;
	at(i) = rand_int_breath;
	x = (at(i) <= deltaT*tics+1);
	mbintscalar(x);
	while (x)
		i = i+1;
		at(i) = at(i-1) + rand_int_breath;
		x = (at(i) <= deltaT*tics+1);
		mbintscalar(x);
	end
	at = at(1:i);

   % at(1) = 0 & length(at) == 1: Qlu constant
   elseif (breathing == 0) 

	at = 0;

   % at(1) = 3 (length(at) == 2): step Qlu (cannot implement dncm)
   elseif (breathing == 4)

	at = [3; th(36)];

   % at(1) = 4 (length(at) == 3): impulse Qlu
   elseif (breathing == 5)

	at = [4; th(36); (((deltaT*tics)+1)+1)];

   % Exiting program with error if breathing is not valid.
   else

	error('breathing is not a valid option');

   end	

% Pre-calculating Qlu, Pth, dPth, and Ppac (Palv) over the entire
% integration period from the generated arrival time vector.
thbreathe = resp_act(th,at,deltaT/2,0,(tics-1)*deltaT);

% Adding a step change to Qlu, Pth, dPth, and Ppac, if mandated via Qfrs.
if ((breathing == 0 | breathing == 1 | breathing == 2) & th(36) ~= 0)
   
   th(103) = 0;
   thbreathes = resp_act(th,[3; th(36)],deltaT/2,0,(tics-1)*deltaT);
   dthbreathe = [-th(98)*ones(1,length(thbreathes)); -th(23)*ones(1,length(thbreathes)); zeros(2,length(thbreathes))]; 
   thbreathe = thbreathe+thbreathes+dthbreathe;

end
   

% Creating the filter characterizing direct neural coupling mechanism
% betweeen Qlu and F at a sampling period of Sgran.  Decimating the 
% above Qlu to a sampling period of Sgran in order to implement this filter.
% (not applicable to step input of Qlu)
if (dncm == 1)

	hilvhr = dncm_filt(th,[Sgran; Slength; Sratio]);
	ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran);
	% Taking care of filtering edge effects.
	ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];
	rsastep = 1;

% Exiting program with error if dncm is not valid.
elseif (dncm ~= 0)

        error('dncm is not a valid option');

end

% Autoregulation of local vascular beds -- exogenous disturbance to Ra.

   % Generating sinusoidal Ra disturbance at a sampling period of Sgran
   % for entire integration period.
   if (dra == 3)

	amp = th(61);
	freq = th(62);
	n = floor(((deltaT*tics)+1)/Sgran)+(Sratio-1);
	i = 1:n;
	tr = (i-1)*Sgran;
	Ralvrp = (amp/(Sgran*2*pi*freq))*(cos(2*pi*freq*tr)-cos(2*pi*freq*(tr+Sgran)));
	Ralvr = zeros(1,length(Ralvrp)-Sratio+1);
	for kt = Sratio:length(Ralvrp)
		Ralvr(kt-(Sratio-1)) = sum(Ralvrp(kt-(Sratio-1):kt))/Sratio;
	end
	Ralvr = Ralvr+((0.01*amp)/sqrt(Sgran))*randn(1,length(Ralvr));

	lvrstep = 2;

   % Generating bandlimited 1/f Ra disturbance at a sampling period of Sgran
   % for entire integration period.
   elseif (dra == 1)

	hra = bl_filt(th(44),Sgran);
	n = floor(((deltaT*tics)+1)/Sgran)+length(hra);
	whitenoise = sqrt(1/Sgran)*th(47)*randn(1,n);
	dmin = -4;
	dmax = 0;
        honeoverfra = oneoverf_filt(th,Sgran,dmin,dmax);
	Ralvrw = conv(honeoverfra,whitenoise);
	Ralvrw = Ralvrw(length(Ralvrw)-length(n)+1:length(Ralvrw));	

	% Filtering 1/f disturbance to Ra.
	Ralvr = conv(hra,Ralvrw);
	Ralvr = Ralvr(length(hra)+1:length(Ralvrw));

	lvrstep = 2;

   % Generating white noise at a sampling period of Sgran for bandlimited 
   % disturbance to Ra.
   elseif (dra == 2)

	wra = sqrt(1/Sgran)*th(47)*randn(length(-60:Sgran:60),1);

   % Exiting program with error if dra is not valid.
   elseif (dra ~= 0)

	error('dra is not a valid option');

   end

% Generating sinusoidal Pasp disturbance at a sampling period of Sgran
% for entire integration period.
if (dpasp == 1);

	amp = th(63);
	freq = th(64);
	n = floor(((deltaT*tics)+1)/Sgran)+(Sratio-1)+Slength;
	i = 1:n;
	tr = (i-1)*Sgran;
	tempPsp = (amp/(Sgran*2*pi*freq))*(cos(2*pi*freq*tr)-cos(2*pi*freq*(tr+Sgran)));
	deltaPsp = zeros(length(tempPsp)-Sratio+1,1);
	for kt = Sratio:length(tempPsp)
		deltaPsp(kt-(Sratio-1)) = sum(tempPsp(kt-(Sratio-1):kt))/Sratio;
	end
	deltaPsp = deltaPsp+((0.1*amp)/sqrt(Sgran))*randn(length(deltaPsp),1);

	ceostep = 2;

elseif (dpasp == 0)

	deltaPsp = 0;

% Exiting program with error if dpasp is not valid.
else
	error('dpasp is not valid');
end

% Creating 1/f filter and white noise at a sampling period of Sgran
% for the exogenous disturbance to F. 
if (df == 1)

        dmin = -4; % Desired minimum frequency decade with 1/f character 
	dmax = 0; % Desired maximum frequency decade with 1/f character
	hans = ans_filt(th,[Sgran; Slength; Sratio]);
	honeoverf = oneoverf_filt(th,Sgran,dmin,dmax);
	hf = conv(honeoverf,hans);
	wf = sqrt(1/Sgran)*th(48)*randn(length(hf),1);

% Exiting program with error if df is not valid.
elseif (df ~= 0)

       error('df is not a valid option');

end

% Generating bandlimited, white disturbance to Qvo at a sampling period of Sgran
% for entire integration period.
if (dqvo == 1)

	hvvo = bl_filt(th(65),Sgran);	
	n = floor(((deltaT*tics)+1)/Sgran)+length(hvvo);
	Vvow = sqrt(1/Sgran)*th(66)*randn(1,n);
	Vvo = conv(hvvo,Vvow);
	Vvo = Vvo(length(hvvo)+1:length(Vvow));
	nvvostep = 2;

% Exiting program with error if dqvo is not valid.
elseif (dqvo ~= 0)

       error('dqvo is not a valid option');

end

% Pre-allocating memory for function output waveforms which is
% necessary for dramatic improvement in execution speed.
P = zeros(11,tics);
Q = zeros(9,tics);
q = zeros(8,tics);
t = zeros(1,tics);
ap = zeros(9,tics);
ve = zeros(4,tics);

% Pre-allocating memory and initializing variables for 
% cardiac function/venous return curve generation.

   if (preparation == 1)

      % Pre-allocating memory.
      Crd = ((th(27)-th(12))/th(32))*th(6);
      final = round((th(27)-th(12)-20)/Crd)+3;
      rarange = [final:-th(88):0, th(23)+1];
      lengthra = length(rarange);
      if (lengthra == 2);

         lengthra = 1;
         rarange = final;

      end
      lengthsa = length(30:th(90):th(31)-10);
      lengtht = lengthra*lengthsa;
      num = zeros(3,lengtht);

      % Initializing variables for family of cardiac function curves.
      if (lengthra > 1 & lengthsa > 1)

         countn = 0;
         countj = 0;
         th(30) = final;
         th(29) = 30;
         tmp = 0;

      % Initializing variables for a single cardiac function curve.
      elseif (lengthra > 1 & lengthsa == 1)

         countn = 0;
         countj = 0;
         th(30) = final;
         tmp = 0;

      % Initializing variables for curve of average ql versus Pa.
      elseif (lengthra == 1 & lengthsa > 1)

          final = th(30);
          rarange = final;
          th(29) = 30;
          tmp = 0;
   
      end
   
      % Initializing variables for averaging Pra and ql over each cardiac cycle.
      summrap = 0;
      summco = 0;
   
      % Initializing variable representing number of beats until voltage
      % source(s) is varied.
      fb_int = 10;
      
   elseif (preparation == 2)

      % Pre-allocating memory.
      Crs = ((th(27)-th(12))/th(32))*th(5);
      lengthra = length([Crs:th(89):60]);
      num = zeros(2,lengthra);
      if (length(num(1,:)) > 1)

         countn = 0;
         Crd = Crs;
         th(6) = (th(32)/(th(27)-th(12)))*Crd;

      end

      % Initializing variables for integrating Pra and qv over each cardiac cycle.
      summrap = 0;
      summco = 0;

      % Initializing variable representing number of beats until Crd
      % is varied.
      fb_int = 10;

   else

      num = zeros(2,1);       

   end

% Resetting ventricular compliance values for the intact circulatory
% preparation with only linear elements.
if (preparation == 5)
   th(1) = th(1)*((th(26)-th(9))/th(31));
   th(2) = th(2)*((th(26)-th(9))/th(31));
   th(5) = th(5)*((th(27)-th(12))/th(32));
   th(6) = th(6)*((th(27)-th(12))/th(32));
   th(31) = th(26)-th(9);
   th(32) = th(27)-th(12);
end

% Initializing the function output waveforms.

   % Adjustable parameters.
   ap(1:6,1) = [th(1)*((th(26)-th(9))/th(31)) th(5)*((th(27)-th(12))/th(32)) th(11) th(16) th(22) .3*sqrt(th(24))]';

   % Pulsatile heart and circulation preparation.
   if (preparation == 0 | preparation == 4 | preparation == 5 | preparation == 7)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = intact_init_cond(th);

   elseif (preparation == 1)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = hlu_init_cond(th);

   elseif (preparation == 2)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = sc_init_cond(th);

   elseif (preparation == 3)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = lv_init_cond(th);

   elseif (preparation == 5)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = init_cond(th);

   elseif (preparation == 6)

	[Ptemp,Q(1:7,1),q(1:7,1),ve(1:2,1)] = third_init_cond(th);
	P(1:6,1) = Ptemp(1:6);
	P(10:11,1) = Ptemp(7:8);
	Q(8,1) = Q(7,1); 

   elseif (preparation == 7)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = nac_init_cond(th);

   elseif (preparation == 8)

	[Ptemp,Q(1:7,1),q(1:7,1),ve(1:2,1)] = third_nac_init_cond(th);
	P(1:6,1) = Ptemp(1:6);
	P(10:11,1) = Ptemp(7:8);
	Q(8,1) = Q(7,1); 

   elseif (preparation == 9)

	[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = apr_init_cond(th);

   elseif (preparation == 10)

        [tempp,tempv,q(1:8,1),ve(:,1)] = a_init_cond(th);
	P(1:6,1) = tempp(1:6);
	P(9,1) = tempp(7);
	P(10:11,1) = tempp(7:8);
	Q(1:6,1) = tempv(1:6);
	Q(8:9,1) = tempv(7:8);

   % Exiting program with error if preparation is not valid,
   else

	error('preparation is not a valid option');

   end

   if (preparation ~= 10)

      P(9,1) = P(4,1);

   end

   % Respiratory-related waveforms.
   P(7,1) = thbreathe(2,1);
   P(8,1) = thbreathe(4,1);
   Q(7,1) = thbreathe(1,1);


% Assigning setpoint variables

  % Setpoint pressure for arterial baroreflex (mmHg).
  Pasp = th(40); 

  % Setpoint pressure for cardiopulmonary baroreflex (mmHg).
  Pratrsp = th(41);

  % F setpoint value (Hz).
  Fsp = th(22); 

  % Ra setpoint value (mmHg-s/ml).
  Rasp = th(16); 

  % Cls, Crs setpoint value (ml/mmHg).
  Clssp = th(1); 
  Crssp = th(5);

  % Qvo setpoint value (ml). 
  Qvosp = th(11); 

  % 7x1 vector of setpoint values.
  thsp = [Pasp; Fsp; Rasp; Clssp; Crssp; Qvosp; Pratrsp];

% Initializing baroreflex-related variables.
if (baro > 0 & baro < 4)

   % Initializing input pressures.  (The lowest
   % index indicates the most recent pressure.)
   bP = zeros(Slength,1);
   bPratr = zeros(Slength,1);
   for i = 1:Slength;
      bP(i) = P(2,1);
      bPratr(i) = P(9,1)-P(7,1);
   end		  
   bPtemp = zeros(Sratio,1);
   bPratrtemp = zeros(Sratio,1);
   for i = 1:Sratio
      bPtemp(i) = P(2,1);
      bPratrtemp(i) = P(9,1)-P(7,1);
   end

   % Initializing variables for integrating Pa, Pratr,
   % and t over each Sgran interval.
   sumP = 0; 
   sumPratr = 0;
   sumt = 0; 

% Exiting program with error if baro is not valid,
elseif (baro ~= 0)

   error('baro is not a valid option');

end

% Initializing parameter vector and other variables in order to
% implement linear interpolation of the controllable parameters.
thn = th;
thp = th;
tn = 0;
tp = 0;

% Initializing parameters for writing data and annotations to file
% in MIT format.  Calling the WAVE window.
if (strcmp(signals,'-1') == 0)

   % Initializing waveform file variables.
   kstart = 1;
   qrs_index_start = 1;
   fid = fopen([outputfile '.dat'],'wa');

   % Initializing annotation file.
   fid2 = fopen([outputfile '.qrs'],'wa');
   fwrite(fid2,1,'bit10');
   fwrite(fid2,1,'ubit6');
   fwrite(fid2,-1,'bit10');
   fwrite(fid2,62,'ubit6');
   Nbytes = length(parameterfile)+2;
   fwrite(fid2,Nbytes,'bit10');
   fwrite(fid2,63,'ubit6');
   fwrite(fid2,[parameterfile '.' num2str(0)],[num2str(Nbytes) '*uchar']);

   % Assigning annotation flag.
   if (th(105) == 0)
      annotator = '-1';
   else
      annotator = 'qrs';
   end

   % Initializing parameter update counter.
   count_pc = 0;

   % Calling the WAVE window.
   wave_remote( outputfile, annotator, 0, signals );

   % Initializing variables which permit annotation flag to be
   % updated on-line.
   wrflag = 1;
   outputfile2 = ['./' outputfile];
	       
end

% Generating the function output waveforms for each integration step.
for k = 2:tics


	% Generating the pressure waveforms and associated time for 
	% the current integration step -- fourth-order Runge-Kutta
	% integration.
	Ptempn = rk4([P(1:6,k-1); P(10:11,k-1)],[Q(1,k-1);Q(4,k-1)],thbreathe(:,2*(k-1)+1:2*(k-1)+3),th,ipfm,sumdeltaT,deltaT,preparation);
	P(1:6,k) = Ptempn(1:6);
	if (preparation == 10)
		P(9,k) = Ptempn(7);
	end
	P(10,k) = Ptempn(7);
	P(11,k) = Ptempn(8);

	% Updating time.
	t(k) = t(k-1)+deltaT;

	% Computing fraction of cardiac cycle that was surpassed
	% during the integration step.
	frac = deltaT*th(22);

	% Obtaining Pth, Ppac (Palv), and Qlu at the current integration step
	P(7,k) = thbreathe(2,2*(k-1)+1);
	P(8,k) = thbreathe(4,2*(k-1)+1);
	Q(7,k) = thbreathe(1,2*(k-1)+1);

	% Integrating Pa data using the Backward Euler method and 
	% denoting the time into the integration.
	sumt = sumt+deltaT;
	sumP = sumP+P(2,k)*deltaT;

	% Computing Pra and integrating Pratr using the Backward
	% Euler method.
	if (preparation ~= 10)
		x = (P(4,k) >= P(3,k));
		mbintscalar(x);
		if (x)
			P(9,k) = P(3,k);
		else
			P(9,k) = P(4,k);
		end
	end
	sumPratr = sumPratr+(P(9,k)-P(7,k))*deltaT;

	% Implementing the short-term regulatory system and/or resting
	% physiologic perturbations each Sgran s.
	x = (sumt > Sgran);
	mbintscalar(x);
	if (x)

		% Setting final th values of previous Sgran to initial th
		% value of current Sgran. 
		thp = thn;
		tp = tn;

		% Initializing thc -- a vector which contains the mandated
		% adjustments to controllable parameters of sampling period
		% Sgran.
		thc = zeros(5,1);

		% Implementing the arterial baroreflex.
		if (baro == 1 | baro == 3)

			% Calculating the integral of Pa precisely
			% over the Sgran interval.
			sumP = sumP-P(2,k)*deltaT;
			sumP = sumP+P(2,k)*(deltaT-(sumt-Sgran));

			% Placing the averaged Pa over the Sgran*Sratio 
			% interval in the most recent slot of the arterial 
			% baroreflex input.
			bPtemp = [(sumP/Sgran); bPtemp(1:Sratio-1)];
			avebP = sum(bPtemp)/Sratio;
			bP = [avebP; bP(1:Slength-1)];

			% Calculating the new baroreflex adjusted parameters 
			% for the next Sgran interval.
			
			   % For constant Pasp.
			   if (length(deltaPsp) == 1)
				thc = thc+abreflex([Sgran; Slength; Sratio],bP,th,thsp,deltaPsp);
			   % For sinusoidally varying Pasp.
			   else 
				thc = thc+abreflex([Sgran; Slength; Sratio],bP,th,thsp,deltaPsp(dpasp:dpasp+length(bP)-1));
				dpasp=dpasp+1;	
			   end
		end

		% Implementing the cardiopulmonary baroreflex.
		if (baro == 2 | baro == 3)

			% Calculating the integral of Pratr precisely
			% over the Sgran interval.
			sumPratr = sumPratr-(P(9,k)-P(7,k))*deltaT;
			sumPratr = sumPratr+(P(9,k)-P(7,k))*(deltaT-(sumt-Sgran));

			% Placing the averaged Pratr over the Sgran*Sratio interval 
			% in the most recent slot of the cardiopulmonary baroreflex 
			% input.
			bPratrtemp = [(sumPratr/Sgran); bPratrtemp(1:Sratio-1)];
			avebPratr = sum(bPratrtemp)/Sratio;
			bPratr = [avebPratr; bPratr(1:Slength-1)];

			% Calculating the new baroreflex adusted parameters 
			% for the next Sgran interval.
			thc = thc+cpreflex([Sgran; Slength; Sratio],bPratr,th,thsp);
		
		end

		% Implementing the direct neural coupling mechanism between
		% Qlu and F.
		if (dncm == 1)

                        % Creating the dncm filter.
			hilvhr = dncm_filt(th,[Sgran; Slength; Sratio]);
                        % Creating the input Qlu vector to be convolved with the filter.
			tempilv = ilv(rsastep:rsastep+length(hilvhr)-1);
			tempilv = tempilv-th(98);
                        % Compute the mandated change in F from the convolution.
			thc(1) = thc(1)+sum(tempilv.*hilvhr(length(hilvhr):-1:1));
			rsastep = rsastep+1;

		end

		% Implementing the exogenous disturbance to Ra.
		if (dra ~= 0)

                   % Bandlimited, white disturbance
		   if (dra == 2)
                        
                        % Creating the bandlimited filter.
		        hra = bl_filt(th(44),Sgran);
                        % Creating white noise input to be convolved with the filter.
		        wra = [sqrt(1/Sgran)*th(47)*randn(1); wra(1:length(wra)-1)];
                        % Computing the mandated change to Ra.
		        thc(2) = thc(2)+sum(wra.*hra');

                   % Other Ra disturbance.
		   else

                        % Assigning mandated change to Ra which was pre-computed.
			thc(2) = thc(2)+Ralvr(lvrstep);
			lvrstep = lvrstep+1;

		   end

		end

		% Implementing the 1/f disturbance to F.
		if (df == 1)

                        % Creating a filter which is a cascade combination of 
                        % an autonomic filter and 1/f filter.
			hans = ans_filt(th,[Sgran; Slength; Sratio]);
			hf = conv(honeoverf,hans);
                        % Creating white noise input to be convolved with the filter.
		        wf = [sqrt(1/Sgran)*th(48)*randn(1); wf(1:length(wf)-1)];
                        % Computing mandated change to F.
			thc(1) = thc(1)+sum(wf.*hf);

		end

		% Implementing disturbance to Qvo.
		if (dqvo == 1)

                        % Assigning mandated change to Qvo which was pre-computed.
			thc(5) = thc(5)+Vvo(nvvostep);
			nvvostep = nvvostep+1;

		end

		% Calculating the values of the adjustable parameters at the 
                % end of the next Sgran interval and assigning these values
                % to the thn vector.
		thn = th;		
		thn(1) = thsp(4)+thc(3);
		thn(5) = thsp(5)+thc(4);
		thn(11) = thsp(6)+thc(5);
		thn(16) = thsp(3)+thc(2);
		thn(22) = thsp(2)+thc(1);

		% Denoting the time at the start of the next Sgran interval.
		tn = t(k)-(sumt-Sgran);	

		% Resetting the initial values for the integration 
		% of Pa, Pratr, and t over the next Sgran interval.
		sumP = P(2,k)*(sumt-Sgran);
		sumPratr = (P(9,k)-P(7,k))*(sumt-Sgran);
		sumt = sumt-Sgran;

	end

	% Calculating the new values of the adjustable parameters
	% at the current time step via linear interpolation and assigning
        % the values to the vector tha.  (Note that the mandated changes
        % to Cls and Crs are implemented at the start of each beat.)
	tha = th;
	tha(11) = ((thn(11)-thp(11))/Sgran)*(t(k)-(tp+Sgran))+thp(11);
	tha(16) = ((thn(16)-thp(16))/Sgran)*(t(k)-(tp+Sgran))+thp(16);
	tha(22) = ((thn(22)-thp(22))/Sgran)*(t(k)-(tp+Sgran))+thp(22);

	% Adjusting Pv in response to any mandated change in Qvo -- for 
        % constant volume considerations.
	P(3,k) = P(3,k)+(th(11)-tha(11))/th(4);
	
	% Setting the th vector to its current values.
	th = tha;

	% Computing Pra and qpv or qv over each cardiac cycle in order
	% to establish cardiac function/venous return curves.
        % (Note that ql is not computed, because it is relatively
        % inaccurate for large integration steps.)
	if (preparation == 1)

		summrap = summrap+0.5*deltaT*(P(9,k-1)+P(9,k));
		x = (P(6,k) > P(1,k));
		mbintscalar(x);
		if (x)
		        qpv = (P(6,k)-P(1,k))/th(20);
		else
		        qpv = 0;
		end
		summco = summco+0.5*deltaT*(q(1,k-1)+qpv);
		
	elseif (preparation == 2)

		summrap = summrap+0.5*deltaT*(P(9,k-1)+P(9,k));
		x = (P(3,k) > P(4,k) & P(4,k) >= th(91));
		y = (P(3,k) > th(91) & th(91) > P(4,k));
		mbintscalar(x);
		mbintscalar(y);
		if (x)
			qv = (P(3,k)-P(4,k))/th(17);
		elseif (y)
			qv = (P(3,k)-th(91))/th(17);
		else
			qv = 0;
		end
		summco = summco+0.5*deltaT*(q(4,k-1)+qv);
	end

	% Recalculating the fraction advance in the cardiac cycle to account 
	% for the mandated change to F. 
	frac = ((deltaT-(sumt-Sgran))*th(22)) + ((sumt-Sgran)*tha(22));

	% Updating time in and fraction of current cardiac cycle.
	% These values are associated with the newly calculated pressure
	% and time data.  The integration of the ipfm parameter is achieved
	% with the Backward Euler method.
	
           sumdeltaT = sumdeltaT + deltaT;

	   % Stopping the initiation of ventricular contractions after a 
           % desired number of beats.
	   if (preparation == 4 & qrs_index > th(43))

		sumdelta = 0;
		frac = 0;

	   end
	   ipfm = ipfm + frac;

	% Completing a cardiac cycle when ipfm is greater than or equal to one.
	ipfmcond = (ipfm >= 1);
	mbintscalar(ipfmcond)
	if (ipfmcond)
	
		% Resetting the ipfm model once the current cardiac contraction
		% is complete.
		sumdeltaT = (ipfm - 1)/th(22);
		ipfm = ipfm - 1;
		qrs_index = qrs_index+1; % Number of beats.
		ap(7,qrs_index) = t(k)-sumdeltaT; % Time of each beat.
		ap(8,qrs_index) = k; % Sample number of each beat.

		% Setting next systolic interval based on previous
		% cardiac cycle length.
		th(24) = t(k)-sumdeltaT-ap(7,qrs_index-1);

		% Implementing mandated changes to Cls and Crs.
			th(1) = thn(1);
		        th(5) = thn(5);

                % Varying parameters and computing averages for generation 
                % of cardiac function/venous return curves
		if (preparation == 1)

			% Calculating the integral of Pra and qpv
			% precisely over the previous cardiac cycle. 	
			summrap = summrap-0.5*deltaT*(P(9,k-1)+P(9,k));
			summrap = summrap+0.5*(deltaT-sumdeltaT)*(P(9,k-1)+P(9,k));
			summco = summco-0.5*deltaT*(q(1,k-1)+qpv);
			summco = summco+0.5*(deltaT-sumdeltaT)*(q(1,k-1)+qpv);

			% Computing the average Pra and qpv over the
			% previous cardiac cycle. 			.
			mrap = summrap/(t(k)-sumdeltaT-ap(7,qrs_index-1));
			mco = summco/(t(k)-sumdeltaT-ap(7,qrs_index-1));

			% Recording mean Pra, Pa, and qpv to the num matrix
			% and varying Pv and/or Pa every fb_int beats until
                        % all of the pre-determined variations have been 
                        % implemented.
			if (qrs_index == fb_int & countn < lengtht)
			       countn = countn+1;
			       countj = countj+1;
			       num(1,countn) = mrap;
			       num(2,countn) = P(2,k);
			       num(3,countn) = mco;
			       if (countj ~= lengthra)
				  th(30) = rarange(countj+1);
			       end
			       if (countj == lengthra & countn ~= lengtht)
			          th(29) = th(29)+th(90);
				  P(2,k) = th(29);
				  % Adjusting range of Pv for large Pa values.
				  if (th(29) > 180) 
				     tmp = tmp+1;
				     final = final-0.5*tmp;
				     rarangei = final:-th(88):0;
				     while (length(rarangei) < lengthra-1)
					   rarangei = [rarangei, rarangei(length(rarangei))-th(88)];
				     end
				     rarange = [rarangei, th(23)+1];
				     if (lengthra == 2);
					lengthra = 1;
					rarange = final;
				     end
				     th(30) = final;
				  end
				  fb_int = fb_int+5;
				  countj = 0;
				  P(3,k) = rarange(countj+1);
			       elseif (countn ~= lengtht)
				  P(3,k) = th(30);
				  fb_int = fb_int+5;
			       end	  
			end

			% Resetting the initial values for the integration 
			% of Pra and qpv over the next cardiac cycle.
			summrap = 0.5*sumdeltaT*(P(9,k-1)+P(9,k));
			summco = 0.5*sumdeltaT*(q(1,k-1)+qpv);
			
		end

		if (preparation == 2)

			% Calculating the integral of Pra and qv
			% precisely over the previous cardiac cycle. 	
			summrap = summrap-0.5*deltaT*(P(9,k-1)+P(9,k));
			summrap = summrap+0.5*(deltaT-sumdeltaT)*(P(9,k-1)+P(9,k));
			summco = summco-0.5*deltaT*(q(4,k-1)+qv);
			summco = summco+0.5*(deltaT-sumdeltaT)*(q(4,k-1)+qv);

			% Computing mean Pra and ql over the previous cardiac cycle.
			mrap = summrap/(t(k)-sumdeltaT-ap(7,qrs_index-1));
			mco = summco/(t(k)-sumdeltaT-ap(7,qrs_index-1));

			% Recording mean Pra and ql to the num matrix
			% and varying Crd after every fb_int beats 
                        % until all the pre-determined variations have
                        % been implemented.
			if (qrs_index == fb_int & Crd <= 60 & lengthra ~= 1)
			   countn = countn+1;		   
			   num(1,countn)=mrap;
			   num(2,countn)=mco;
			   Crdo = Crd;
			   Crd = Crd+th(89);

			   % Adjustment of Pr in response to variation in Crd -- for
                           % constant volume considerations.
			   th(6) = (th(32)/(th(27)-th(12)))*Crd;
			   ypres = (P(4,k)-P(7,k))/th(32);
			   xvol = vent_vol(Q(4,k-1),ypres,1/Crdo);
			   alpha = 1/th(6);
			   x0 = 1/(alpha+1);
			   kscale = 50;
			   c1 = 1+exp(-kscale*(1-x0));
			   d1 = 1+exp(kscale*x0);
			   denb1 = (1+(1/kscale)*log(c1/d1));
			   beta = (1-alpha)/denb1;
		  	   c = 1+exp(-kscale*(xvol-x0));
			   d = 1+exp(kscale*x0);
			   denb = (xvol+(1/kscale)*log(c/d));
			   P(4,k) = th(32)*(alpha*xvol+beta*denb)+P(7,k);

			   % Incrementing fb_int.
			   fb_int = fb_int+5;

			end

			% Resetting the initial values for the integration 
			% of Pra and qv over the next cardiac cycle.
			summrap = 0.5*sumdeltaT*(P(9,k-1)+P(9,k));
			summco = 0.5*sumdeltaT*(q(4,k-1)+qv);

		end


	end

	% Recording the adjustable parameter values at the current integration step.
	ap(1:6,k) = [th(1)*((th(26)-th(9))/th(31)) th(5)*((th(27)-th(12))/th(32)) th(11) th(16) th(22) .3*sqrt(th(24))]';

	% Recording the variable elastance values at the current integration step.
	if (preparation ~= 10)
		[El,dEl] = var_cap(th(1),th(2),th(24),sumdeltaT);
		[Er,dEr] = var_cap(th(5),th(6),th(24),sumdeltaT);
		ve(1:2,k) = [El*(th(31)/(th(26)-th(9))); Er*(th(32)/(th(27)-th(12)))];
	else
		[El,dEl] = var_vcap(th(1),th(2),th(24),sumdeltaT);
		[Er,dEr] = var_vcap(th(5),th(6),th(24),sumdeltaT);
		[Ela,dEla] = var_acap(th(80),th(81),th(24),sumdeltaT);
		[Era,dEra] = var_acap(th(84),th(85),th(24),sumdeltaT);
		ve(:,k) = [El*(th(31)/(th(26)-th(9))); Er*(th(32)/(th(27)-th(12))); Ela; Era];
	end		

	% Obtaining the volume values at the current integration step.  Equations depend
	% on the preparation being implemented.
	if (preparation ~= 10)
		Q(1:6,k) = [(1/El) th(3) th(4) (1/Er) th(7) th(8)]'.*(P(1:6,k)-P(7,k)*[1 0 0 1 1 1]') + [th(9:14)];
	else
		Q(1:6,k) = [(1/El) th(3) th(4) (1/Er) th(7) th(8)]'.*([P(1:6,k)]-P(7,k)*[1 0 0 1 1 1]')+[th(9:14)];
		Q(8:9,k) = [(1/Era) (1/Ela)]'.*([P(10:11,k)]-P(7,k)*[1 1]')+[th(86); th(82)];
	end

	if (preparation == 6)
		Q(2,k) = th(72)*P(2,k)+th(74);
		Q(8,k) = th(73)*P(10,k)+th(75);
	elseif (preparation == 1)
		Q(2,k) = 0;
		Q(3,k) = 0;
	elseif (preparation == 3)
		Q(2,k) = 0;
		Q(6,k) = 0;
	elseif (preparation == 7)
		alphatau = -th(70)/th(3);
		Qamax = th(71)-th(3)^2/th(70);
		Qao = th(71)+(th(3)^2/th(70))*exp(-th(70)*th(40)/th(3))*(1-exp(th(70)*th(40)/th(3)));
		Q(2,k) = (Qamax-Qao)*(1-exp(-alphatau*P(2,k)))+Qao;
	elseif (preparation == 8)
		alphatau = -th(70)/th(72);
		Qamax = th(76)-th(72)^2/th(70);
		Qao = th(76)+(th(72)^2/th(70))*exp(-th(70)*th(40)/th(72))*(1-exp(th(70)*th(40)/th(72)));
		Q(2,k) = (Qamax-Qao)*(1-exp(-alphatau*P(2,k)))+Qao;
		Q(8,k) = th(73)*P(10,k)+th(75);
	end

	if (preparation ~= 5)
		x = (P(1,k) < P(7,k));
		mbintscalar(x);
		if (x)
			Q(1,k) = th(9)*(2/pi)*atan(((P(1,k)-P(7,k))/((2/pi)*th(9)*El))) + th(9);
		else
			yl = (P(1,k)-P(7,k))/th(31);
			xl = vent_vol(Q(1,k-1),yl,El);
			Q(1,k) = (th(26)-th(9))*xl+th(9);
		end

		x = (P(4,k) < P(7,k));
		mbintscalar(x);
		if (x)
			Q(4,k) = th(12)*(2/pi)*atan(((P(4,k)-P(7,k))/((2/pi)*th(12)*Er))) + th(12);
		else
			yr = (P(4,k)-P(7,k))/th(32);
			xr = vent_vol(Q(4,k-1),yr,Er);
			Q(4,k) = (th(27)-th(12))*xr+th(12);
		end
	end

	% Correcting Qv and Pv in order to keep total blood volume in the
        % intact circulatory preparations, which may vary due to integration
	% error, constant.
	if (preparation == 0 | preparation == 4 | preparation == 5 | preparation == 6 | preparation == 7 | preparation == 8 | preparation == 9 | preparation == 10)

		if (preparation == 6 | preparation == 8)
		   Qtot = sum(Q([1:6 8],k));
		elseif (preparation == 10)
		   Qtot = sum(Q([1:6 8 9],k));
		else
		   Qtot = sum(Q(1:6,k));
		end
		Q(3,k) = Q(3,k) - (Qtot-th(21));
		P(3,k) = (Q(3,k)-th(11))/th(4);

	end

	% Checking for a request to pause simulation after the completion
	% of every beat provided that waveforms are being viewed on-line.
	if (strcmp(signals,'-1') == 0 & ipfmcond)

	   flagnew = read_key;

	   if (flagnew == 1)

	       % Re-reading parameter file, in the event of a possible update.
	       thspold = thsp;
	       thold = th;
	       fidp = fopen(parameterfile,'r');
	       [th,signals] = read_param(fidp);
	       fclose(fidp);

	       % Precluding adjustment of unchangeable parameters.
	       th(98) = thold(98); % Qfr
	       th(33) = thold(33); % Pms
	       th(88) = thold(88); % Pvs
	       th(89) = thold(89); % Crds
	       th(90) = thold(90); % Pas	       

	       % Precluding adjustment of Pv, Pa, and Crd if they are being
	       % automatically varied.
	       if (preparation == 1 & lengthra > 1)
		  th(30) = thold(30);
	       end
	       if (preparation == 1 & lengthsa > 1)
		  th(29) = thold(29);
	       end
	       if (preparation == 2 & lengthra > 1)
		  th(6) = thold(6);
	       end

	       % Setting annotator flag. 
	       if (th(105) == 0)
		  annotator = 'wqrs';
	       else
		  annotator = 'qrs';
	       end
	       
	       % Determining if parameter update is relevant to
	       % status of current simulation.
	       parameter_change = param_change(th,thold,flag);	       

	       mbintscalar(parameter_change == 1);
	       if (parameter_change == 1)

	          % Annotating and documenting the relevant parameter update.
		  count_pc = count_pc+1;
		  ap(9,qrs_index) = count_pc;
		  write_param(parameterfile,num2str(count_pc));

		  % Adjusting setpoint pressures.
		  thsp(1) = th(40);
		  thsp(7) = th(41);

		  % Conserving Qfr by altering Pth at the functional reserve
		  % volume of the lungs.
	          th(23) = th(25)-((th(98)-th(102))/th(101));

		  % Updating respiratory-related pressures and volumes
		  % over the remainder of the integration period, if necessary.

		  % No breathing.
		  if (breathing == 0 & th(23) ~= thold(23))

		     thbreathe = resp_act(th,at,deltaT/2,0,(tics-1)*deltaT);

		  % Fixed-rate breathing.
		  elseif (breathing == 1 & (th(42) ~= thold(42) | th(77) ~= thold(77) | th(23) ~= thold(23) | th(100) ~= thold(100)))

		     atcount = 2;
		     x = (at(atcount) < t(k)+1.5);
		     mbintscalar(x);
		     while (x)
			atcount=atcount+1;
			x = (at(atcount) < t(k)+1.5);
			mbintscalar(x);
		     end

		     start_index = k+ceil((at(atcount)-t(k))/deltaT);
		     start_time = t(k)+deltaT*ceil((at(atcount)-t(k))/deltaT)-at(atcount);
		     end_time = (length(thbreathe)-(2*(start_index-1)+1))*(deltaT/2)+start_time;

		     if (start_index < tics)

			at = [at(1:atcount); zeros(length(start_index:tics)+1,1)];      
			x = (at(atcount) <= deltaT*tics+1);
			mbintscalar(x);
			while (x)
			   atcount=atcount+1;
			   at(atcount) = at(atcount-1)+th(42);
			   x = (at(atcount) <= deltaT*tics+1);
			   mbintscalar(x);
		        end	
			at = at(1:atcount); 

			thbreathe(:,2*(start_index-1)+1:end) = resp_act(th,zeros(4,1),deltaT/2,start_time,end_time);

			ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran);
			ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];

		    end

                  % Random-interval breathing.
		  elseif (breathing == 2 & (th(42) ~= thold(42) | th(77) ~= thold(77) | th(99) ~= thold(99) | th(23) ~= thold(23) | th(100) ~= thold(100)))

		       atcount = 2;
		       x = (at(atcount) < t(k)+1.5);
		       mbintscalar(x);
		       while (x)
			  atcount=atcount+1;
			  x = (at(atcount) < t(k)+1.5);
			  mbintscalar(x);
		       end

		       start_index = k+ceil((at(atcount)-t(k))/deltaT);
		       start_time = t(k)+deltaT*ceil((at(atcount)-t(k))/deltaT)-at(atcount);
		       end_time = (length(thbreathe)-(2*(start_index-1)+1))*(deltaT/2)+start_time;

		       if (start_index < tics)

		          thbreathe(:,2*(start_index-1)+1:end) = resp_act(th,[1; (at(atcount+1:end)-at(atcount))],deltaT/2,start_time,end_time);
 		          ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran); 
			  ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];

		       end

		  end

		  % Implementing a step change to Qlu, if mandated. 
		  if (th(36) ~= 0)
   
		    if (breathing == 1 | breathing == 2)

			 atcount = 2;
			 x = (at(atcount) < t(k)+1.5);
			 mbintscalar(x);
			 while (x)
			    atcount=atcount+1;
			    x = (at(atcount) < t(k)+1.5);
			    mbintscalar(x);
		         end

			 start_index = k+ceil((at(atcount)-t(k))/deltaT);
			 start_time = (start_index-1)*deltaT;
			 th(103) = start_time+1.5;

		    elseif (breathing == 0)

			    th(103) = t(k)+1.5;

		    end			
		  
		    if (th(103) < (tics-1)*deltaT)

			 thbreathes = resp_act(th,[3; th(36)],deltaT/2,0,(tics-1)*deltaT);
			 dthbreathe = [-th(98)*ones(1,length(thbreathes)); -th(23)*ones(1,length(thbreathes)); zeros(2,length(thbreathes))]; 
			 thbreathe = thbreathe+thbreathes+dthbreathe;
			 ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran); 
			 ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];

		    end

                 end

                 % Assigning values to Pth, Ppac (Palv), and Qlu
		 P(7,k) = thbreathe(2,2*(k-1)+1);
		 P(8,k) = thbreathe(4,2*(k-1)+1);
		 Q(7,k) = thbreathe(1,2*(k-1)+1);

                 % Conserving volume in each compartment of the pulsatile
                 % heart and circulation through pressure adjustments.
		 P(1:6,k) = conserve_vol(th,sumdeltaT,Q(1:6,k),P(7,k));
		 if (preparation == 1)
		    P(2,k) = th(29);
		    P(3,k) = th(30);
		 end

		 % Updating adjustable parameters via instantaneous values
                 % and setpoint values.
		   if (thspold(3) ~= th(16))
		      thn(16) = th(16);
		      thp(16) = th(16);
		      thsp(3) = th(16);
		      ap(4,k) = th(16);
		   else
		      th(16) = thold(16);
		   end

		   if (thspold(2) ~= th(22))
		      thn(22) = th(22);
		      thp(22) = th(22);
		      thsp(2) = th(22);
		      ap(5,k) = th(22);	
		  else
		      th(22) = thold(22);
		  end

		  if (thspold(6) ~= th(11))
		     thn(11) = th(11);
		     thp(11) = th(11);
		     thsp(6) = th(11);
		     ap(3,k) = th(11);
		  else
		     th(11) = thold(11);
		     if (preparation == 0)
			P(3,k) = (Q(3,k)-th(11))/th(4);
		     end	
		  end

		  thn([1:10 12:15 17:21 23:end]) = th([1:10 12:15 17:21 23:end]);
		  thp([1:10 12:15 17:21 23:end]) = th([1:10 12:15 17:21 23:end]);


		  if (thspold(4) ~= th(1))
		     thsp(4) = th(1);
		  else
		     th(1) = thold(1);
		  end

		  if (thspold(5) ~= th(5))
		     thsp(5) = th(5);
		  else
		     th(5) = thold(5);
		  end

	       end

	       end

	end

	% Obtaining the flow rate values at the current integration step.  Equations
        % are dependent on the preparation being implemented.
	if (preparation == 0 | preparation == 1 | preparation == 2 | preparation == 4 | preparation == 5 | preparation == 6 | preparation == 7 | preparation == 8)
		x = (P(6,k) > P(1,k));
		mbintscalar(x);
		if (x)
		        q(1,k) = (P(6,k)-P(1,k))/th(20);
		else
		        q(1,k) = 0;
		end
	elseif (preparation == 3)
		x = (th(34) > P(1,k));
		mbintscalar(x);
		if (x)
		        q(1,k) = (th(34)-P(1,k))/th(20);
		else
		        q(1,k) = 0;
		end
	elseif (preparation == 10)
	        x = (P(11,k) > P(1,k));
		mbintscalar(x);
		if (x)
			q(1,k) = (P(11,k)-P(1,k))/th(83);
		else
		        q(1,k) = 0;
		end

	end
	
	if (preparation == 0 | preparation == 1 | preparation == 4 | preparation == 5 | preparation == 6 | preparation == 7 | preparation == 8 | preparation == 9 | preparation == 10) 
		x = (P(1,k) > P(2,k));
		mbintscalar(x);
		if (x)
		        q(2,k) = (P(1,k)-P(2,k))/th(15);
		else
		        q(2,k) = 0;	
		end
	elseif (preparation == 2)
		x = (P(1,k) > th(33));
		mbintscalar(x);
		if (x)
		        q(2,k) = (P(1,k)-th(33))/th(15);
		else
		        q(2,k) = 0;
		end
	elseif (preparation == 3)
		x = (P(1,k) > th(29));
		mbintscalar(x);
		if (x)
		        q(2,k) = (P(1,k)-th(29))/th(15);
		else
		        q(2,k) = 0;	
		end
	end

	if (preparation == 0 | preparation == 2 | preparation == 4 | preparation == 5 | preparation == 7 | preparation == 9 | preparation == 10)
		q(3,k) = (P(2,k)-P(3,k))/th(16);
	elseif (preparation == 1 | preparation == 3)
		q(3,k) = 0;
	elseif (preparation == 6 | preparation == 8)
		q(3,k) = (P(10,k)-P(3,k))/th(16);
	end

	if (preparation ~= 10 & preparation ~= 5)
		x = (P(3,k) > P(4,k) & P(4,k) >= th(91));
		y = (P(3,k) > th(91) & th(91) > P(4,k));
		mbintscalar(x);
		mbintscalar(y);
		if (x)
			q(4,k) = (P(3,k)-P(4,k))/th(17);
		elseif (y)
			q(4,k) = (P(3,k)-th(91))/th(17);
		else
			q(4,k) = 0;
		end
	elseif (preparation == 5)
	       x = (P(3,k) >= P(4,k));
	       mbintscalar(x);
	       if (x)
			q(4,k) = (P(3,k)-P(4,k))/th(17);
	       else
			q(4,k) = 0;
	       end
	else
		x = (P(9,k) >= P(4,k));
		mbintscalar(x);
		if (x)
			   q(4,k) = (P(9,k)-P(4,k))/th(87);
		else
			   q(4,k) = 0;
		end

	end

	if (preparation ~= 2)
		x = (P(4,k) > P(5,k));
		mbintscalar(x);
		if (x)
		        q(5,k) = (P(4,k)-P(5,k))/th(18);
		else
		        q(5,k) = 0;
		end
	elseif (preparation == 2)
		x = (P(4,k) > th(28));
		mbintscalar(x);
		if (x)
		        q(5,k) = (P(4,k)-th(28))/th(18);
		else
		        q(5,k) = 0;
		end
	end

	if (preparation == 2 | preparation == 5)
		q(6,k) = (P(5,k)-P(6,k))/th(19);
	elseif (preparation == 0 | preparation == 1 | preparation == 4 | preparation == 6 |preparation == 7 | preparation == 8 | preparation == 9 | preparation == 10)
		h1 = (1330/(1.06*980))*(P(6,k)-P(8,k));
		x = h1 < -10;
		mbintscalar(x);
		y = h1 > 20;
		mbintscalar(y);
		if (x)
			h1 = -10;
		elseif (y)	
			h1 = 20;
		end

		h2 = (1330/(1.06*980))*(P(5,k)-P(8,k));
		x = h2 < -10;
		mbintscalar(x);
		y = h2 > 20;
		mbintscalar(y);
		if (x)
			h2 = -10;
		elseif (y)	
			h2 = 20;
		end

		q(6,k) = (((P(5,k)-P(6,k))/(30*th(19)))*(h1+10)) + (((P(5,k)-P(8,k))/(30*th(19)))*(h2-h1)) - (((1.06*980)/(2660*(30*th(19))))*(h2^2-h1^2));

	elseif (preparation == 3)
		q(6,k) = 0;
	end

	if (preparation == 10)
		x = (P(3,k) >= th(91) & P(9,k) >= th(91));
		y = (P(3,k) > th(91) & th(91) > P(9,k));
		z = (P(9,k) > th(91) & th(91) > P(3,k));
		mbintscalar(x);
		mbintscalar(y);
		mbintscalar(z);
		if (x)
			 q(7,k) = (P(3,k)-P(9,k))/th(17);
		elseif (y)	 
			 q(7,k) = (P(3,k)-th(91))/th(17);
		elseif (z)	 
			 q(7,k) = (th(91)-P(9,k))/th(17);
		else
			 q(7,k) = 0;
		end
	elseif (preparation == 6 | preparation == 8)
		q(7,k) = P(11,k);
	end

        if (preparation == 10)
		q(8,k) = (P(6,k)-P(11,k))/th(20);
	end	


	% Writing waveforms and annotations to files in MIT format every step parameter
        % seconds once a window parameter duration of data have been calculated.
	if (strcmp(signals,'-1') == 0 & t(k) > th(104) & (t(k)-t(kstart)) > th(78))

           % Writing newly generated waveform data to MIT format files.
	   X = [P(1:9,kstart:k); Q(1:6,kstart:k); Q(7,kstart:k)*0.1; q(1:6,kstart:k); ap(1:2,kstart:k)*100; ap(3,kstart:k); ap(4,kstart:k)*100; ap(5,kstart:k)*600; ve(1:2,kstart:k)*10]*10;
	   fwrite(fid,X,'short');

           % Updating variables necessary for writing waveforms.
	   kstart = k+1;
	   start_time = t(k)-th(104);

           % Writing newly generated annotations to MIT format files.
	   qrso = [];
	   qrso = diff(ap(8,qrs_index_start:qrs_index))';
	   for i = 1:length(qrso)
	       fwrite(fid2,qrso(i),'bit10');
	       fwrite(fid2,1,'ubit6');

	       % If parameter update has been made, writing auxillary information.
	       x = (ap(9,qrs_index_start+i) ~= 0);
	       mbintscalar(x)
	       if (x)
		  suffix = num2str(ap(9,qrs_index_start+i));
		  fwrite(fid2,Nbytes+(length(suffix)-1),'bit10');
		  fwrite(fid2,63,'ubit6');
		  fwrite(fid2,[parameterfile '.' suffix],[num2str(Nbytes+length(suffix)-1) '*uchar']);
	       end
	   end
	   qrs_index_start = qrs_index;

           % Alternating outputfile signaled to WAVE in order to permit
	   % on-line updating of the annotation flag.
	   if (wrflag == 0)

	      wave_remote( outputfile, annotator, start_time, signals );
	      wrflag = 1;

	   else
	      
	      wave_remote( outputfile2, annotator, start_time, signals );
	      wrflag = 0;

	   end

	end

	% If cardiac function/venous return curves are complete, there
	% is no need for further integration.  That is, for-loop is exited.
	if (preparation == 1 | preparation == 2)
	   if (countn >= length(num'))
	      tics = k;
	      break;
	   end
	end

end

% Writing remaining portion of simulated data to the MIT format files.
if (strcmp(signals,'-1') == 0)

        % Writing last portion of generated waveform data.
	X = [P(1:9,kstart:tics); Q(1:6,kstart:tics); Q(7,kstart:tics)*0.1; q(1:6,kstart:tics); ap(1:2,kstart:tics)*100; ap(3,kstart:tics); ap(4,kstart:tics)*100; ap(5,kstart:tics)*600; ve(1:2,kstart:tics)*10]*10;
	fwrite(fid,X,'short');
     
        % Writing last annotations.
	qrso = [];
	qrso = diff(ap(8,qrs_index_start:qrs_index))';
	for i = 1:length(qrso)

	    fwrite(fid2,qrso(i),'bit10');
	    fwrite(fid2,1,'ubit6');

	    % If parameter updates have been made, writing auxillary info. 
	    x = (ap(9,qrs_index_start+i) ~= 0);
	    mbintscalar(x)
	    if (x)

		suffix = num2str(ap(9,qrs_index_start+i));
		fwrite(fid2,Nbytes+(length(suffix)-1),'bit10');
		fwrite(fid2,63,'ubit6');
		fwrite(fid2,[parameterfile '.' suffix],[num2str(Nbytes+length(suffix)-1) '*uchar']);

	    end

	end
	fwrite(fid2,0,'bit10');
	fwrite(fid2,0,'ubit6');
	fclose(fid2);
	fclose(fid);

end

% If MATLAB execution, writing qrs vector as times of ventricular contractions.
% If Linux execution, writing qrs vector as samples of ventricular contractions.
if (nargin == 2)
   qrs = ap(7,1:qrs_index);
else
   qrs = ap(8,1:qrs_index);
end

% Only supplying as output the pressures that are generated by
% the model, for MATLAB execution only.
if (preparation == 6 | preparation == 8)

	P = P(1:10,:);
	q = q(1:7,:);
	Q = Q(1:8,:);
	ve = ve(1:2,:);

elseif (preparation == 10 | preparation == 13)	

        P = P([1:9 11],:);

else

	P = P(1:9,1:tics);
	q = q(1:6,1:tics);
	Q = Q(1:7,1:tics);
	ve = ve(1:2,1:tics);
	ap = ap(:,1:tics);
	t = t(1:tics);

end