function [E_up,E_down] = evalPBO(PBO,Att_mtr,start_node,stop_node,Ref_Topology,Ref_Noise);
%% ===========================================================================
% evalPBO -     Calculates PBO for up (and downstream)
%
%        INPUT
%        -----
%        PBOmethod :    PBO struct (method  and param)
%        tot_Att   :    matrix with Attenuation in dB/km
%        start_node:    Start node for PBO modem
%        stop_node :    Stop node for PBO modem
%        Ref_Topology :    Reference Topology for PBO calculation
%        Ref_Noise :    Reference Noise for PBO method
%        Tx :           Transmitted signal from far end without PBO
%        Rx :           Recieved signal at our end without PBO (not available?)
%
%        OUTPUT
%        ------
%        E_up      :    Scaling of transmit upstream PSD
%        E_down    :    Scaling of transmit downstream PSD
%
% Algorithmic details:
%
% Reference:
%     [1] ETSI Document TS 101 270-2, Transmission and Multiplexing
%         (TM); Access transmission systems on metal access cables; Very
%         high-speed Digital Subscriber Line (VDSL); Part 1: Tranceiver
%         Specification; V0.1.2 (2000-03)
%     [2] FSAN xDSL simulation tool manual
%     [3] FSAN VDSL working group, "Power-backoff methods for VDSL",
%         ETSI TM6 983T17A0, Lulea, Sweden June 1998
%     [4] ITU Q4/15 FI-102, Fiji 2000
%     [5] ITU Q4/15 FI-114, Fiji 2000
%     [6] ITU Q4/15 D.815, Geneva 2000
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1998-1999 by Telia Research AB, Lulea, Sweden;                
%       2000-2001 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: evalPBO.m,v 1.12 2001/05/31 10:11:58 tono Exp $
%% ===========================================================================
% Change History
%      1998-12-16 (DaB)  Created     
%      1998-12-22 (ToNo) Added and updated comments
%      1999-02-04 (DaB)  Modified to allowmultiple cables 
%      1999-03-10 (DaB)  Added new PBO methods
%      1999-03-11 (DaB)  New input arguments
%      1999-08-11 (DaB)  Rewritten for new structure
%      1999-09-30 (DaB)  Rewritten for new ex struct
%      1999-11-16 (ToNo) Implementing the multiple reference length method
%      2000-03-23 (ToNo) Implementing the multiple reference freq method
%      2000-04-04 (ToNo) Implementing a number of multiple constant PBO methods
%      2001-05-30 (ToNo) Fixed a bug in how to calculate the
%                        reference attenuation for various lengths
%% ===========================================================================
global ex;

Zterm=ex.param.Zterm;            
f=ex.param.frequency.f; 	% Frequency axis

switch PBO.method

%% ===========================================================================

 case 'None', 	
  % No PBO method
  E_up=ones(1,length(f));
  
%% ===========================================================================
 
 case 'RefLen',	
  % Reference length method
  att=Att_mtr{start_node,stop_node};
  Ref_Topology{2,1} = PBO.param.len(1);
  Ref_Att=getAtt(ex,1,2,Ref_Topology,f,Zterm);
  tmp=Ref_Att-att;
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*10.^(tmp/10);
  
 case 'RefLenM',	
  % Reference length method with multiple references
  att=Att_mtr{start_node,stop_node};
  high=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  for band=1:length(PBO.param.len),
      Ref_Topology{2,1} = PBO.param.len(band);
      Ref_Att=getAtt(ex,1,2,Ref_Topology,f,Zterm);
      tmp=Ref_Att-att;
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;
   
%% ===========================================================================

 case 'RefFEXT',	
  % Reference FEXT method
  l_ref=PBO.param.len(1);
  att=Att_mtr{start_node,stop_node};
  len_att=sumLen(start_node,stop_node);
  Ref_Topology{2,1} = l_ref;
  Ref_Att=getAtt(ex,1,2,Ref_Topology,f,Zterm);
  tmp=Ref_Att+10*log10(l_ref) ...
      -att-10*log10(len_att);
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*10.^(tmp/10);

 case 'ConstantRefFEXTM'	
  % Flat spectrum reference FEXT (ITU FI-114)
  att=Att_mtr{start_node,stop_node};
  len_att=sumLen(start_node,stop_node);
  high=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  
  for band=1:length(PBO.param.len),
      l_ref=PBO.param.len(band);
      f_ref=PBO.param.freq(band); % find the level at the reference frequency 
      f_index=max(find(f<=f_ref));
      Ref_Topology{2,1} = l_ref;
      Ref_Att=getAtt(ex,1,2,Ref_Topology,f,Zterm);
      tmp=Ref_Att(f_index)+10*log10(l_ref) ...
	  -att(f_index)-10*log10(len_att);
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;

