%  Program "nyq_resp.m": Stability analysis of the chemoreflex
%    control of respiration in normals vs congestive heart failure
%     using Nyquist plots

% Common Parameter values
Vlung=2.5; Kco2=0.0065; Gp=0.02; Gc=0.04;
taup=20; tauc=120; 
Vdote=0.12; Vdotd=0.03; PaCO2=40;PICO2=0;
f = [0.01:0.001:0.1]';
w  = 2*pi*f;

% Construct Loop Transfer Function for normal(N) human 
QN=0.1; TpN=6.1; TcN=7.1;
GlungN=(PaCO2-PICO2)/(Vdote-Vdotd+863*QN*Kco2);
taulungN=Vlung/(Vdote-Vdotd+863*QN*Kco2);
num1N=[GlungN*Gp];
den1N=[taulungN*taup  (taulungN+taup)  1];
Hs1N = tf(num1N, den1N);
num2N=[GlungN*Gc];
den2N=[(taulungN*tauc)  (taulungN+tauc)  1];
Hs2N = tf(num2N, den2N);
% Compute Nyquist plot without delay
[R1N,I1N] = nyquist(Hs1N,w); R1N=squeeze(R1N); I1N=squeeze(I1N);
[R2N,I2N] = nyquist(Hs2N,w); R2N=squeeze(R2N); I2N=squeeze(I2N);
% Add delay to results
R1delN=real((R1N+j*I1N).*exp(-j*w*TpN));
I1delN = imag((R1N+j*I1N).*exp(-j*w*TpN));
R2delN=real((R2N+j*I2N).*exp(-j*w*TcN));
I2delN=imag((R2N+j*I2N).*exp(-j*w*TcN));
RdelN=R1delN+R2delN;
IdelN=I1delN+I2delN;
% Plot Nyquist diagram of overall frequency response
input(' Display Nyquist Plot for Normal Human: Hit <Enter> to continue >>')
close all; figure;
axis square; plot(RdelN,IdelN); grid;
title(' Normal') 

% Construct Loop Transfer Function for Congestive Heart Failure (C)  
QC=0.05; TpC=12.2; TcC=14.2;
GlungC=(PaCO2-PICO2)/(Vdote-Vdotd+863*QC*Kco2);
taulungC=Vlung/(Vdote-Vdotd+863*QC*Kco2);
num1C=[GlungC*Gp];
den1C=[taulungC*taup  (taulungC+taup)  1];
Hs1C = tf(num1C, den1C);
num2C=[GlungC*Gc];
den2C=[(taulungC*tauc)  (taulungC+tauc)  1];
Hs2C = tf(num2C, den2C);
% Compute Nyquist plot without delay
[R1C,I1C] = nyquist(Hs1C,w); R1C=squeeze(R1C); I1C=squeeze(I1C);
[R2C,I2C] = nyquist(Hs2C,w); R2C=squeeze(R2C); I2C=squeeze(I2C);
% Add delay to results
R1delC=real((R1C+j*I1C).*exp(-j*w*TpC));
I1delC=imag((R1C+j*I1C).*exp(-j*w*TpC));
R2delC=real((R2C+j*I2C).*exp(-j*w*TcC));
I2delC=imag((R2C+j*I2C).*exp(-j*w*TcC));
RdelC=R1delC+R2delC;
IdelC=I1delC+I2delC;
% Plot Nyquist diagram of overall frequency response
input('Display Nyquist Plot for CHF: Hit <Enter> to continue >>')
figure;
axis square; plot(RdelC,IdelC); grid;
title(' Congestive Heart Failure')