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 FSAN 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: 	Calulate tot_noise from analysed 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:
%    FSAN xDSL simulation tool manual
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 1998 by Telia Research AB, Lulea, Sweden; All rights reserved. 
% Project       : FSAN duplex model
% Author(s)     : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%
% CVS:       $Id: evalExperiment.m,v 1.20 2000/07/20 13:43:12 tono Exp $
%% ===========================================================================
% 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 backgorund 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 possiblity to control XTlevel
%      1999-02-04 (DaB)  Added input paramter 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 whitenoise 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) Generlized list (tfplist) handling
%      1999-03-01 (DaB)  Fixed xtalk + bg addition bug, renamed xtalk to tot
%      1999-03-01 (DaB)  Rewritten HAM band defintion 
%      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 getAtt
%      2000-04-06 (UvAn) Using ex.param.expNoise instead of "-12"
%      2000-07-19 (ToNo) Added hook to test the "full" next function
%% ===========================================================================
global ex;

% Need to change
% ---------------
% getMinFreqAxis


% ---------------------------------------
% Test ex structure for consistency
% ---------------------------------------

[ok errstr]=checkEx;
if ~ok,
    message=sprintf('Ex struct not correct:\n%s', errstr);
    error(message);
else 
    if length(errstr)~=0,
        fprintf([errstr '\n'])
    end;
end;

% Extract disturbed modems
dm=getDM(ex.tt.traffic,ex.param.modemlist);

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

temp=size(tt.topology);
no_nodes=temp(1);

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


% ----------------------------------
% Step 1: Calculation of cable model
% ----------------------------------

len(1)=0;
tot_len(1)=0;
for index=2:no_nodes,  
    len(index)=tt.topology{index,1}./1000;
    tot_len(index)=sum(len(1:index));
end; % for 

temp=size(tt.traffic);
for id=1:no_nodes,
   Att_mtr{id,id}=zeros(size(f));
end;

for id=1:no_nodes,
   for i=(id+1):no_nodes,
       start=id;
       stop=i;    
       Att_mtr{start,stop}=getAtt(ex,start,stop,tt.topology,f,Zterm);
    end;
end;

% Set reference cable for powerbackoff 
ind=find(max(len)==len);
start=ind(1)-1;
stop=ind(1);
Ref_Att=getAtt(ex,start,stop,tt.topology,f,Zterm)./len(ind(1));

% Set Reference noise for PBO method
Ref_Noise=ones(size(f)).*10.^(ex.param.backgroundNoise/10);

%------------------------
% calculate HAM band mask
% ------------------------

tfplan = getList(ex.tfplist,ex.param.HAMBandName);
HAM_up=eval(tfplan.PSD.upstream);
HAM_down=eval(tfplan.PSD.downstream); 

%-----------------------------------------------------
% Step 2: Calculate Alien noise or use FSAN models
%----------------------------------------------------

