function PSD=modelPSD_SDSL(lcname,f,dir)
%% ===========================================================================
% modelPSD_SDSL - Calculates PSDs for SDSL coded 16 PAM 
%
% Parameter:    lcname    Linecode name
% 		f         Frequency
% 		dir       True (1) if it is the downstream (LTU)
% Returns:      PSD       Resulting SDSL PSD
%
%
% References:
%   ETSI TR 101 830-2,  Spectral management on metallic access
%   networks; Part 2: Technical  methods and performance
%   evaluations, m01p20a9
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 2003-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                     All rights reserved. 
% Project       : FTW's xDSLsimu
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: modelPSD_SDSL_sym_SpM.m,v 1.3 2005/01/04 10:32:28 tono Exp $
%% ===========================================================================
% Change History
%      2003-11-11 (ToNo) Created based on an earlier version
%      2004-01-14 (ToNo) First update according to ETSI SpM part 2 complete
%      2005-01-11 (ToNo) Renaming and integration into main code
%      2005-01-13 (ToNo) Made it generic for both sym and asym
%      2008-08-28 (ToNo) Introduced consistent naming (pam -> mod)
%% ===========================================================================
global ex

if nargin<3,
    dir=1;
end;

lc=getList(ex.lclist,lcname);

bits_per_sym = lc.param.mod.bpsym;           % Payload bits per PAM symbol
brate        = lc.param.mod.brate.rate;      % Payload bit rate 
overh        = lc.param.mod.brate.ohead;     % Overhead
fc           = lc.param.mod.highpass3dBf;    % Coupling transformer cutoff fq
R            = ex.param.Zterm; 	             % Line termination impedance
% N            = lc.param.mod.filterorder;     % Order of Tx butterworth filter
floorf       = lc.param.mod.floorf;          % When the PSD start being flat

% Some parameters that are constant (for now)
Kx           = 0.5683e-4;                    % [W/Hz]

% Get the other parameters for the transmit PSD
stype = ~isempty(findstr(lcname,'-sym'));
[fxfactor,fhfactor,N,K]=etsi_SDSLmodelParameters(stype,brate,dir);

% Calculate the values for the variables used in the formula
fsym         = (brate+overh)/bits_per_sym;   % Transmitted symbol rate
fx           = fsym*fxfactor;
fh           = fx*fhfactor;

% The formulas
% Main PSD in mW/Hz
P1 = 1e3 * K/(fx*R) * sinc(f/fx).^2 ./ (1+(fc./f).^2) ./ (1+(f./fh).^(2*N));

% Floor PSD in mW/Hz
P2 = 1e3 * Kx * f.^(-1.5);                   

% Find the best intersect point in frequency vector between P1 and P2
indint=find((P2>=P1));
dind=diff(indint);
tmp1=find(dind~=1)+1;
intersect=indint(tmp1);

% Resulting PSD:
PSD(1:intersect-1)       = P1(1:intersect-1);
PSD(intersect:length(f)) = P2(intersect:length(f));
PSD(f>floorf)=10^-11;    % Fix "P3"