function result = evalExperiment
%% ===========================================================================
%evalExperiment - Run the experiment description in the global structure ex
%
% Returns :      result		Structure describing the result
%
% Example(s):
%   res = evalExperiment;
%
% Algorithmic details:
% Step 1:	Cable attenuation calculation
% Step 2:	Calculate Alien noise or use predetermined noise models
% Step 2.1: 	  Find disturbed modem in traffic routes (Main loop start)
% Step 2.1.1: 	    PSD mask with PBO (only upstream implemented) for 
%                   disturbed modem
% Step 2.1.2 	    Calculate noise in LT and NT end for disturbed modem
% Step 2.1.3 	    Writing Alien_noise for disturbed modem n in result struct
% Step 3: 	Calculate tot_noise from analyzed systems
% Step 3.1.1: 	    PSD mask with PBO (only upstream PBO implemented) for
%                   disturbed modem
% Step 3.1.2: 	    Calculate tot-noise for disturbed modem
% Step 3.2 	  Write tot_noise and Tx, Rx_signals for dist. modem n in 
%                 res. struct
%
%
% Reference:
%    FTW's xDSLsimu manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1998-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)
%               : Bo Engstrom (bosse@upzide.com)
%
% CVS:       $Id: evalExperiment.m 752 2009-01-02 13:03:52Z tono $
%% ===========================================================================
% Change History
%      1998-10-19 (ToNo) Created
%      1998-11-09 (DaB)  Editing to get a first simulation going
%      1998-11-12 (ToNo) Restructured the experiment structure
%      1998-11-13 (DaB)  Noise Calculation method added with new structure
%      1998-11-15 (DaB)  Changed to fit new structure
%      1998-12-04 (DaB)  update
%      1998-12-15 (DaB)	 Saving of background noise implemented
%      1998-12-16 (DaB)  Splitted into 3 functions evalExperiment, 
%                        evalNoise, evalPBO
%      1998-12-17 (DaB)  Fixed bug in attenuation 
%      1998-12-21 (ToNo) Fixed bug in saving ex structure
%      1998-12-30 (DaB)  Rewritten cable model
%      1999-01-05 (DaB/ToNo) Added FSAN noise models
%      1999-01-07 (DaB)  Added choice of disturbed modem 
%      1999-01-08 (DaB)  fixed num_dist bug
%      1999-01-09 (DaB)  fixed num_dist bug part 2
%      1999-01-14 (ToNo) Added check for ex struct
%      1999-01-17 (ToNo) Tested setting of max frequency depending on tfplan
%      1999-02-04 (DaB)  Allows multiple cable types
%      1999-02-04 (DaB)  Added possibility to control XTlevel
%      1999-02-04 (DaB)  Added input parameter with disturbed modemlist
%      1999-02-05 (DaB)  Fixed bug in cable calculation
%      1999-02-08 (DaB)  Added modem in result struct
%      1999-02-08 (DaB)  Removed maxfrequency override by vdslduplex
%      1999-02-09 (DaB)  Fixed white noise add bug
%      1999-02-12 (DaB)  Added ex.param.PBOmaxlen as max length for PBO
%      1999-02-23 (DaB)  Corrected dmax bug
%      1999-02-25 (ToNo) Generalized list (tfplist) handling
%      1999-03-01 (DaB)  Fixed xtalk + bg addition bug, renamed xtalk to tot
%      1999-03-01 (DaB)  Rewritten HAM band definition 
%      1999-03-02 (DaB)  Renamed bg to alien
%      1999-03-03 (DaB)  Added new way of calculation noise models 
%      1999-03-11 (DaB)  Updated to work with new PBO methods
%      1999-08-11 (DaB)  Updated for Bridge tap cables
%      1999-08-19 (DaB)  Updated for new tfplan struct
%      1999-09-29 (DaB)  BugFix in Alien_noise
%      1999-09-30 (DaB)  Rewritten for new ex struct
%      2000-04-03 (UvAn) Changed according to new interface of getIL
%      2000-04-06 (UvAn) Using ex.param.expNoise instead of "-12"
%      2000-07-19 (ToNo) Added hook to test the "full" next function
%      2001-02-23 (ToNo) Added support to have different power of f (slant)
%                        for NEXT and FEXT.
%      2001-05-30 (ToNo) Fixed a bug in how to calculate the
%                        reference attenuation for various lengths
%      2001-06-15 (Bosse) Octave port
%      2001-07-05 (ToNo) Now we support two-sided colored background noise
%      2001-10-08 (ToNo) Added support for NEXT with fourth order terms
%      2001-12-05 (ToNo) Fixed 4-order term to be |H|^2^2 instead of |H|^2^4
%      2002-07-01 (PeKa) Octave to Matlab compatibility port
%      2003-01-02 (ToNo) Changed parameters for evalPBO,
%                        instead of a scaling we now use the actual PSD
%      2003-01-03 (ToNo) Split up upstream and downstream PBO functions
%      2003-11-03 (ToNo) Lists now use cell arrays for both Octave and Matlab
%      2004-10-05 (ToNo) Fixed bug in FEXT calculation for nodes
%                        with only one modem present
%      2005-01-04 (ToNo) Removed the need for Hamband definitions if not used
%      2005-04-01 (ToNo) Added the return of total FEXT noise
%      2008-08-25 (ToNo) Added check for "no modem under test" situations
%      2008-08-28 (ToNo) Again fixed 4-order term, for given Noise model
%% ===========================================================================
global ex

