function system = ftld(Dn,sigma,v,Ra,Rb,varargin)
%FTLD generates transfer function approximations for G1, G2 and G3 
%     for the following causality
%
%                       |ZoQa| |G1  -G2||Ps|
%                       |    |=|       ||  |
%                       |ZoQb| |G2  -G3||Pc|
%
%              SYSTEM = FTLD(Dn,sigma,gamma,Ra,Rb,w,r,vis,'G1')
%              SYSTEM = FTLD(Dn,sigma,gamma,Ra,Rb,w,r,vis,'G2')
%              SYSTEM = FTLD(Dn,sigma,gamma,Ra,Rb,w,r,vis,'G3')
%
%Where,     "SYSTEM"     denotes "s-function" for TF object
%           "Dn"         Dimensionless dissipation number
%           "sigma"      Prandtl number
%           "gamma"      Specific heat ratio
%           "Ra"         Impedance ratio 'Ra/Zo' at the left end of the line
%           "Rb"         Impedance ratio 'Zo/Rb' at the right end of the line
%           "w"          The maximum desired frequency, rad/sec
%           "r"          Inside radius of the line, meter
%           "vis"        Kinematic viscosity, m^2/sec
%  'G1','G2' and 'G3'    Strings denoting transfer functions
%
% ***Note that the transfer function numerator will be of order 2*m+3 and 
%the transfer function denominator of order 2*m+4. This means that the  
%transfer function consists of 2 first-order modes and 'm+1' second-order modes.
%   A second option for using FTLD is to specify the number of second  
%order modes to be matched in the approximation.
%
%              SYSTEM = FTLD(Dn,sigma,gamma,Ra,Rb,m,'G1')
%              SYSTEM = FTLD(Dn,sigma,gamma,Ra,Rb,m,'G2')
%              SYSTEM = FTLD(Dn,sigma,gamma,Ra,Rb,m,'G3')
%
%where,"m" is the desired number of second order modes.
%
%    Version 1.101, 08/30/01

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(varargin) == 4
   [input,r,viz,strng]=deal(varargin{:});
   w=input;
   wn=r^2*w/viz;
   if v == 1
      i=0.37832*Dn*wn^0.9825;
   else
      i=0.42612*Dn*wn^0.966;
   end
   dummy=int32(i);
   new_i=double(dummy);
   if new_i == 0
      i=1;
   elseif i-new_i >= 0.5
      i=new_i+1;
   end
   m=i;
   if v == 1
      wn=(m/(0.37832*Dn))^(1/0.9825);
   else
      wn=(m/(0.42612*Dn))^(1/0.966);
   end
elseif length(varargin) == 2
   [input,strng]=deal(varargin{:});
   m=input;
   if v == 1
      wn=(m/(0.37832*Dn))^(1/0.9825);
   else
      wn=(m/(0.42612*Dn))^(1/0.966);
   end   
else
   error('Not enough arguments!! Please read the help statement')
end

%     selecting the transfer functions to perform curve-fitting
if (strcmp(strng,'G1'))
   k=1;
elseif (strcmp(strng,'G2'))
   k=2;
elseif (strcmp(strng,'G3'))
   k=3;
else
   error('String variables need to be G1 G2 or G3')
end

e=wn;
if e <= 11
   Eomega=[1:100];
elseif e <= 101
   Eomega=[1:1000];
elseif e <= 1001
   Eomega=[1:10000];
else
   Eomega=[10:10000];
end
% select set of sampling frequencies
if e <=11
   omega=[(1:0.005:e)];
elseif e <=101
   omega=[(1:0.01:10),(10.01:0.01:10.99),(11:0.1:e)];
elseif e <=1001
   omega=[(1:0.01:10),(10.01:0.01:10.99),(11:0.1:100),(100.01:0.1:100.9),...
         (101:0.5:e)];
else
   omega=[(1:0.01:10),(10.01:0.01:10.99),(11:0.1:100),(100.01:0.1:100.9),...
         (101:0.5:1000),(1000.01:0.01:1000.09),(1001:1:e)];
