function [Z0,gamma,R0]=getTwoPortModel(cable,f)
%getTwoPortModel - get the two port model [Z0,gamma] for cable 'name'
%
% Parameter:    cable         Cable model
% Parameter:    f	      frequency axis
% Returns:      Z0	      Characteristic Impedance [Ohm]
% Returns:      gamma	      Propagation constant [S/km]
% Returns:      R0 (optional) Cable parameter R [Ohms/km] at f==0
%
% Example(s):
%   cable = getList(ex.clist,'DTAG-40');
%   [Z0,G] = getTwoPortModel(cable,ex.param.frequency.f); 
%
% Reference:
%   ETSI VDSL TM6(97)02 "Cable Reference models for simulating
%   metallic access networks" (970p02r3, version 01-07-98)
%   Musson, J. "Maximum likelihood estimation of the primary parameters of
%   twisted pair cables", ETSI Contribution 981t08a1, Madrid, 1998
%   W. Henkel, et al. "Update of cable reference models", 
%   ETSI Contribution 985t04c0, Sopha Antipolis, France, Nov 24-27, 1998
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1999 by Telia Research AB, Lulea, Sweden;                
%       2000-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Gernot Schmid (Gernot.Schmid@arcs.ac.at)
%		: Petr Kadlec (kadlec@ftw.at)
%
% CVS:       $Id: getTwoPortModel.m 752 2009-01-02 13:03:52Z tono $
%% ===========================================================================
% Change History
%      1999-08-11 (DaB)  Created
%      2000-03-14 (ToNo) Changed to reflect the new list structure of cables 
%      2000-03-16 (GS)   Added model for SDSL Testloop Cables RLC interp.
%      2000-04-03 (UvAn) Changed the interface; now the cable is
%                        passed instead of name
%      2000-05-04 (ToNo) Cleaned up the RLC section
%      2000-06-26 (ToNo) Use linear interpolation for the RLC parameters
%      2000-07-18 (ToNo) Corrected the tan_fi values for KPN cables
%      2000-09-13 (ToNo) Added the MAR model
%      2001-12-15 (ToNo) Noted the problem with matlab interpolation
%      2001-12-20 (ToNo) Added the modified DTAG model, cf 985t04c0
%      2002-01-30 (ToNo) Added the MAR2 model
%      2002-07-02 (Peka) Octave to Matlab port
%      2002-07-30 (ToNo) Speed things up by reducing duplicated calculations
%      2002-12-27 (ToNo) Workaround for a bug in matlab 6.5 in the
%                        BT model, idea by David Ferguson (videonetworks.com)
%      2004-01-11 (R. Stolle, Infineon) Returns cable parameter R for f==0,
%                        useful for impulse response and continuity purposes 
%      2005-10-31 (ToNo) Added the Swisscom model (Not Fully Functional Yet!)
%% ===========================================================================

if isempty(cable)
   estr=sprintf('Unknown cable model %s', name);
   error(estr);       
end;