% Extract disturbed modems
dm=getDM(ex.tt.traffic,ex.param.modemlist);
if(length(dm)<1)
    error('No modems under test (ex.param.modemlist) was found in ex.tt.traffic?');
end;

% Use internal tt struct to change generic xDSL names to specific
tt=setTT(ex.tt,ex.param.xDSLlist);

% Decide what frequency axis to use 
if isfield(ex.param.frequency,'set') && ~ex.param.frequency.set
    if ex.param.frequency.fastrecalc,
        tmp=[];
    else
        tmp=ex.param.frequency.min:ex.param.frequency.granularity:...
             ex.param.frequency.max;
    end;  
    ex.param.frequency.f = union(getMinFreqaxis(ex.param.frequency.max),tmp);
end

no_nodes=length(tt.topology);

% ----------------------
% Locally used variabels
% ----------------------
f=ex.param.frequency.f;       % Frequency axis for simulation
Zterm=ex.param.Zterm;

% ----------------------------------
% Step 1: Calculation of cable model
% ----------------------------------
len(1)=0;
tot_len=zeros(1,no_nodes);
for index=2:no_nodes,
    len(index)=tt.topology{index}{1}/1000;
    tot_len(index)=sum(len(1:index));
end

tt_topology=tt.topology;
Att_mtr=buildAttMtr(ex,f,no_nodes,tt_topology,Zterm);

%-----------------------------------------------------
% Set up reference cable for powerbackoff
% Uses the longest cable section as reference cable
%-----------------------------------------------------
ind=find(max(len)==len); % Find the longest sections
start=ind(1)-1;          % Use the cable for the first max length section
% stop=ind(1);
% Set up a Reference Topology for us in PBO calcs
Ref_Topology={{0,'','CO',''},... 
              {1000,ex.tt.topology{start+1}{2},'NT','TBD_len'}};

% Set Reference noise for PBO method
Ref_Noise=ones(size(f)).*10.^(eval([ex.param.bgNoise.LT, ';'])/10);

%------------------------
% Calculate HAM band mask
% ------------------------
if isfield(ex.param,'HAMBandName') && (length(ex.param.HAMBandName)>0)
    tfplan = getList(ex.tfplist,ex.param.HAMBandName);
    HAM_up=eval(sprintf('%s;',tfplan.PSD.upstream));
    HAM_down=eval(sprintf('%s;',tfplan.PSD.downstream)); 
