function [bulk_parameters]=generate_bulk_par(scmpar,linkpar,antpar)
%GENERATE_BULK_PAR Generation of SCME bulk parameters
%   [BULK_PAR]=GENERATE_BULK_PAR(SCMPAR,LINKPAR,ANTPAR) generates the
%   "bulk" parameters according to 3GPP TR 25.996. For explanation of
%   the input structs, see SCMPARSET, LINKPARSET, and ANTPARSET. 
%   Denoting with K the number of links, N the number of paths, 
%   M the number of subpaths, the fields BULK_PAR are as follows: 
%
%   When scmpar.ScmOptions is 'none' or 'urban_canyon':
%   delays          - path delays in seconds [KxN]
%   path_powers     - relative path powers [KxN]
%   aods            - angles of departure in degrees over (-180,180) [KxNxM]
%   aoas            - angles of arrival in degrees over (-180,180) [KxNxM]
%   subpath_phases  - random phases for subpaths in degrees over (0,360) [KxNxM]
%   path_losses     - path losses in linear scale [Kx1]
%   shadow_fading   - shadow fading losses in linear scale [Kx1]
% 
%   In addition, when scmpar.ScmOptions is 'los' (in addition to the above):
%   K_factors       - K factors for all links [Kx1]
%   Phi_LOS         - random phases for LOS paths in degrees over (-180,180) [Kx1]
%
%   When scmpar.ScmOptions is 'polarized' (in addition to scmpar.ScmOptions='none'):
%   subpath_phases  - random phases for subpaths in degrees over (0,360)
%                     [Kx4xNxM], where the second dimension are the [VV VH HV HH]
%                     components (iid).
%   xpd             - cross-polarization ratios in linear scale [Kx2xN],
%                     where the (:,1,:)th dimension is the V-to-H power coupling, 
%                     and (:,2,:)th dimension is the H-to-V power coupling.
%
%   In addition, when scmpar.IntraClusterDsUsed='yes', the following fields
%   are included [L=LENGTH(NumSubPathsPerMidpath)]:
%   NumSubPathsPerMidpath   - sum over this vector is M  [1xL]
%   MidPathOrder            - order in which midpaths are summed [1xM]
%   MidPathPowers           - powers of midpaths relative to path power [1xL]
%   MidPathDelays           - delays of midpaths relative to path delay [1xL]
%   path_powers_all         - powers of all midpaths for all links [Kx(L*N)]
%   
%   Also delays has size [Kx(L*N)]. The number of midpaths per path and the
%   delays and powers of midpaths are constant for all links but may vary in
%   different scenarios. 
%
%   Ref. [1]: 3GPP TR 25.996 v6.1.0 (2003-09)
%        [2]: D. Baum et al, VTC spring 2005
%
%   See also SCM.

%   Authors: Jari Salo (HUT), Daniela Laselva (EBIT), Giovanni Del Galdo (TUI),
%   Marko Milojevic (TUI), Pekka Ky鰏ti (EBIT), Christian Schneider (TUI)
%   $Revision: 0.41$  $Date: May 20, 2005$


% Input parameter validity checking is done in the main function.


% extract certain parameters from the input structs
Scenario=scmpar.Scenario;


switch lower(Scenario)
    
    % SUBURBAN MACRO AND URBAN MACRO, [1, Sec. 5.3.1]
    case {'suburban_macro','urban_macro'}
        
        bulk_parameters=macro(scmpar,linkpar,antpar);
        
        
        % URBAN MICRO, [1, Sec. 5.3.2]    
    case {'urban_micro'}        
        
        bulk_parameters=micro(scmpar,linkpar,antpar);
        
end     % end of user parameter generation main program








% FUNCTION DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%
% Fixed AoDs for different scenarios
% Needed when scmpar.FixedAnglesUsed='yes'
function aods=fixedAods(scmpar)
% NOTE: one must set linkpar.ThetaBs=0 and linkpar.ThetaMs=0 to repeat
% Table 5 in [2].


N=6;    % number of paths/clusters
M=20;   % number of subpaths per path


