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 duplex simulation memo, Telia 1998
%% ===========================================================================

%% ===========================================================================
% Copyright (C) 1998 by Telia Research AB, Lulea, Sweden; All rights reserved. 
% Project       : FSAN duplex model
% Author(s)     : Tomas Nordstrom (Tomas.B.Nordstrom@Telia.se)
%               : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%
% CVS:       $Id: evalExperiment.m,v 1.5 1999/03/31 15:16:41 dab 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
%% ===========================================================================
global ex;

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

ex.tt.dm=getDM(ex.tt.traffic,ex.param.modemlist);

if ex.param.fastrecalc,
    tmp=[];
else
    tmp=[ex.param.minfrequency:ex.param.granularity:ex.param.maxfrequency];
end;  
ex.param.f = union(getMinFreqaxis(ex.param.maxfrequency),tmp);

temp=size(ex.tt.topology);
no_nodes=temp(1);
f=ex.param.f; 		% Frequency axis for simulation 
                        % // not right when loading ex?
filename=ex.tt.name;

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

len(1)=0;
tot_len(1)=0;
Att(1,:)=zeros(size(ex.param.f));
tot_Att(1,:)=Att(1,:);
for index=2:no_nodes,  
    len(index)=ex.tt.topology{index,1}./1000;
    tot_len(index)=sum(len(1:index));
    cable_model=ex.tt.topology{index,2};    
    Att(index,:) = getCableModel(cable_model,ex.param.f).*len(index);
    tot_Att(index,:) = tot_Att(index-1,:)+Att(index,:);
end; % for 