else
    % Even if no hamband is defined some tfp's might set it active
    HAM_up=ones(size(f))*1000;
    HAM_down = HAM_up;
end;

%-----------------------------------------------------
% Step 2: Calculate Alien noise or use predetermined noise models
%----------------------------------------------------
l_result={};
Alien.NT={};
Alien.LT={};

if ~strcmpi(ex.param.NoiseModel,'Calculate'),
    tfplan = getList(ex.tfplist,ex.param.NoiseModel);
    PSD_upstream=eval([tfplan.PSD.upstream, ';']);
    PSD_downstream=eval([tfplan.PSD.downstream, ';']);
    if ~ex.param.XTlevel.NEXTx4,
        % Typically we do not take into account the factor sqrt(1-att^4)
        dB_xtalk=ex.param.XTlevel.NEXT+ex.param.XTlevel.NEXTslant*log10(f./1e6);
    else
        % Add the factor sqrt(1-att^4) to the next function
        % This will not make any difference unless the cable is very short
        % A difference of 0.1dB (att. of 4.1dB) at 3kHz   => short is < 900m
        % A difference of 0.1dB (att. of 4.1dB) at 100kHz => short is < 600m
        % A difference of 0.1dB (att. of 4.1dB) at 1MHz   => short is < 200m
        dist_route=dm(1);              % Get the attenuation (beginning to end)
        LT_node=tt.traffic{dist_route}{1};
        NT_node=tt.traffic{dist_route}{2};
        att=Att_mtr{LT_node,NT_node};

        dB_xtalk=ex.param.XTlevel.NEXT+...
                 ex.param.XTlevel.NEXTslant*log10(f./1e6)+ ...
                 10*log10(1-10.^(2*att/10));
    end;

    next=10.^(dB_xtalk/10);
    nr_dist_routes=length(dm);
    for n=1:nr_dist_routes,
        dist_route=dm(n);
        LT_node=tt.traffic{dist_route}{1};
        NT_node=tt.traffic{dist_route}{2};
%        dist_modem=tt.traffic{dist_route}{3};
        len_fext=tot_len(NT_node)-tot_len(LT_node);
        att=Att_mtr{LT_node,NT_node};
        dB_xtalk=ex.param.XTlevel.FEXT+ ...
                 ex.param.XTlevel.FEXTslant.*log10(f/1e6)+...
                 10.*log10(len_fext)+att;
        fext=10.^(dB_xtalk/10);
        noise_NT.NEXT=next.*PSD_upstream;
        noise_NT.FEXT=fext.*PSD_downstream;
        noise_NT.thirdCXT=zeros(size(f));
        noise_LT.NEXT=next.*PSD_downstream;
        noise_LT.FEXT=fext.*PSD_upstream;
        noise_LT.thirdCXT=zeros(size(f));

        Alien.NT{n}=noise_NT;
        Alien.LT{n}=noise_LT;
        res.NT.Alien_noise=noise_NT.NEXT+noise_NT.FEXT+noise_NT.thirdCXT;
        res.LT.Alien_noise=noise_LT.NEXT+noise_LT.FEXT+noise_LT.thirdCXT;
        l_result{n}=res;
    end;
    result=l_result;