%% ===========================================================================

 case {'RefFreq' , 'Constant'} 
  % Reference frequency method (Constant power)
  att=Att_mtr{start_node,stop_node};

  f_ref=PBO.param.freq(1); % find the constant
  f_index=max(find(f<=f_ref));
  Ref_Topology{2,1} = PBO.param.maxlen;
  Ref_Att=getAtt(ex,1,2,Ref_Topology,f,Zterm);
  tmp=Ref_Att(f_index)-att(f_index);
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*(10.^(tmp/10));
  
 case {'RefFreqM', 'ConstantM'}	
  % Reference frequency (Constant) method with multiple references
  att=Att_mtr{start_node,stop_node};
  high=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  Ref_Topology{2,1} = PBO.param.maxlen;
  Ref_Att=getAtt(ex,1,2,Ref_Topology,f,Zterm);
  for band=1:length(PBO.param.freq),
      f_ref=PBO.param.freq(band); % find the constant
      f_index=max(find(f<=f_ref));
      tmp=Ref_Att(f_index)-att(f_index);
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;
  
%% ===========================================================================

 case 'RefNoise'	
  % Reference Noise method	
  att=Att_mtr{start_node,stop_node};
  len_att=sumLen(start_node,stop_node)./1000;
  tmp=10*log10(Ref_Noise)-att-20*log10(f./1e6)-10*log10(len_att)-ex.param.XTlevel.FEXT;
  tmp=min(tmp,zeros(size(f)));
  E_up=ones(1,length(f)).*(10.^(tmp/10));

%% ===========================================================================
  
 case 'ConstantAvgLogM'	% NOTE/WARNING - This is not ready yet!
  % Flat spectrum constant average log method (ITU Q4/15 D.815, Geneva 2000)
  Tx=-60*ones(1,length(f)); % ToDo: Get correct values for Tx and Rx!
  Rx=-60*ones(1,length(f)); 
  high=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  for band=1:length(PBO.param.k),
      k_ref=PBO.param.k(band);
      ix=find(f<high(band) & f>low(band));
      k=trapz(f(ix),log2(Rx(ix)));
      level=2^((k_ref-k)/(high(band)-low(band))); %%% Should this be 10*log10 ?
      E_up(ix)=level*Tx(ix);
  end;
  
 case 'ConstantRefNoiseM'	
  % Flat spectrum reference noise (ITU Q4/15 FI-114, Fiji 2000)
  att=Att_mtr{start_node,stop_node};
  high=PBO.param.plan.up.fc+PBO.param.plan.up.bw/2;
  low =PBO.param.plan.up.fc-PBO.param.plan.up.bw/2;
  E_up=ones(1,length(f));
  len_att=sumLen(start_node,stop_node)./1000;
  
  for band=1:length(PBO.param.freq),
      q_ref=PBO.param.q(band);    % get the boosting parameter
      f_ref=PBO.param.freq(band); % find the level at the reference frequency 
      f_index=max(find(f<=f_ref));
      tmp=10*log10(q_ref)+10*log10(Ref_Noise) ...
	  -att(f_index)-20*log10(f_ref./1e6)-10*log10(len_att)-...
          ex.param.XTlevel.FEXT;
      tmp=min(tmp,zeros(size(f)));
      ix=find(f<high(band) & f>low(band));
      E_up(ix)=10.^(tmp(ix)/10);
  end;

%% ===========================================================================

 otherwise
  sprintf('Unknown PBO method: %s',PBO.method)
end; %switch

E_down=ones(1,length(f));	% Downstream PBO not implemented !!!!!!!


function len=sumLen(startno,stopno)
global ex;    
len=0;
for nodeno=startno+1:stopno,
    len=len+ex.tt.topology{nodeno,1};
end;