clear all
close all
clc
% Main_Predictor_CrossingBridges_BJTU is the main function that is able to
% predict the propagation (received power and loss) in the area with
% presence of semi-closed obstacles, in this case, the crossing bridges, 
% by using empirical models based on multi-edge diffraction_CBGfraction models, 
% in this case, the Deygout model with weight coefficients estimated by 
% measurement data from high-speed railway at 930 MHz in China.   
% INPUT:
% 'Case 1'or 'Case 2' - system setup information including many parameters
% 'RasterDatabase_map' - digital elevation model (DEM) raster database of 
% the map in the related area
% 'RasterDatabase_map' - Measurement data (received power) that can be used
% to estimate weight coefficients in the model or to validate the model
% 'gain_M45' - antenna pattern that provides antenna gain at every direction
%
% OUTPUT:
% Propagation loss and received power of measurement and various models,
% such as free space model, Hata model, and Hata model plus Deygout model,
% ICB model and CBG model.
% 
% Author: Dr. Ke Guan
% Date: 15.11.2016

%%
clc
clear
load('M45'); % load transmitting (Tx) antenna pattern

scenario = 'Case 2'; 
% 'Case 1' -- 10 crossing bridges, corresponding to the case of TABLE V, 
% Fig. 12 and Fig. 13 in K. Guan, Z. Zhong, B. Ai and T. Kuerner, 
% "Propagation Measurements and Modeling of Crossing Bridges on High-Speed 
% Railway at 930 MHz," in IEEE Transactions on Vehicular Technology, 
% vol. 63, no. 2, pp. 502-517, Feb. 2014.  
% 'Case 2' -- 4 independent crossing bridges, corresponding to the "Case 1" of TABLE I, 
% Fig. 4 in K. Guan, Z. Zhong, B. Ai and T. Kuerner, 
% "Propagation Measurements and Modeling of Crossing Bridges on High-Speed 
% Railway at 930 MHz," in IEEE Transactions on Vehicular Technology, 
% vol. 63, no. 2, pp. 502-517, Feb. 2014.  

switch scenario
    case 'Case 1'
        load('measurementData_BJTU') % load the measured result
        load('RasterDatabase_map') % load map information
        f=932.4; % carrier frequency [MHz] 
        Tz=34; %Tx height (relative to the ground) [m]
        theta_Tx =3; % elevation (tilt) angle of transmitting antenna [degree], 
        % 0 degree means that the antenna beam is towards the horizontal direction; 
        % 3 degrees means the beam is 3 degrees lower than the horizontal direction, 
        % which is a common setup in cell communications. 
        phi_Tx = 109; % azimuth angle of transmitting antenna [degree], 
        % 0 degree norward, clockwise rotation
        br_height = [15.15,17.93,14.04,16.00,17.11,...
                     10.50,11.95,11.92,11.00, 8.85];
                 % heights of ten crossing bridges relative to the track 
                 % in this scenario [m]
       % weight coefficients in the semi-closed obstacle models
       A_A = 0.51; % alpha in zone A
       B_A = 1-A_A; % beta in Zone A
       A_B = 0.68; % alpha in zone B
       B_B = 1-A_B; % beta in Zone B
       A_C = 0.55; % alpha in zone C
       B_C = 1-A_C; % beta in Zone C
       A_D = 0; % alpha in zone D
       B_D = 1-A_D; % beta in zone D
    case 'Case 2'
        load('measurementData_BJTU_case2') % load the measured result
        load('RasterDatabase_map_case2') % load map information
        f=930.2; % carrier frequency [MHz] 
        Tz=33; %Tx height (relative to the ground) [m]
        theta_Tx =4; % elevation (tilt) angle of transmitting antenna [degree], 
        % 0 degree means that the antenna beam is towards the horizontal direction; 
        % 3 degrees means the beam is 3 degrees lower than the horizontal direction, 
        % which is a common setup in cell communications. 
        phi_Tx = 110; % azimuth angle of transmitting antenna [degree], 
        % 0 degree norward, clockwise rotation
        br_height = [1,    1,    1,    1,    1,...
                     1,18.87,17.11,12.36,12.24];
                 % heights of ten crossing bridges relative to the track 
                 % in this scenario [m]
       % weight coefficients in the semi-closed obstacle models
       A_A = 0.51; % alpha in zone A
       B_A = 1-A_A; % beta in Zone A
       A_B = 0.68; % alpha in zone B
       B_B = 1-A_B; % beta in Zone B
       A_C = 0.55; % alpha in zone C
       B_C = 1-A_C; % beta in Zone C
       A_D = 0; % alpha in zone D
       B_D = 1-A_D; % beta in zone D