switch lower(scmpar.Scenario)
    
    case {'urban_macro'}
        
        aod_2deg     = [0.0894 0.2826 0.4984 0.7431 1.0257 1.3594 1.7688 2.2961 3.0389 4.3101];     % [1, Table 5.2]
        delta_nm_aod = [aod_2deg; -aod_2deg];
        delta_nm_aod = delta_nm_aod(:);   % this (M x 1) vector is the same for all users and paths
        delta_nm_aod = repmat(delta_nm_aod,1,N);  % a (M x N) matrix
        
        
        if strcmpi(scmpar.BsUrbanMacroAS,'eight')
            
            mean_aods=[81.9720 79.6210 80.5354 98.6319 102.1308 107.0643];
            aods=repmat(mean_aods,M,1)+delta_nm_aod;
            
        else    % AS is 'fifteen'
            
            mean_aods=[-127.2788  -136.8071 -129.9678 -96.2155 -159.5999 173.1860];
            aods=repmat(mean_aods,M,1)+delta_nm_aod;
            
        end
        
    case {'suburban_macro'}  % suburban macro [2]
        
        aod_2deg     = [0.0894 0.2826 0.4984 0.7431 1.0257 1.3594 1.7688 2.2961 3.0389 4.3101];     % [1, Table 5.2]
        delta_nm_aod = [aod_2deg; -aod_2deg];
        delta_nm_aod = delta_nm_aod(:);   % this (M x 1) vector is the same for all users and paths
        delta_nm_aod = repmat(delta_nm_aod,1,N);  % a (M x N) matrix
        
        mean_aods=[-101.3376 -110.9587 -100.8629 -112.9888 -115.5088 -118.0681];
        aods=repmat(mean_aods,M,1)+delta_nm_aod;
        
        
    case {'urban_micro'}
        
        aod_5deg=[0.2236 0.7064 1.2461 1.8578 2.5642 3.3986 4.4220 5.7403 7.5974 10.7753]; % [1, Table 5.2]
        delta_nm_aod = [aod_5deg; -aod_5deg];
        delta_nm_aod=delta_nm_aod(:);   % this (M x 1) vector is the same for all users and paths
        delta_nm_aod=repmat(delta_nm_aod,1,N);  % a (M x N) matrix
        
        mean_aods=[6.6100 50.8297 14.1360 38.3972 6.6690 40.2849];
        aods=repmat(mean_aods,M,1)+delta_nm_aod;
        
        
end % switch


% Fixed angles for AoAs
function aoas=fixedAoas(scmpar)
% NOTE: one must set linkpar.ThetaBs=0 and linkpar.ThetaMs=0 to repeat
% Table 5 in [2].

N=6;    % number of paths/clusters
M=20;   % number of subpaths per path

aoa_35deg    = [1.5679 4.9447 8.7224 13.0045 17.9492 23.7899 30.9538 40.1824 53.1816 75.4274];      % [1, Table 5.2]
delta_nm_aoa = [aoa_35deg; -aoa_35deg];
delta_nm_aoa = delta_nm_aoa(:);       % these are the same for all users and paths
delta_nm_aoa = repmat(delta_nm_aoa,1,N); % a (M x N) matrix


switch lower(scmpar.Scenario)

    case {'urban_macro'}
        if strcmpi(scmpar.BsUrbanMacroAS,'eight')
            
            mean_aoas=[65.7489 143.1863 45.6454 32.5131 -91.0551 -19.1657];
            aoas=repmat(mean_aoas,M,1)+delta_nm_aoa;
            
            
        else    % AS is 'fifteen'
            
            mean_aoas=[76.4750 -14.5707 -11.8704 17.7089 167.6567 139.0774];
            aoas=repmat(mean_aoas,M,1)+delta_nm_aoa;
            
        end

    case {'suburban_macro'}

        mean_aoas=[156.1507 39.3383 -137.2020 115.1626 91.1897 4.6769];
        aoas=repmat(mean_aoas,M,1)+delta_nm_aoa;
        
        
    case {'urban_micro'}
        
        mean_aoas=[0.6966 146.0669 -13.2268  -30.5485 -11.4412 -1.0587];
        aoas=repmat(mean_aoas,M,1)+delta_nm_aoa;
        
        
