function [Rate_LT, Rate_NT, Margin_LT, Margin_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:
%    FSAN xDSL simulation tool manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 1999,2000 by Telia Research AB, Lulea, Sweden;
%               2000 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                     All rights reserved. 
% Project       : FSAN duplex model
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%
% CVS:       $Id: calcResultTheo.m
%% ===========================================================================
% 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    
%% ===========================================================================

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

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;
Px             = lc.param.Px;

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

% Pick up frequency plan info
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, 
% and 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 (NOT USED?)
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) resultshould be 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 should be adjusted \n'],Px_NT,Px);
else
   scale_NT=1;
end;

SNR_LT=Rx_LT./noise_LT;
SNR_NT=Rx_NT./noise_NT;

% SNR_LT=scale_LT.*Rx_LT./noise_LT;
% SNR_NT=scale_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);

% Select over which frequencies to integrate
index_LT=find((f>=f_start_LT)&(f<=f_stop_LT));
index_NT=find((f>=f_start_NT)&(f<=f_stop_NT));

if tfplan.fixBitrate.active, 
    % If the bitrate is fixed we should calculate the margin
    % NOTE: Experimental code
    reqSNR         = lc.param.reqSNR;

    nani = isnan(SNR_LT); % Fix any NaN
    SNR_LT(nani) = 1e-20;
    nani = isnan(SNR_NT);
    SNR_NT(nani) = 1e-20;
    
    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-reqSNR+codingGain-SNRloss-signal_margin;

    BW_NT=(f_stop_NT-f_start_NT)*2;
    SumSNR_NT=10*log10(trapz(f(index_LT),SNR_NT(index_LT))./BW_LT)
    Margin_NT=SumSNR_NT-reqSNR+codingGain-SNRloss-signal_margin;

    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;
    Margin_NT=xtalk_margin+signal_margin+SNRloss;
end;