end

%% simulation parameters
c = 3*10^8;     % light speed [m/s]
w = c/(f*10^6); % wavelength [m]
res = [10 10];  % resolution of the raster database, which is 10m x 10m 
Rx = 4.8;       % Rx height [m]
Faktor = 1;     % a factor that can be used to decrease the resolution, for possible faster calculation
[row1,col1] = size(HEIGHT); % absolute height of the ground and track, in [m]
height = HEIGHT(1:Faktor:row1,1:Faktor:col1);
clutter = CLUTTER(1:Faktor:row1,1:Faktor:col1);
[row,col] = size(height);

data_mea = data_cal(1:end,3); % data_cal(longitude, latitude, received power [dBm])

%% calculation of the propagation loss by using Hata model plus CBD model, ICB model, Deygout model or no diffraction 
% transmitting antenna index
Xindex = bs_y;
Yindex = bs_x;
a=[Yindex/Faktor Xindex/Faktor];

tic;    
bar = waitbar(0,'Prediction at every point on the track...');
for i=1:length(Rx_x);

    [dis(i,1)] = ...
        Hata_calculator_without_diffraction_BJTU(a,[Rx_x(i) Rx_y(i)], height, clutter, res.*Faktor, Tz, Rx, f, br_height);

    waitbar(i/size(Rx_x,1),bar) 
end
close(bar)   
toc;

%% add the antenna pattern to get the Tx gain at any direction
theta_Rx = atan2(abs(Tz-Rx),(dis*1000))*180/pi;

for i = 1:length(theta_Rx)
    if theta_Rx(i,1) > theta_Tx
        theta_Gain(i,1) = round(90+(theta_Rx(i,1) - theta_Tx))+1;
    else
        theta_Gain(i,1) = round(90-(theta_Tx - theta_Rx(i,1)))+1;
    end
    if phi_Rx(i,1) > phi_Tx
        phi_Gain(i,1) = round(phi_Rx(i,1) - phi_Tx)+1;
    else
        phi_Gain(i,1) = round(360-(phi_Tx - phi_Rx(i,1)))+1;
    end
    Gain_Tx(i,1) = gain(theta_Gain(i,1),phi_Gain(i,1));
end

%% measurement parametersment campaign
% Transmitting power = 40 dBm; transmitting antenna gain =17 dBi; 
sysFactor = 40+17-3.3-3.5-3.3;  % system gain and loss in the measure
% loss of feeder: 3.3 dB + 3.5 dB; 
% loss of coupling: 3.3 dB;
% gain of the receiving anttenna: 0 dB (omni-directional antenna);
mea_PL = sysFactor+Gain_Tx-data_cal(:,3); % measured propagation loss


%% Free space propagation model
PL_free = -(-32.4-20*log10(f)-10.*log10(dis))+Gain_Tx;% path loss of using free space propagation model
Pred_power_free = - PL_free+sysFactor+Gain_Tx; % predicted power of using free space propagation model

%% Plot results

figure % figure for propagation loss
plot(dis*1000,mea_PL,'r');
hold on
plot(dis*1000,PL_free(:,1),'g');
hold on
legend('Změřený útlum šířením','Model útlumu šířením pro volný prostor')
    
title('Srovnání měřených hodnot útlumu s modelem pro volný prostor')
xlabel('Vzdálenost vysílače a přijímače [m]')
ylabel('Útlum šířením [dB]')
grid on

clear qwe coef1 coef2 coef3
coef3=0;
for qwe = 1:length(dis)
    if (dis(qwe) > 0.45 & dis(qwe) < 1.2)
        coef3=coef3+1;
    end
end
clear qwe

dis_2 = zeros(coef3,1);
data_cal_2 = zeros(coef3,1);
coef1 = 0;
coef2 = 0;

for qwe = 1:length(dis)
    coef2=coef2+1;
    if (dis(qwe) > 0.45 & dis(qwe) < 1.2)
        coef1 = coef1+1;
        dis_2(coef1) = dis(coef2) ;
        data_cal_2(coef1) = data_cal(coef2,3);
    end
end

figure; % figure for received power
plot(dis_2*1000,data_cal_2,'r');

hold on

legend('Měřený přijatý výkon')
title('Měření přijatého výkonu')
xlabel('Vzdálenost vysílače a přijímače [m]')
ylabel('Přijatý výkon [dBm]')
grid on