end % switch






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A function to generate sigma_as, sigma_ds and sigma_sf for all links
% Step 3 in [1, Sec 5.3.1], see also [1, Sec 5.6].
% Here Section 5.6 in [1] is interpreted so that it describes the channel matrix
% generation for a single MS only. Hence, there is inter-site correlation
% only between radio links between a single MS and multiple BSs.
% There is no correlation in shadow fading between different MSs. 
function sigmas=step3(scmpar,linkpar)

% extract certain parameters from the input structs
MsNumber=linkpar.MsNumber(:);  
Scenario=scmpar.Scenario;
ScmOptions=scmpar.ScmOptions;

NumLinks=length(MsNumber);

% matrices from [1,Sec. 5.6]
%alfa_beta=0.5; 
%gamma_beta=-0.6; 
%gamma_alfa=-0.6;
%A=[1 alfa_beta gamma_alfa; alfa_beta 1 gamma_beta;gamma_alfa gamma_beta 1];
%B=[0 0 0; 0 0 0; 0 0 0.5];
Bsq=[0 0 0; 0 0 0; 0 0 0.7071];
bsq=0.7071;     % Bsq(3,3)
% pre-computed value: C=sqrtm(A-B)
C = [0.8997 0.1926 -0.3917; 0.1926 0.8997 -0.3917; -0.3917 -0.3917  0.4395];

% the number of different MS
NumOfMs= max(MsNumber); % MsNumber is a vector! 
if (NumOfMs>10*NumLinks)
    warning('MATLAB:SparseMsNumberVector','Max index of linkpar.MsNumber is large compared to number of links!')
end

switch lower(Scenario)
    
    case {'suburban_macro'}
        
        % general environment parameters for suburban macro [1, Table 5.1]
        mu_as      =  0.69 ;
        epsilon_as =  0.13 ;
        mu_ds      = -6.80 ;
        epsilon_ds =  0.288;
        
        % shadow fading 
        if ~strcmpi(scmpar.AlternativePathloss,'yes')
            switch lower(scmpar.ScmOptions)
                case {'none','polarized','urban_canyon'}
                    sigma_sf_ave   =  8    ; % in dB
                case{'los'}
                    sigma_sf_ave   =  4    ; % in dB
                    warning('no SCM LOS standard deviation of shadowing is defined for suburban macro - urban micro LOS standard deviation of shadowing is used')
            end
        else
            switch lower(scmpar.ScmOptions)
                case {'none','polarized','urban_canyon'}
                    sigma_sf_ave   =  10    ; % in dB
                case{'los'}
                    sigma_sf_ave   =  4    ; % in dB
            end 
        end
        
        % generate alphas, betas and gammas for all links
        abc = C*randn(3,NumLinks);
        
        % inter-site correlation terms for all different MSs 
        gamma= bsq*randn(1,NumOfMs);    % bsq*ksi_3 for all different MSs
        gammas=gamma(MsNumber); gammas=gammas(:).';  % so that works also when NumOfMs==1
        abc(3,:)=abc(3,:) + gammas;    % add inter-site correlation term
      
        sigma_ds = 10.^(epsilon_ds*abc(1,:).' + mu_ds);       
        sigma_as = 10.^(epsilon_as*abc(2,:).' + mu_as);
        sigma_sf = 10.^(0.1*sigma_sf_ave*abc(3,:).');
        
        % output
        sigmas=[sigma_ds sigma_as sigma_sf];
        
        
    case {'urban_macro'}    
        
        % general environment parameters for urban macro [1, Table 5.1]
        if strcmp(scmpar.BsUrbanMacroAS,'fifteen')
            mu_as      =  1.18 ;
            epsilon_as =  0.210;
        else     % Note: 8 degree angle spread is set automatically if no match to 'fifteen'
            mu_as      =  0.810;
            epsilon_as =  0.34 ;
        end    
        
        mu_ds      = -6.18 ;      
            
        % shadow fading    
        if ~strcmpi(scmpar.AlternativePathloss,'yes')
            switch lower(scmpar.ScmOptions)
                case {'none','polarized','urban_canyon'}
                    sigma_sf_ave   =  8    ; % in dB