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