close all
clear all
rng(0);
% if exist('viewer','var') && isvalid(viewer) % viewer handle exists and viewer window is open
%     clearMap(viewer);
% else
%     viewer = siteviewer("Basemap","openstreetmap","Buildings","hongkong.osm");    
% end
simPar = setupSimParameters();
mx_number = simPar.N;
% mz_number = 50;
% my_number = 50;
% threshold = 0.2;
% threhold = 1e-3;
% mulength = 40;
% mz_total = -mz_number:mz_number;
% my_total = -my_number:my_number;
lam = simPar.lambda;
dy = lam/8; dz =  dy;
% iter_test = 30;
% R = zeros(2*my_number+1,2*mz_number+1);
% N = (2*mz_number+1)*(2*my_number+1);

bsPosition = [22.287495, 114.140706]; % lat, lon
bsAntSize = [mx_number 1];                    % number of rows and columns in rectangular array (base station)
bsArrayOrientation = [0 0].';       % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards) in deg
uePosition = [22.287521+rand(1)*0.000008,114.140696+rand(1)*20*0.000001];  % lat, lon
ueAntSize = [1 1];                    % number of rows and columns in rectangular array (UE).
ueArrayOrientation = [0 0].';      % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards)  in deg
reflectionsOrder = 2;                 % number of reflections for ray tracing analysis (0 for LOS)
bsSite = txsite("Name","Base station", ...
    "Latitude",bsPosition(1),"Longitude",bsPosition(2),...
    "AntennaAngle",bsArrayOrientation(1:2),...
    "AntennaHeight",1,...  % in m
    "TransmitterFrequency",3e8/lam);

ueSite = rxsite("Name","UE", ...
    "Latitude",uePosition(1),"Longitude",uePosition(2),...
    "AntennaHeight",1,... % in m
    "AntennaAngle",ueArrayOrientation(1:2));

% bsSite = txsite("Name","Base station", ...
%     "Latitude",bsPosition(1),"Longitude",bsPosition(2),...
%     "AntennaHeight",0,...  % in m
%     "TransmitterFrequency",3e8/lam);
% 
% ueSite = rxsite("Name","UE", ...
%     "Latitude",uePosition(1),"Longitude",uePosition(2),...
%     "AntennaHeight",0);

% show(bsSite);
% show(ueSite);

pm = propagationModel("raytracing","Method","sbr","MaxNumReflections",reflectionsOrder);
rays = raytrace(bsSite,ueSite,pm,"Type","pathloss");