end
n=0;
switch k
case 1  
   for s=j*omega;
      n=n+1;
      Bsigma=2*besselj(1,j*sqrt(sigma*s))/(j*sqrt(sigma*s)...
         *besselj(0,j*sqrt(sigma*s)));
      B=2*besselj(1,j*sqrt(s))/(j*sqrt(s)*besselj(0,j*sqrt(s)));
      Gamma=Dn*s*sqrt((1+(v-1)*Bsigma)/(1-B));
      Z=1/(sqrt(1-B)*sqrt(1+(v-1)*Bsigma)); % already in term Z/Zo
      C1=cosh(Gamma)/(Z*sinh(Gamma));
      Gone(n,1)=(Rb*C1+(1/Z)^2)/(Rb+C1*(Rb*Ra+1)+Ra*(1/Z)^2);
   end
   mag=20*log10(abs(Gone));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(Gone,omega',na,nd,[],100);
   system=tf(Numer,Denom);
   Modal=freqresp(system,j*Eomega);
   Modal=Modal(:);
   Mmag=20*log10(abs(Modal));
   Momega=Eomega;
   for s=j*Eomega;
      Bsigma=2*besselj(1,j*sqrt(sigma*s))/(j*sqrt(sigma*s)...
         *besselj(0,j*sqrt(sigma*s)));
      B=2*besselj(1,j*sqrt(s))/(j*sqrt(s)*besselj(0,j*sqrt(s)));
      Gamma=Dn*s*sqrt((1+(v-1)*Bsigma)/(1-B));
      Z=1/(sqrt(1-B)*sqrt(1+(v-1)*Bsigma));
      C1=cosh(Gamma)/(Z*sinh(Gamma));
      Exact(abs(s),1)=(Rb*C1+(1/Z)^2)/(Rb+C1*(Rb*Ra+1)+Ra*(1/Z)^2);
   end
case 2
   for s=j*omega;
      n=n+1;
      Bsigma=2*besselj(1,j*sqrt(sigma*s))/(j*sqrt(sigma*s)...
         *besselj(0,j*sqrt(sigma*s)));
      B=2*besselj(1,j*sqrt(s))/(j*sqrt(s)*besselj(0,j*sqrt(s)));
      Gamma=Dn*s*sqrt((1+(v-1)*Bsigma)/(1-B));
      Z=1/(sqrt(1-B)*sqrt(1+(v-1)*Bsigma)); % already in term Z/Zo
      C1=cosh(Gamma)/(Z*sinh(Gamma));
      C2=1/(Z*sinh(Gamma));
      Gtwo(n,1)=Rb*C2/(Rb+C1*(Rb*Ra+1)+Ra*(1/Z)^2); 
   end
   mag=20*log10(abs(Gtwo));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(Gtwo,omega',na,nd,[],100);
   system=tf(Numer,Denom);
   Modal=freqresp(system,j*Eomega);
   Modal=Modal(:);
   Mmag=20*log10(abs(Modal));
   Momega=Eomega;
   for s=j*Eomega;
      Bsigma=2*besselj(1,j*sqrt(sigma*s))/(j*sqrt(sigma*s)...
         *besselj(0,j*sqrt(sigma*s)));
      B=2*besselj(1,j*sqrt(s))/(j*sqrt(s)*besselj(0,j*sqrt(s)));
      Gamma=Dn*s*sqrt((1+(v-1)*Bsigma)/(1-B));
      Z=1/(sqrt(1-B)*sqrt(1+(v-1)*Bsigma));
      C1=cosh(Gamma)/(Z*sinh(Gamma));
      C2=1/(Z*sinh(Gamma));
      Exact(abs(s),1)=Rb*C2/(Rb+C1*(Rb*Ra+1)+Ra*(1/Z)^2);
   end
case 3
   for s=j*omega;
      n=n+1;
      Bsigma=2*besselj(1,j*sqrt(sigma*s))/(j*sqrt(sigma*s)...
         *besselj(0,j*sqrt(sigma*s)));
      B=2*besselj(1,j*sqrt(s))/(j*sqrt(s)*besselj(0,j*sqrt(s)));
      Gamma=Dn*s*sqrt((1+(v-1)*Bsigma)/(1-B));
      Z=1/(sqrt(1-B)*sqrt(1+(v-1)*Bsigma)); % already in term Z/Zo
      C1=cosh(Gamma)/(Z*sinh(Gamma));
      Gthree(n,1)=((Rb*C1)+Ra*Rb*(1/Z)^2)/(Rb+C1*(Rb*Ra+1)+Ra*(1/Z)^2);
   end
   mag=20*log10(abs(Gthree));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(Gthree,omega',na,nd,[],100);
   system=tf(Numer,Denom);
   Modal=freqresp(system,j*Eomega);
   Modal=Modal(:);
   Mmag=20*log10(abs(Modal));
   Momega=Eomega;
   for s=j*Eomega;
      Bsigma=2*besselj(1,j*sqrt(sigma*s))/(j*sqrt(sigma*s)...
         *besselj(0,j*sqrt(sigma*s)));
      B=2*besselj(1,j*sqrt(s))/(j*sqrt(s)*besselj(0,j*sqrt(s)));
      Gamma=Dn*s*sqrt((1+(v-1)*Bsigma)/(1-B));
      Z=1/(sqrt(1-B)*sqrt(1+(v-1)*Bsigma));
      C1=cosh(Gamma)/(Z*sinh(Gamma));
      Exact(abs(s),1)=((Rb*C1)+Ra*Rb*(1/Z)^2)/...
         (Rb+C1*(Rb*Ra+1)+Ra*(1/Z)^2);
   end
end

% plotting graphs
if e > 1000
   Emag=20*log10(abs(Exact(10:10000)));
else
   Emag=20*log10(abs(Exact));
end
semilogx(Momega,Mmag,'r',Eomega,Emag);grid;
title('Magnitude plot');
xlabel('normalized frequency');ylabel('dB');
figure
if e > 1000
   AngleExact=(90/pi)*angle(Exact(10:10000));
else
   AngleExact=(90/pi)*angle(Exact);
end
semilogx(Momega,AngleExact,Eomega,(90/pi)*angle(Modal),'r');grid;
title('Phase plot');
xlabel('normalized frequency');ylabel('degree');
% end of code