if ~strcmp(ex.param.FSANNoiseModel,'Calculate'),
   tfplan = getList(ex.tfplist,ex.param.FSANNoiseModel);
   PSD_upstream=eval(tfplan.PSD.upstream);
   PSD_downstream=eval(tfplan.PSD.downstream);
   if 1,
       % Typically we do not take into account the factor sqrt(1-att^4)
       dB_xtalk=ex.param.XTlevel.NEXT+15*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)
       tmp=tt.traffic(dist_route,:);
       LT_node=tmp{1};
       NT_node=tmp{2};
       att=Att_mtr{LT_node,NT_node};  
       dB_xtalk=ex.param.XTlevel.NEXT+15*log10(f./1e6)+ ...
                10*log10(1-10.^(4*att/10));
   end;
   next=10.^(dB_xtalk/10);
   nr_dist_routes=length(dm);
   for n=1:nr_dist_routes,
      dist_route=dm(n);
      tmp=tt.traffic(dist_route,:);
      LT_node=tmp{1};
      NT_node=tmp{2};
      dist_modem=tmp{3};  
      len_fext=tot_len(NT_node)-tot_len(LT_node);
      att=Att_mtr{LT_node,NT_node};  
      dB_xtalk=ex.param.XTlevel.FEXT+20.*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(n).NT=noise_NT;
      Alien(n).LT=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);
   	tmp=tt.traffic(dist_route,:);
   	LT_node=tmp{1};
   	NT_node=tmp{2};
   	dist_modem=tmp{3};        
      % ----------------------------------------------------------------------------  
      % Step 2.1.1: PSD mask with PBO (only upstream implemented for disturbed modem
      % ----------------------------------------------------------------------------   
        
      tfplan = getList(ex.tfplist,dist_modem);
      PSD_up=eval(tfplan.PSD.upstream);
      PSD_down=eval(tfplan.PSD.downstream);  
      HAM=tfplan.PSD.HAM.active;

      minValue_up=min(PSD_up);
      minValue_down=min(PSD_down);

      [E_up, E_down]=evalPBO(tfplan.PSD.PBO,Att_mtr,LT_node,NT_node,Ref_Att,Ref_Noise./PSD_up);
      pbo_index_up=find(PSD_up>10.^expNoise);
      pbo_index_down=find(PSD_down>10.^expNoise);
      PSD_up(pbo_index_up)=PSD_up(pbo_index_up).*E_up(pbo_index_up);
      PSD_down(pbo_index_down)=PSD_down(pbo_index_down).*E_down(pbo_index_down);
      PSD_up(find(PSD_up<minValue_up))=minValue_up*ones(size(find(PSD_up<minValue_up)));
      PSD_down(find(PSD_down<minValue_down))=minValue_down*ones(size(find(PSD_down<minValue_down)));
      if HAM,
         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);
      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 index=1:traffic_routes,
          traffic=tt.traffic(index,:);
          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);
          HAM=tfplan.PSD.HAM.active;
          PSD_up=eval(tfplan.PSD.upstream);
          PSD_down=eval(tfplan.PSD.downstream); 
          minValue_up=min(PSD_up);
          minValue_down=min(PSD_down);
          [E_up, E_down]=evalPBO(tfplan.PSD.PBO,Att_mtr,start_node,stop_node,Ref_Att,Ref_Noise./PSD_up);
          pbo_index_up=find(PSD_up>10.^expNoise);
          pbo_index_down=find(PSD_down>10.^expNoise);
          PSD_up(pbo_index_up)=PSD_up(pbo_index_up).*E_up(pbo_index_up);
          PSD_down(pbo_index_down)=PSD_down(pbo_index_down).*E_down(pbo_index_down);
          PSD_up(find(PSD_up<minValue_up))=minValue_up*ones(size(find(PSD_up<minValue_up)));
          PSD_down(find(PSD_down<minValue_down))=minValue_down*ones(size(find(PSD_down<minValue_down)));
            if HAM,
                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(n).NT=noise_NT;
        Alien(n).LT=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; % fsan_noise

% --------------------------------------------------
% Step 3: Calculate tot_noise from analysed systems
% --------------------------------------------------