switch cable.model
   case 'DTAGorg'
   Ka1=cable.param.Ka1; Ka2=cable.param.Ka2; Ka3=cable.param.Ka3; 
   Kb1=cable.param.Kb1; Kb2=cable.param.Kb2;
   Kz1=cable.param.Kz1; Kz2=cable.param.Kz2; Kz3=cable.param.Kz3;
   Kx1=cable.param.Kx1; Kx2=cable.param.Kx2; Kx3=cable.param.Kx3;
   freq=f;
   gamma=[];
   start=1;
   stop=max(find(freq<0.5e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(1); K2=Ka2(1); K3=Ka3(1);
   gamma=(K1+K2.*(f.^K3))/20*log(10)+j*(Kb1.*f+Kb2.*sqrt(f));
   start=stop+1;
   stop=max(find(freq<5e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(2); K2=Ka2(2); K3=Ka3(2);
   x=(K1+K2.*(f.^K3))/20*log(10)+j*(Kb1.*f+Kb2*sqrt(f));
   gamma=[gamma x];
   start=stop+1;
   stop=max(find(freq<=30e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(3); K2=Ka2(3); K3=Ka3(3);
   x=(K1+K2.*(f.^K3))/20*log(10)+j*(Kb1.*f+Kb2*sqrt(f));
   gamma=[gamma x];
   f=freq;
   Z0=(Kz1+Kz2./((f/1e6).^Kz3)).*exp(-j*Kx1./((Kx2+f/1e6).^Kx3));
   if size(Z0)==size(gamma),
   else
      fprintf('function not defined over 30MHz');
   end;
   f=75e-3; % lowest frequency/MHz where model is said to be valid
   g=(K1+K2.*(f.^K3))/20*log(10)+j*(Kb1.*f+Kb2.*sqrt(f));
   z=(Kz1+Kz2./(f.^Kz3)).*exp(-j*Kx1./((Kx2+f).^Kx3));
   R0=real(g*z);

   case 'DTAG'
   Ka1=cable.param.Ka1; Ka2=cable.param.Ka2; Ka3=cable.param.Ka3; 
   Kb1=cable.param.Kb1; Kb2=cable.param.Kb2; Kb3=cable.param.Kb3;
   Kz1=cable.param.Kz1; Kz2=cable.param.Kz2; Kz3=cable.param.Kz3; 
   Kx1=cable.param.Kx1; Kx2=cable.param.Kx2; Kx3=cable.param.Kx3;
   freq=f;
   gamma=[];
   start=1;
   stop=max(find(freq<0.5e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(1); K2=Ka2(1); K3=Ka3(1);
   L1=Kb1(1); L2=Kb2(1); L3=Kb3(1);
   gamma=(K1+K2.*(f.^K3))/20*log(10)+j*(L1.*f+L2.*(f.^L3));
   start=stop+1;
   stop=max(find(freq<5e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(2); K2=Ka2(2); K3=Ka3(2);
   L1=Kb1(2); L2=Kb2(2); L3=Kb3(2);
   gamma=[gamma (K1+K2.*(f.^K3))/20*log(10)+j*(L1.*f+L2.*(f.^L3))];
   start=stop+1;
   stop=max(find(freq<=30e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(3); K2=Ka2(3); K3=Ka3(3);
   L1=Kb1(3); L2=Kb2(3); L3=Kb3(3);
   gamma=[gamma (K1+K2.*(f.^K3))/20*log(10)+j*(L1.*f+L2.*(f.^L3))];
   f=freq;
   Z0=(Kz1+Kz2./((f/1e6).^Kz3)).*exp(-j*Kx1./((Kx2+f/1e6).^Kx3));
   
   if size(Z0)==size(gamma),
   else
      fprintf('function not defined over 30MHz');
   end;
   f=75e-3; % lowest frequency/MHz where model is said to be valid
   g=(K1+K2.*(f.^K3))/20*log(10)+j*(L1.*f+L2.*(f.^L3));
   z=(Kz1+Kz2./(f.^Kz3)).*exp(-j*Kx1./((Kx2+f).^Kx3));
   R0=real(g*z);

case 'BT'
   r0c=cable.param.r0c; r0s=cable.param.r0s; 
   ac=cable.param.ac; as=cable.param.as; 
   l0=cable.param.l0; loo=cable.param.loo; b=cable.param.b;
   fm=cable.param.fm; 
   coo=cable.param.coo; c0=cable.param.c0; ce=cable.param.ce; 
   g0=cable.param.g0; ge=cable.param.ge;
   f2 = f.^2;
   if (as==0) && (r0s==0),
       % No steel core, fall back to the simplified model
       R=sqrt(sqrt(r0c.^4+ac.*f2)); % That is, r0s = inf;
       R0=r0c;
   else
       R=1./((1./sqrt(sqrt(r0c.^4+ac.*f2))) + (1./sqrt(sqrt(r0s.^4+ as.*f2))));
       R0=r0c*r0s/(r0c+r0s);
   end;
   ffmb=(f./fm).^b;
   L=(l0+loo*ffmb)./(1+ffmb);
   C=coo+c0.*f.^(-ce);
   G=g0.*f.^(ge);
   Z=R+j.*2*pi.*f.*L;
   Y=G+j*2*pi.*f.*C;
   gamma=sqrt(Z.*Y);
   Z0=sqrt(Z./Y);
    
case 'SWC'
   % Not Fully Functional Yet!
   Z00 = cable.param.Z00; f1 = cable.param.f1; f2 = cable.param.f2; 
   f3 = cable.param.f3; f4 = cable.param.f4; f5 = cable.param.f5;
   Ne1 = cable.param.Ne1; Ne2 = cable.param.Ne2;
   Ne3 = cable.param.Ne3; Ne4 = cable.param.Ne4;
   c1 = cable.param.C1; c2 = cable.param.C2; c3 = cable.param.C3;

   Z0 = Z00 *  (1 + f1./f ).^Ne1 .* exp(j*(-pi/4 + c1*atan(f./f2)));
   gamma = c2 * log(10)/20 * ((1+ f./f4)./( 1 + f3./f)).^Ne2 + ...
           j*pi/180 * c3 .* (f./f5).^Ne3 .* (1 + f./f5).^Ne4;
   R0=Z0(1); % This is not the correct value,
             % but this model is not good below 10kHz anyway
 
case 'KPN'
   Z0oo=cable.param.Z0oo; c=cable.param.c; Rss00=cable.param.Rss00; 
   tan_fi=cable.param.tan_fi; Kf=cable.param.Kf; K1=cable.param.K1;
   Kn=cable.param.Kn; Kc=cable.param.Kc; N=cable.param.N; 
   fc0=cable.param.fc0; M=cable.param.M;
   tji=(1+j).*sqrt(f.*4*pi*1e-7/Rss00/(Kn*Kf));
   Z=Rss00.*(1+K1*Kf*(tji.*coth(4/3.*tji)-3/4))+j*2*pi.*f.*Z0oo./c;
   Y=j*2*pi.*f/Z0oo/c.*(1+(Kc-1)./(1+(f./fc0).^N))+tan_fi./(Z0oo.*c).*(2*pi*f).^M;
   Z0=sqrt(Z./Y);
   gamma=sqrt(Z.*Y).*1e3;
   R0=Rss00*1e3;   

case 'MAR'
   % Evaluate Zs, Yp, for MAR1 cable model in z, y
   R0=cable.param.R0;           % Roo/km
   Linf=cable.param.Linf;       % Linf/km
   a=cable.param.a;             % (proximity factor)
   b=cable.param.b; 
   c=cable.param.c; 
   C1M=cable.param.C1M;         % Cp/km(@1MHz)
   delta=cable.param.delta;     % Delta(constant)
   
   s=j/9*16*4e-4*pi/R0*f;
   jwL=j*2*pi*f*Linf;              
   x=1+a.*s.*(s+b)./(s+c);      
   Z=jwL + R0*(.25 + .75*sqrt(x)); % Z = jwLinf + R0(.25 + .75sqrt( z ))
   jwC=j*2*pi*C1M; 
   yy=delta*2/pi; 
   Y=jwC*f./(1e-6*j*f).^yy;        % Y = jw C1M (j f1M )^-2d/pi
   gamma=sqrt(Z.*Y);
   Z0=sqrt(Z./Y);

case 'MAR2'
   % Evaluate Zs, Yp, for MAR1 cable model in z, y
   R0=cable.param.R0;           % Roo/km
   Linf=cable.param.Linf;       % Linf/km
   a=cable.param.a;             % (proximity factor)
   b=2; 
   c=2.765; 
   C1M=cable.param.C1M;         % Cp/km(@1MHz)
   delta=cable.param.delta;     % Delta(constant)
   
   s=j/9*16*4e-4*pi/R0*f;
   jwL=j*2*pi*f*Linf;              
   x=1+a.*s.*(s+b)./(s+c);      
   Z=jwL + R0*(.25 + .75*sqrt(x)); % Z = jwLinf + R0(.25 + .75sqrt( z ))
   jwC=j*2*pi*C1M; 
   yy=delta*2/pi; 
   Y=jwC*f./(1e-6*j*f).^yy;        % Y = jw C1M (j f1M )^-2d/pi
   gamma=sqrt(Z.*Y);
   Z0=sqrt(Z./Y);
   
case 'RLC'   
   % for f > max(cable.param.ftmpl) not defined !!! 
   % for f<max(cable.param.ftmpl)
   %     linear interpolation of template values is used,
  
   ftmpl=cable.param.ftmpl; 
   
   if max(f)>max(ftmpl);
       estr=sprintf(['WARNING: the used cable type (%s)'...
                     'is only defined up to %gkHz!\n'], ...
		     cable.name, max(ftmpl)/1e3);
   end
   
   % interpolation according to used f-values:
   %method='linear';
   R=interpolate(ftmpl, cable.param.Rtmpl, f);
   L=interpolate(ftmpl, cable.param.Ltmpl, f);
   C=interpolate(ftmpl, cable.param.Ctmpl, f);
   if isfield(cable.param,'Gtmpl')
       G=interpolate(ftmpl, cable.param.Gtmpl, f);
   else
       G=zeros(size(f));
   end;
   
   Z=R+j*2*pi.*f.*L; 
   Y=G+j*2*pi.*f.*C;
   Z0=sqrt(Z./Y);
   gamma=sqrt(Z.*Y)*1e3;      
   R0=cable.param.Rtmpl(1)*1e3; % usually the value at 0 or 1 Hz

end;