A Discrete-Event Network Simulator
API
three-gpp-channel-model.cc
Go to the documentation of this file.
1 /* -*- Mode: C++; c-file-style: "gnu"; indent-tabs-mode:nil; -*- */
2 /*
3  * Copyright (c) 2019 SIGNET Lab, Department of Information Engineering,
4  * University of Padova
5  * Copyright (c) 2015, NYU WIRELESS, Tandon School of Engineering,
6  * New York University
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License version 2 as
10  * published by the Free Software Foundation;
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
24 #include "ns3/log.h"
25 #include "ns3/three-gpp-antenna-array-model.h"
26 #include "ns3/node.h"
27 #include "ns3/double.h"
28 #include "ns3/string.h"
29 #include "ns3/integer.h"
30 #include <algorithm>
31 #include <random>
32 #include "ns3/log.h"
33 #include <ns3/simulator.h>
34 #include "ns3/mobility-model.h"
35 #include "ns3/pointer.h"
36 
37 namespace ns3 {
38 
39 NS_LOG_COMPONENT_DEFINE ("ThreeGppChannelModel");
40 
41 NS_OBJECT_ENSURE_REGISTERED (ThreeGppChannelModel);
42 
43 //Table 7.5-3: Ray offset angles within a cluster, given for rms angle spread normalized to 1.
44 static const double offSetAlpha[20] = {
45  0.0447,-0.0447,0.1413,-0.1413,0.2492,-0.2492,0.3715,-0.3715,0.5129,-0.5129,0.6797,-0.6797,0.8844,-0.8844,1.1481,-1.1481,1.5195,-1.5195,2.1551,-2.1551
46 };
47 
48 /*
49  * The cross correlation matrix is constructed according to table 7.5-6.
50  * All the square root matrix is being generated using the Cholesky decomposition
51  * and following the order of [SF,K,DS,ASD,ASA,ZSD,ZSA].
52  * The parameter K is ignored in NLOS.
53  *
54  * The Matlab file to generate the matrices can be found in
55  * https://github.com/nyuwireless-unipd/ns3-mmwave/blob/master/src/mmwave/model/BeamFormingMatrix/SqrtMatrix.m
56  *
57  */
58 static const double sqrtC_RMa_LOS[7][7] = {
59  {1, 0, 0, 0, 0, 0, 0},
60  {0, 1, 0, 0, 0, 0, 0},
61  {-0.5, 0, 0.866025, 0, 0, 0, 0},
62  {0, 0, 0, 1, 0, 0, 0},
63  {0, 0, 0, 0, 1, 0, 0},
64  {0.01, 0, -0.0519615, 0.73, -0.2, 0.651383, 0},
65  {-0.17, -0.02, 0.21362, -0.14, 0.24, 0.142773, 0.909661},
66 };
67 
68 static const double sqrtC_RMa_NLOS[6][6] = {
69  {1, 0, 0, 0, 0, 0},
70  {-0.5, 0.866025, 0, 0, 0, 0},
71  {0.6, -0.11547, 0.791623, 0, 0, 0},
72  {0, 0, 0, 1, 0, 0},
73  {-0.04, -0.138564, 0.540662, -0.18, 0.809003, 0},
74  {-0.25, -0.606218, -0.240013, 0.26, -0.231685, 0.625392},
75 };
76 
77 static const double sqrtC_RMa_O2I[6][6] = {
78  {1, 0, 0, 0, 0, 0},
79  {0, 1, 0, 0, 0, 0},
80  {0, 0, 1, 0, 0, 0},
81  {0, 0, -0.7, 0.714143, 0, 0},
82  {0, 0, 0.66, -0.123225, 0.741091, 0},
83  {0, 0, 0.47, 0.152631, -0.393194, 0.775373},
84 };
85 
86 static const double sqrtC_UMa_LOS[7][7] = {
87  {1, 0, 0, 0, 0, 0, 0},
88  {0, 1, 0, 0, 0, 0, 0},
89  {-0.4, -0.4, 0.824621, 0, 0, 0, 0},
90  {-0.5, 0, 0.242536, 0.83137, 0, 0, 0},
91  {-0.5, -0.2, 0.630593, -0.484671, 0.278293, 0, 0},
92  {0, 0, -0.242536, 0.672172, 0.642214, 0.27735, 0},
93  {-0.8, 0, -0.388057, -0.367926, 0.238537, -3.58949e-15, 0.130931},
94 };
95 
96 
97 static const double sqrtC_UMa_NLOS[6][6] = {
98  {1, 0, 0, 0, 0, 0},
99  {-0.4, 0.916515, 0, 0, 0, 0},
100  {-0.6, 0.174574, 0.78072, 0, 0, 0},
101  {0, 0.654654, 0.365963, 0.661438, 0, 0},
102  {0, -0.545545, 0.762422, 0.118114, 0.327327, 0},
103  {-0.4, -0.174574, -0.396459, 0.392138, 0.49099, 0.507445},
104 };
105 
106 static const double sqrtC_UMa_O2I[6][6] = {
107  {1, 0, 0, 0, 0, 0},
108  {-0.5, 0.866025, 0, 0, 0, 0},
109  {0.2, 0.57735, 0.791623, 0, 0, 0},
110  {0, 0.46188, -0.336861, 0.820482, 0, 0},
111  {0, -0.69282, 0.252646, 0.493742, 0.460857, 0},
112  {0, -0.23094, 0.16843, 0.808554, -0.220827, 0.464515},
113 
114 };
115 
116 static const double sqrtC_UMi_LOS[7][7] = {
117  {1, 0, 0, 0, 0, 0, 0},
118  {0.5, 0.866025, 0, 0, 0, 0, 0},
119  {-0.4, -0.57735, 0.711805, 0, 0, 0, 0},
120  {-0.5, 0.057735, 0.468293, 0.726201, 0, 0, 0},
121  {-0.4, -0.11547, 0.805464, -0.23482, 0.350363, 0, 0},
122  {0, 0, 0, 0.688514, 0.461454, 0.559471, 0},
123  {0, 0, 0.280976, 0.231921, -0.490509, 0.11916, 0.782603},
124 };
125 
126 static const double sqrtC_UMi_NLOS[6][6] = {
127  {1, 0, 0, 0, 0, 0},
128  {-0.7, 0.714143, 0, 0, 0, 0},
129  {0, 0, 1, 0, 0, 0},
130  {-0.4, 0.168034, 0, 0.90098, 0, 0},
131  {0, -0.70014, 0.5, 0.130577, 0.4927, 0},
132  {0, 0, 0.5, 0.221981, -0.566238, 0.616522},
133 };
134 
135 static const double sqrtC_UMi_O2I[6][6] = {
136  {1, 0, 0, 0, 0, 0},
137  {-0.5, 0.866025, 0, 0, 0, 0},
138  {0.2, 0.57735, 0.791623, 0, 0, 0},
139  {0, 0.46188, -0.336861, 0.820482, 0, 0},
140  {0, -0.69282, 0.252646, 0.493742, 0.460857, 0},
141  {0, -0.23094, 0.16843, 0.808554, -0.220827, 0.464515},
142 };
143 
144 static const double sqrtC_office_LOS[7][7] = {
145  {1, 0, 0, 0, 0, 0, 0},
146  {0.5, 0.866025, 0, 0, 0, 0, 0},
147  {-0.8, -0.11547, 0.588784, 0, 0, 0, 0},
148  {-0.4, 0.23094, 0.520847, 0.717903, 0, 0, 0},
149  {-0.5, 0.288675, 0.73598, -0.348236, 0.0610847, 0, 0},
150  {0.2, -0.11547, 0.418943, 0.541106, 0.219905, 0.655744, 0},
151  {0.3, -0.057735, 0.73598, -0.348236, 0.0610847, -0.304997, 0.383375},
152 };
153 
154 static const double sqrtC_office_NLOS[6][6] = {
155  {1, 0, 0, 0, 0, 0},
156  {-0.5, 0.866025, 0, 0, 0, 0},
157  {0, 0.46188, 0.886942, 0, 0, 0},
158  {-0.4, -0.23094, 0.120263, 0.878751, 0, 0},
159  {0, -0.311769, 0.55697, -0.249198, 0.728344, 0},
160  {0, -0.069282, 0.295397, 0.430696, 0.468462, 0.709214},
161 };
162 
164 {
165  NS_LOG_FUNCTION (this);
166  m_uniformRv = CreateObject<UniformRandomVariable> ();
167 
168  m_normalRv = CreateObject<NormalRandomVariable> ();
169  m_normalRv->SetAttribute ("Mean", DoubleValue (0.0));
170  m_normalRv->SetAttribute ("Variance", DoubleValue (1.0));
171 }
172 
174 {
175  NS_LOG_FUNCTION (this);
176 }
177 
178 TypeId
180 {
181  static TypeId tid = TypeId ("ns3::ThreeGppChannelModel")
182  .SetParent<Object> ()
183  .SetGroupName ("Spectrum")
184  .AddConstructor<ThreeGppChannelModel> ()
185  .AddAttribute ("Frequency",
186  "The operating Frequency in Hz",
187  DoubleValue (500.0e6),
190  MakeDoubleChecker<double> ())
191  .AddAttribute ("Scenario",
192  "The 3GPP scenario (RMa, UMa, UMi-StreetCanyon, InH-OfficeOpen, InH-OfficeMixed)",
193  StringValue ("UMa"),
197  .AddAttribute ("ChannelConditionModel",
198  "Pointer to the channel condition model",
199  PointerValue (),
202  MakePointerChecker<ChannelConditionModel> ())
203  .AddAttribute ("UpdatePeriod",
204  "Specify the channel coherence time",
205  TimeValue (MilliSeconds (0)),
207  MakeTimeChecker ())
208  // attributes for the blockage model
209  .AddAttribute ("Blockage",
210  "Enable blockage model A (sec 7.6.4.1)",
211  BooleanValue (false),
214  .AddAttribute ("NumNonselfBlocking",
215  "number of non-self-blocking regions",
216  IntegerValue (4),
218  MakeIntegerChecker<uint16_t> ())
219  .AddAttribute ("PortraitMode",
220  "true for portrait mode, false for landscape mode",
221  BooleanValue (true),
224  .AddAttribute ("BlockerSpeed",
225  "The speed of moving blockers, the unit is m/s",
226  DoubleValue (1),
228  MakeDoubleChecker<double> ())
229  ;
230  return tid;
231 }
232 
233 void
235 {
236  NS_LOG_FUNCTION (this);
237  m_channelConditionModel = model;
238 }
239 
242 {
243  NS_LOG_FUNCTION (this);
245 }
246 
247 void
249 {
250  NS_LOG_FUNCTION (this);
251  NS_ASSERT_MSG (f >= 500.0e6 && f <= 100.0e9, "Frequency should be between 0.5 and 100 GHz but is " << f);
252  m_frequency = f;
253 }
254 
255 double
257 {
258  NS_LOG_FUNCTION (this);
259  return m_frequency;
260 }
261 
262 void
263 ThreeGppChannelModel::SetScenario (const std::string &scenario)
264 {
265  NS_LOG_FUNCTION (this);
266  NS_ASSERT_MSG (scenario == "RMa" || scenario == "UMa" || scenario == "UMi-StreetCanyon" ||
267  scenario == "InH-OfficeOpen" || scenario == "InH-OfficeMixed",
268  "Unknown scenario, choose between RMa, UMa, UMi-StreetCanyon, InH-OfficeOpen or InH-OfficeMixed");
269  m_scenario = scenario;
270 }
271 
272 std::string
274 {
275  NS_LOG_FUNCTION (this);
276  return m_scenario;
277 }
278 
280 ThreeGppChannelModel::GetThreeGppTable (bool los, bool o2i, double hBS, double hUT, double distance2D) const
281 {
282  NS_LOG_FUNCTION (this);
283 
284  double fcGHz = m_frequency / 1e9;
285  Ptr<ParamsTable> table3gpp = Create<ParamsTable> ();
286  // table3gpp includes the following parameters:
287  // numOfCluster, raysPerCluster, uLgDS, sigLgDS, uLgASD, sigLgASD,
288  // uLgASA, sigLgASA, uLgZSA, sigLgZSA, uLgZSD, sigLgZSD, offsetZOD,
289  // cDS, cASD, cASA, cZSA, uK, sigK, rTau, uXpr, sigXpr, shadowingStd
290 
291  // In NLOS case, parameter uK and sigK are not used and they are set to 0
292  if (m_scenario == "RMa")
293  {
294  if (los && !o2i)
295  {
296  // 3GPP mentioned that 3.91 ns should be used when the Cluster DS (cDS)
297  // entry is N/A.
298  table3gpp->m_numOfCluster = 11;
299  table3gpp->m_raysPerCluster = 20;
300  table3gpp->m_uLgDS = -7.49;
301  table3gpp->m_sigLgDS = 0.55;
302  table3gpp->m_uLgASD = 0.90;
303  table3gpp->m_sigLgASD = 0.38;
304  table3gpp->m_uLgASA = 1.52;
305  table3gpp->m_sigLgASA = 0.24;
306  table3gpp->m_uLgZSA = 0.47;
307  table3gpp->m_sigLgZSA = 0.40;
308  table3gpp->m_uLgZSD = 0.34;
309  table3gpp->m_sigLgZSD = std::max (-1.0, -0.17 * (distance2D / 1000) - 0.01 * (hUT - 1.5) + 0.22);
310  table3gpp->m_offsetZOD = 0;
311  table3gpp->m_cDS = 3.91e-9;
312  table3gpp->m_cASD = 2;
313  table3gpp->m_cASA = 3;
314  table3gpp->m_cZSA = 3;
315  table3gpp->m_uK = 7;
316  table3gpp->m_sigK = 4;
317  table3gpp->m_rTau = 3.8;
318  table3gpp->m_uXpr = 12;
319  table3gpp->m_sigXpr = 4;
320  table3gpp->m_perClusterShadowingStd = 3;
321 
322  for (uint8_t row = 0; row < 7; row++)
323  {
324  for (uint8_t column = 0; column < 7; column++)
325  {
326  table3gpp->m_sqrtC[row][column] = sqrtC_RMa_LOS[row][column];
327  }
328  }
329  }
330  else if (!los && !o2i)
331  {
332  table3gpp->m_numOfCluster = 10;
333  table3gpp->m_raysPerCluster = 20;
334  table3gpp->m_uLgDS = -7.43;
335  table3gpp->m_sigLgDS = 0.48;
336  table3gpp->m_uLgASD = 0.95;
337  table3gpp->m_sigLgASD = 0.45;
338  table3gpp->m_uLgASA = 1.52;
339  table3gpp->m_sigLgASA = 0.13;
340  table3gpp->m_uLgZSA = 0.58,
341  table3gpp->m_sigLgZSA = 0.37;
342  table3gpp->m_uLgZSD = std::max (-1.0, -0.19 * (distance2D / 1000) - 0.01 * (hUT - 1.5) + 0.28);
343  table3gpp->m_sigLgZSD = 0.30;
344  table3gpp->m_offsetZOD = atan ((35 - 3.5) / distance2D) - atan ((35 - 1.5) / distance2D);
345  table3gpp->m_cDS = 3.91e-9;
346  table3gpp->m_cASD = 2;
347  table3gpp->m_cASA = 3;
348  table3gpp->m_cZSA = 3;
349  table3gpp->m_uK = 0;
350  table3gpp->m_sigK = 0;
351  table3gpp->m_rTau = 1.7;
352  table3gpp->m_uXpr = 7;
353  table3gpp->m_sigXpr = 3;
354  table3gpp->m_perClusterShadowingStd = 3;
355 
356  for (uint8_t row = 0; row < 6; row++)
357  {
358  for (uint8_t column = 0; column < 6; column++)
359  {
360  table3gpp->m_sqrtC[row][column] = sqrtC_RMa_NLOS[row][column];
361  }
362  }
363  }
364  else // o2i
365  {
366  table3gpp->m_numOfCluster = 10;
367  table3gpp->m_raysPerCluster = 20;
368  table3gpp->m_uLgDS = -7.47;
369  table3gpp->m_sigLgDS = 0.24;
370  table3gpp->m_uLgASD = 0.67;
371  table3gpp->m_sigLgASD = 0.18;
372  table3gpp->m_uLgASA = 1.66;
373  table3gpp->m_sigLgASA = 0.21;
374  table3gpp->m_uLgZSA = 0.93,
375  table3gpp->m_sigLgZSA = 0.22;
376  table3gpp->m_uLgZSD = std::max (-1.0, -0.19 * (distance2D / 1000) - 0.01 * (hUT - 1.5) + 0.28);
377  table3gpp->m_sigLgZSD = 0.30;
378  table3gpp->m_offsetZOD = atan ((35 - 3.5) / distance2D) - atan ((35 - 1.5) / distance2D);
379  table3gpp->m_cDS = 3.91e-9;
380  table3gpp->m_cASD = 2;
381  table3gpp->m_cASA = 3;
382  table3gpp->m_cZSA = 3;
383  table3gpp->m_uK = 0;
384  table3gpp->m_sigK = 0;
385  table3gpp->m_rTau = 1.7;
386  table3gpp->m_uXpr = 7;
387  table3gpp->m_sigXpr = 3;
388  table3gpp->m_perClusterShadowingStd = 3;
389 
390  for (uint8_t row = 0; row < 6; row++)
391  {
392  for (uint8_t column = 0; column < 6; column++)
393  {
394  table3gpp->m_sqrtC[row][column] = sqrtC_RMa_O2I[row][column];
395  }
396  }
397  }
398  }
399  else if (m_scenario == "UMa")
400  {
401  if (los && !o2i)
402  {
403  table3gpp->m_numOfCluster = 12;
404  table3gpp->m_raysPerCluster = 20;
405  table3gpp->m_uLgDS = -6.955 - 0.0963 * log10 (fcGHz);
406  table3gpp->m_sigLgDS = 0.66;
407  table3gpp->m_uLgASD = 1.06 + 0.1114 * log10 (fcGHz);
408  table3gpp->m_sigLgASD = 0.28;
409  table3gpp->m_uLgASA = 1.81;
410  table3gpp->m_sigLgASA = 0.20;
411  table3gpp->m_uLgZSA = 0.95;
412  table3gpp->m_sigLgZSA = 0.16;
413  table3gpp->m_uLgZSD = std::max (-0.5, -2.1 * distance2D / 1000 - 0.01 * (hUT - 1.5) + 0.75);
414  table3gpp->m_sigLgZSD = 0.40;
415  table3gpp->m_offsetZOD = 0;
416  table3gpp->m_cDS = std::max (0.25, -3.4084 * log10 (fcGHz) + 6.5622) * 1e-9;
417  table3gpp->m_cASD = 5;
418  table3gpp->m_cASA = 11;
419  table3gpp->m_cZSA = 7;
420  table3gpp->m_uK = 9;
421  table3gpp->m_sigK = 3.5;
422  table3gpp->m_rTau = 2.5;
423  table3gpp->m_uXpr = 8;
424  table3gpp->m_sigXpr = 4;
425  table3gpp->m_perClusterShadowingStd = 3;
426 
427  for (uint8_t row = 0; row < 7; row++)
428  {
429  for (uint8_t column = 0; column < 7; column++)
430  {
431  table3gpp->m_sqrtC[row][column] = sqrtC_UMa_LOS[row][column];
432  }
433  }
434  }
435  else
436  {
437  double uLgZSD = std::max (-0.5, -2.1 * distance2D / 1000 - 0.01 * (hUT - 1.5) + 0.9);
438 
439  double afc = 0.208 * log10 (fcGHz) - 0.782;
440  double bfc = 25;
441  double cfc = -0.13 * log10 (fcGHz) + 2.03;
442  double efc = 7.66 * log10 (fcGHz) - 5.96;
443 
444  double offsetZOD = efc - std::pow (10, afc * log10 (std::max (bfc,distance2D)) + cfc);
445 
446  if (!los && !o2i)
447  {
448  table3gpp->m_numOfCluster = 20;
449  table3gpp->m_raysPerCluster = 20;
450  table3gpp->m_uLgDS = -6.28 - 0.204 * log10 (fcGHz);
451  table3gpp->m_sigLgDS = 0.39;
452  table3gpp->m_uLgASD = 1.5 - 0.1144 * log10 (fcGHz);
453  table3gpp->m_sigLgASD = 0.28;
454  table3gpp->m_uLgASA = 2.08 - 0.27 * log10 (fcGHz);
455  table3gpp->m_sigLgASA = 0.11;
456  table3gpp->m_uLgZSA = -0.3236 * log10 (fcGHz) + 1.512;
457  table3gpp->m_sigLgZSA = 0.16;
458  table3gpp->m_uLgZSD = uLgZSD;
459  table3gpp->m_sigLgZSD = 0.49;
460  table3gpp->m_offsetZOD = offsetZOD;
461  table3gpp->m_cDS = std::max (0.25, -3.4084 * log10 (fcGHz) + 6.5622) * 1e-9;;
462  table3gpp->m_cASD = 2;
463  table3gpp->m_cASA = 15;
464  table3gpp->m_cZSA = 7;
465  table3gpp->m_uK = 0;
466  table3gpp->m_sigK = 0;
467  table3gpp->m_rTau = 2.3;
468  table3gpp->m_uXpr = 7;
469  table3gpp->m_sigXpr = 3;
470  table3gpp->m_perClusterShadowingStd = 3;
471 
472  for (uint8_t row = 0; row < 6; row++)
473  {
474  for (uint8_t column = 0; column < 6; column++)
475  {
476  table3gpp->m_sqrtC[row][column] = sqrtC_UMa_NLOS[row][column];
477  }
478  }
479  }
480  else //(o2i)
481  {
482  table3gpp->m_numOfCluster = 12;
483  table3gpp->m_raysPerCluster = 20;
484  table3gpp->m_uLgDS = -6.62;
485  table3gpp->m_sigLgDS = 0.32;
486  table3gpp->m_uLgASD = 1.25;
487  table3gpp->m_sigLgASD = 0.42;
488  table3gpp->m_uLgASA = 1.76;
489  table3gpp->m_sigLgASA = 0.16;
490  table3gpp->m_uLgZSA = 1.01;
491  table3gpp->m_sigLgZSA = 0.43;
492  table3gpp->m_uLgZSD = uLgZSD;
493  table3gpp->m_sigLgZSD = 0.49;
494  table3gpp->m_offsetZOD = offsetZOD;
495  table3gpp->m_cDS = 11e-9;
496  table3gpp->m_cASD = 5;
497  table3gpp->m_cASA = 8;
498  table3gpp->m_cZSA = 3;
499  table3gpp->m_uK = 0;
500  table3gpp->m_sigK = 0;
501  table3gpp->m_rTau = 2.2;
502  table3gpp->m_uXpr = 9;
503  table3gpp->m_sigXpr = 5;
504  table3gpp->m_perClusterShadowingStd = 4;
505 
506  for (uint8_t row = 0; row < 6; row++)
507  {
508  for (uint8_t column = 0; column < 6; column++)
509  {
510  table3gpp->m_sqrtC[row][column] = sqrtC_UMa_O2I[row][column];
511  }
512  }
513 
514  }
515 
516  }
517 
518  }
519  else if (m_scenario == "UMi-StreetCanyon")
520  {
521  if (los && !o2i)
522  {
523  table3gpp->m_numOfCluster = 12;
524  table3gpp->m_raysPerCluster = 20;
525  table3gpp->m_uLgDS = -0.24 * log10 (1 + fcGHz) - 7.14;
526  table3gpp->m_sigLgDS = 0.38;
527  table3gpp->m_uLgASD = -0.05 * log10 (1 + fcGHz) + 1.21;
528  table3gpp->m_sigLgASD = 0.41;
529  table3gpp->m_uLgASA = -0.08 * log10 (1 + fcGHz) + 1.73;
530  table3gpp->m_sigLgASA = 0.014 * log10 (1 + fcGHz) + 0.28;
531  table3gpp->m_uLgZSA = -0.1 * log10 (1 + fcGHz) + 0.73;
532  table3gpp->m_sigLgZSA = -0.04 * log10 (1 + fcGHz) + 0.34;
533  table3gpp->m_uLgZSD = std::max (-0.21, -14.8 * distance2D / 1000 + 0.01 * std::abs (hUT - hBS) + 0.83);
534  table3gpp->m_sigLgZSD = 0.35;
535  table3gpp->m_offsetZOD = 0;
536  table3gpp->m_cDS = 5e-9;
537  table3gpp->m_cASD = 3;
538  table3gpp->m_cASA = 17;
539  table3gpp->m_cZSA = 7;
540  table3gpp->m_uK = 9;
541  table3gpp->m_sigK = 5;
542  table3gpp->m_rTau = 3;
543  table3gpp->m_uXpr = 9;
544  table3gpp->m_sigXpr = 3;
545  table3gpp->m_perClusterShadowingStd = 3;
546 
547  for (uint8_t row = 0; row < 7; row++)
548  {
549  for (uint8_t column = 0; column < 7; column++)
550  {
551  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_LOS[row][column];
552  }
553  }
554  }
555  else
556  {
557  double uLgZSD = std::max (-0.5, -3.1 * distance2D / 1000 + 0.01 * std::max (hUT - hBS,0.0) + 0.2);
558  double offsetZOD = -1 * std::pow (10, -1.5 * log10 (std::max (10.0, distance2D)) + 3.3);
559  if (!los && !o2i)
560  {
561  table3gpp->m_numOfCluster = 19;
562  table3gpp->m_raysPerCluster = 20;
563  table3gpp->m_uLgDS = -0.24 * log10 (1 + fcGHz) - 6.83;
564  table3gpp->m_sigLgDS = 0.16 * log10 (1 + fcGHz) + 0.28;
565  table3gpp->m_uLgASD = -0.23 * log10 (1 + fcGHz) + 1.53;
566  table3gpp->m_sigLgASD = 0.11 * log10 (1 + fcGHz) + 0.33;
567  table3gpp->m_uLgASA = -0.08 * log10 (1 + fcGHz) + 1.81;
568  table3gpp->m_sigLgASA = 0.05 * log10 (1 + fcGHz) + 0.3;
569  table3gpp->m_uLgZSA = -0.04 * log10 (1 + fcGHz) + 0.92;
570  table3gpp->m_sigLgZSA = -0.07 * log10 (1 + fcGHz) + 0.41;
571  table3gpp->m_uLgZSD = uLgZSD;
572  table3gpp->m_sigLgZSD = 0.35;
573  table3gpp->m_offsetZOD = offsetZOD;
574  table3gpp->m_cDS = 11e-9;
575  table3gpp->m_cASD = 10;
576  table3gpp->m_cASA = 22;
577  table3gpp->m_cZSA = 7;
578  table3gpp->m_uK = 0;
579  table3gpp->m_sigK = 0;
580  table3gpp->m_rTau = 2.1;
581  table3gpp->m_uXpr = 8;
582  table3gpp->m_sigXpr = 3;
583  table3gpp->m_perClusterShadowingStd = 3;
584 
585  for (uint8_t row = 0; row < 6; row++)
586  {
587  for (uint8_t column = 0; column < 6; column++)
588  {
589  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_NLOS[row][column];
590  }
591  }
592  }
593  else //(o2i)
594  {
595  table3gpp->m_numOfCluster = 12;
596  table3gpp->m_raysPerCluster = 20;
597  table3gpp->m_uLgDS = -6.62;
598  table3gpp->m_sigLgDS = 0.32;
599  table3gpp->m_uLgASD = 1.25;
600  table3gpp->m_sigLgASD = 0.42;
601  table3gpp->m_uLgASA = 1.76;
602  table3gpp->m_sigLgASA = 0.16;
603  table3gpp->m_uLgZSA = 1.01;
604  table3gpp->m_sigLgZSA = 0.43;
605  table3gpp->m_uLgZSD = uLgZSD;
606  table3gpp->m_sigLgZSD = 0.35;
607  table3gpp->m_offsetZOD = offsetZOD;
608  table3gpp->m_cDS = 11e-9;
609  table3gpp->m_cASD = 5;
610  table3gpp->m_cASA = 8;
611  table3gpp->m_cZSA = 3;
612  table3gpp->m_uK = 0;
613  table3gpp->m_sigK = 0;
614  table3gpp->m_rTau = 2.2;
615  table3gpp->m_uXpr = 9;
616  table3gpp->m_sigXpr = 5;
617  table3gpp->m_perClusterShadowingStd = 4;
618 
619  for (uint8_t row = 0; row < 6; row++)
620  {
621  for (uint8_t column = 0; column < 6; column++)
622  {
623  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_O2I[row][column];
624  }
625  }
626  }
627  }
628  }
629  else if (m_scenario == "InH-OfficeMixed"||m_scenario == "InH-OfficeOpen")
630  {
631  NS_ASSERT_MSG (!o2i, "The indoor scenario does out support outdoor to indoor");
632  if (los)
633  {
634  table3gpp->m_numOfCluster = 15;
635  table3gpp->m_raysPerCluster = 20;
636  table3gpp->m_uLgDS = -0.01 * log10 (1 + fcGHz) - 7.692;
637  table3gpp->m_sigLgDS = 0.18;
638  table3gpp->m_uLgASD = 1.60;
639  table3gpp->m_sigLgASD = 0.18;
640  table3gpp->m_uLgASA = -0.19 * log10 (1 + fcGHz) + 1.781;
641  table3gpp->m_sigLgASA = 0.12 * log10 (1 + fcGHz) + 0.119;
642  table3gpp->m_uLgZSA = -0.26 * log10 (1 + fcGHz) + 1.44;
643  table3gpp->m_sigLgZSA = -0.04 * log10 (1 + fcGHz) + 0.264;
644  table3gpp->m_uLgZSD = -1.43 * log10 (1 + fcGHz) + 2.228;
645  table3gpp->m_sigLgZSD = 0.13 * log10 (1 + fcGHz) + 0.30;
646  table3gpp->m_offsetZOD = 0;
647  table3gpp->m_cDS = 3.91e-9;
648  table3gpp->m_cASD = 5;
649  table3gpp->m_cASA = 8;
650  table3gpp->m_cZSA = 9;
651  table3gpp->m_uK = 7;
652  table3gpp->m_sigK = 4;
653  table3gpp->m_rTau = 3.6;
654  table3gpp->m_uXpr = 11;
655  table3gpp->m_sigXpr = 4;
656  table3gpp->m_perClusterShadowingStd = 6;
657 
658  for (uint8_t row = 0; row < 7; row++)
659  {
660  for (uint8_t column = 0; column < 7; column++)
661  {
662  table3gpp->m_sqrtC[row][column] = sqrtC_office_LOS[row][column];
663  }
664  }
665  }
666  else
667  {
668  table3gpp->m_numOfCluster = 19;
669  table3gpp->m_raysPerCluster = 20;
670  table3gpp->m_uLgDS = -0.28 * log10 (1 + fcGHz) - 7.173;
671  table3gpp->m_sigLgDS = 0.1 * log10 (1 + fcGHz) + 0.055;
672  table3gpp->m_uLgASD = 1.62;
673  table3gpp->m_sigLgASD = 0.25;
674  table3gpp->m_uLgASA = -0.11 * log10 (1 + fcGHz) + 1.863;
675  table3gpp->m_sigLgASA = 0.12 * log10 (1 + fcGHz) + 0.059;
676  table3gpp->m_uLgZSA = -0.15 * log10 (1 + fcGHz) + 1.387;
677  table3gpp->m_sigLgZSA = -0.09 * log10 (1 + fcGHz) + 0.746;
678  table3gpp->m_uLgZSD = 1.08;
679  table3gpp->m_sigLgZSD = 0.36;
680  table3gpp->m_offsetZOD = 0;
681  table3gpp->m_cDS = 3.91e-9;
682  table3gpp->m_cASD = 5;
683  table3gpp->m_cASA = 11;
684  table3gpp->m_cZSA = 9;
685  table3gpp->m_uK = 0;
686  table3gpp->m_sigK = 0;
687  table3gpp->m_rTau = 3;
688  table3gpp->m_uXpr = 10;
689  table3gpp->m_sigXpr = 4;
690  table3gpp->m_perClusterShadowingStd = 3;
691 
692  for (uint8_t row = 0; row < 6; row++)
693  {
694  for (uint8_t column = 0; column < 6; column++)
695  {
696  table3gpp->m_sqrtC[row][column] = sqrtC_office_NLOS[row][column];
697  }
698  }
699  }
700  }
701  else
702  {
703  NS_FATAL_ERROR ("unkonw scenarios");
704  }
705 
706  return table3gpp;
707 }
708 
709 bool
711 {
712  NS_LOG_FUNCTION (this);
713 
714  bool update = false;
715 
716  // if the los condition is different the channel has to be updated
717  if (channelMatrix->m_los != los)
718  {
719  NS_LOG_DEBUG ("old los condition " << channelMatrix->m_los << " new los condition " << los);
720  update = true;
721  }
722 
723  // if the coherence time is over the channel has to be updated
724  if (!m_updatePeriod.IsZero () && Simulator::Now() - channelMatrix->m_generatedTime > m_updatePeriod)
725  {
726  NS_LOG_DEBUG ("Generation time " << channelMatrix->m_generatedTime.GetNanoSeconds () << " now " << Simulator::Now ().GetNanoSeconds ());
727  update = true;
728  }
729 
730  return update;
731 }
732 
738 {
739  NS_LOG_FUNCTION (this);
740 
741  // Compute the channel key. The key is reciprocal, i.e., key (a, b) = key (b, a)
742  uint32_t x1 = std::min (aMob->GetObject<Node> ()->GetId (), bMob->GetObject<Node> ()->GetId ());
743  uint32_t x2 = std::max (aMob->GetObject<Node> ()->GetId (), bMob->GetObject<Node> ()->GetId ());
744  uint32_t channelId = GetKey (x1, x2);
745 
746  // retrieve the channel condition
747  Ptr<const ChannelCondition> condition = m_channelConditionModel->GetChannelCondition (aMob, bMob);
748  bool los = (condition->GetLosCondition () == ChannelCondition::LosConditionValue::LOS);
749  bool o2i = false; // TODO include the o2i condition in the channel condition model
750 
751  // Check if the channel is present in the map and return it, otherwise
752  // generate a new channel
753  bool update = false;
754  bool notFound = false;
755  Ptr<ThreeGppChannelMatrix> channelMatrix;
756  if (m_channelMap.find (channelId) != m_channelMap.end ())
757  {
758  // channel matrix present in the map
759  NS_LOG_DEBUG ("channel matrix present in the map");
760  channelMatrix = m_channelMap[channelId];
761 
762  // check if it has to be updated
763  update = ChannelMatrixNeedsUpdate (channelMatrix, los);
764  }
765  else
766  {
767  NS_LOG_DEBUG ("channel matrix not found");
768  notFound = true;
769  }
770 
771  // If the channel is not present in the map or if it has to be updated
772  // generate a new realization
773  if (notFound || update)
774  {
775  // channel matrix not found or has to be updated, generate a new one
776  Angles txAngle (bMob->GetPosition (), aMob->GetPosition ());
777  Angles rxAngle (aMob->GetPosition (), bMob->GetPosition ());
778 
779  double x = aMob->GetPosition ().x - bMob->GetPosition ().x;
780  double y = aMob->GetPosition ().y - bMob->GetPosition ().y;
781  double distance2D = sqrt (x * x + y * y);
782 
783  // NOTE we assume hUT = min (height(a), height(b)) and
784  // hBS = max (height (a), height (b))
785  double hUt = std::min (aMob->GetPosition ().z, bMob->GetPosition ().z);
786  double hBs = std::max (aMob->GetPosition ().z, bMob->GetPosition ().z);
787 
788  // TODO this is not currently used, it is needed for the computation of the
789  // additional blockage in case of spatial consistent update
790  // I do not know who is the UT, I can use the relative distance between
791  // tx and rx instead
792  Vector locUt = Vector (0.0, 0.0, 0.0);
793 
794  channelMatrix = GetNewChannel (locUt, los, o2i, aAntenna, bAntenna, rxAngle, txAngle, distance2D, hBs, hUt);
795  channelMatrix->m_nodeIds = std::make_pair (aMob->GetObject<Node> ()->GetId (), bMob->GetObject<Node> ()->GetId ());
796 
797  // store or replace the channel matrix in the channel map
798  m_channelMap[channelId] = channelMatrix;
799  }
800 
801  return channelMatrix;
802 }
803 
805 ThreeGppChannelModel::GetNewChannel (Vector locUT, bool los, bool o2i,
808  Angles &uAngle, Angles &sAngle,
809  double dis2D, double hBS, double hUT) const
810 {
811  NS_LOG_FUNCTION (this);
812 
813  NS_ASSERT_MSG (m_frequency > 0.0, "Set the operating frequency first!");
814 
815  // get the 3GPP parameters
816  Ptr<const ParamsTable> table3gpp = GetThreeGppTable (los, o2i, hBS, hUT, dis2D);
817 
818  // get the number of clusters and the number of rays per cluster
819  uint8_t numOfCluster = table3gpp->m_numOfCluster;
820  uint8_t raysPerCluster = table3gpp->m_raysPerCluster;
821 
822  // create a channel matrix instance
823  Ptr<ThreeGppChannelMatrix> channelParams = Create<ThreeGppChannelMatrix> ();
824  channelParams->m_los = los; // set the LOS condition
825  channelParams->m_o2i = o2i; // set the O2I condition
826  channelParams->m_generatedTime = Simulator::Now ();
827 
828  // compute the 3D distance using eq. 7.4-1
829  double dis3D = std::sqrt (dis2D * dis2D + (hBS - hUT) * (hBS - hUT));
830 
831  //Step 4: Generate large scale parameters. All LSPS are uncorrelated.
832  DoubleVector LSPsIndep, LSPs;
833  uint8_t paramNum;
834  if (los)
835  {
836  paramNum = 7;
837  }
838  else
839  {
840  paramNum = 6;
841  }
842  //Generate paramNum independent LSPs.
843  for (uint8_t iter = 0; iter < paramNum; iter++)
844  {
845  LSPsIndep.push_back (m_normalRv->GetValue ());
846  }
847  for (uint8_t row = 0; row < paramNum; row++)
848  {
849  double temp = 0;
850  for (uint8_t column = 0; column < paramNum; column++)
851  {
852  temp += table3gpp->m_sqrtC[row][column] * LSPsIndep[column];
853  }
854  LSPs.push_back (temp);
855  }
856 
857  // NOTE the shadowing is generated in the propagation loss model
858 
859  double DS,ASD,ASA,ZSA,ZSD,K_factor = 0;
860  if (los)
861  {
862  K_factor = LSPs[1] * table3gpp->m_sigK + table3gpp->m_uK;
863  DS = pow (10, LSPs[2] * table3gpp->m_sigLgDS + table3gpp->m_uLgDS);
864  ASD = pow (10, LSPs[3] * table3gpp->m_sigLgASD + table3gpp->m_uLgASD);
865  ASA = pow (10, LSPs[4] * table3gpp->m_sigLgASA + table3gpp->m_uLgASA);
866  ZSD = pow (10, LSPs[5] * table3gpp->m_sigLgZSD + table3gpp->m_uLgZSD);
867  ZSA = pow (10, LSPs[6] * table3gpp->m_sigLgZSA + table3gpp->m_uLgZSA);
868  }
869  else
870  {
871  DS = pow (10, LSPs[1] * table3gpp->m_sigLgDS + table3gpp->m_uLgDS);
872  ASD = pow (10, LSPs[2] * table3gpp->m_sigLgASD + table3gpp->m_uLgASD);
873  ASA = pow (10, LSPs[3] * table3gpp->m_sigLgASA + table3gpp->m_uLgASA);
874  ZSD = pow (10, LSPs[4] * table3gpp->m_sigLgZSD + table3gpp->m_uLgZSD);
875  ZSA = pow (10, LSPs[5] * table3gpp->m_sigLgZSA + table3gpp->m_uLgZSA);
876 
877  }
878  ASD = std::min (ASD, 104.0);
879  ASA = std::min (ASA, 104.0);
880  ZSD = std::min (ZSD, 52.0);
881  ZSA = std::min (ZSA, 52.0);
882 
883  channelParams->m_DS = DS;
884  channelParams->m_K = K_factor;
885 
886  NS_LOG_INFO ("K-factor=" << K_factor << ",DS=" << DS << ", ASD=" << ASD << ", ASA=" << ASA << ", ZSD=" << ZSD << ", ZSA=" << ZSA);
887 
888  //Step 5: Generate Delays.
889  DoubleVector clusterDelay;
890  double minTau = 100.0;
891  for (uint8_t cIndex = 0; cIndex < numOfCluster; cIndex++)
892  {
893  double tau = -1*table3gpp->m_rTau*DS*log (m_uniformRv->GetValue (0,1)); //(7.5-1)
894  if (minTau > tau)
895  {
896  minTau = tau;
897  }
898  clusterDelay.push_back (tau);
899  }
900 
901  for (uint8_t cIndex = 0; cIndex < numOfCluster; cIndex++)
902  {
903  clusterDelay[cIndex] -= minTau;
904  }
905  std::sort (clusterDelay.begin (), clusterDelay.end ()); //(7.5-2)
906 
907  /* since the scaled Los delays are not to be used in cluster power generation,
908  * we will generate cluster power first and resume to compute Los cluster delay later.*/
909 
910  //Step 6: Generate cluster powers.
911  DoubleVector clusterPower;
912  double powerSum = 0;
913  for (uint8_t cIndex = 0; cIndex < numOfCluster; cIndex++)
914  {
915  double power = exp (-1 * clusterDelay[cIndex] * (table3gpp->m_rTau - 1) / table3gpp->m_rTau / DS) *
916  pow (10,-1 * m_normalRv->GetValue () * table3gpp->m_perClusterShadowingStd / 10); //(7.5-5)
917  powerSum += power;
918  clusterPower.push_back (power);
919  }
920  double powerMax = 0;
921 
922  for (uint8_t cIndex = 0; cIndex < numOfCluster; cIndex++)
923  {
924  clusterPower[cIndex] = clusterPower[cIndex] / powerSum; //(7.5-6)
925  }
926 
927  DoubleVector clusterPowerForAngles; // this power is only for equation (7.5-9) and (7.5-14), not for (7.5-22)
928  if (los)
929  {
930  double K_linear = pow (10,K_factor / 10);
931 
932  for (uint8_t cIndex = 0; cIndex < numOfCluster; cIndex++)
933  {
934  if (cIndex == 0)
935  {
936  clusterPowerForAngles.push_back (clusterPower[cIndex] / (1 + K_linear) + K_linear / (1 + K_linear)); //(7.5-8)
937  }
938  else
939  {
940  clusterPowerForAngles.push_back (clusterPower[cIndex] / (1 + K_linear)); //(7.5-8)
941  }
942  if (powerMax < clusterPowerForAngles[cIndex])
943  {
944  powerMax = clusterPowerForAngles[cIndex];
945  }
946  }
947  }
948  else
949  {
950  for (uint8_t cIndex = 0; cIndex < numOfCluster; cIndex++)
951  {
952  clusterPowerForAngles.push_back (clusterPower[cIndex]); //(7.5-6)
953  if (powerMax < clusterPowerForAngles[cIndex])
954  {
955  powerMax = clusterPowerForAngles[cIndex];
956  }
957  }
958  }
959 
960  //remove clusters with less than -25 dB power compared to the maxim cluster power;
961  //double thresh = pow(10,-2.5);
962  double thresh = 0.0032;
963  for (uint8_t cIndex = numOfCluster; cIndex > 0; cIndex--)
964  {
965  if (clusterPowerForAngles[cIndex - 1] < thresh * powerMax )
966  {
967  clusterPowerForAngles.erase (clusterPowerForAngles.begin () + cIndex - 1);
968  clusterPower.erase (clusterPower.begin () + cIndex - 1);
969  clusterDelay.erase (clusterDelay.begin () + cIndex - 1);
970  }
971  }
972  uint8_t numReducedCluster = clusterPower.size ();
973 
974  channelParams->m_numCluster = numReducedCluster;
975  // Resume step 5 to compute the delay for LoS condition.
976  if (los)
977  {
978  double C_tau = 0.7705 - 0.0433 * K_factor + 2e-4 * pow (K_factor,2) + 17e-6 * pow (K_factor,3); //(7.5-3)
979  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
980  {
981  clusterDelay[cIndex] = clusterDelay[cIndex] / C_tau; //(7.5-4)
982  }
983  }
984 
985  //Step 7: Generate arrival and departure angles for both azimuth and elevation.
986 
987  double C_NLOS, C_phi;
988  //According to table 7.5-6, only cluster number equals to 8, 10, 11, 12, 19 and 20 is valid.
989  //Not sure why the other cases are in Table 7.5-2.
990  switch (numOfCluster) // Table 7.5-2
991  {
992  case 4:
993  C_NLOS = 0.779;
994  break;
995  case 5:
996  C_NLOS = 0.860;
997  break;
998  case 8:
999  C_NLOS = 1.018;
1000  break;
1001  case 10:
1002  C_NLOS = 1.090;
1003  break;
1004  case 11:
1005  C_NLOS = 1.123;
1006  break;
1007  case 12:
1008  C_NLOS = 1.146;
1009  break;
1010  case 14:
1011  C_NLOS = 1.190;
1012  break;
1013  case 15:
1014  C_NLOS = 1.221;
1015  break;
1016  case 16:
1017  C_NLOS = 1.226;
1018  break;
1019  case 19:
1020  C_NLOS = 1.273;
1021  break;
1022  case 20:
1023  C_NLOS = 1.289;
1024  break;
1025  default:
1026  NS_FATAL_ERROR ("Invalide cluster number");
1027  }
1028 
1029  if (los)
1030  {
1031  C_phi = C_NLOS * (1.1035 - 0.028 * K_factor - 2e-3 * pow (K_factor,2) + 1e-4 * pow (K_factor,3)); //(7.5-10))
1032  }
1033  else
1034  {
1035  C_phi = C_NLOS; //(7.5-10)
1036  }
1037 
1038  double C_theta;
1039  switch (numOfCluster) //Table 7.5-4
1040  {
1041  case 8:
1042  C_NLOS = 0.889;
1043  break;
1044  case 10:
1045  C_NLOS = 0.957;
1046  break;
1047  case 11:
1048  C_NLOS = 1.031;
1049  break;
1050  case 12:
1051  C_NLOS = 1.104;
1052  break;
1053  case 15:
1054  C_NLOS = 1.1088;
1055  break;
1056  case 19:
1057  C_NLOS = 1.184;
1058  break;
1059  case 20:
1060  C_NLOS = 1.178;
1061  break;
1062  default:
1063  NS_FATAL_ERROR ("Invalide cluster number");
1064  }
1065 
1066  if (los)
1067  {
1068  C_theta = C_NLOS * (1.3086 + 0.0339 * K_factor - 0.0077 * pow (K_factor,2) + 2e-4 * pow (K_factor,3)); //(7.5-15)
1069  }
1070  else
1071  {
1072  C_theta = C_NLOS;
1073  }
1074 
1075 
1076  DoubleVector clusterAoa, clusterAod, clusterZoa, clusterZod;
1077  double angle;
1078  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
1079  {
1080  angle = 2*ASA*sqrt (-1 * log (clusterPowerForAngles[cIndex] / powerMax)) / 1.4 / C_phi; //(7.5-9)
1081  clusterAoa.push_back (angle);
1082  angle = 2*ASD*sqrt (-1 * log (clusterPowerForAngles[cIndex] / powerMax)) / 1.4 / C_phi; //(7.5-9)
1083  clusterAod.push_back (angle);
1084  angle = -1*ZSA*log (clusterPowerForAngles[cIndex] / powerMax) / C_theta; //(7.5-14)
1085  clusterZoa.push_back (angle);
1086  angle = -1*ZSD*log (clusterPowerForAngles[cIndex] / powerMax) / C_theta;
1087  clusterZod.push_back (angle);
1088  }
1089 
1090  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
1091  {
1092  int Xn = 1;
1093  if (m_uniformRv->GetValue (0,1) < 0.5)
1094  {
1095  Xn = -1;
1096  }
1097  clusterAoa[cIndex] = clusterAoa[cIndex] * Xn + (m_normalRv->GetValue () * ASA / 7) + uAngle.phi * 180 / M_PI; //(7.5-11)
1098  clusterAod[cIndex] = clusterAod[cIndex] * Xn + (m_normalRv->GetValue () * ASD / 7) + sAngle.phi * 180 / M_PI;
1099  if (o2i)
1100  {
1101  clusterZoa[cIndex] = clusterZoa[cIndex] * Xn + (m_normalRv->GetValue () * ZSA / 7) + 90; //(7.5-16)
1102  }
1103  else
1104  {
1105  clusterZoa[cIndex] = clusterZoa[cIndex] * Xn + (m_normalRv->GetValue () * ZSA / 7) + uAngle.theta * 180 / M_PI; //(7.5-16)
1106  }
1107  clusterZod[cIndex] = clusterZod[cIndex] * Xn + (m_normalRv->GetValue () * ZSD / 7) + sAngle.theta * 180 / M_PI + table3gpp->m_offsetZOD; //(7.5-19)
1108 
1109  }
1110 
1111  if (los)
1112  {
1113  //The 7.5-12 can be rewrite as Theta_n,ZOA = Theta_n,ZOA - (Theta_1,ZOA - Theta_LOS,ZOA) = Theta_n,ZOA - diffZOA,
1114  //Similar as AOD, ZSA and ZSD.
1115  double diffAoa = clusterAoa[0] - uAngle.phi * 180 / M_PI;
1116  double diffAod = clusterAod[0] - sAngle.phi * 180 / M_PI;
1117  double diffZsa = clusterZoa[0] - uAngle.theta * 180 / M_PI;
1118  double diffZsd = clusterZod[0] - sAngle.theta * 180 / M_PI;
1119 
1120  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
1121  {
1122  clusterAoa[cIndex] -= diffAoa; //(7.5-12)
1123  clusterAod[cIndex] -= diffAod;
1124  clusterZoa[cIndex] -= diffZsa; //(7.5-17)
1125  clusterZod[cIndex] -= diffZsd;
1126 
1127  }
1128  }
1129 
1130  double rayAoa_radian[numReducedCluster][raysPerCluster]; //rayAoa_radian[n][m], where n is cluster index, m is ray index
1131  double rayAod_radian[numReducedCluster][raysPerCluster]; //rayAod_radian[n][m], where n is cluster index, m is ray index
1132  double rayZoa_radian[numReducedCluster][raysPerCluster]; //rayZoa_radian[n][m], where n is cluster index, m is ray index
1133  double rayZod_radian[numReducedCluster][raysPerCluster]; //rayZod_radian[n][m], where n is cluster index, m is ray index
1134 
1135  for (uint8_t nInd = 0; nInd < numReducedCluster; nInd++)
1136  {
1137  for (uint8_t mInd = 0; mInd < raysPerCluster; mInd++)
1138  {
1139  double tempAoa = clusterAoa[nInd] + table3gpp->m_cASA * offSetAlpha[mInd]; //(7.5-13)
1140  while (tempAoa > 360)
1141  {
1142  tempAoa -= 360;
1143  }
1144 
1145  while (tempAoa < 0)
1146  {
1147  tempAoa += 360;
1148 
1149  }
1150  NS_ASSERT_MSG (tempAoa >= 0 && tempAoa <= 360, "the AOA should be the range of [0,360]");
1151  rayAoa_radian[nInd][mInd] = tempAoa * M_PI / 180;
1152 
1153  double tempAod = clusterAod[nInd] + table3gpp->m_cASD * offSetAlpha[mInd];
1154  while (tempAod > 360)
1155  {
1156  tempAod -= 360;
1157  }
1158 
1159  while (tempAod < 0)
1160  {
1161  tempAod += 360;
1162  }
1163  NS_ASSERT_MSG (tempAod >= 0 && tempAod <= 360, "the AOD should be the range of [0,360]");
1164  rayAod_radian[nInd][mInd] = tempAod * M_PI / 180;
1165 
1166  double tempZoa = clusterZoa[nInd] + table3gpp->m_cZSA * offSetAlpha[mInd]; //(7.5-18)
1167 
1168  while (tempZoa > 360)
1169  {
1170  tempZoa -= 360;
1171  }
1172 
1173  while (tempZoa < 0)
1174  {
1175  tempZoa += 360;
1176  }
1177 
1178  if (tempZoa > 180)
1179  {
1180  tempZoa = 360 - tempZoa;
1181  }
1182 
1183  NS_ASSERT_MSG (tempZoa >= 0&&tempZoa <= 180, "the ZOA should be the range of [0,180]");
1184  rayZoa_radian[nInd][mInd] = tempZoa * M_PI / 180;
1185 
1186  double tempZod = clusterZod[nInd] + 0.375 * pow (10,table3gpp->m_uLgZSD) * offSetAlpha[mInd]; //(7.5-20)
1187 
1188  while (tempZod > 360)
1189  {
1190  tempZod -= 360;
1191  }
1192 
1193  while (tempZod < 0)
1194  {
1195  tempZod += 360;
1196  }
1197  if (tempZod > 180)
1198  {
1199  tempZod = 360 - tempZod;
1200  }
1201  NS_ASSERT_MSG (tempZod >= 0&&tempZod <= 180, "the ZOD should be the range of [0,180]");
1202  rayZod_radian[nInd][mInd] = tempZod * M_PI / 180;
1203  }
1204  }
1205  DoubleVector angle_degree;
1206  double sizeTemp = clusterZoa.size ();
1207  for (uint8_t ind = 0; ind < 4; ind++)
1208  {
1209  switch (ind)
1210  {
1211  case 0:
1212  angle_degree = clusterAoa;
1213  break;
1214  case 1:
1215  angle_degree = clusterZoa;
1216  break;
1217  case 2:
1218  angle_degree = clusterAod;
1219  break;
1220  case 3:
1221  angle_degree = clusterZod;
1222  break;
1223  default:
1224  NS_FATAL_ERROR ("Programming Error");
1225  }
1226 
1227  for (uint8_t nIndex = 0; nIndex < sizeTemp; nIndex++)
1228  {
1229  while (angle_degree[nIndex] > 360)
1230  {
1231  angle_degree[nIndex] -= 360;
1232  }
1233 
1234  while (angle_degree[nIndex] < 0)
1235  {
1236  angle_degree[nIndex] += 360;
1237  }
1238 
1239  if (ind == 1 || ind == 3)
1240  {
1241  if (angle_degree[nIndex] > 180)
1242  {
1243  angle_degree[nIndex] = 360 - angle_degree[nIndex];
1244  }
1245  }
1246  }
1247  switch (ind)
1248  {
1249  case 0:
1250  clusterAoa = angle_degree;
1251  break;
1252  case 1:
1253  clusterZoa = angle_degree;
1254  break;
1255  case 2:
1256  clusterAod = angle_degree;
1257  break;
1258  case 3:
1259  clusterZod = angle_degree;
1260  break;
1261  default:
1262  NS_FATAL_ERROR ("Programming Error");
1263  }
1264  }
1265 
1266  DoubleVector attenuation_dB;
1267  if (m_blockage)
1268  {
1269  attenuation_dB = CalcAttenuationOfBlockage (channelParams, clusterAoa, clusterZoa);
1270  for (uint8_t cInd = 0; cInd < numReducedCluster; cInd++)
1271  {
1272  clusterPower[cInd] = clusterPower[cInd] / pow (10,attenuation_dB[cInd] / 10);
1273  }
1274  }
1275  else
1276  {
1277  attenuation_dB.push_back (0);
1278  }
1279 
1280  //Step 8: Coupling of rays within a cluster for both azimuth and elevation
1281  //shuffle all the arrays to perform random coupling
1282  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
1283  {
1284  std::shuffle (&rayAod_radian[cIndex][0],&rayAod_radian[cIndex][raysPerCluster],std::default_random_engine (cIndex * 1000 + 100));
1285  std::shuffle (&rayAoa_radian[cIndex][0],&rayAoa_radian[cIndex][raysPerCluster],std::default_random_engine (cIndex * 1000 + 200));
1286  std::shuffle (&rayZod_radian[cIndex][0],&rayZod_radian[cIndex][raysPerCluster],std::default_random_engine (cIndex * 1000 + 300));
1287  std::shuffle (&rayZoa_radian[cIndex][0],&rayZoa_radian[cIndex][raysPerCluster],std::default_random_engine (cIndex * 1000 + 400));
1288  }
1289 
1290  //Step 9: Generate the cross polarization power ratios
1291  //Step 10: Draw initial phases
1292  Double2DVector crossPolarizationPowerRatios; // vector containing the cross polarization power ratios, as defined by 7.5-21
1293  Double3DVector clusterPhase; //rayAoa_radian[n][m], where n is cluster index, m is ray index
1294  for (uint8_t nInd = 0; nInd < numReducedCluster; nInd++)
1295  {
1296  DoubleVector temp; // used to store the XPR values
1297  Double2DVector temp2; // used to store the PHI values for all the possible combination of polarization
1298  for (uint8_t mInd = 0; mInd < raysPerCluster; mInd++)
1299  {
1300  double uXprLinear = pow (10, table3gpp->m_uXpr / 10); // convert to linear
1301  double sigXprLinear = pow (10, table3gpp->m_sigXpr / 10); // convert to linear
1302 
1303  temp.push_back (std::pow (10, (m_normalRv->GetValue () * sigXprLinear + uXprLinear) / 10));
1304  DoubleVector temp3; // used to store the PHI valuse
1305  for (uint8_t pInd = 0; pInd < 4; pInd++)
1306  {
1307  temp3.push_back (m_uniformRv->GetValue (-1 * M_PI, M_PI));
1308  }
1309  temp2.push_back (temp3);
1310  }
1311  crossPolarizationPowerRatios.push_back (temp);
1312  clusterPhase.push_back (temp2);
1313  }
1314  channelParams->m_clusterPhase = clusterPhase;
1315 
1316  //Step 11: Generate channel coefficients for each cluster n and each receiver
1317  // and transmitter element pair u,s.
1318 
1319  Complex3DVector H_NLOS; // channel coefficients H_NLOS [u][s][n],
1320  // where u and s are receive and transmit antenna element, n is cluster index.
1321  uint64_t uSize = uAntenna->GetNumberOfElements ();
1322  uint64_t sSize = sAntenna->GetNumberOfElements ();
1323 
1324  uint8_t cluster1st = 0, cluster2nd = 0; // first and second strongest cluster;
1325  double maxPower = 0;
1326  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
1327  {
1328  if (maxPower < clusterPower[cIndex])
1329  {
1330  maxPower = clusterPower[cIndex];
1331  cluster1st = cIndex;
1332  }
1333  }
1334  maxPower = 0;
1335  for (uint8_t cIndex = 0; cIndex < numReducedCluster; cIndex++)
1336  {
1337  if (maxPower < clusterPower[cIndex] && cluster1st != cIndex)
1338  {
1339  maxPower = clusterPower[cIndex];
1340  cluster2nd = cIndex;
1341  }
1342  }
1343 
1344  NS_LOG_INFO ("1st strongest cluster:" << (int)cluster1st << ", 2nd strongest cluster:" << (int)cluster2nd);
1345 
1346  Complex3DVector H_usn; //channel coffecient H_usn[u][s][n];
1347  // NOTE Since each of the strongest 2 clusters are divided into 3 sub-clusters,
1348  // the total cluster will be numReducedCLuster + 4.
1349 
1350  H_usn.resize (uSize);
1351  for (uint64_t uIndex = 0; uIndex < uSize; uIndex++)
1352  {
1353  H_usn[uIndex].resize (sSize);
1354  for (uint64_t sIndex = 0; sIndex < sSize; sIndex++)
1355  {
1356  H_usn[uIndex][sIndex].resize (numReducedCluster);
1357  }
1358  }
1359 
1360  // The following for loops computes the channel coefficients
1361  for (uint64_t uIndex = 0; uIndex < uSize; uIndex++)
1362  {
1363  Vector uLoc = uAntenna->GetElementLocation (uIndex);
1364 
1365  for (uint64_t sIndex = 0; sIndex < sSize; sIndex++)
1366  {
1367 
1368  Vector sLoc = sAntenna->GetElementLocation (sIndex);
1369 
1370  for (uint8_t nIndex = 0; nIndex < numReducedCluster; nIndex++)
1371  {
1372  //Compute the N-2 weakest cluster, only vertical polarization. (7.5-22)
1373  if (nIndex != cluster1st && nIndex != cluster2nd)
1374  {
1375  std::complex<double> rays (0,0);
1376  for (uint8_t mIndex = 0; mIndex < raysPerCluster; mIndex++)
1377  {
1378  DoubleVector initialPhase = clusterPhase[nIndex][mIndex];
1379  double k = crossPolarizationPowerRatios[nIndex][mIndex];
1380  //lambda_0 is accounted in the antenna spacing uLoc and sLoc.
1381  double rxPhaseDiff = 2 * M_PI * (sin (rayZoa_radian[nIndex][mIndex]) * cos (rayAoa_radian[nIndex][mIndex]) * uLoc.x
1382  + sin (rayZoa_radian[nIndex][mIndex]) * sin (rayAoa_radian[nIndex][mIndex]) * uLoc.y
1383  + cos (rayZoa_radian[nIndex][mIndex]) * uLoc.z);
1384 
1385  double txPhaseDiff = 2 * M_PI * (sin (rayZod_radian[nIndex][mIndex]) * cos (rayAod_radian[nIndex][mIndex]) * sLoc.x
1386  + sin (rayZod_radian[nIndex][mIndex]) * sin (rayAod_radian[nIndex][mIndex]) * sLoc.y
1387  + cos (rayZod_radian[nIndex][mIndex]) * sLoc.z);
1388  // NOTE Doppler is computed in the CalcBeamformingGain function and is simplified to only account for the center anngle of each cluster.
1389 
1390  double rxFieldPatternPhi, rxFieldPatternTheta, txFieldPatternPhi, txFieldPatternTheta;
1391  std::tie (rxFieldPatternPhi, rxFieldPatternTheta) = uAntenna->GetElementFieldPattern (Angles (rayAoa_radian[nIndex][mIndex], rayZoa_radian[nIndex][mIndex]));
1392  std::tie (txFieldPatternPhi, txFieldPatternTheta) = sAntenna->GetElementFieldPattern (Angles (rayAod_radian[nIndex][mIndex], rayZod_radian[nIndex][mIndex]));
1393 
1394  rays += (exp (std::complex<double> (0, initialPhase[0])) * rxFieldPatternTheta * txFieldPatternTheta +
1395  +exp (std::complex<double> (0, initialPhase[1])) * std::sqrt (1 / k) * rxFieldPatternTheta * txFieldPatternPhi +
1396  +exp (std::complex<double> (0, initialPhase[2])) * std::sqrt (1 / k) * rxFieldPatternPhi * txFieldPatternTheta +
1397  +exp (std::complex<double> (0, initialPhase[3])) * rxFieldPatternPhi * txFieldPatternPhi)
1398  * exp (std::complex<double> (0, rxPhaseDiff))
1399  * exp (std::complex<double> (0, txPhaseDiff));
1400  }
1401  rays *= sqrt (clusterPower[nIndex] / raysPerCluster);
1402  H_usn[uIndex][sIndex][nIndex] = rays;
1403  }
1404  else //(7.5-28)
1405  {
1406  std::complex<double> raysSub1 (0,0);
1407  std::complex<double> raysSub2 (0,0);
1408  std::complex<double> raysSub3 (0,0);
1409 
1410  for (uint8_t mIndex = 0; mIndex < raysPerCluster; mIndex++)
1411  {
1412  double k = crossPolarizationPowerRatios[nIndex][mIndex];
1413 
1414  //ZML:Just remind me that the angle offsets for the 3 subclusters were not generated correctly.
1415 
1416  DoubleVector initialPhase = clusterPhase[nIndex][mIndex];
1417  double rxPhaseDiff = 2 * M_PI * (sin (rayZoa_radian[nIndex][mIndex]) * cos (rayAoa_radian[nIndex][mIndex]) * uLoc.x
1418  + sin (rayZoa_radian[nIndex][mIndex]) * sin (rayAoa_radian[nIndex][mIndex]) * uLoc.y
1419  + cos (rayZoa_radian[nIndex][mIndex]) * uLoc.z);
1420  double txPhaseDiff = 2 * M_PI * (sin (rayZod_radian[nIndex][mIndex]) * cos (rayAod_radian[nIndex][mIndex]) * sLoc.x
1421  + sin (rayZod_radian[nIndex][mIndex]) * sin (rayAod_radian[nIndex][mIndex]) * sLoc.y
1422  + cos (rayZod_radian[nIndex][mIndex]) * sLoc.z);
1423 
1424  double rxFieldPatternPhi, rxFieldPatternTheta, txFieldPatternPhi, txFieldPatternTheta;
1425  std::tie (rxFieldPatternPhi, rxFieldPatternTheta) = uAntenna->GetElementFieldPattern (Angles (rayAoa_radian[nIndex][mIndex], rayZoa_radian[nIndex][mIndex]));
1426  std::tie (txFieldPatternPhi, txFieldPatternTheta) = sAntenna->GetElementFieldPattern (Angles (rayAod_radian[nIndex][mIndex], rayZod_radian[nIndex][mIndex]));
1427 
1428  switch (mIndex)
1429  {
1430  case 9:
1431  case 10:
1432  case 11:
1433  case 12:
1434  case 17:
1435  case 18:
1436  raysSub2 += (exp (std::complex<double> (0, initialPhase[0])) * rxFieldPatternTheta * txFieldPatternTheta +
1437  +exp (std::complex<double> (0, initialPhase[1])) * std::sqrt (1 / k) * rxFieldPatternTheta * txFieldPatternPhi +
1438  +exp (std::complex<double> (0, initialPhase[2])) * std::sqrt (1 / k) * rxFieldPatternPhi * txFieldPatternTheta +
1439  +exp (std::complex<double> (0, initialPhase[3])) * rxFieldPatternPhi * txFieldPatternPhi)
1440  * exp (std::complex<double> (0, rxPhaseDiff))
1441  * exp (std::complex<double> (0, txPhaseDiff));
1442  break;
1443  case 13:
1444  case 14:
1445  case 15:
1446  case 16:
1447  raysSub3 += (exp (std::complex<double> (0, initialPhase[0])) * rxFieldPatternTheta * txFieldPatternTheta +
1448  +exp (std::complex<double> (0, initialPhase[1])) * std::sqrt (1 / k) * rxFieldPatternTheta * txFieldPatternPhi +
1449  +exp (std::complex<double> (0, initialPhase[2])) * std::sqrt (1 / k) * rxFieldPatternPhi * txFieldPatternTheta +
1450  +exp (std::complex<double> (0, initialPhase[3])) * rxFieldPatternPhi * txFieldPatternPhi)
1451  * exp (std::complex<double> (0, rxPhaseDiff))
1452  * exp (std::complex<double> (0, txPhaseDiff));
1453  break;
1454  default: //case 1,2,3,4,5,6,7,8,19,20
1455  raysSub1 += (exp (std::complex<double> (0, initialPhase[0])) * rxFieldPatternTheta * txFieldPatternTheta +
1456  +exp (std::complex<double> (0, initialPhase[1])) * std::sqrt (1 / k) * rxFieldPatternTheta * txFieldPatternPhi +
1457  +exp (std::complex<double> (0, initialPhase[2])) * std::sqrt (1 / k) * rxFieldPatternPhi * txFieldPatternTheta +
1458  +exp (std::complex<double> (0, initialPhase[3])) * rxFieldPatternPhi * txFieldPatternPhi)
1459  * exp (std::complex<double> (0, rxPhaseDiff))
1460  * exp (std::complex<double> (0, txPhaseDiff));
1461  break;
1462  }
1463  }
1464  raysSub1 *= sqrt (clusterPower[nIndex] / raysPerCluster);
1465  raysSub2 *= sqrt (clusterPower[nIndex] / raysPerCluster);
1466  raysSub3 *= sqrt (clusterPower[nIndex] / raysPerCluster);
1467  H_usn[uIndex][sIndex][nIndex] = raysSub1;
1468  H_usn[uIndex][sIndex].push_back (raysSub2);
1469  H_usn[uIndex][sIndex].push_back (raysSub3);
1470 
1471  }
1472  }
1473  if (los) //(7.5-29) && (7.5-30)
1474  {
1475  std::complex<double> ray (0,0);
1476  double rxPhaseDiff = 2 * M_PI * (sin (uAngle.theta) * cos (uAngle.phi) * uLoc.x
1477  + sin (uAngle.theta) * sin (uAngle.phi) * uLoc.y
1478  + cos (uAngle.theta) * uLoc.z);
1479  double txPhaseDiff = 2 * M_PI * (sin (sAngle.theta) * cos (sAngle.phi) * sLoc.x
1480  + sin (sAngle.theta) * sin (sAngle.phi) * sLoc.y
1481  + cos (sAngle.theta) * sLoc.z);
1482 
1483  double rxFieldPatternPhi, rxFieldPatternTheta, txFieldPatternPhi, txFieldPatternTheta;
1484  std::tie (rxFieldPatternPhi, rxFieldPatternTheta) = uAntenna->GetElementFieldPattern (Angles (uAngle.phi, uAngle.theta));
1485  std::tie (txFieldPatternPhi, txFieldPatternTheta) = sAntenna->GetElementFieldPattern (Angles (sAngle.phi, sAngle.theta));
1486 
1487  double lambda = 3e8 / m_frequency; // the wavelength of the carrier frequency
1488 
1489  ray = (rxFieldPatternTheta * txFieldPatternTheta - rxFieldPatternPhi * txFieldPatternPhi)
1490  * exp (std::complex<double> (0, - 2 * M_PI * dis3D / lambda))
1491  * exp (std::complex<double> (0, rxPhaseDiff))
1492  * exp (std::complex<double> (0, txPhaseDiff));
1493 
1494  double K_linear = pow (10,K_factor / 10);
1495  // the LOS path should be attenuated if blockage is enabled.
1496  H_usn[uIndex][sIndex][0] = sqrt (1 / (K_linear + 1)) * H_usn[uIndex][sIndex][0] + sqrt (K_linear / (1 + K_linear)) * ray / pow (10,attenuation_dB[0] / 10); //(7.5-30) for tau = tau1
1497  double tempSize = H_usn[uIndex][sIndex].size ();
1498  for (uint8_t nIndex = 1; nIndex < tempSize; nIndex++)
1499  {
1500  H_usn[uIndex][sIndex][nIndex] *= sqrt (1 / (K_linear + 1)); //(7.5-30) for tau = tau2...taunN
1501  }
1502 
1503  }
1504  }
1505  }
1506 
1507  // store the delays and the angles for the subclusters
1508  if (cluster1st == cluster2nd)
1509  {
1510  clusterDelay.push_back (clusterDelay[cluster1st] + 1.28 * table3gpp->m_cDS);
1511  clusterDelay.push_back (clusterDelay[cluster1st] + 2.56 * table3gpp->m_cDS);
1512 
1513  clusterAoa.push_back (clusterAoa[cluster1st]);
1514  clusterAoa.push_back (clusterAoa[cluster1st]);
1515 
1516  clusterZoa.push_back (clusterZoa[cluster1st]);
1517  clusterZoa.push_back (clusterZoa[cluster1st]);
1518 
1519  clusterAod.push_back (clusterAod[cluster1st]);
1520  clusterAod.push_back (clusterAod[cluster1st]);
1521 
1522  clusterZod.push_back (clusterZod[cluster1st]);
1523  clusterZod.push_back (clusterZod[cluster1st]);
1524  }
1525  else
1526  {
1527  double min, max;
1528  if (cluster1st < cluster2nd)
1529  {
1530  min = cluster1st;
1531  max = cluster2nd;
1532  }
1533  else
1534  {
1535  min = cluster2nd;
1536  max = cluster1st;
1537  }
1538  clusterDelay.push_back (clusterDelay[min] + 1.28 * table3gpp->m_cDS);
1539  clusterDelay.push_back (clusterDelay[min] + 2.56 * table3gpp->m_cDS);
1540  clusterDelay.push_back (clusterDelay[max] + 1.28 * table3gpp->m_cDS);
1541  clusterDelay.push_back (clusterDelay[max] + 2.56 * table3gpp->m_cDS);
1542 
1543  clusterAoa.push_back (clusterAoa[min]);
1544  clusterAoa.push_back (clusterAoa[min]);
1545  clusterAoa.push_back (clusterAoa[max]);
1546  clusterAoa.push_back (clusterAoa[max]);
1547 
1548  clusterZoa.push_back (clusterZoa[min]);
1549  clusterZoa.push_back (clusterZoa[min]);
1550  clusterZoa.push_back (clusterZoa[max]);
1551  clusterZoa.push_back (clusterZoa[max]);
1552 
1553  clusterAod.push_back (clusterAod[min]);
1554  clusterAod.push_back (clusterAod[min]);
1555  clusterAod.push_back (clusterAod[max]);
1556  clusterAod.push_back (clusterAod[max]);
1557 
1558  clusterZod.push_back (clusterZod[min]);
1559  clusterZod.push_back (clusterZod[min]);
1560  clusterZod.push_back (clusterZod[max]);
1561  clusterZod.push_back (clusterZod[max]);
1562 
1563 
1564  }
1565 
1566  NS_LOG_INFO ("size of coefficient matrix =[" << H_usn.size () << "][" << H_usn[0].size () << "][" << H_usn[0][0].size () << "]");
1567 
1568  channelParams->m_channel = H_usn;
1569  channelParams->m_delay = clusterDelay;
1570 
1571  channelParams->m_angle.clear ();
1572  channelParams->m_angle.push_back (clusterAoa);
1573  channelParams->m_angle.push_back (clusterZoa);
1574  channelParams->m_angle.push_back (clusterAod);
1575  channelParams->m_angle.push_back (clusterZod);
1576 
1577  return channelParams;
1578 }
1579 
1582  const DoubleVector &clusterAOA,
1583  const DoubleVector &clusterZOA) const
1584 {
1585  NS_LOG_FUNCTION (this);
1586 
1587  DoubleVector powerAttenuation;
1588  uint8_t clusterNum = clusterAOA.size ();
1589  for (uint8_t cInd = 0; cInd < clusterNum; cInd++)
1590  {
1591  powerAttenuation.push_back (0); //Initial power attenuation for all clusters to be 0 dB;
1592  }
1593  //step a: the number of non-self blocking blockers is stored in m_numNonSelfBlocking.
1594 
1595  //step b:Generate the size and location of each blocker
1596  //generate self blocking (i.e., for blockage from the human body)
1597  double phi_sb, x_sb, theta_sb, y_sb;
1598  //table 7.6.4.1-1 Self-blocking region parameters.
1599  if (m_portraitMode)
1600  {
1601  phi_sb = 260;
1602  x_sb = 120;
1603  theta_sb = 100;
1604  y_sb = 80;
1605  }
1606  else // landscape mode
1607  {
1608  phi_sb = 40;
1609  x_sb = 160;
1610  theta_sb = 110;
1611  y_sb = 75;
1612  }
1613 
1614  //generate or update non-self blocking
1615  if (params->m_nonSelfBlocking.size () == 0) //generate new blocking regions
1616  {
1617  for (uint16_t blockInd = 0; blockInd < m_numNonSelfBlocking; blockInd++)
1618  {
1619  //draw value from table 7.6.4.1-2 Blocking region parameters
1620  DoubleVector table;
1621  table.push_back (m_normalRv->GetValue ()); //phi_k: store the normal RV that will be mapped to uniform (0,360) later.
1622  if (m_scenario == "InH-OfficeMixed" || m_scenario == "InH-OfficeOpen")
1623  {
1624  table.push_back (m_uniformRv->GetValue (15, 45)); //x_k
1625  table.push_back (90); //Theta_k
1626  table.push_back (m_uniformRv->GetValue (5, 15)); //y_k
1627  table.push_back (2); //r
1628  }
1629  else
1630  {
1631  table.push_back (m_uniformRv->GetValue (5, 15)); //x_k
1632  table.push_back (90); //Theta_k
1633  table.push_back (5); //y_k
1634  table.push_back (10); //r
1635  }
1636  params->m_nonSelfBlocking.push_back (table);
1637  }
1638  }
1639  else
1640  {
1641  double deltaX = sqrt (pow (params->m_preLocUT.x - params->m_locUT.x, 2) + pow (params->m_preLocUT.y - params->m_locUT.y, 2));
1642  //if deltaX and speed are both 0, the autocorrelation is 1, skip updating
1643  if (deltaX > 1e-6 || m_blockerSpeed > 1e-6)
1644  {
1645  double corrDis;
1646  //draw value from table 7.6.4.1-4: Spatial correlation distance for different m_scenarios.
1647  if (m_scenario == "InH-OfficeMixed" || m_scenario == "InH-OfficeOpen")
1648  {
1649  //InH, correlation distance = 5;
1650  corrDis = 5;
1651  }
1652  else
1653  {
1654  if (params->m_o2i) // outdoor to indoor
1655  {
1656  corrDis = 5;
1657  }
1658  else //LOS or NLOS
1659  {
1660  corrDis = 10;
1661  }
1662  }
1663  double R;
1664  if (m_blockerSpeed > 1e-6) // speed not equal to 0
1665  {
1666  double corrT = corrDis / m_blockerSpeed;
1667  R = exp (-1 * (deltaX / corrDis + (Now ().GetSeconds () - params->m_generatedTime.GetSeconds ()) / corrT));
1668  }
1669  else
1670  {
1671  R = exp (-1 * (deltaX / corrDis));
1672  }
1673 
1674  NS_LOG_INFO ("Distance change:" << deltaX << " Speed:" << m_blockerSpeed
1675  << " Time difference:" << Now ().GetSeconds () - params->m_generatedTime.GetSeconds ()
1676  << " correlation:" << R);
1677 
1678  //In order to generate correlated uniform random variables, we first generate correlated normal random variables and map the normal RV to uniform RV.
1679  //Notice the correlation will change if the RV is transformed from normal to uniform.
1680  //To compensate the distortion, the correlation of the normal RV is computed
1681  //such that the uniform RV would have the desired correlation when transformed from normal RV.
1682 
1683  //The following formula was obtained from MATLAB numerical simulation.
1684 
1685  if (R * R * (-0.069) + R * 1.074 - 0.002 < 1) //transform only when the correlation of normal RV is smaller than 1
1686  {
1687  R = R * R * (-0.069) + R * 1.074 - 0.002;
1688  }
1689  for (uint16_t blockInd = 0; blockInd < m_numNonSelfBlocking; blockInd++)
1690  {
1691 
1692  //Generate a new correlated normal RV with the following formula
1693  params->m_nonSelfBlocking[blockInd][PHI_INDEX] =
1694  R * params->m_nonSelfBlocking[blockInd][PHI_INDEX] + sqrt (1 - R * R) * m_normalRv->GetValue ();
1695  }
1696  }
1697 
1698  }
1699 
1700  //step c: Determine the attenuation of each blocker due to blockers
1701  for (uint8_t cInd = 0; cInd < clusterNum; cInd++)
1702  {
1703  NS_ASSERT_MSG (clusterAOA[cInd] >= 0 && clusterAOA[cInd] <= 360, "the AOA should be the range of [0,360]");
1704  NS_ASSERT_MSG (clusterZOA[cInd] >= 0 && clusterZOA[cInd] <= 180, "the ZOA should be the range of [0,180]");
1705 
1706  //check self blocking
1707  NS_LOG_INFO ("AOA=" << clusterAOA[cInd] << " Block Region[" << phi_sb - x_sb / 2 << "," << phi_sb + x_sb / 2 << "]");
1708  NS_LOG_INFO ("ZOA=" << clusterZOA[cInd] << " Block Region[" << theta_sb - y_sb / 2 << "," << theta_sb + y_sb / 2 << "]");
1709  if ( std::abs (clusterAOA[cInd] - phi_sb) < (x_sb / 2) && std::abs (clusterZOA[cInd] - theta_sb) < (y_sb / 2))
1710  {
1711  powerAttenuation[cInd] += 30; //anttenuate by 30 dB.
1712  NS_LOG_INFO ("Cluster[" << (int)cInd << "] is blocked by self blocking region and reduce 30 dB power,"
1713  "the attenuation is [" << powerAttenuation[cInd] << " dB]");
1714  }
1715 
1716  //check non-self blocking
1717  double phiK, xK, thetaK, yK;
1718  for (uint16_t blockInd = 0; blockInd < m_numNonSelfBlocking; blockInd++)
1719  {
1720  //The normal RV is transformed to uniform RV with the desired correlation.
1721  phiK = (0.5 * erfc (-1 * params->m_nonSelfBlocking[blockInd][PHI_INDEX] / sqrt (2))) * 360;
1722  while (phiK > 360)
1723  {
1724  phiK -= 360;
1725  }
1726 
1727  while (phiK < 0)
1728  {
1729  phiK += 360;
1730  }
1731 
1732  xK = params->m_nonSelfBlocking[blockInd][X_INDEX];
1733  thetaK = params->m_nonSelfBlocking[blockInd][THETA_INDEX];
1734  yK = params->m_nonSelfBlocking[blockInd][Y_INDEX];
1735  NS_LOG_INFO ("AOA=" << clusterAOA[cInd] << " Block Region[" << phiK - xK << "," << phiK + xK << "]");
1736  NS_LOG_INFO ("ZOA=" << clusterZOA[cInd] << " Block Region[" << thetaK - yK << "," << thetaK + yK << "]");
1737 
1738  if ( std::abs (clusterAOA[cInd] - phiK) < (xK)
1739  && std::abs (clusterZOA[cInd] - thetaK) < (yK))
1740  {
1741  double A1 = clusterAOA[cInd] - (phiK + xK / 2); //(7.6-24)
1742  double A2 = clusterAOA[cInd] - (phiK - xK / 2); //(7.6-25)
1743  double Z1 = clusterZOA[cInd] - (thetaK + yK / 2); //(7.6-26)
1744  double Z2 = clusterZOA[cInd] - (thetaK - yK / 2); //(7.6-27)
1745  int signA1, signA2, signZ1, signZ2;
1746  //draw sign for the above parameters according to table 7.6.4.1-3 Description of signs
1747  if (xK / 2 < clusterAOA[cInd] - phiK && clusterAOA[cInd] - phiK <= xK)
1748  {
1749  signA1 = -1;
1750  }
1751  else
1752  {
1753  signA1 = 1;
1754  }
1755  if (-1 * xK < clusterAOA[cInd] - phiK && clusterAOA[cInd] - phiK <= -1 * xK / 2)
1756  {
1757  signA2 = -1;
1758  }
1759  else
1760  {
1761  signA2 = 1;
1762  }
1763 
1764  if (yK / 2 < clusterZOA[cInd] - thetaK && clusterZOA[cInd] - thetaK <= yK)
1765  {
1766  signZ1 = -1;
1767  }
1768  else
1769  {
1770  signZ1 = 1;
1771  }
1772  if (-1 * yK < clusterZOA[cInd] - thetaK && clusterZOA[cInd] - thetaK <= -1 * yK / 2)
1773  {
1774  signZ2 = -1;
1775  }
1776  else
1777  {
1778  signZ2 = 1;
1779  }
1780  double lambda = 3e8 / m_frequency;
1781  double F_A1 = atan (signA1 * M_PI / 2 * sqrt (M_PI / lambda *
1782  params->m_nonSelfBlocking[blockInd][R_INDEX] * (1 / cos (A1 * M_PI / 180) - 1))) / M_PI; //(7.6-23)
1783  double F_A2 = atan (signA2 * M_PI / 2 * sqrt (M_PI / lambda *
1784  params->m_nonSelfBlocking[blockInd][R_INDEX] * (1 / cos (A2 * M_PI / 180) - 1))) / M_PI;
1785  double F_Z1 = atan (signZ1 * M_PI / 2 * sqrt (M_PI / lambda *
1786  params->m_nonSelfBlocking[blockInd][R_INDEX] * (1 / cos (Z1 * M_PI / 180) - 1))) / M_PI;
1787  double F_Z2 = atan (signZ2 * M_PI / 2 * sqrt (M_PI / lambda *
1788  params->m_nonSelfBlocking[blockInd][R_INDEX] * (1 / cos (Z2 * M_PI / 180) - 1))) / M_PI;
1789  double L_dB = -20 * log10 (1 - (F_A1 + F_A2) * (F_Z1 + F_Z2)); //(7.6-22)
1790  powerAttenuation[cInd] += L_dB;
1791  NS_LOG_INFO ("Cluster[" << (int)cInd << "] is blocked by no-self blocking, "
1792  "the loss is [" << L_dB << "]" << " dB");
1793 
1794  }
1795  }
1796  }
1797  return powerAttenuation;
1798 }
1799 
1800 int64_t
1802 {
1803  NS_LOG_FUNCTION (this << stream);
1804  m_normalRv->SetStream (stream);
1805  m_uniformRv->SetStream (stream + 1);
1806  return 2;
1807 }
1808 
1809 } // namespace ns3
static const double offSetAlpha[20]
static const double sqrtC_UMa_LOS[7][7]
Ptr< const AttributeChecker > MakeStringChecker(void)
Definition: string.cc:30
Double3DVector m_clusterPhase
the initial random phases
std::unordered_map< uint32_t, Ptr< ThreeGppChannelMatrix > > m_channelMap
map containing the channel realizations
Smart pointer class similar to boost::intrusive_ptr.
Definition: ptr.h:73
#define NS_LOG_FUNCTION(parameters)
If log level LOG_FUNCTION is enabled, this macro will output all input parameters separated by "...
void SetStream(int64_t stream)
Specifies the stream number for the RngStream.
AttributeValue implementation for Boolean.
Definition: boolean.h:36
static const double sqrtC_office_NLOS[6][6]
Vector m_preLocUT
location of UT when generating the previous channel
uint32_t GetId(void) const
Definition: node.cc:107
std::vector< Double2DVector > Double3DVector
type definition for 3D matrices of doubles
#define NS_OBJECT_ENSURE_REGISTERED(type)
Register an Object subclass with the TypeId system.
Definition: object-base.h:45
NS_ASSERT_MSG(false, "Ipv4AddressGenerator::MaskToIndex(): Impossible")
Hold variables of type string.
Definition: string.h:41
#define min(a, b)
Definition: 80211b.c:42
static const double sqrtC_UMi_O2I[6][6]
DoubleVector CalcAttenuationOfBlockage(Ptr< ThreeGppChannelMatrix > params, const DoubleVector &clusterAOA, const DoubleVector &clusterZOA) const
Applies the blockage model A described in 3GPP TR 38.901.
std::vector< Complex2DVector > Complex3DVector
type definition for complex 3D matrices
static const uint8_t THETA_INDEX
index of the THETA value in the m_nonSelfBlocking array
Ptr< const AttributeAccessor > MakeBooleanAccessor(T1 a1)
Create an AttributeAccessor for a class data member, or a lone class get functor or set method...
Definition: boolean.h:85
Ptr< NormalRandomVariable > m_normalRv
normal random variable
double GetSeconds(void) const
Get an approximation of the time stored in this instance in the indicated unit.
Definition: nstime.h:361
std::string GetScenario(void) const
Returns the propagation scenario.
Hold a signed integer type.
Definition: integer.h:44
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition: log.h:205
Time MilliSeconds(uint64_t value)
Construct a Time in the indicated unit.
Definition: nstime.h:1078
std::vector< DoubleVector > Double2DVector
type definition for matrices of doubles
#define NS_LOG_INFO(msg)
Use NS_LOG to output a message of level LOG_INFO.
Definition: log.h:281
#define NS_FATAL_ERROR(msg)
Report a fatal error with a message and terminate.
Definition: fatal-error.h:162
double theta
the inclination angle in radians
Definition: angles.h:117
static const uint8_t R_INDEX
index of the R value in the m_nonSelfBlocking array
void SetChannelConditionModel(Ptr< ChannelConditionModel > model)
Set the channel condition model.
Ptr< ChannelConditionModel > m_channelConditionModel
the channel condition model
std::string m_scenario
the 3GPP scenario
static const double sqrtC_UMa_O2I[6][6]
Ptr< const ThreeGppChannelMatrix > GetChannel(Ptr< const MobilityModel > aMob, Ptr< const MobilityModel > bMob, Ptr< const ThreeGppAntennaArrayModel > aAntenna, Ptr< const ThreeGppAntennaArrayModel > bAntenna)
Looks for the channel matrix associated to the aMob and bMob pair in m_channelMap.
Complex3DVector m_channel
channel matrix H[u][s][n].
static const uint8_t X_INDEX
index of the X value in the m_nonSelfBlocking array
Ptr< const AttributeChecker > MakeTimeChecker(const Time min, const Time max)
Helper to make a Time checker with bounded range.
Definition: time.cc:449
Double2DVector m_nonSelfBlocking
store the blockages
Ptr< const AttributeAccessor > MakeIntegerAccessor(T1 a1)
Create an AttributeAccessor for a class data member, or a lone class get functor or set method...
Definition: integer.h:45
Ptr< const ParamsTable > GetThreeGppTable(bool los, bool o2i, double hBS, double hUT, double distance2D) const
Get the parameters needed to apply the channel generation procedure.
static const double sqrtC_office_LOS[7][7]
void SetFrequency(double f)
Sets the center frequency of the model.
double GetFrequency(void) const
Returns the center frequency.
Ptr< const AttributeAccessor > MakePointerAccessor(T1 a1)
Create an AttributeAccessor for a class data member, or a lone class get functor or set method...
Definition: pointer.h:227
#define max(a, b)
Definition: 80211b.c:43
bool IsZero(void) const
Definition: nstime.h:294
AttributeValue implementation for Time.
Definition: nstime.h:1132
uint16_t m_numNonSelfBlocking
number of non-self-blocking regions
static const double sqrtC_RMa_LOS[7][7]
Time m_updatePeriod
the channel update period
int64_t AssignStreams(int64_t stream)
Assign a fixed random variable stream number to the random variables used by this model...
Channel Matrix Generation following 3GPP TR 38.901.
bool m_portraitMode
true if potrait mode, false if landscape
bool ChannelMatrixNeedsUpdate(Ptr< const ThreeGppChannelMatrix > channelMatrix, bool isLos) const
Check if the channel matrix has to be updated.
double m_blockerSpeed
the blocker speed
Double2DVector m_angle
cluster angle angle[direction][n], where direction = 0(aoa), 1(zoa), 2(aod), 3(zod) in degree...
double f(double x, void *params)
Definition: 80211b.c:70
Every class exported by the ns3 library is enclosed in the ns3 namespace.
Hold objects of type Ptr<T>.
Definition: pointer.h:36
std::pair< uint32_t, uint32_t > m_nodeIds
the first element is the s-node ID, the second element is the u-node ID
int64_t GetNanoSeconds(void) const
Get an approximation of the time stored in this instance in the indicated unit.
Definition: nstime.h:373
Ptr< const AttributeChecker > MakeBooleanChecker(void)
Definition: boolean.cc:121
double GetValue(double min, double max)
Get the next random value, as a double in the specified range .
Ptr< UniformRandomVariable > m_uniformRv
uniform random variable
Ptr< const AttributeAccessor > MakeTimeAccessor(T1 a1)
Create an AttributeAccessor for a class data member, or a lone class get functor or set method...
Definition: nstime.h:1133
static Time Now(void)
Return the current simulation virtual time.
Definition: simulator.cc:193
static const double sqrtC_RMa_NLOS[6][6]
static const double sqrtC_RMa_O2I[6][6]
Ptr< const AttributeAccessor > MakeDoubleAccessor(T1 a1)
Create an AttributeAccessor for a class data member, or a lone class get functor or set method...
Definition: double.h:42
double max(double x, double y)
void SetScenario(const std::string &scenario)
Sets the propagation scenario.
std::vector< double > DoubleVector
type definition for vectors of doubles
static const uint8_t PHI_INDEX
index of the PHI value in the m_nonSelfBlocking array
if(desigRtr==addrLocal)
bool m_blockage
enables the blockage model A
static TypeId GetTypeId()
Get the type ID.
A network Node.
Definition: node.h:56
#define NS_LOG_DEBUG(msg)
Use NS_LOG to output a message of level LOG_DEBUG.
Definition: log.h:273
double phi
the azimuth angle in radians
Definition: angles.h:111
double min(double x, double y)
Time Now(void)
create an ns3::Time instance which contains the current simulation time.
Definition: simulator.cc:309
struct holding the azimuth and inclination angles of spherical coordinates.
Definition: angles.h:71
A base class which provides memory management and object aggregation.
Definition: object.h:87
Ptr< const AttributeAccessor > MakeStringAccessor(T1 a1)
Create an AttributeAccessor for a class data member, or a lone class get functor or set method...
Definition: string.h:42
static const double sqrtC_UMa_NLOS[6][6]
This class can be used to hold variables of floating point type such as &#39;double&#39; or &#39;float&#39;...
Definition: double.h:41
Ptr< ChannelConditionModel > GetChannelConditionModel() const
Get the associated channel condition model.
static const uint8_t Y_INDEX
index of the Y value in the m_nonSelfBlocking array
a unique identifier for an interface.
Definition: type-id.h:58
static const double sqrtC_UMi_LOS[7][7]
Ptr< ThreeGppChannelMatrix > GetNewChannel(Vector locUT, bool los, bool o2i, Ptr< const ThreeGppAntennaArrayModel > sAntenna, Ptr< const ThreeGppAntennaArrayModel > uAntenna, Angles &uAngle, Angles &sAngle, double dis2D, double hBS, double hUT) const
Compute the channel matrix between two devices using the procedure described in 3GPP TR 38...
TypeId SetParent(TypeId tid)
Set the parent TypeId.
Definition: type-id.cc:923
static const double sqrtC_UMi_NLOS[6][6]
static constexpr uint32_t GetKey(uint32_t x1, uint32_t x2)
Calculate the channel key using the Cantor function.
double m_frequency
the operating frequency