function [H, delays, full_output]=scm(scmpar,linkpar,antpar,initvalues)
%SCM 3GPP Spatial Channel Model (3GPP TR 25.996)
%   H=SCM(SCMPAR,LINKPAR,ANTPAR) is a 5D-array of channel coefficients. For
%   explanation of the input parameter structs, see SCMPARSET, LINKPARSET,
%   and ANTPARSET. SIZE(H)=[U S N T K], where U is the number of MS (RX) 
%   elements, S is the number of BS (TX) elements, N is the number of
%   paths,  T is the number of time samples, and K is the number of links.
%   If K=1, the final dimension will be dropped, i.e. H is a 4D-array.
%
%   [H DELAYS]=SCM(...) outputs also a [KxN] matrix of path delays (in seconds). 
%
%   [H DELAYS BULKPAR]=SCM(...) outputs also the struct BULKPAR, whose fields
%   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  - final 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]
%   delta_t         - time sampling intervals for all links [Kx1]
%
%   In addition, when scmpar.ScmOptions is 'los' (in addition to the above):
%
%   K_factors       - K factors for all links [Kx1]
%   Phi_LOS         - final phases for LOS paths in degrees over (-180,180) [Kx1]
%
%   When scmpar.ScmOptions is 'polarized' (in addition to scmpar.ScmOptions='none'):
%
%   subpath_phases  - final 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.
                      %%%%%%%這里假設對于每一條path，其所有20條subpath的XPD是相同的%%%%%%
%
%
%   [H ...]=SCM(...,INIT_VALUES) uses initial values given in the struct
%   INIT_VALUES, instead of random parameter generation. INIT_VALUES has
%   the same format as BULKPAR, except that SUBPATH_PHASES are now the
%   initial phases. Also, time sampling intervals (delta_t) are not used
%   (they are recalculated for every call of SCM).
%
%   The 'far scatterer clusters' option [1, Sec. 5.5.2] is not currently
%   supported. The SCM options are mutually exclusive, i.e. one cannot, for
%   instance, choose 'polarized' and 'los' simultaneously.
%
%   Examples:
%       % to generate matrices for 10 links with default parameters
%       H=scm(scmparset,linkparset(10),antparset);
%       % to generate matrices for 'urban_macro' scenario
%       scmpar=scmparset;scmpar.Scenario='urban_macro';
%       H=scm(scmpar,linkparset(10),antparset);
%
%   Ref. [1]: 3GPP TR 25.996 v6.1.0 (2003-09)
%
%   See also SCMPARSET, LINKPARSET, ANTPARSET

%   Authors: Jari Salo (HUT), Giovanni Del Galdo (TUI), Pekka Ky鰏ti (EBIT), 
%   Daniela Laselva (EBIT), Marko Milojevic (TUI), Christian Schneider (TUI)
%   $Revision: 0.34$  $Date: Dec 12, 2004$


% Note: all units are in degrees, meters, Hertz (1/s) and meters/second (m/s)





ni=nargin;
if (ni<3 || ni>4)
    error('SCM requires three or four input arguments !')
end



% SCM parameters, common to all links
Scenario=scmpar.Scenario;
SampleDensity=scmpar.SampleDensity;
NumTimeSamples=scmpar.NumTimeSamples;
N=scmpar.NumPaths;
M=scmpar.NumSubPathsPerPath;
CenterFrequency=scmpar.CenterFrequency;
ScmOptions=scmpar.ScmOptions;
DelaySamplingInterval=scmpar.DelaySamplingInterval;
PathLossModel=scmpar.PathLossModel;
RandomSeed=scmpar.RandomSeed;
UniformTimeSampling=scmpar.UniformTimeSampling;
PathLossModelUsed=scmpar.PathLossModelUsed;
ShadowingModelUsed=scmpar.ShadowingModelUsed;
AnsiC_core=scmpar.AnsiC_core;
LookUpTable=scmpar.LookUpTable;

% antenna parameters 
BsGainPattern=antpar.BsGainPattern;
BsGainAnglesAz=antpar.BsGainAnglesAz;
BsElementPosition=antpar.BsElementPosition;
MsGainPattern=antpar.MsGainPattern;
MsGainAnglesAz=antpar.MsGainAnglesAz;
MsElementPosition=antpar.MsElementPosition;
InterpFunction=antpar.InterpFunction;
InterpMethod=antpar.InterpMethod;   

