function [Rate_LT, Rate_NT, Margin_LT, Margin_NT,TPSD_LT,TPSD_NT] = calcResultTheo(tfplan,res,lc,f)
%% ===========================================================================
%calcResultTheo - Calculates the total theoretical bit rate
%
% Parameter:    tfplan	        The TF plan used for calculating the bit rate
% Parameter:    res	        The result structure
% Parameter:    lc		The linecode  structure
% 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       Resulting margin SNR at LT side (in dB)
% Returns:      Margin_NT       Resulting margin SNR at NT side (in dB)
%
% Example(s):
%
% Algorithmic details:
%
% 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)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%
% CVS:       $Id: calcResultTheo.m 752 2009-01-02 13:03:52Z tono $
%% ===========================================================================
% Change History
%      1999-11-25 (DaB)  Created
%      2000-01-12 (DaB)  Added power check
%      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-06-28 (ToNo) Added (partial) support for fixed bitrate systems
%      2002-10-09 (ToNo) Implemented two-sided power constraints
%      2003-03-20 (ToNo) Implemented disjoint active frequency regions
%      2003-12-01 (ToNo) Made it possible to define active freq with eval
%      2003-12-02 (ToNo) Added support for lc dependent modem noise
%      2008-04-18 (ToNo) Added support for return of actual true Tx PSDs
%% ===========================================================================

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

% Generic parameters defined in lc list
efficiencyLoss = lc.param.efficiencyLoss;
SNRmax	       = lc.param.SNRMax;
codingGain     = lc.param.codingGain;
SNRloss	       = lc.param.SNRloss;
refSNR         = lc.param.refSNR;
xtalk_margin   = lc.param.xtalk_margin;
signal_margin  = lc.param.signal_margin;
PxDown         = lc.param.PxDown;
PxUp           = lc.param.PxUp;

ShannonGap     = refSNR-codingGain+SNRloss+signal_margin;
% Note: The xtalk_margin is not included here!! Instead the total noise 
% is increased by xtalk_margin dB:s

% Get some time and frequency parameters
TimeDivision_up         = tfplan.timeDivision.up;
TimeDivision_down       = tfplan.timeDivision.down;
sync                    = tfplan.timeDivision.sync;
f_start_NT              = tfplan.PSD.active.downstream(1);
f_stop_NT               = tfplan.PSD.active.downstream(2);
f_start_LT              = tfplan.PSD.active.upstream(1);
f_stop_LT               = tfplan.PSD.active.upstream(2);

% Pick up results from the result structure
Rx_NT=res.NT.Rx_signal;
Tx_NT=res.NT.Tx_signal;
Rx_LT=res.LT.Rx_signal;
Tx_LT=res.LT.Tx_signal;
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;

% The total noise is increased by xtalk_margin dB:s
% (which is not in the ShannonGap).
noise_LT=Tot_LT.*10.^(xtalk_margin/10);
noise_NT=Tot_NT.*10.^(xtalk_margin/10);

% Check if Total power is OK, if not set up a scaling factor
% At the moment it only issues warnings and do not compensate for it
Px_LT=10*log10(trapz(f,Tx_LT));
Px_NT=10*log10(trapz(f,Tx_NT));
if Px_LT>PxDown,
   scale_LT=10.^((PxDown-Px_LT)/10);
   fprintf(['Warning power exceeded at LT side: %g dBm (%g dBm' ...
            ' allowed) result should be adjusted \n'],Px_LT,PxDown);
else
   scale_LT=1;
end;
if Px_NT>PxUp,
   scale_NT=10.^((PxUp-Px_NT)/10);
   fprintf(['Warning power exceeded at NT side: %g dBm (%g dBm' ...
            ' allowed) result should be adjusted \n'],Px_NT,PxUp);
else
   scale_NT=1;
end;
% SNR_LT=scale_LT.*Rx_LT./noise_LT;
% SNR_NT=scale_NT.*Rx_NT./noise_NT;

% 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;

% Calculate the SNR(f)
SNR_LT=Rx_LT./noise_LT;
SNR_NT=Rx_NT./noise_NT;

%% ===========================================================================
% Calculate theoretical capacity based on Shannon gap formula

% Clip at SNRmax
tmp=find(SNR_LT>=10.^(SNRmax/10));
SNR_LT(tmp)=ones(size(tmp)).*10.^(SNRmax/10);
tmp=find(SNR_NT>=10.^(SNRmax/10));
SNR_NT(tmp)=ones(size(tmp)).*10.^(SNRmax/10);

% Clear SNR for non-active frequencies (if defined)
if isfield(tfplan.PSD.active,'actmaskds') &&  ~isempty(tfplan.PSD.active.actmaskds)
    % Clear nonactive freq
    if ~ischar(tfplan.PSD.active.actmaskds)
        SNR_LT(find(tfplan.PSD.active.actmaskus<1))=1e-20; 
        SNR_NT(find(tfplan.PSD.active.actmaskds<1))=1e-20;
    else
        SNR_LT(find(eval(tfplan.PSD.active.actmaskus)<1))=1e-20; 
        SNR_NT(find(eval(tfplan.PSD.active.actmaskds)<1))=1e-20;
    end;
end
% Set over which frequencies we should integrate
index_LT=find((f>=f_start_LT)&(f<=f_stop_LT));
index_NT=find((f>=f_start_NT)&(f<=f_stop_NT));

nani = isnan(SNR_LT); % Fix any NaN
SNR_LT(nani) = 1e-20;
nani = isnan(SNR_NT);
SNR_NT(nani) = 1e-20;

if tfplan.fixBitrate.active, 
    % If the bitrate is fixed we should calculate the margin
    fprintf('Warning (calcResultTheo): fixBitrate is experimental code\n');

    % Note the difference between refSNR and reqSNR (is it handled correctly)?
    targetSNR = lc.param.reqSNR-codingGain+SNRloss+signal_margin;
    
    % Incorrect
    BW_LT=(f_stop_LT-f_start_LT)*2; % What is the correct value/formula?
    SumSNR_LT=10*log10(trapz(f(index_LT),SNR_LT(index_LT))./BW_LT);
    Margin_LT=SumSNR_LT-targetSNR;

    BW_NT=(f_stop_NT-f_start_NT)*2;
    SumSNR_NT=10*log10(trapz(f(index_NT),SNR_NT(index_NT))./BW_NT);
    Margin_NT=SumSNR_NT-targetSNR;

    Rate_LT=0; % How to get the rate used?
    Rate_NT=0;

else
    % If the margin is fixed we should calculate the bitrate
    b_LT=log2(1+SNR_LT./10.^(ShannonGap/10));   
    b_NT=log2(1+SNR_NT./10.^(ShannonGap/10));
    Rate_LT=trapz(f(index_LT),b_LT(index_LT))./1e6.*...
	TimeDivision_up.*(1-efficiencyLoss);
    Rate_NT=trapz(f(index_NT),b_NT(index_NT))./1e6.*...
	TimeDivision_down.*(1-efficiencyLoss);
    
    Margin_LT=xtalk_margin+signal_margin+SNRloss; % Fixed margin!
    Margin_NT=xtalk_margin+signal_margin+SNRloss;

end;

TPSD_NT=res.NT.Tx_signal*scale_NT;
TPSD_LT=res.LT.Tx_signal*scale_LT;