function system = transdyn(Dn,sigma,v,varargin)
%TRANSDYN generates a transfer function approximation forseven unique combinations
%of the hyperbolic functions associated with the four possible causalities associated  
%with fluid transmission lines
% 
%    String names assigned to the hyperbolic functions are as follows:
%
%                'C1' = Zo*cosh(Gamma)/Z*sinh(Gamma)
%                'C2' = Zo/Z*sinh(Gamma)
%                'C3' = Zo*sinh(Gamma)/Z*cosh(Gamma)
%
%                'C4' = 1/cosh(Gamma)
%
%                'C5' = Z*cosh(Gamma)/Zo*sinh(Gamma)
%                'C6' = Z/Zo*sinh(Gamma)
%                'C7' = Z*sinh(Gamma)/Zo*cosh(Gamma)
%
%    The function arguments are,
%
%              SYSTEM = TRANSDYN(Dn,sigma,gamma,w,r,vis,'C*')
%
%Where,     "SYSTEM"     denotes "s-function" for TF object
%           "Dn"         Dimensionless dissipation number
%           "sigma"      Prandtl number
%           "gamma"      Specific heat ratio
%           "w"          The maximum desired frequency
%           "r"          Inside radius of the line
%           "vis"        Kinematic viscosity
%           'C*'         The string name of desired 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 TRANSDYN is to specify the number of second  
%order modes to be matched in the approximation.
%
%              SYSTEM = TRANSDYN(Dn,sigma,gamma,m,'C*')
%
%where,"m" is the desired number of second order modes.
%
%    Version 1.01, 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,'C1'))
   k=1;
elseif (strcmp(strng,'C2'))
   k=2;
elseif (strcmp(strng,'C3'))
   k=3;
elseif (strcmp(strng,'C4'))
   k=4;
elseif (strcmp(strng,'C5'))
   k=5;
elseif (strcmp(strng,'C6'))
   k=6;
elseif (strcmp(strng,'C7'))
   k=7;
else
   error('String variables need to be only C1 to C7')
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(n,1)=cosh(Gamma)/(Z*sinh(Gamma));
   end
   mag=20*log10(abs(C1));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C1,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));
      Exact(abs(s),1)=cosh(Gamma)/(Z*sinh(Gamma));
   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
      C2(n,1)=1/(Z*sinh(Gamma));
   end
   mag=20*log10(abs(C2));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C2,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));
      Exact(abs(s),1)=1/(Z*sinh(Gamma));
   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
      C3(n,1)=sinh(Gamma)/(Z*cosh(Gamma));
   end
   mag=20*log10(abs(C3));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C3,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));
      Exact(abs(s),1)=sinh(Gamma)/(Z*cosh(Gamma));
   end
case 4  
   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
      C4(n,1)=1/cosh(Gamma);
   end
   mag=20*log10(abs(C4));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C4,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));
      Exact(abs(s),1)=1/cosh(Gamma);
   end  
case 5  
   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
      C5(n,1)=Z*cosh(Gamma)/sinh(Gamma);
   end
   mag=20*log10(abs(C5));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C5,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));
      Exact(abs(s),1)=Z*cosh(Gamma)/sinh(Gamma);
   end
case 6
   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
      C6(n,1)=Z/sinh(Gamma);
   end
   mag=20*log10(abs(C6));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C6,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));
      Exact(abs(s),1)=Z/sinh(Gamma);
   end
case 7
   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
      C7(n,1)=Z*sinh(Gamma)/cosh(Gamma);
   end
   mag=20*log10(abs(C7));
   na=3+2+2*(m-1);
   nd=4+2+2*(m-1);
   [Numer,Denom]=invfreqs(C7,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));
      Exact(abs(s),1)=Z*sinh(Gamma)/cosh(Gamma);
   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