else

    % ------------------------------------------------------------------
    % Step 2.1: Find disturbed modem in traffic routes (Main loop start)
    % ------------------------------------------------------------------
    nr_dist_routes=length(dm);
    for n=1:nr_dist_routes,
        dist_route=dm(n);
        LT_node=tt.traffic{dist_route}{1};
        NT_node=tt.traffic{dist_route}{2};
        dist_modem=tt.traffic{dist_route}{3};
        % -----------------------------------------------------------------
        % Step 2.1.1: PSD mask with PBO 
        % (Note: only upstream implemented for disturbed modem)
        % -----------------------------------------------------------------
        
        tfplan = getList(ex.tfplist,dist_modem);
        PSD_down=eval(sprintf('%s;',tfplan.PSD.downstream));
        PSD_upi=eval(sprintf('%s;',tfplan.PSD.upstream));
        PSD_up=calcUPBO(PSD_upi,tfplan.PSD.UPBO,Att_mtr,LT_node,NT_node,...
                        Ref_Topology,Ref_Noise./PSD_upi);

        if tfplan.PSD.HAM.active,
            PSD_up=min(PSD_up,HAM_up);
            PSD_down=min(PSD_down,HAM_down);
        end;

        % ---------------------------------------------------------------
        % Step 2.1.2 Calculate noise in LT and NT end for disturbed modem
        % ---------------------------------------------------------------
        %tmp=size(tt.traffic);
        %traffic_routes=tmp(1);
        traffic_routes=length(tt.traffic);
        noise_LT.NEXT=10.^(-300/10).*ones(size(f)); % approx zero (not set to zero due to log calc)
        noise_LT.FEXT=noise_LT.NEXT;
        noise_LT.thirdCXT=noise_LT.NEXT;
        noise_NT=noise_LT;

        for trix=1:traffic_routes,
            traffic=tt.traffic{trix};
            start_node=traffic{1};
            stop_node=traffic{2};
            dist=traffic{3};
            num_dist=traffic{4};
            xtalk_noise=strcmp(dist,dist_modem);
            tfplan = getList(ex.tfplist,dist);
            PSD_down=eval(sprintf('%s;',tfplan.PSD.downstream)); 
            PSD_upi=eval(sprintf('%s;',tfplan.PSD.upstream));
            PSD_up=calcUPBO(PSD_upi,tfplan.PSD.UPBO,...
                            Att_mtr,start_node, stop_node,...
                            Ref_Topology,Ref_Noise./PSD_upi);

            if tfplan.PSD.HAM.active,
                PSD_up=min(PSD_up,HAM_up);
                PSD_down=min(PSD_down,HAM_down);
            end;
            
            if num_dist>0,
                PSD_up=PSD_up.*(num_dist).^0.6;
                PSD_down=PSD_down.*(num_dist).^0.6;   
            else
                PSD_up=zeros(size(f));
                PSD_down=zeros(size(f));
            end;
            if ~xtalk_noise,
                [add_noise_LT,add_noise_NT]=...
                    evalNoise(f,NT_node,LT_node,start_node,stop_node,...
                              tot_len,Att_mtr,PSD_up,PSD_down,ex.param.XTlevel);
                noise_NT.NEXT=addFSAN(noise_NT.NEXT,add_noise_NT.NEXT);
                noise_NT.FEXT=addFSAN(noise_NT.FEXT,add_noise_NT.FEXT);
                noise_NT.thirdCXT=addFSAN(noise_NT.thirdCXT,...
                                          add_noise_NT.thirdCXT);
                
                noise_LT.NEXT=addFSAN(noise_LT.NEXT,add_noise_LT.NEXT);
                noise_LT.FEXT=addFSAN(noise_LT.FEXT,add_noise_LT.FEXT);
                noise_LT.thirdCXT=addFSAN(noise_LT.thirdCXT,...
                                          add_noise_LT.thirdCXT);   
            end; % xtalk_noise  
        end; % traffic_route
        
        tot_noise_NT=noise_NT.NEXT+noise_NT.FEXT+noise_NT.thirdCXT;
        tot_noise_LT=noise_LT.NEXT+noise_LT.FEXT+noise_LT.thirdCXT;
        
        %----------------------------------------------------------------------
        % Step 2.1.3 Writing Alien_noise for disturbed modem n in result struct
        %----------------------------------------------------------------------  
        Alien.NT{n}=noise_NT;
        Alien.LT{n}=noise_LT;
        res.NT.Alien_noise=tot_noise_NT;
        res.LT.Alien_noise=tot_noise_LT;
        l_result{n}=res;
        
    end; % dist_routes
    result=l_result;
