function [Rate_LT, Rate_NT, Margin_LT, Margin_NT] = calcResultDMT(tfplan,res,lc,f);
%% ===========================================================================
%calcDMTRate - Calculates the total rate for a DMT 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 
% Returns:      Margin_NT
% Example(s):
%
% 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:
%    FSAN xDSL simulation tool manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1999-2000 by Telia Research AB, Lulea, Sweden;                
%       2001 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : xDSL simulation tool
% Author(s)     : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: calcResultDMT.m,v 1.9 2001/12/18 09:44:17 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
%      2000-01-12 (DaB)  Using dmt function as bitloading function
%      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
%      2001-07-03 (ToNo) Calculate bit-rates for active tones
%      2001-09-12 (ToNo) Add support for "Maximum aggregate transmit power"
%      2001-12-14 (ToNo) Cleaning up the code (remove unused things)
%% ===========================================================================

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

% 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;
constellation   = lc.param.constellation;
xtalk_margin    = lc.param.xtalk_margin;
signal_margin   = lc.param.signal_margin;
PxUp            = lc.param.PxUp;
PxDown          = lc.param.PxDown;

df              = lc.param.dmt.df;
carriers        = lc.param.dmt.carriers;
ustones         = lc.param.dmt.tonesus;
dstones         = lc.param.dmt.tonesds;
srate           = lc.param.dmt.srate;

newf= df:df:df*carriers; % The frequency of the DMT tones

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

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

% Increase the noise with xtalk_margin 
noise_LT=Tot_LT.*10.^(xtalk_margin/10);
noise_NT=Tot_NT.*10.^(xtalk_margin/10);

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

H_LT=interp1(f,Rx_LT./Tx_NT,newf);
N_LT=interp1(f,noise_LT,newf);
H_NT=interp1(f,Rx_NT./Tx_LT,newf);
N_NT=interp1(f,noise_NT,newf);
PSD_NT=interp1(f,Tx_NT,newf);
PSD_LT=interp1(f,Tx_LT,newf);
Gap=10.^((refSNR+SNRloss-codingGain)/10);

if lc.name(1) == 'A'
    % Limit the total power (this is a formula for ADSL)
    mDown=ceil(10*(-3.65+10*log10(length(dstones)+1)))/10;
    mDown=min(PxDown,mDown);
    mDown=10.^(mDown./10)./df;
    mUp=  ceil(10*(-1.65+10*log10(length(ustones)+1)))/10;
    mUp=min(PxUp,mUp);
    mUp=10.^(mUp./10)./df;
else
    mDown=10.^(lc.param.PxDown./10); % Not check yet!
    mUp  =10.^(lc.param.PxUp./10);
end;

% Do DMT water-filling
[b_LT, E_LT, g_LT, M_LT] =dmt(sqrt(abs(H_LT)),N_LT,mUp,Gap,PSD_NT,constellation,ustones);
[b_NT, E_NT, g_NT, M_NT] =dmt(sqrt(abs(H_NT)),N_NT,mDown,Gap,PSD_LT,constellation,dstones);

%Base bitrate on the active tones
Rate_LT=sum(b_LT(ustones)).*srate./1e6.*TimeDivision_up;
Rate_NT=sum(b_NT(dstones)).*srate./1e6.*TimeDivision_down;

Margin_LT=mean(M_LT)+xtalk_margin+signal_margin;
Margin_NT=mean(M_NT)+xtalk_margin+signal_margin;