% % plot(rays{1})
% 
pathToAs = [rays{1}.PropagationDelay]-min([rays{1}.PropagationDelay]);  % Time of arrival of each ray (normalized to 0 sec)
avgPathGains  = -[rays{1}.PathLoss];                                    % Average path gains of each ray
pathAoDs = [rays{1}.AngleOfDeparture];                                  % AoD of each ray
pathAoAs = [rays{1}.AngleOfArrival];                                    % AoA of each ray
isLOS = any([rays{1}.LineOfSight]);                                     % Line of sight flag
% 
% % for i = 1:40
cdl_model = nrCDLChannel('ChannelFiltering', false);
cdl_model.DelayProfile = 'Custom';
cdl_model.PathDelays = pathToAs;
cdl_model.AveragePathGains = avgPathGains;%!
cdl_model.AnglesAoD = pathAoDs(1,:);       % azimuth of departure
cdl_model.AnglesZoD = 90-pathAoDs(2,:);    % channel uses zenith angle, rays use elevation
cdl_model.AnglesAoA = pathAoAs(1,:);       % azimuth of arrival
cdl_model.AnglesZoA = 90-pathAoAs(2,:);    % channel uses zenith angle, rays use elevation
cdl_model.HasLOSCluster = isLOS;
% % cdl_model.MaximumDopplerShift = 300;
% % cdl_model.CarrierFrequency = 3e8/lam;
% % cdl_model.SampleRate = 10e2;
% %cdl_model.XPR = 0;
% % cdl_model.SampleDensity = Inf;
% % cdl_model.NumTimeSamples = 1;
% % cdl_model.TransmitAntennaArray.Size = [1, 1, 1, 1, 1]; % [M N P Mg Ng]
% % cdl_model.ReceiveAntennaArray.Size = [2*my_number+1, 2*mz_number+1, 1, 1, 1];
% % cdl_model.ReceiveAntennaArray.PolarizationAngles = [45,0];%?
% % cdl_model.ReceiveAntennaArray.ElementSpacing = [dy/lam,dz/lam,1,1];%?
% %cdl_model.ChannelFiltering = false;
% % cdl_model.Seed = randi(100*i);
% % [pathGains, ~] = cdl_model();
% % h = reshape(sum(pathGains(1, :, :, :), 2), [2*my_number+1, 2*mz_number+1]);
%     cdl_model.MaximumDopplerShift       = 2*simPar.v/simPar.lambda;
%     cdl_model.UTDirectionOfTravel       = [0;0];
%     cdl_model.CarrierFrequency          = simPar.fc;
%     cdl_model.SampleRate                = 1/simPar.Ts;
%     cdl_model.SampleDensity             = Inf;
%     cdl_model.NumTimeSamples            = simPar.T+simPar.T_pre;
%     cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
%     cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
%     cdl_model.Seed                      = randi(10000);
%     [pathGains, ~] = cdl_model();
%     for idx = 1:simPar.T+simPar.T_pre
%          h1(:,idx) = reshape(sum(pathGains(idx, :, :, :), 2), [simPar.N, 1]);
%     end
%     h = h1;
     cdl_model = nrCDLChannel('ChannelFiltering', false);
    cdl_model.MaximumDopplerShift       = 2*simPar.v/simPar.lambda;
    cdl_model.UTDirectionOfTravel       = [0;0];
    cdl_model.CarrierFrequency          = simPar.fc;
    cdl_model.SampleRate                = 1/simPar.Ts;
    cdl_model.SampleDensity             = Inf;
    cdl_model.NumTimeSamples            = simPar.T+simPar.T_pre;
    cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
    cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
    cdl_model.Seed                      = randi(10000);
    [pathGains, ~] = cdl_model();
    for idx = 1:simPar.T+simPar.T_pre
         h1(:,idx) = reshape(sum(pathGains(idx, :, :, :), 2), [simPar.N, 1]);
    end
    h = h1;
figure('color', [1 1 1]);

% C = linspecer(5);
C = [   0.0000    0.0000    0.0000
        0.2000    0.4000    0.8494
        0.3718    0.7176    0.4000
        1.0000    0.5482    0.1000
        0.9650    0.1110    0.1330
        0.8650    0.1110    0.8330
        0.0000    0.0000    0.9000
    ]; 
plot(real(h1(3,:)),       'Marker', 'diamond', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(1, :));
% K1_temp = dftmtx(2*my_number+1);
% K2_temp = dftmtx(2*mz_number+1);
% h_fft = K1_temp*h*K2_temp';
% h_fft = abs(h_fft);
%R = R+h*h(my_number+1,mz_number+1)';
% R = R+h*h(my_number+1,1)';
% end
%[x,y]=meshgrid(1:(2*my_number+1),1:(2*mz_number+1));
% [x,y]=meshgrid(my_total,mz_total);
%displayChannel(cdl_model,'LinkEnd','Rx')
% [x,y]=meshgrid(1:(2*my_number+1),1:(2*mz_number+1));
% pcolor(x,y,abs(R));
% figure
% DFT_Rela = abs(K1_temp*h*K2_temp');
% pcolor(x,y,abs(R));
% shading interp; 
% xlabel('y axis','Interpreter','latex','FontSize',14);
% ylabel('z axis','Interpreter','latex','FontSize',14);
% Rela_cdl = R;
% save('rela_ray_tracing','Rela_cdl');