end; % etsi_noise

% --------------------------------------------------
% Step 3: Calculate tot_noise from analysed systems
% --------------------------------------------------
nr_dist_routes=length(dm);

for n=1:nr_dist_routes,
    dist_route=dm(n);
    LT_node=tt.traffic{dist_route}{1};
    NT_node=tt.traffic{dist_route}{2};
    dist_modem=tt.traffic{dist_route}{3};
    % ---------------------------------------------------
    % Step 3.1.1: PSD mask with PBO for disturbed modem
    % (Note: only upstream PBO implemented)
    % ---------------------------------------------------
    tfplan = getList(ex.tfplist,dist_modem);
    PSD_down=eval(sprintf('%s;',tfplan.PSD.downstream));
    x=result{n};
    Ref_Noise=x.LT.Alien_noise+10.^(eval([ex.param.bgNoise.LT, ';'])/10);
    PSD_upi=eval(sprintf('%s;',tfplan.PSD.upstream));
    PSD_up=calcUPBO(PSD_upi,tfplan.PSD.UPBO,Att_mtr,LT_node,NT_node,...
                    Ref_Topology,Ref_Noise./PSD_upi);
    
    if tfplan.PSD.HAM.active,
        PSD_up=min(PSD_up,HAM_up);
        PSD_down=min(PSD_down,HAM_down);
    end;
    Att_mtr_temp=Att_mtr{LT_node,NT_node};
    H=10.^(Att_mtr_temp./20);
    tr_signal_NT=PSD_up;
    tr_signal_LT=PSD_down;
    
    rec_signal_NT=PSD_down.*(H.^2);
    rec_signal_LT=PSD_up.*(H.^2);

    % -----------------------------------------------------
    % Step 3.1.2: Calculate tot-noise for disturbed modem
    % -----------------------------------------------------
    traffic_routes=length(tt.traffic);
    x=Alien.NT{n};
    tot_NT.NEXT=ones(2,1) * x.NEXT;
    tot_NT.FEXT=ones(2,1) * x.FEXT;
    tot_NT.thirdCXT=ones(2,1) * x.thirdCXT;
    x=Alien.LT{n};
    tot_LT.NEXT=ones(2,1) * x.NEXT;
    tot_LT.FEXT=ones(2,1) * x.FEXT;
    tot_LT.thirdCXT=ones(2,1) * x.thirdCXT;
    
    for trix=1:traffic_routes,
        traffic=tt.traffic{trix};
        start_node=traffic{1};
        stop_node=traffic{2};
        dist=traffic{3};
        num_dist=traffic{4};
        if strcmp(dist,dist_modem)
            xtalk_noise=1;         
            tfplan = getList(ex.tfplist,dist);
            PSD_down=eval(sprintf('%s;',tfplan.PSD.downstream));
            PSD_upi=eval(sprintf('%s;',tfplan.PSD.upstream));
            PSD_up=calcUPBO(PSD_upi,tfplan.PSD.UPBO,...
                            Att_mtr,start_node,stop_node,...
                            Ref_Topology,Ref_Noise./PSD_upi);
            
            if tfplan.PSD.HAM.active,
                PSD_up=min(PSD_up,HAM_up);
                PSD_down=min(PSD_down,HAM_down);
            end
            
            if (dist_route==trix), % We should not disturb ourself
                num_dist=num_dist-1;
            end
            
            if num_dist>0,
                PSD_up=PSD_up.*(num_dist).^0.6;
                PSD_down=PSD_down.*(num_dist).^0.6;
            else
                PSD_up=zeros(size(f));
                PSD_down=zeros(size(f));
            end

            if tfplan.timeDivision.up==1,
                TDD=0;
            else
                TDD=1;
            end 
            
            loops=2;
            for loop=1:loops,         
                if TDD,
                    if loop==1,
                        PSD_down2=PSD_down;
                        PSD_down=zeros(size(f));
                    else
