function [Rate_LT, Rate_NT, Margin_LT, Margin_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.
%       FIXME: The local function definition should be moved into a
%       seperate file for octave support
%
% 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-2003 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,v 3.3 2003/12/17 10:15:02 tono Exp $
%% ===========================================================================
% 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
%% ===========================================================================

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 exist('lc.param.modemNoise.LT')~=1
    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]=carrierRate(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]=carrierRate(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]=carrierRate(f,lc,alpha,SNRmax,reqSNR,SNR_LT,...
                             tfplan.PSD.fc.up,tfplan.PSD.fs.up);	
    [R_NT, M_NT]=carrierRate(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_up.*(1-efficiencyLoss);
Margin_LT=M_LT+xtalk_margin+signal_margin;
Margin_NT=M_NT+xtalk_margin+signal_margin;

%
% Local function carrierRate
%
function [rate, margin]=carrierRate(f,lc,alpha,SNRmax,reqSNR,SNR,fc,fs,useFixBitrate,fixConstellation)

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

if nargin<9,
	useFixBitrate = 0;
   fixConstellation = 0;
end

FILT=sqrcPSD(f,alpha,fs,fc);
Rx=SNR.*FILT;
Rx(find(Rx>=10.^(SNRmax/10)))=10.^(SNRmax/10).*ones(size(find(Rx>=10.^(SNRmax/10))));
for a=1:length(fs),
   T=1./fs(a);
   f_min=-(1/(2*T).*(1+alpha))+fc(a);
   f_max=(1/(2*T).*(1+alpha))+fc(a);
   f_ind=find((f>=f_min)&(f<=f_max));
   
   totSNR=10*log10(trapz(f(f_ind),Rx(f_ind))./(f_max-f_min));
   ind=find((totSNR-reqSNR)>=0);
 
   if isempty(ind),
      b(a)=0;
      M(a)=0;
   else
      % Check if the constellation is fixed? 
      if useFixBitrate,
         newIndex=find(fixConstellation(a) >= lc.param.constellation);
         % Change index vector if the fixed one point out smaller constellations
         if max(newIndex) < max(ind),
            ind = newIndex;      
         end;
      end;
      b(a)=lc.param.constellation(max(ind));
      M(a)=totSNR-reqSNR(max(ind));	
   end;
   R(a)=b(a).*fs(a);
end;      
rate=sum(R);
margin=mean(M);