% Set reference cable for powerbackoff 
ind=find(max(len)==len); 
Ref_Att=Att(ind(1),:)./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.upstream);
HAM_down=eval(tfplan.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.upstream);
   PSD_downstream=eval(tfplan.downstream);
   dB_xtalk=ex.param.XTlevel.NEXT+15*log10(ex.param.f./1e6);
   next=10.^(dB_xtalk/10);
   nr_dist_routes=length(ex.tt.dm);
   for n=1:nr_dist_routes,
      dist_route=ex.tt.dm(n);
      tmp=ex.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=tot_Att(NT_node,:)-tot_Att(LT_node,:);
      dB_xtalk=ex.param.XTlevel.FEXT+20.*log10(ex.param.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(ex.param.f));
      noise_LT.NEXT=next.*PSD_downstream;
      noise_LT.FEXT=fext.*PSD_upstream;
      noise_LT.thirdCXT=zeros(size(ex.param.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(ex.tt.dm);
    for n=1:nr_dist_routes,
   	dist_route=ex.tt.dm(n);
   	tmp=ex.tt.traffic(dist_route,:);
   	LT_node=tmp{1};
   	NT_node=tmp{2};
   	dist_modem=tmp{3};   
   	if strncmp(dist_modem,'VDSL',4)
            dist_modem=ex.param.vdslDuplex;
            HAM=ex.param.HAM;
        else
            HAM=0;
   	end; %if
        
      % ----------------------------------------------------------------------------  
   	% Step 2.1.1: PSD mask with PBO (only upstream implemented for disturbed modem
      % ----------------------------------------------------------------------------   
        
      tfplan = getList(ex.tfplist,dist_modem);
      PSD_up=eval(tfplan.upstream);
      PSD_down=eval(tfplan.downstream);  
      minValue_up=min(PSD_up);
      minValue_down=min(PSD_down);

      d=tot_len(NT_node)-tot_len(LT_node);
      dmax=ex.param.PBOmaxlen./1000;
      PBOmethod=tfplan.PBOmethod;
      PBOparam=tfplan.PBOparam;
      [E_up, E_down]=evalPBO(PBOmethod,PBOparam,tot_Att,LT_node,NT_node,Ref_Att,Ref_Noise,dmax);
      pbo_index_up=find(PSD_up>10.^-12);
      pbo_index_down=find(PSD_down>10.^-12);
      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(ex.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=ex.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);
            index;
            if strncmp(dist,'VDSL',4),
                   dist=ex.param.vdslDuplex;
                   HAM=ex.param.HAM;
            end; 
            if strncmp(dist,dist_modem,4)      
                xtalk_noise=1;
            else
                xtalk_noise=0;
                HAM=0;
            end; %if
            tfplan = getList(ex.tfplist,dist);
            PSD_up=eval(tfplan.upstream);
            PSD_down=eval(tfplan.downstream); 
            minValue_up=min(PSD_up);
            minValue_down=min(PSD_down);
            PBOmethod=tfplan.PBOmethod;
            PBOparam=tfplan.PBOparam;
            [E_up, E_down]=evalPBO(PBOmethod,PBOparam,tot_Att,start_node,stop_node,Ref_Att,Ref_Noise,dmax);
            pbo_index_up=find(PSD_up>10.^-12);
            pbo_index_down=find(PSD_down>10.^-12);
            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,tot_Att,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
                  
         
        % Add White noise to Alien_noise
		  tot_noise_NT=noise_NT.NEXT+noise_NT.FEXT+noise_NT.thirdCXT+10.^(ex.param.backgroundNoise/10);
        tot_noise_LT=noise_LT.NEXT+noise_LT.FEXT+noise_LT.thirdCXT+10.^(ex.param.backgroundNoise/10);
 
        %------------------------------------------------------------------- 
        % 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(ex.tt.dm);
for n=1:nr_dist_routes,
    dist_route=ex.tt.dm(n);
    tmp=ex.tt.traffic(dist_route,:);
    LT_node=tmp{1};
    NT_node=tmp{2};
    dist_modem=tmp{3};   
    if strncmp(dist_modem,'VDSL',4)
        dist_modem=ex.param.vdslDuplex;
        HAM=ex.param.HAM;
    else
        HAM=0;
    end; %if
    
    % ---------------------------------------------------------------------------------  
    % Step 3.1.1: PSD mask with PBO (only upstream PBO implemented) for disturbed modem
    % ---------------------------------------------------------------------------------
    tfplan = getList(ex.tfplist,dist_modem);
    PSD_up=eval(tfplan.upstream);
    PSD_down=eval(tfplan.downstream);  
    minValue_up=min(PSD_up);
    minValue_down=min(PSD_down);
    d=tot_len(NT_node)-tot_len(LT_node);
    dmax=ex.param.PBOmaxlen./1000; 
    PBOmethod=tfplan.PBOmethod;
    PBOparam=tfplan.PBOparam;
    Ref_Noise=result(n).LT.Alien_noise+10.^(ex.param.backgroundNoise/10);
    [E_up, E_down]=evalPBO(PBOmethod,PBOparam,tot_Att,LT_node,NT_node,Ref_Att,Ref_Noise,dmax);
    pbo_index_up=find(PSD_up>10.^-12);
    pbo_index_down=find(PSD_down>10.^-12);
    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.^((tot_Att(NT_node,:)-tot_Att(LT_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(ex.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=ex.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);
        dmax=ex.param.PBOmaxlen./1000;
        if strncmp(dist,'VDSL',4),
               dist=ex.param.vdslDuplex;
        end;
        if strncmp(dist,dist_modem,4)
            xtalk_noise=1;         
            tfplan = getList(ex.tfplist,dist);
            PSD_up=eval(tfplan.upstream);
            PSD_down=eval(tfplan.downstream); 
            minValue_up=min(PSD_up);
            minValue_down=min(PSD_down);
            PBOmethod=tfplan.PBOmethod;
            PBOparam=tfplan.PBOparam;
            [E_up, E_down]=evalPBO(PBOmethod,PBOparam,tot_Att,start_node,stop_node,Ref_Att,Ref_Noise,dmax);            
            pbo_index_up=find(PSD_up>10.^-12);
            pbo_index_down=find(PSD_down>10.^-12);
            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 ex.param.HAM,
         		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,tot_Att,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=dist_modem;
    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