%                        PSD_up2=PSD_up;
                        PSD_down=PSD_down2;
                        PSD_up=zeros(size(f));
                    end
                end
                [add_noise_LT,add_noise_NT]=...
                    evalNoise(f,NT_node,LT_node,start_node,...
                              stop_node,tot_len,Att_mtr,PSD_up,...
                              PSD_down,ex.param.XTlevel);
                tot_LT.NEXT(loop,:)=addFSAN(tot_LT.NEXT(loop,:),...
                                            add_noise_LT.NEXT); 
                tot_LT.FEXT(loop,:)=addFSAN(tot_LT.FEXT(loop,:),...
                                            add_noise_LT.FEXT);
                tot_LT.thirdCXT(loop,:)=addFSAN(tot_LT.thirdCXT(loop,:),...
                                                add_noise_LT.thirdCXT);
                tot_NT.NEXT(loop,:)=addFSAN(tot_NT.NEXT(loop,:),...
                                            add_noise_NT.NEXT); 
                tot_NT.FEXT(loop,:)=addFSAN(tot_NT.FEXT(loop,:),...
                                            add_noise_NT.FEXT);     
                tot_NT.thirdCXT(loop,:)=addFSAN(tot_NT.thirdCXT(loop,:),...
                                                add_noise_NT.thirdCXT);     
                
            end
        else
            xtalk_noise=0;
        end; % iff strcmp(dist,dist_modem)
    end; % traffic_routes
    
    tot_noise_NT(1,:)=tot_NT.NEXT(1,:)+tot_NT.FEXT(1,:)+...
        tot_NT.thirdCXT(1,:)+10.^(eval([ex.param.bgNoise.NT, ';'])/10);
    tot_noise_NT(2,:)=tot_NT.NEXT(2,:)+tot_NT.FEXT(2,:)+...
        tot_NT.thirdCXT(2,:)+10.^(eval([ex.param.bgNoise.NT, ';'])/10);
    tot_noise_LT(1,:)=tot_LT.NEXT(1,:)+tot_LT.FEXT(1,:)+...
        tot_LT.thirdCXT(1,:)+10.^(eval([ex.param.bgNoise.LT, ';'])/10);
    tot_noise_LT(2,:)=tot_LT.NEXT(2,:)+tot_LT.FEXT(2,:)+...
        tot_LT.thirdCXT(2,:)+10.^(eval([ex.param.bgNoise.LT, ';'])/10);
    %-----------------------------------------------------------
    % Step 3.2 Write tot_noise and Tx, Rx_signals for disturbed
    % modem n in result structure
    %-----------------------------------------------------------

    res=result{n};
    res.Modem.Name=dist_modem;
    res.Modem.LT_Node=LT_node;
    res.Modem.NT_Node=NT_node;
    res.NT.Rx_signal=rec_signal_NT;
    res.NT.Tx_signal=tr_signal_NT;
    res.NT.Tot_noise.up=tot_noise_NT(1,:);
    res.NT.Tot_noise.down=tot_noise_NT(2,:);
    % Save total FEXT noise (useful for PBO parameter optimization)
    res.NT.FEXT_noise.up=tot_NT.FEXT(1,:);
    res.NT.FEXT_noise.down=tot_NT.FEXT(2,:);    
    res.LT.FEXT_noise.up=tot_LT.FEXT(1,:);    
    res.LT.FEXT_noise.down=tot_LT.FEXT(2,:);    
    %
    res.LT.Rx_signal=rec_signal_LT;
    res.LT.Tx_signal=tr_signal_LT;
    res.LT.Tot_noise.up=tot_noise_LT(1,:);
    res.LT.Tot_noise.down=tot_noise_LT(2,:);
    result{n}=res;

end; % dist_routes