function [Rate_LT, Rate_NT, Margin_LT, Margin_NT,TPSD_LT,TPSD_NT] = calcResultSCM(tfplan,res,lc,f)
%% ===========================================================================
%calcSCMRate - Calculates the total rate for a SCM system
% Parameter:    tfplan	     The TF plan used for calculating the bit rate
% Parameter:    res  	     Result structure
% Parameter:    lc           The linecode 
% Parameter:    f            Frequency axis
%
% Returns:      Rate__LT     Resulting bitrates at LT side (in Mbit)
% Returns:      Rate__NT     Resulting bitrates at NT side (in Mbit)
% Returns:      Margin_LT    Margin at LT side (in dB)    
% Returns:      Margin_NT    Margin at NT side (in dB)    
% Example(s):
%
% NOTE: This routine is not up-to-date and a new implementation
%       which uses folded spectra and filters between various bands
%       is needed.
%
% Algorithmic details:
%	The calculation of the rate is basically an approximative
%	calculation of an integral, with b and f as the 
%	y- and x-axis. The algorithm approximates the integral by 
%	taking the average value between two points multiplied 
%  	with the frequency difference of the two points. 
%
% Reference:
%    FTW's xDSLsimu manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1999-2000 by Telia Research AB, Lulea, Sweden;                
%       2000-2009 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: calcResultSCM.m 752 2009-01-02 13:03:52Z tono $
%% ===========================================================================
% Change History
%      1999-08-13 (DaB)  Created
%      1999-10-11 (ToNo) Reimplemented, now based on the new calcTheoRate
%      1999-11-22 (DaB)  Margin Output added
%      1999-11-25 (DaB)  xtalk_margin defined
%      2000-01-12 (DaB)  power control added
%      2000-04-03 (DaB)  CalcRate removed
%      2000-04-04 (DaB)  Signal_margin added
%      2000-04-06 (UvAn) Passing ftplan, lc and f into the
%                        function; removed global ex
%      2000-04-11 (DaB)  added sync switch
%      2000-05-08 (UvAn) Removed f=ex.freq... in the code 
%                        and using the PSDfs and PSD.fc 
%                        defined in the ftplan instead of those in lc.
%      2000-05-26 (UvAn) Added possibillities to fix the bitrate 
%      2000-05-29 (UvAn) Fixed error when accessing non defined 
%                        tfplan.fixBitrate.param...  
%      2003-12-02 (ToNo) Added support for lc dependent modem noise
%      2005-05-13 (ToNo) Split out the local funcion carrierRate
%      2008-04-18 (ToNo) Added support for return of actual true Tx PSDs
%% ===========================================================================

if nargin<4,
   error('Not enough input arguments.');
end

% Parameters defined in lc list 
alpha	          = lc.param.alpha;
refSNR	          = lc.param.refSNR;
SNRloss           = lc.param.SNRloss;
codingGain        = lc.param.codingGain;
constellation     = lc.param.constellation;
efficiencyLoss    = lc.param.efficiencyLoss;
xtalk_margin      = lc.param.xtalk_margin;
signal_margin     = lc.param.signal_margin;
SNRmax            = lc.param.SNRMax;
Px                = lc.param.Px;

TimeDivision_up	  = tfplan.timeDivision.up;
TimeDivision_down = tfplan.timeDivision.down;
sync              = tfplan.timeDivision.sync;

reqSNR = 10*log10(10.^((refSNR+SNRloss-codingGain+signal_margin)/10).*(2.^(constellation)-1)); %

% Calculate Theoretical SNR
Rx_NT=res.NT.Rx_signal;
Tx_NT=res.NT.Tx_signal;
tmp=find(Tx_NT>10.^-14);
%LT_mask=10.^14.*ones(size(f));
%LT_mask(tmp)=ones(size(tmp));
if (sync==1),
   Tot_NT=res.NT.Tot_noise.down;
   Tot_LT=res.LT.Tot_noise.up;
else
   Tot_NT=res.NT.Tot_noise.up;
   Tot_LT=res.LT.Tot_noise.down;
end;
Rx_LT=res.LT.Rx_signal;
Tx_LT=res.LT.Tx_signal;
tmp=find(Tx_LT>10.^-14);
%NT_mask=10.^14.*ones(size(f));
%NT_mask(tmp)=ones(size(tmp));

noise_LT=Tot_LT.*10.^(xtalk_margin/10);
noise_NT=Tot_NT.*10.^(xtalk_margin/10);

% Check if Total power is OK
Px_LT=10*log10(trapz(f,Tx_LT));
Px_NT=10*log10(trapz(f,Tx_NT));
if Px_LT>Px,
   scale_LT=10.^((Px-Px_LT)/10);
   fprintf('Warning power exceeded at LT side: %g dBm (%g dBm allowed) result is adjusted \n',Px_LT,Px);
else
   scale_LT=1;
end;
if Px_NT>Px,
   scale_NT=10.^((Px-Px_NT)/10);
   fprintf('Warning power exceeded at NT side: %g dBm (%g dBm allowed) result is adjusted \n',Px_NT,Px);
else
   scale_NT=1;
end;

% Add modem noise
if isfield(lc.param,'modemNoise')
    noise_LT = noise_LT + 10.^(eval([lc.param.modemNoise.LT ';'])/10);
    noise_NT = noise_NT + 10.^(eval([lc.param.modemNoise.NT ';'])/10);
end;

SNR_LT=scale_LT.*Rx_LT./noise_LT;
SNR_NT=scale_NT.*Rx_NT./noise_NT;


if tfplan.fixBitrate.active,
    [R_LT, M_LT]=calcResultSCM_cRate(f,lc,alpha,SNRmax,reqSNR,SNR_LT,...
                             tfplan.PSD.fc.up,tfplan.PSD.fs.up,...
                             tfplan.fixBitrate.active,...
                             tfplan.fixBitrate.param.constellation.up);
    [R_NT, M_NT]=calcResultSCM_cRate(f,lc,alpha,SNRmax,reqSNR,SNR_NT,...
                             tfplan.PSD.fc.down,tfplan.PSD.fs.down,...
                             tfplan.fixBitrate.active,...
                             tfplan.fixBitrate.param.constellation.down);
else
    [R_LT, M_LT]=calcResultSCM_cRate(f,lc,alpha,SNRmax,reqSNR,SNR_LT,...
                             tfplan.PSD.fc.up,tfplan.PSD.fs.up);	
    [R_NT, M_NT]=calcResultSCM_cRate(f,lc,alpha,SNRmax,reqSNR,SNR_NT,...
                             tfplan.PSD.fc.down,tfplan.PSD.fs.down);
end;

Rate_LT=R_LT./1e6.*TimeDivision_up.*(1-efficiencyLoss);
Rate_NT=R_NT./1e6.*TimeDivision_down.*(1-efficiencyLoss);
Margin_LT=M_LT+xtalk_margin+signal_margin;
Margin_NT=M_NT+xtalk_margin+signal_margin;
TPSD_NT=res.NT.Tx_signal*scale_NT;
TPSD_LT=res.LT.Tx_signal*scale_LT;