% link parameters
MsBsDistance=linkpar.MsBsDistance;
ThetaBs=linkpar.ThetaBs;
ThetaMs=linkpar.ThetaMs;
OmegaMs=linkpar.OmegaMs; 
MsVelocity=linkpar.MsVelocity; 
MsDirection=linkpar.MsDirection;
MsHeight=linkpar.MsHeight; 
BsHeight=linkpar.BsHeight; 
MsNumber=linkpar.MsNumber;


% check that the scenario is a valid string
if(any(strcmpi(Scenario,{'suburban_macro','urban_macro','urban_micro'}))==0) 
    error('scmpar.Scenario must be ''suburban_macro'', ''urban_macro'', or ''urban_micro''')
end

% check that the ScmOptions is a valid string
if(any(strcmpi(ScmOptions,{'none','polarized','los','urban_canyon'}))==0) 
    error('scmpar.Scmoptions must be ''none'', ''polarized'', ''los'', or ''urban_canyon'' ')
end


% check that SCM options comply with the selected scenario
if (strcmpi(ScmOptions,'urban_canyon')==1 && strcmpi(Scenario,'suburban_macro')==1 )
    scmpar.Scenario='urban_macro';
    warning('MATLAB:UrbanCanyonWrongScenario','Urban canyon option cannot be selected with "suburban_macro" -> scenario changed to "urban_macro"')
end

if (strcmp(ScmOptions,'los')==1 && strcmp(Scenario,'urban_micro')==0 )
    scmpar.Scenario='urban_micro';
    warning('MATLAB:LineOfSightWrongScenario','LOS option can only be selected with "urban_micro" -> scenario changed to "urban_micro"')
end


% extract the number of links links的數目必須和LINKPAR中給出的各參數的長度相等
NumLinks=length(MsBsDistance);

% Check that the struct linkpar has the same number of parameters in
% each of its fields. This is also the number of links/users.
if (    NumLinks ~= length(ThetaBs)     ||... 
        NumLinks ~= length(ThetaMs)     ||...
        NumLinks ~= length(OmegaMs)     ||...
        NumLinks ~= length(MsVelocity)  ||...
        NumLinks ~= length(MsDirection) ||...
        NumLinks ~= length(MsHeight)    ||...
        NumLinks ~= length(BsHeight)    ||...
        NumLinks ~= length(MsNumber))
    error('All fields in input struct LINKPAR must be of same size!')
end



% Set random seeds if given 
if (isempty(RandomSeed)==0)
    rand('state',RandomSeed);
    randn('state',RandomSeed);
end



% determine the size of the MIMO system
% S - number of BS array antenna elements
if (numel(BsGainPattern)==1)
    S=scmpar.NumBsElements;
else
    S=size(BsGainPattern,1);
end

% U - number of MS array antenna elements
if (numel(MsGainPattern)==1)
    U=scmpar.NumMsElements;
else
    U=size(MsGainPattern,1);  %%%return the number of rows%%%
end

% check that element displacement vector is of right size
if (length(BsElementPosition)~=S && length(BsElementPosition)~=1)
    error('antpar.BsElementPosition has wrong size!')
end

if (length(MsElementPosition)~=U && length(MsElementPosition)~=1)
    error('antpar.MsElementPosition has wrong size!')
end


% check that LUT size is a power-of-two     %%%Look Up Table，用于ANSI-C中計算余弦函數查表%%%
%   
if (strcmpi(AnsiC_core,'yes')==1)
    if (LookUpTable>0)
        if (2^nextpow2(LookUpTable)-LookUpTable~=0)
            scmpar.LookUpTable=2^nextpow2(LookUpTable);
            warning('MATLAB:LUTSizeChanged',['scmpar.LookUpTable is not a power-of-2: size changed to ' num2str(scmpar.LookUpTable) '.'])
        end
    end
end


% These features are not included in this version, so they are fixed
FixedPdpUsed='no'; FixedAnglesUsed='no';
if (strcmpi(FixedPdpUsed,'yes')==1 && N~=6)
    scmpar.NumPaths=6; N=6;