nr_dist_routes=length(dm);
for n=1:nr_dist_routes,
    dist_route=dm(n);
    tmp=tt.traffic(dist_route,:);
    LT_node=tmp{1};
    NT_node=tmp{2};
    dist_modem=tmp{3};   
    % ---------------------------------------------------------------------------------  
    % Step 3.1.1: PSD mask with PBO (only upstream PBO implemented) for disturbed modem
    % ---------------------------------------------------------------------------------
    tfplan = getList(ex.tfplist,dist_modem);
    HAM=tfplan.PSD.HAM.active;
    PSD_up=eval(tfplan.PSD.upstream);
    PSD_down=eval(tfplan.PSD.downstream);  
    minValue_up=min(PSD_up);
    minValue_down=min(PSD_down);
    d=tot_len(NT_node)-tot_len(LT_node);
    Ref_Noise=result(n).LT.Alien_noise+10.^(ex.param.backgroundNoise/10);
    [E_up, E_down]=evalPBO(tfplan.PSD.PBO,Att_mtr,LT_node,NT_node,Ref_Att,Ref_Noise./PSD_up);
    pbo_index_up=find(PSD_up>10.^expNoise);
    pbo_index_down=find(PSD_down>10.^expNoise);
    PSD_up(pbo_index_up)=PSD_up(pbo_index_up).*E_up(pbo_index_up);
    PSD_down(pbo_index_down)=PSD_down(pbo_index_down).*E_down(pbo_index_down);
    PSD_up(find(PSD_up<minValue_up))=minValue_up*ones(size(find(PSD_up<minValue_up)));
    PSD_down(find(PSD_down<minValue_down))=minValue_down*ones(size(find(PSD_down<minValue_down)));

    if HAM,
   	PSD_up=min(PSD_up,HAM_up);
      PSD_down=min(PSD_down,HAM_down);
    end;
    H=10.^(Att_mtr{LT_node,NT_node}./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
    % -----------------------------------------------------
    tmp=size(tt.traffic);
    traffic_routes=tmp(1);   
    tot_NT.NEXT=ones(2,1)*Alien(n).NT.NEXT; 
    tot_NT.FEXT=ones(2,1)*Alien(n).NT.FEXT;
    tot_NT.thirdCXT=ones(2,1)*Alien(n).NT.thirdCXT;
    tot_LT.NEXT=ones(2,1)*Alien(n).LT.NEXT; 
    tot_LT.FEXT=ones(2,1)*Alien(n).LT.FEXT;
    tot_LT.thirdCXT=ones(2,1)*Alien(n).LT.thirdCXT;
    
    for index=1:traffic_routes,
        index;
        traffic=tt.traffic(index,:);
        start_node=traffic{1};
        stop_node=traffic{2};
        dist=traffic{3};
        num_dist=traffic{4};
        d=tot_len(stop_node)-tot_len(start_node);
        if strcmp(dist,dist_modem)
            xtalk_noise=1;         
            tfplan = getList(ex.tfplist,dist);
            PSD_up=eval(tfplan.PSD.upstream);
            PSD_down=eval(tfplan.PSD.downstream); 
            minValue_up=min(PSD_up);
            minValue_down=min(PSD_down);
            [E_up, E_down]=evalPBO(tfplan.PSD.PBO,Att_mtr,start_node,stop_node,Ref_Att,Ref_Noise./PSD_up);            
            pbo_index_up=find(PSD_up>10.^expNoise);
            pbo_index_down=find(PSD_down>10.^expNoise);
            PSD_up(pbo_index_up)=PSD_up(pbo_index_up).*E_up(pbo_index_up);
            PSD_down(pbo_index_down)=PSD_down(pbo_index_down).*E_down(pbo_index_down);
            PSD_up(find(PSD_up<minValue_up))=minValue_up*ones(size(find(PSD_up<minValue_up)));
            PSD_down(find(PSD_down<minValue_down))=minValue_down*ones(size(find(PSD_down<minValue_down)));

            if tfplan.PSD.HAM.active,
         		PSD_up=min(PSD_up,HAM_up);
         		PSD_down=min(PSD_down,HAM_down);
            end;
            
            if (start_node==LT_node)&(stop_node==NT_node),
               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; % if 
            
            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; %if
         	end; %if            
                [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; % for looop     
         else,
            xtalk_noise=0;
     	end; % if VDSL 
    end; % traffic_routes
    
    tot_noise_NT(1,:)=tot_NT.NEXT(1,:)+tot_NT.FEXT(1,:)+tot_NT.thirdCXT(1,:)+10.^(ex.param.backgroundNoise/10);
    tot_noise_NT(2,:)=tot_NT.NEXT(2,:)+tot_NT.FEXT(2,:)+tot_NT.thirdCXT(2,:)+10.^(ex.param.backgroundNoise/10);
    tot_noise_LT(1,:)=tot_LT.NEXT(1,:)+tot_LT.FEXT(1,:)+tot_LT.thirdCXT(1,:)+10.^(ex.param.backgroundNoise/10);
    tot_noise_LT(2,:)=tot_LT.NEXT(2,:)+tot_LT.FEXT(2,:)+tot_LT.thirdCXT(2,:)+10.^(ex.param.backgroundNoise/10);
    %-----------------------------------------------------------------------------------
    % Step 3.2 Write tot_noise and Tx, Rx_signals for disturbed modem n in res. struct
    %-----------------------------------------------------------------------------------
    result(n).Modem.Name=dist_modem;
    result(n).Modem.LT_Node=LT_node;
    result(n).Modem.NT_Node=NT_node;
    result(n).NT.Rx_signal=rec_signal_NT;
    result(n).NT.Tx_signal=tr_signal_NT;
    result(n).NT.Tot_noise.up=tot_noise_NT(1,:);
    result(n).NT.Tot_noise.down=tot_noise_NT(2,:);
    result(n).LT.Rx_signal=rec_signal_LT;
    result(n).LT.Tx_signal=tr_signal_LT;
    result(n).LT.Tot_noise.up=tot_noise_LT(1,:);
    result(n).LT.Tot_noise.down=tot_noise_LT(2,:);
    end; % dist_routes