function [Z0,gamma]=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]
%
% Example(s):
%   cable = getList(ex.clist,'DTAG-40');
%   [Z0,G] = getTwoPortModel(cable,ex.param.frequency.f); 
%
% Reference:
%   ETSI VDSL TM6(97)02 (version 01-07-98)
%   Musson, J. "Maximum likelihood estimation of the primary parameters of
%   twisted pair cables", ETSI Contribution 981t08a1, Madrid, 1998
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1999 by Telia Research AB, Lulea, Sweden;                
%       2000 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FSAN xDSL simulation tool
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Gernot Schmid (Gernot.Schmid@arcs.ac.at)
%
% CVS:       $Id: getTwoPortModel.m,v 1.8 2000/09/14 10:42:01 tono Exp $
%% ===========================================================================
% 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
%% ===========================================================================

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

switch cable.model
case 'DTAG'
   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);
   gamma=[gamma (K1+K2.*(f.^K3))/20*log(10)+j*(Kb1.*f+Kb2*sqrt(f))];
   start=stop+1;
   stop=max(find(freq<=30e6));
   f=freq(start:stop)/1e6;
   K1=Ka1(3); K2=Ka2(3); K3=Ka3(3);
   gamma=[gamma (K1+K2.*(f.^K3))/20*log(10)+j*(Kb1.*f+Kb2*sqrt(f))];
   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;
   
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;
   if (as==0) & (r0s==0),
       % no steel core, fall back to the simplified model
       r0s = inf;
   end;
   R=1./( (1./sqrt(sqrt(r0c.^4+ac.*f.^2))) + (1./sqrt(sqrt(r0s.^4+as.*f.^2))));
   L=(l0+loo*(f./fm).^b)./(1+(f./fm).^b);
   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 '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;
   
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 '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; 
   Rtmpl=cable.param.Rtmpl; 
   Ltmpl=cable.param.Ltmpl; 
   Ctmpl=cable.param.Ctmpl;
   
   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=interp1(ftmpl, cable.param.Rtmpl, f, method);
   L=interp1(ftmpl, cable.param.Ltmpl, f, method);
   C=interp1(ftmpl, cable.param.Ctmpl, f, method);
   if exist('cable.param.Gtmpl')
       G=interp1(ftmpl, cable.param.Gtmpl, f, method);
   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;      

end;