A Discrete-Event Network Simulator
API
three-gpp-channel-model.cc
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2019 SIGNET Lab, Department of Information Engineering,
3  * University of Padova
4  * Copyright (c) 2015, NYU WIRELESS, Tandon School of Engineering,
5  * New York University
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License version 2 as
9  * published by the Free Software Foundation;
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
23 
24 #include "ns3/double.h"
25 #include "ns3/integer.h"
26 #include "ns3/log.h"
27 #include "ns3/mobility-model.h"
28 #include "ns3/node.h"
29 #include "ns3/phased-array-model.h"
30 #include "ns3/pointer.h"
31 #include "ns3/string.h"
32 #include <ns3/simulator.h>
33 
34 #include <algorithm>
35 #include <random>
36 
37 namespace ns3
38 {
39 
40 NS_LOG_COMPONENT_DEFINE("ThreeGppChannelModel");
41 
42 NS_OBJECT_ENSURE_REGISTERED(ThreeGppChannelModel);
43 
46 static const double offSetAlpha[20] = {
47  0.0447, -0.0447, 0.1413, -0.1413, 0.2492, -0.2492, 0.3715, -0.3715, 0.5129, -0.5129,
48  0.6797, -0.6797, 0.8844, -0.8844, 1.1481, -1.1481, 1.5195, -1.5195, 2.1551, -2.1551,
49 };
50 
59 static const double sqrtC_RMa_LOS[7][7] = {
60  {1, 0, 0, 0, 0, 0, 0},
61  {0, 1, 0, 0, 0, 0, 0},
62  {-0.5, 0, 0.866025, 0, 0, 0, 0},
63  {0, 0, 0, 1, 0, 0, 0},
64  {0, 0, 0, 0, 1, 0, 0},
65  {0.01, 0, -0.0519615, 0.73, -0.2, 0.651383, 0},
66  {-0.17, -0.02, 0.21362, -0.14, 0.24, 0.142773, 0.909661},
67 };
68 
78 static const double sqrtC_RMa_NLOS[6][6] = {
79  {1, 0, 0, 0, 0, 0},
80  {-0.5, 0.866025, 0, 0, 0, 0},
81  {0.6, -0.11547, 0.791623, 0, 0, 0},
82  {0, 0, 0, 1, 0, 0},
83  {-0.04, -0.138564, 0.540662, -0.18, 0.809003, 0},
84  {-0.25, -0.606218, -0.240013, 0.26, -0.231685, 0.625392},
85 };
86 
95 static const double sqrtC_RMa_O2I[6][6] = {
96  {1, 0, 0, 0, 0, 0},
97  {0, 1, 0, 0, 0, 0},
98  {0, 0, 1, 0, 0, 0},
99  {0, 0, -0.7, 0.714143, 0, 0},
100  {0, 0, 0.66, -0.123225, 0.741091, 0},
101  {0, 0, 0.47, 0.152631, -0.393194, 0.775373},
102 };
103 
112 static const double sqrtC_UMa_LOS[7][7] = {
113  {1, 0, 0, 0, 0, 0, 0},
114  {0, 1, 0, 0, 0, 0, 0},
115  {-0.4, -0.4, 0.824621, 0, 0, 0, 0},
116  {-0.5, 0, 0.242536, 0.83137, 0, 0, 0},
117  {-0.5, -0.2, 0.630593, -0.484671, 0.278293, 0, 0},
118  {0, 0, -0.242536, 0.672172, 0.642214, 0.27735, 0},
119  {-0.8, 0, -0.388057, -0.367926, 0.238537, -3.58949e-15, 0.130931},
120 };
121 
131 static const double sqrtC_UMa_NLOS[6][6] = {
132  {1, 0, 0, 0, 0, 0},
133  {-0.4, 0.916515, 0, 0, 0, 0},
134  {-0.6, 0.174574, 0.78072, 0, 0, 0},
135  {0, 0.654654, 0.365963, 0.661438, 0, 0},
136  {0, -0.545545, 0.762422, 0.118114, 0.327327, 0},
137  {-0.4, -0.174574, -0.396459, 0.392138, 0.49099, 0.507445},
138 };
139 
148 static const double sqrtC_UMa_O2I[6][6] = {
149  {1, 0, 0, 0, 0, 0},
150  {-0.5, 0.866025, 0, 0, 0, 0},
151  {0.2, 0.57735, 0.791623, 0, 0, 0},
152  {0, 0.46188, -0.336861, 0.820482, 0, 0},
153  {0, -0.69282, 0.252646, 0.493742, 0.460857, 0},
154  {0, -0.23094, 0.16843, 0.808554, -0.220827, 0.464515},
155 
156 };
157 
166 static const double sqrtC_UMi_LOS[7][7] = {
167  {1, 0, 0, 0, 0, 0, 0},
168  {0.5, 0.866025, 0, 0, 0, 0, 0},
169  {-0.4, -0.57735, 0.711805, 0, 0, 0, 0},
170  {-0.5, 0.057735, 0.468293, 0.726201, 0, 0, 0},
171  {-0.4, -0.11547, 0.805464, -0.23482, 0.350363, 0, 0},
172  {0, 0, 0, 0.688514, 0.461454, 0.559471, 0},
173  {0, 0, 0.280976, 0.231921, -0.490509, 0.11916, 0.782603},
174 };
175 
185 static const double sqrtC_UMi_NLOS[6][6] = {
186  {1, 0, 0, 0, 0, 0},
187  {-0.7, 0.714143, 0, 0, 0, 0},
188  {0, 0, 1, 0, 0, 0},
189  {-0.4, 0.168034, 0, 0.90098, 0, 0},
190  {0, -0.70014, 0.5, 0.130577, 0.4927, 0},
191  {0, 0, 0.5, 0.221981, -0.566238, 0.616522},
192 };
193 
202 static const double sqrtC_UMi_O2I[6][6] = {
203  {1, 0, 0, 0, 0, 0},
204  {-0.5, 0.866025, 0, 0, 0, 0},
205  {0.2, 0.57735, 0.791623, 0, 0, 0},
206  {0, 0.46188, -0.336861, 0.820482, 0, 0},
207  {0, -0.69282, 0.252646, 0.493742, 0.460857, 0},
208  {0, -0.23094, 0.16843, 0.808554, -0.220827, 0.464515},
209 };
210 
219 static const double sqrtC_office_LOS[7][7] = {
220  {1, 0, 0, 0, 0, 0, 0},
221  {0.5, 0.866025, 0, 0, 0, 0, 0},
222  {-0.8, -0.11547, 0.588784, 0, 0, 0, 0},
223  {-0.4, 0.23094, 0.520847, 0.717903, 0, 0, 0},
224  {-0.5, 0.288675, 0.73598, -0.348236, 0.0610847, 0, 0},
225  {0.2, -0.11547, 0.418943, 0.541106, 0.219905, 0.655744, 0},
226  {0.3, -0.057735, 0.73598, -0.348236, 0.0610847, -0.304997, 0.383375},
227 };
228 
238 static const double sqrtC_office_NLOS[6][6] = {
239  {1, 0, 0, 0, 0, 0},
240  {-0.5, 0.866025, 0, 0, 0, 0},
241  {0, 0.46188, 0.886942, 0, 0, 0},
242  {-0.4, -0.23094, 0.120263, 0.878751, 0, 0},
243  {0, -0.311769, 0.55697, -0.249198, 0.728344, 0},
244  {0, -0.069282, 0.295397, 0.430696, 0.468462, 0.709214},
245 };
246 
248 {
249  NS_LOG_FUNCTION(this);
250  m_uniformRv = CreateObject<UniformRandomVariable>();
251  m_uniformRvShuffle = CreateObject<UniformRandomVariable>();
252  m_uniformRvDoppler = CreateObject<UniformRandomVariable>();
253 
254  m_normalRv = CreateObject<NormalRandomVariable>();
255  m_normalRv->SetAttribute("Mean", DoubleValue(0.0));
256  m_normalRv->SetAttribute("Variance", DoubleValue(1.0));
257 }
258 
260 {
261  NS_LOG_FUNCTION(this);
262 }
263 
264 void
266 {
267  NS_LOG_FUNCTION(this);
269  {
270  m_channelConditionModel->Dispose();
271  }
272  m_channelMatrixMap.clear();
273  m_channelParamsMap.clear();
274  m_channelConditionModel = nullptr;
275 }
276 
277 TypeId
279 {
280  static TypeId tid =
281  TypeId("ns3::ThreeGppChannelModel")
282  .SetGroupName("Spectrum")
284  .AddConstructor<ThreeGppChannelModel>()
285  .AddAttribute("Frequency",
286  "The operating Frequency in Hz",
287  DoubleValue(500.0e6),
290  MakeDoubleChecker<double>())
291  .AddAttribute(
292  "Scenario",
293  "The 3GPP scenario (RMa, UMa, UMi-StreetCanyon, InH-OfficeOpen, InH-OfficeMixed)",
294  StringValue("UMa"),
298  .AddAttribute("ChannelConditionModel",
299  "Pointer to the channel condition model",
300  PointerValue(),
303  MakePointerChecker<ChannelConditionModel>())
304  .AddAttribute("UpdatePeriod",
305  "Specify the channel coherence time",
308  MakeTimeChecker())
309  // attributes for the blockage model
310  .AddAttribute("Blockage",
311  "Enable blockage model A (sec 7.6.4.1)",
312  BooleanValue(false),
315  .AddAttribute("NumNonselfBlocking",
316  "number of non-self-blocking regions",
317  IntegerValue(4),
319  MakeIntegerChecker<uint16_t>())
320  .AddAttribute("PortraitMode",
321  "true for portrait mode, false for landscape mode",
322  BooleanValue(true),
325  .AddAttribute("BlockerSpeed",
326  "The speed of moving blockers, the unit is m/s",
327  DoubleValue(1),
329  MakeDoubleChecker<double>())
330  .AddAttribute("vScatt",
331  "Maximum speed of the vehicle in the layout (see 3GPP TR 37.885 v15.3.0, "
332  "Sec. 6.2.3)."
333  "Used to compute the additional contribution for the Doppler of"
334  "delayed (reflected) paths",
335  DoubleValue(0.0),
337  MakeDoubleChecker<double>(0.0))
338 
339  ;
340  return tid;
341 }
342 
343 void
345 {
346  NS_LOG_FUNCTION(this);
347  m_channelConditionModel = model;
348 }
349 
352 {
353  NS_LOG_FUNCTION(this);
355 }
356 
357 void
359 {
360  NS_LOG_FUNCTION(this);
361  NS_ASSERT_MSG(f >= 500.0e6 && f <= 100.0e9,
362  "Frequency should be between 0.5 and 100 GHz but is " << f);
363  m_frequency = f;
364 }
365 
366 double
368 {
369  NS_LOG_FUNCTION(this);
370  return m_frequency;
371 }
372 
373 void
374 ThreeGppChannelModel::SetScenario(const std::string& scenario)
375 {
376  NS_LOG_FUNCTION(this);
377  NS_ASSERT_MSG(scenario == "RMa" || scenario == "UMa" || scenario == "UMi-StreetCanyon" ||
378  scenario == "InH-OfficeOpen" || scenario == "InH-OfficeMixed" ||
379  scenario == "V2V-Urban" || scenario == "V2V-Highway",
380  "Unknown scenario, choose between: RMa, UMa, UMi-StreetCanyon, "
381  "InH-OfficeOpen, InH-OfficeMixed, V2V-Urban or V2V-Highway");
382  m_scenario = scenario;
383 }
384 
385 std::string
387 {
388  NS_LOG_FUNCTION(this);
389  return m_scenario;
390 }
391 
394  double hBS,
395  double hUT,
396  double distance2D) const
397 {
398  NS_LOG_FUNCTION(this);
399 
400  double fcGHz = m_frequency / 1.0e9;
401  Ptr<ParamsTable> table3gpp = Create<ParamsTable>();
402  // table3gpp includes the following parameters:
403  // numOfCluster, raysPerCluster, uLgDS, sigLgDS, uLgASD, sigLgASD,
404  // uLgASA, sigLgASA, uLgZSA, sigLgZSA, uLgZSD, sigLgZSD, offsetZOD,
405  // cDS, cASD, cASA, cZSA, uK, sigK, rTau, uXpr, sigXpr, shadowingStd
406 
407  bool los = channelCondition->IsLos();
408  bool o2i = channelCondition->IsO2i();
409 
410  // In NLOS case, parameter uK and sigK are not used and they are set to 0
411  if (m_scenario == "RMa")
412  {
413  if (los && !o2i)
414  {
415  // 3GPP mentioned that 3.91 ns should be used when the Cluster DS (cDS)
416  // entry is N/A.
417  table3gpp->m_numOfCluster = 11;
418  table3gpp->m_raysPerCluster = 20;
419  table3gpp->m_uLgDS = -7.49;
420  table3gpp->m_sigLgDS = 0.55;
421  table3gpp->m_uLgASD = 0.90;
422  table3gpp->m_sigLgASD = 0.38;
423  table3gpp->m_uLgASA = 1.52;
424  table3gpp->m_sigLgASA = 0.24;
425  table3gpp->m_uLgZSA = 0.47;
426  table3gpp->m_sigLgZSA = 0.40;
427  table3gpp->m_uLgZSD = 0.34;
428  table3gpp->m_sigLgZSD =
429  std::max(-1.0, -0.17 * (distance2D / 1000.0) - 0.01 * (hUT - 1.5) + 0.22);
430  table3gpp->m_offsetZOD = 0;
431  table3gpp->m_cDS = 3.91e-9;
432  table3gpp->m_cASD = 2;
433  table3gpp->m_cASA = 3;
434  table3gpp->m_cZSA = 3;
435  table3gpp->m_uK = 7;
436  table3gpp->m_sigK = 4;
437  table3gpp->m_rTau = 3.8;
438  table3gpp->m_uXpr = 12;
439  table3gpp->m_sigXpr = 4;
440  table3gpp->m_perClusterShadowingStd = 3;
441 
442  for (uint8_t row = 0; row < 7; row++)
443  {
444  for (uint8_t column = 0; column < 7; column++)
445  {
446  table3gpp->m_sqrtC[row][column] = sqrtC_RMa_LOS[row][column];
447  }
448  }
449  }
450  else if (!los && !o2i)
451  {
452  table3gpp->m_numOfCluster = 10;
453  table3gpp->m_raysPerCluster = 20;
454  table3gpp->m_uLgDS = -7.43;
455  table3gpp->m_sigLgDS = 0.48;
456  table3gpp->m_uLgASD = 0.95;
457  table3gpp->m_sigLgASD = 0.45;
458  table3gpp->m_uLgASA = 1.52;
459  table3gpp->m_sigLgASA = 0.13;
460  table3gpp->m_uLgZSA = 0.58;
461  table3gpp->m_sigLgZSA = 0.37;
462  table3gpp->m_uLgZSD =
463  std::max(-1.0, -0.19 * (distance2D / 1000.0) - 0.01 * (hUT - 1.5) + 0.28);
464  table3gpp->m_sigLgZSD = 0.30;
465  table3gpp->m_offsetZOD = atan((35 - 3.5) / distance2D) - atan((35 - 1.5) / distance2D);
466  table3gpp->m_cDS = 3.91e-9;
467  table3gpp->m_cASD = 2;
468  table3gpp->m_cASA = 3;
469  table3gpp->m_cZSA = 3;
470  table3gpp->m_uK = 0;
471  table3gpp->m_sigK = 0;
472  table3gpp->m_rTau = 1.7;
473  table3gpp->m_uXpr = 7;
474  table3gpp->m_sigXpr = 3;
475  table3gpp->m_perClusterShadowingStd = 3;
476 
477  for (uint8_t row = 0; row < 6; row++)
478  {
479  for (uint8_t column = 0; column < 6; column++)
480  {
481  table3gpp->m_sqrtC[row][column] = sqrtC_RMa_NLOS[row][column];
482  }
483  }
484  }
485  else // o2i
486  {
487  table3gpp->m_numOfCluster = 10;
488  table3gpp->m_raysPerCluster = 20;
489  table3gpp->m_uLgDS = -7.47;
490  table3gpp->m_sigLgDS = 0.24;
491  table3gpp->m_uLgASD = 0.67;
492  table3gpp->m_sigLgASD = 0.18;
493  table3gpp->m_uLgASA = 1.66;
494  table3gpp->m_sigLgASA = 0.21;
495  table3gpp->m_uLgZSA = 0.93;
496  table3gpp->m_sigLgZSA = 0.22;
497  table3gpp->m_uLgZSD =
498  std::max(-1.0, -0.19 * (distance2D / 1000.0) - 0.01 * (hUT - 1.5) + 0.28);
499  table3gpp->m_sigLgZSD = 0.30;
500  table3gpp->m_offsetZOD = atan((35 - 3.5) / distance2D) - atan((35 - 1.5) / distance2D);
501  table3gpp->m_cDS = 3.91e-9;
502  table3gpp->m_cASD = 2;
503  table3gpp->m_cASA = 3;
504  table3gpp->m_cZSA = 3;
505  table3gpp->m_uK = 0;
506  table3gpp->m_sigK = 0;
507  table3gpp->m_rTau = 1.7;
508  table3gpp->m_uXpr = 7;
509  table3gpp->m_sigXpr = 3;
510  table3gpp->m_perClusterShadowingStd = 3;
511 
512  for (uint8_t row = 0; row < 6; row++)
513  {
514  for (uint8_t column = 0; column < 6; column++)
515  {
516  table3gpp->m_sqrtC[row][column] = sqrtC_RMa_O2I[row][column];
517  }
518  }
519  }
520  }
521  else if (m_scenario == "UMa")
522  {
523  if (los && !o2i)
524  {
525  table3gpp->m_numOfCluster = 12;
526  table3gpp->m_raysPerCluster = 20;
527  table3gpp->m_uLgDS = -6.955 - 0.0963 * log10(fcGHz);
528  table3gpp->m_sigLgDS = 0.66;
529  table3gpp->m_uLgASD = 1.06 + 0.1114 * log10(fcGHz);
530  table3gpp->m_sigLgASD = 0.28;
531  table3gpp->m_uLgASA = 1.81;
532  table3gpp->m_sigLgASA = 0.20;
533  table3gpp->m_uLgZSA = 0.95;
534  table3gpp->m_sigLgZSA = 0.16;
535  table3gpp->m_uLgZSD =
536  std::max(-0.5, -2.1 * distance2D / 1000.0 - 0.01 * (hUT - 1.5) + 0.75);
537  table3gpp->m_sigLgZSD = 0.40;
538  table3gpp->m_offsetZOD = 0;
539  table3gpp->m_cDS = std::max(0.25, -3.4084 * log10(fcGHz) + 6.5622) * 1e-9;
540  table3gpp->m_cASD = 5;
541  table3gpp->m_cASA = 11;
542  table3gpp->m_cZSA = 7;
543  table3gpp->m_uK = 9;
544  table3gpp->m_sigK = 3.5;
545  table3gpp->m_rTau = 2.5;
546  table3gpp->m_uXpr = 8;
547  table3gpp->m_sigXpr = 4;
548  table3gpp->m_perClusterShadowingStd = 3;
549 
550  for (uint8_t row = 0; row < 7; row++)
551  {
552  for (uint8_t column = 0; column < 7; column++)
553  {
554  table3gpp->m_sqrtC[row][column] = sqrtC_UMa_LOS[row][column];
555  }
556  }
557  }
558  else
559  {
560  double uLgZSD = std::max(-0.5, -2.1 * distance2D / 1000.0 - 0.01 * (hUT - 1.5) + 0.9);
561 
562  double afc = 0.208 * log10(fcGHz) - 0.782;
563  double bfc = 25;
564  double cfc = -0.13 * log10(fcGHz) + 2.03;
565  double efc = 7.66 * log10(fcGHz) - 5.96;
566 
567  double offsetZOD = efc - std::pow(10, afc * log10(std::max(bfc, distance2D)) + cfc);
568 
569  if (!los && !o2i)
570  {
571  table3gpp->m_numOfCluster = 20;
572  table3gpp->m_raysPerCluster = 20;
573  table3gpp->m_uLgDS = -6.28 - 0.204 * log10(fcGHz);
574  table3gpp->m_sigLgDS = 0.39;
575  table3gpp->m_uLgASD = 1.5 - 0.1144 * log10(fcGHz);
576  table3gpp->m_sigLgASD = 0.28;
577  table3gpp->m_uLgASA = 2.08 - 0.27 * log10(fcGHz);
578  table3gpp->m_sigLgASA = 0.11;
579  table3gpp->m_uLgZSA = -0.3236 * log10(fcGHz) + 1.512;
580  table3gpp->m_sigLgZSA = 0.16;
581  table3gpp->m_uLgZSD = uLgZSD;
582  table3gpp->m_sigLgZSD = 0.49;
583  table3gpp->m_offsetZOD = offsetZOD;
584  table3gpp->m_cDS = std::max(0.25, -3.4084 * log10(fcGHz) + 6.5622) * 1e-9;
585  table3gpp->m_cASD = 2;
586  table3gpp->m_cASA = 15;
587  table3gpp->m_cZSA = 7;
588  table3gpp->m_uK = 0;
589  table3gpp->m_sigK = 0;
590  table3gpp->m_rTau = 2.3;
591  table3gpp->m_uXpr = 7;
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_UMa_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.49;
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_UMa_O2I[row][column];
634  }
635  }
636  }
637  }
638  }
639  else if (m_scenario == "UMi-StreetCanyon")
640  {
641  if (los && !o2i)
642  {
643  table3gpp->m_numOfCluster = 12;
644  table3gpp->m_raysPerCluster = 20;
645  table3gpp->m_uLgDS = -0.24 * log10(1 + fcGHz) - 7.14;
646  table3gpp->m_sigLgDS = 0.38;
647  table3gpp->m_uLgASD = -0.05 * log10(1 + fcGHz) + 1.21;
648  table3gpp->m_sigLgASD = 0.41;
649  table3gpp->m_uLgASA = -0.08 * log10(1 + fcGHz) + 1.73;
650  table3gpp->m_sigLgASA = 0.014 * log10(1 + fcGHz) + 0.28;
651  table3gpp->m_uLgZSA = -0.1 * log10(1 + fcGHz) + 0.73;
652  table3gpp->m_sigLgZSA = -0.04 * log10(1 + fcGHz) + 0.34;
653  table3gpp->m_uLgZSD =
654  std::max(-0.21, -14.8 * distance2D / 1000.0 + 0.01 * std::abs(hUT - hBS) + 0.83);
655  table3gpp->m_sigLgZSD = 0.35;
656  table3gpp->m_offsetZOD = 0;
657  table3gpp->m_cDS = 5e-9;
658  table3gpp->m_cASD = 3;
659  table3gpp->m_cASA = 17;
660  table3gpp->m_cZSA = 7;
661  table3gpp->m_uK = 9;
662  table3gpp->m_sigK = 5;
663  table3gpp->m_rTau = 3;
664  table3gpp->m_uXpr = 9;
665  table3gpp->m_sigXpr = 3;
666  table3gpp->m_perClusterShadowingStd = 3;
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_UMi_LOS[row][column];
673  }
674  }
675  }
676  else
677  {
678  double uLgZSD =
679  std::max(-0.5, -3.1 * distance2D / 1000.0 + 0.01 * std::max(hUT - hBS, 0.0) + 0.2);
680  double offsetZOD = -1 * std::pow(10, -1.5 * log10(std::max(10.0, distance2D)) + 3.3);
681  if (!los && !o2i)
682  {
683  table3gpp->m_numOfCluster = 19;
684  table3gpp->m_raysPerCluster = 20;
685  table3gpp->m_uLgDS = -0.24 * log10(1 + fcGHz) - 6.83;
686  table3gpp->m_sigLgDS = 0.16 * log10(1 + fcGHz) + 0.28;
687  table3gpp->m_uLgASD = -0.23 * log10(1 + fcGHz) + 1.53;
688  table3gpp->m_sigLgASD = 0.11 * log10(1 + fcGHz) + 0.33;
689  table3gpp->m_uLgASA = -0.08 * log10(1 + fcGHz) + 1.81;
690  table3gpp->m_sigLgASA = 0.05 * log10(1 + fcGHz) + 0.3;
691  table3gpp->m_uLgZSA = -0.04 * log10(1 + fcGHz) + 0.92;
692  table3gpp->m_sigLgZSA = -0.07 * log10(1 + fcGHz) + 0.41;
693  table3gpp->m_uLgZSD = uLgZSD;
694  table3gpp->m_sigLgZSD = 0.35;
695  table3gpp->m_offsetZOD = offsetZOD;
696  table3gpp->m_cDS = 11e-9;
697  table3gpp->m_cASD = 10;
698  table3gpp->m_cASA = 22;
699  table3gpp->m_cZSA = 7;
700  table3gpp->m_uK = 0;
701  table3gpp->m_sigK = 0;
702  table3gpp->m_rTau = 2.1;
703  table3gpp->m_uXpr = 8;
704  table3gpp->m_sigXpr = 3;
705  table3gpp->m_perClusterShadowingStd = 3;
706 
707  for (uint8_t row = 0; row < 6; row++)
708  {
709  for (uint8_t column = 0; column < 6; column++)
710  {
711  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_NLOS[row][column];
712  }
713  }
714  }
715  else //(o2i)
716  {
717  table3gpp->m_numOfCluster = 12;
718  table3gpp->m_raysPerCluster = 20;
719  table3gpp->m_uLgDS = -6.62;
720  table3gpp->m_sigLgDS = 0.32;
721  table3gpp->m_uLgASD = 1.25;
722  table3gpp->m_sigLgASD = 0.42;
723  table3gpp->m_uLgASA = 1.76;
724  table3gpp->m_sigLgASA = 0.16;
725  table3gpp->m_uLgZSA = 1.01;
726  table3gpp->m_sigLgZSA = 0.43;
727  table3gpp->m_uLgZSD = uLgZSD;
728  table3gpp->m_sigLgZSD = 0.35;
729  table3gpp->m_offsetZOD = offsetZOD;
730  table3gpp->m_cDS = 11e-9;
731  table3gpp->m_cASD = 5;
732  table3gpp->m_cASA = 8;
733  table3gpp->m_cZSA = 3;
734  table3gpp->m_uK = 0;
735  table3gpp->m_sigK = 0;
736  table3gpp->m_rTau = 2.2;
737  table3gpp->m_uXpr = 9;
738  table3gpp->m_sigXpr = 5;
739  table3gpp->m_perClusterShadowingStd = 4;
740 
741  for (uint8_t row = 0; row < 6; row++)
742  {
743  for (uint8_t column = 0; column < 6; column++)
744  {
745  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_O2I[row][column];
746  }
747  }
748  }
749  }
750  }
751  else if (m_scenario == "InH-OfficeMixed" || m_scenario == "InH-OfficeOpen")
752  {
753  NS_ASSERT_MSG(!o2i, "The indoor scenario does out support outdoor to indoor");
754  if (los)
755  {
756  table3gpp->m_numOfCluster = 15;
757  table3gpp->m_raysPerCluster = 20;
758  table3gpp->m_uLgDS = -0.01 * log10(1 + fcGHz) - 7.692;
759  table3gpp->m_sigLgDS = 0.18;
760  table3gpp->m_uLgASD = 1.60;
761  table3gpp->m_sigLgASD = 0.18;
762  table3gpp->m_uLgASA = -0.19 * log10(1 + fcGHz) + 1.781;
763  table3gpp->m_sigLgASA = 0.12 * log10(1 + fcGHz) + 0.119;
764  table3gpp->m_uLgZSA = -0.26 * log10(1 + fcGHz) + 1.44;
765  table3gpp->m_sigLgZSA = -0.04 * log10(1 + fcGHz) + 0.264;
766  table3gpp->m_uLgZSD = -1.43 * log10(1 + fcGHz) + 2.228;
767  table3gpp->m_sigLgZSD = 0.13 * log10(1 + fcGHz) + 0.30;
768  table3gpp->m_offsetZOD = 0;
769  table3gpp->m_cDS = 3.91e-9;
770  table3gpp->m_cASD = 5;
771  table3gpp->m_cASA = 8;
772  table3gpp->m_cZSA = 9;
773  table3gpp->m_uK = 7;
774  table3gpp->m_sigK = 4;
775  table3gpp->m_rTau = 3.6;
776  table3gpp->m_uXpr = 11;
777  table3gpp->m_sigXpr = 4;
778  table3gpp->m_perClusterShadowingStd = 6;
779 
780  for (uint8_t row = 0; row < 7; row++)
781  {
782  for (uint8_t column = 0; column < 7; column++)
783  {
784  table3gpp->m_sqrtC[row][column] = sqrtC_office_LOS[row][column];
785  }
786  }
787  }
788  else
789  {
790  table3gpp->m_numOfCluster = 19;
791  table3gpp->m_raysPerCluster = 20;
792  table3gpp->m_uLgDS = -0.28 * log10(1 + fcGHz) - 7.173;
793  table3gpp->m_sigLgDS = 0.1 * log10(1 + fcGHz) + 0.055;
794  table3gpp->m_uLgASD = 1.62;
795  table3gpp->m_sigLgASD = 0.25;
796  table3gpp->m_uLgASA = -0.11 * log10(1 + fcGHz) + 1.863;
797  table3gpp->m_sigLgASA = 0.12 * log10(1 + fcGHz) + 0.059;
798  table3gpp->m_uLgZSA = -0.15 * log10(1 + fcGHz) + 1.387;
799  table3gpp->m_sigLgZSA = -0.09 * log10(1 + fcGHz) + 0.746;
800  table3gpp->m_uLgZSD = 1.08;
801  table3gpp->m_sigLgZSD = 0.36;
802  table3gpp->m_offsetZOD = 0;
803  table3gpp->m_cDS = 3.91e-9;
804  table3gpp->m_cASD = 5;
805  table3gpp->m_cASA = 11;
806  table3gpp->m_cZSA = 9;
807  table3gpp->m_uK = 0;
808  table3gpp->m_sigK = 0;
809  table3gpp->m_rTau = 3;
810  table3gpp->m_uXpr = 10;
811  table3gpp->m_sigXpr = 4;
812  table3gpp->m_perClusterShadowingStd = 3;
813 
814  for (uint8_t row = 0; row < 6; row++)
815  {
816  for (uint8_t column = 0; column < 6; column++)
817  {
818  table3gpp->m_sqrtC[row][column] = sqrtC_office_NLOS[row][column];
819  }
820  }
821  }
822  }
823  else if (m_scenario == "V2V-Urban")
824  {
825  if (channelCondition->IsLos())
826  {
827  // 3GPP mentioned that 3.91 ns should be used when the Cluster DS (cDS)
828  // entry is N/A.
829  table3gpp->m_numOfCluster = 12;
830  table3gpp->m_raysPerCluster = 20;
831  table3gpp->m_uLgDS = -0.2 * log10(1 + fcGHz) - 7.5;
832  table3gpp->m_sigLgDS = 0.1;
833  table3gpp->m_uLgASD = -0.1 * log10(1 + fcGHz) + 1.6;
834  table3gpp->m_sigLgASD = 0.1;
835  table3gpp->m_uLgASA = -0.1 * log10(1 + fcGHz) + 1.6;
836  table3gpp->m_sigLgASA = 0.1;
837  table3gpp->m_uLgZSA = -0.1 * log10(1 + fcGHz) + 0.73;
838  table3gpp->m_sigLgZSA = -0.04 * log10(1 + fcGHz) + 0.34;
839  table3gpp->m_uLgZSD = -0.1 * log10(1 + fcGHz) + 0.73;
840  table3gpp->m_sigLgZSD = -0.04 * log10(1 + fcGHz) + 0.34;
841  table3gpp->m_offsetZOD = 0;
842  table3gpp->m_cDS = 5;
843  table3gpp->m_cASD = 17;
844  table3gpp->m_cASA = 17;
845  table3gpp->m_cZSA = 7;
846  table3gpp->m_uK = 3.48;
847  table3gpp->m_sigK = 2;
848  table3gpp->m_rTau = 3;
849  table3gpp->m_uXpr = 9;
850  table3gpp->m_sigXpr = 3;
851  table3gpp->m_perClusterShadowingStd = 4;
852 
853  for (uint8_t row = 0; row < 7; row++)
854  {
855  for (uint8_t column = 0; column < 7; column++)
856  {
857  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_LOS[row][column];
858  }
859  }
860  }
861  else if (channelCondition->IsNlos())
862  {
863  table3gpp->m_numOfCluster = 19;
864  table3gpp->m_raysPerCluster = 20;
865  table3gpp->m_uLgDS = -0.3 * log10(1 + fcGHz) - 7;
866  table3gpp->m_sigLgDS = 0.28;
867  table3gpp->m_uLgASD = -0.08 * log10(1 + fcGHz) + 1.81;
868  table3gpp->m_sigLgASD = 0.05 * log10(1 + fcGHz) + 0.3;
869  table3gpp->m_uLgASA = -0.08 * log10(1 + fcGHz) + 1.81;
870  table3gpp->m_sigLgASA = 0.05 * log10(1 + fcGHz) + 0.3;
871  table3gpp->m_uLgZSA = -0.04 * log10(1 + fcGHz) + 0.92;
872  table3gpp->m_sigLgZSA = -0.07 * log10(1 + fcGHz) + 0.41;
873  table3gpp->m_uLgZSD = -0.04 * log10(1 + fcGHz) + 0.92;
874  table3gpp->m_sigLgZSD = -0.07 * log10(1 + fcGHz) + 0.41;
875  table3gpp->m_offsetZOD = 0;
876  table3gpp->m_cDS = 11;
877  table3gpp->m_cASD = 22;
878  table3gpp->m_cASA = 22;
879  table3gpp->m_cZSA = 7;
880  table3gpp->m_uK = 0; // N/A
881  table3gpp->m_sigK = 0; // N/A
882  table3gpp->m_rTau = 2.1;
883  table3gpp->m_uXpr = 8;
884  table3gpp->m_sigXpr = 3;
885  table3gpp->m_perClusterShadowingStd = 4;
886 
887  for (uint8_t row = 0; row < 6; row++)
888  {
889  for (uint8_t column = 0; column < 6; column++)
890  {
891  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_NLOS[row][column];
892  }
893  }
894  }
895  else if (channelCondition->IsNlosv())
896  {
897  table3gpp->m_numOfCluster = 19;
898  table3gpp->m_raysPerCluster = 20;
899  table3gpp->m_uLgDS = -0.4 * log10(1 + fcGHz) - 7;
900  table3gpp->m_sigLgDS = 0.1;
901  table3gpp->m_uLgASD = -0.1 * log10(1 + fcGHz) + 1.7;
902  table3gpp->m_sigLgASD = 0.1;
903  table3gpp->m_uLgASA = -0.1 * log10(1 + fcGHz) + 1.7;
904  table3gpp->m_sigLgASA = 0.1;
905  table3gpp->m_uLgZSA = -0.04 * log10(1 + fcGHz) + 0.92;
906  table3gpp->m_sigLgZSA = -0.07 * log10(1 + fcGHz) + 0.41;
907  table3gpp->m_uLgZSD = -0.04 * log10(1 + fcGHz) + 0.92;
908  table3gpp->m_sigLgZSD = -0.07 * log10(1 + fcGHz) + 0.41;
909  table3gpp->m_offsetZOD = 0;
910  table3gpp->m_cDS = 11;
911  table3gpp->m_cASD = 22;
912  table3gpp->m_cASA = 22;
913  table3gpp->m_cZSA = 7;
914  table3gpp->m_uK = 0;
915  table3gpp->m_sigK = 4.5;
916  table3gpp->m_rTau = 2.1;
917  table3gpp->m_uXpr = 8;
918  table3gpp->m_sigXpr = 3;
919  table3gpp->m_perClusterShadowingStd = 4;
920 
921  for (uint8_t row = 0; row < 6; row++)
922  {
923  for (uint8_t column = 0; column < 6; column++)
924  {
925  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_LOS[row][column];
926  }
927  }
928  }
929  else
930  {
931  NS_FATAL_ERROR("Unknown channel condition");
932  }
933  }
934  else if (m_scenario == "V2V-Highway")
935  {
936  if (channelCondition->IsLos())
937  {
938  table3gpp->m_numOfCluster = 12;
939  table3gpp->m_raysPerCluster = 20;
940  table3gpp->m_uLgDS = -8.3;
941  table3gpp->m_sigLgDS = 0.2;
942  table3gpp->m_uLgASD = 1.4;
943  table3gpp->m_sigLgASD = 0.1;
944  table3gpp->m_uLgASA = 1.4;
945  table3gpp->m_sigLgASA = 0.1;
946  table3gpp->m_uLgZSA = -0.1 * log10(1 + fcGHz) + 0.73;
947  table3gpp->m_sigLgZSA = -0.04 * log10(1 + fcGHz) + 0.34;
948  table3gpp->m_uLgZSD = -0.1 * log10(1 + fcGHz) + 0.73;
949  table3gpp->m_sigLgZSD = -0.04 * log10(1 + fcGHz) + 0.34;
950  table3gpp->m_offsetZOD = 0;
951  table3gpp->m_cDS = 5;
952  table3gpp->m_cASD = 17;
953  table3gpp->m_cASA = 17;
954  table3gpp->m_cZSA = 7;
955  table3gpp->m_uK = 9;
956  table3gpp->m_sigK = 3.5;
957  table3gpp->m_rTau = 3;
958  table3gpp->m_uXpr = 9;
959  table3gpp->m_sigXpr = 3;
960  table3gpp->m_perClusterShadowingStd = 4;
961 
962  for (uint8_t row = 0; row < 7; row++)
963  {
964  for (uint8_t column = 0; column < 7; column++)
965  {
966  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_LOS[row][column];
967  }
968  }
969  }
970  else if (channelCondition->IsNlosv())
971  {
972  table3gpp->m_numOfCluster = 19;
973  table3gpp->m_raysPerCluster = 20;
974  table3gpp->m_uLgDS = -8.3;
975  table3gpp->m_sigLgDS = 0.3;
976  table3gpp->m_uLgASD = 1.5;
977  table3gpp->m_sigLgASD = 0.1;
978  table3gpp->m_uLgASA = 1.5;
979  table3gpp->m_sigLgASA = 0.1;
980  table3gpp->m_uLgZSA = -0.04 * log10(1 + fcGHz) + 0.92;
981  table3gpp->m_sigLgZSA = -0.07 * log10(1 + fcGHz) + 0.41;
982  table3gpp->m_uLgZSD = -0.04 * log10(1 + fcGHz) + 0.92;
983  table3gpp->m_sigLgZSD = -0.07 * log10(1 + fcGHz) + 0.41;
984  table3gpp->m_offsetZOD = 0;
985  table3gpp->m_cDS = 11;
986  table3gpp->m_cASD = 22;
987  table3gpp->m_cASA = 22;
988  table3gpp->m_cZSA = 7;
989  table3gpp->m_uK = 0;
990  table3gpp->m_sigK = 4.5;
991  table3gpp->m_rTau = 2.1;
992  table3gpp->m_uXpr = 8.0;
993  table3gpp->m_sigXpr = 3;
994  table3gpp->m_perClusterShadowingStd = 4;
995 
996  for (uint8_t row = 0; row < 6; row++)
997  {
998  for (uint8_t column = 0; column < 6; column++)
999  {
1000  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_LOS[row][column];
1001  }
1002  }
1003  }
1004  else if (channelCondition->IsNlos())
1005  {
1006  NS_LOG_WARN(
1007  "The fast fading parameters for the NLOS condition in the Highway scenario are not "
1008  "defined in TR 37.885, use the ones defined in TDoc R1-1803671 instead");
1009 
1010  table3gpp->m_numOfCluster = 19;
1011  table3gpp->m_raysPerCluster = 20;
1012  table3gpp->m_uLgDS = -0.3 * log10(1 + fcGHz) - 7;
1013  table3gpp->m_sigLgDS = 0.28;
1014  table3gpp->m_uLgASD = -0.08 * log10(1 + fcGHz) + 1.81;
1015  table3gpp->m_sigLgASD = 0.05 * log10(1 + fcGHz) + 0.3;
1016  table3gpp->m_uLgASA = -0.08 * log10(1 + fcGHz) + 1.81;
1017  table3gpp->m_sigLgASA = 0.05 * log10(1 + fcGHz) + 0.3;
1018  table3gpp->m_uLgZSA = -0.04 * log10(1 + fcGHz) + 0.92;
1019  table3gpp->m_sigLgZSA = -0.07 * log10(1 + fcGHz) + 0.41;
1020  table3gpp->m_uLgZSD = -0.04 * log10(1 + fcGHz) + 0.92;
1021  table3gpp->m_sigLgZSD = -0.07 * log10(1 + fcGHz) + 0.41;
1022  table3gpp->m_offsetZOD = 0;
1023  table3gpp->m_cDS = 11;
1024  table3gpp->m_cASD = 22;
1025  table3gpp->m_cASA = 22;
1026  table3gpp->m_cZSA = 7;
1027  table3gpp->m_uK = 0; // N/A
1028  table3gpp->m_sigK = 0; // N/A
1029  table3gpp->m_rTau = 2.1;
1030  table3gpp->m_uXpr = 8;
1031  table3gpp->m_sigXpr = 3;
1032  table3gpp->m_perClusterShadowingStd = 4;
1033 
1034  for (uint8_t row = 0; row < 6; row++)
1035  {
1036  for (uint8_t column = 0; column < 6; column++)
1037  {
1038  table3gpp->m_sqrtC[row][column] = sqrtC_UMi_NLOS[row][column];
1039  }
1040  }
1041  }
1042  else
1043  {
1044  NS_FATAL_ERROR("Unknown channel condition");
1045  }
1046  }
1047  else
1048  {
1049  NS_FATAL_ERROR("unknown scenarios");
1050  }
1051 
1052  return table3gpp;
1053 }
1054 
1055 bool
1057  Ptr<const ChannelCondition> channelCondition) const
1058 {
1059  NS_LOG_FUNCTION(this);
1060 
1061  bool update = false;
1062 
1063  // if the channel condition is different the channel has to be updated
1064  if (!channelCondition->IsEqual(channelParams->m_losCondition, channelParams->m_o2iCondition))
1065  {
1066  NS_LOG_DEBUG("Update the channel condition");
1067  update = true;
1068  }
1069 
1070  // if the coherence time is over the channel has to be updated
1071  if (!m_updatePeriod.IsZero() &&
1072  Simulator::Now() - channelParams->m_generatedTime > m_updatePeriod)
1073  {
1074  NS_LOG_DEBUG("Generation time " << channelParams->m_generatedTime.As(Time::NS) << " now "
1075  << Now().As(Time::NS));
1076  update = true;
1077  }
1078 
1079  return update;
1080 }
1081 
1082 bool
1084  Ptr<const ChannelMatrix> channelMatrix)
1085 {
1086  if (channelParams->m_generatedTime > channelMatrix->m_generatedTime)
1087  {
1088  return true;
1089  }
1090  else
1091  {
1092  return false;
1093  }
1094 }
1095 
1099  Ptr<const PhasedArrayModel> aAntenna,
1100  Ptr<const PhasedArrayModel> bAntenna)
1101 {
1102  NS_LOG_FUNCTION(this);
1103 
1104  // Compute the channel params key. The key is reciprocal, i.e., key (a, b) = key (b, a)
1105  uint64_t channelParamsKey =
1106  GetKey(aMob->GetObject<Node>()->GetId(), bMob->GetObject<Node>()->GetId());
1107  // Compute the channel matrix key. The key is reciprocal, i.e., key (a, b) = key (b, a)
1108  uint64_t channelMatrixKey = GetKey(aAntenna->GetId(), bAntenna->GetId());
1109 
1110  // retrieve the channel condition
1111  Ptr<const ChannelCondition> condition =
1112  m_channelConditionModel->GetChannelCondition(aMob, bMob);
1113 
1114  // Check if the channel is present in the map and return it, otherwise
1115  // generate a new channel
1116  bool updateParams = false;
1117  bool updateMatrix = false;
1118  bool notFoundParams = false;
1119  bool notFoundMatrix = false;
1120  Ptr<ChannelMatrix> channelMatrix;
1121  Ptr<ThreeGppChannelParams> channelParams;
1122 
1123  if (m_channelParamsMap.find(channelParamsKey) != m_channelParamsMap.end())
1124  {
1125  channelParams = m_channelParamsMap[channelParamsKey];
1126  // check if it has to be updated
1127  updateParams = ChannelParamsNeedsUpdate(channelParams, condition);
1128  }
1129  else
1130  {
1131  NS_LOG_DEBUG("channel params not found");
1132  notFoundParams = true;
1133  }
1134 
1135  double x = aMob->GetPosition().x - bMob->GetPosition().x;
1136  double y = aMob->GetPosition().y - bMob->GetPosition().y;
1137  double distance2D = sqrt(x * x + y * y);
1138 
1139  // NOTE we assume hUT = min (height(a), height(b)) and
1140  // hBS = max (height (a), height (b))
1141  double hUt = std::min(aMob->GetPosition().z, bMob->GetPosition().z);
1142  double hBs = std::max(aMob->GetPosition().z, bMob->GetPosition().z);
1143 
1144  // get the 3GPP parameters
1145  Ptr<const ParamsTable> table3gpp = GetThreeGppTable(condition, hBs, hUt, distance2D);
1146 
1147  if (notFoundParams || updateParams)
1148  {
1149  // Step 4: Generate large scale parameters. All LSPS are uncorrelated.
1150  // Step 5: Generate Delays.
1151  // Step 6: Generate cluster powers.
1152  // Step 7: Generate arrival and departure angles for both azimuth and elevation.
1153  // Step 8: Coupling of rays within a cluster for both azimuth and elevation
1154  // shuffle all the arrays to perform random coupling
1155  // Step 9: Generate the cross polarization power ratios
1156  // Step 10: Draw initial phases
1157  channelParams = GenerateChannelParameters(condition, table3gpp, aMob, bMob);
1158  // store or replace the channel parameters
1159  m_channelParamsMap[channelParamsKey] = channelParams;
1160  }
1161 
1162  if (m_channelMatrixMap.find(channelMatrixKey) != m_channelMatrixMap.end())
1163  {
1164  // channel matrix present in the map
1165  NS_LOG_DEBUG("channel matrix present in the map");
1166  channelMatrix = m_channelMatrixMap[channelMatrixKey];
1167  updateMatrix = ChannelMatrixNeedsUpdate(channelParams, channelMatrix);
1168  }
1169  else
1170  {
1171  NS_LOG_DEBUG("channel matrix not found");
1172  notFoundMatrix = true;
1173  }
1174 
1175  // If the channel is not present in the map or if it has to be updated
1176  // generate a new realization
1177  if (notFoundMatrix || updateMatrix)
1178  {
1179  // channel matrix not found or has to be updated, generate a new one
1180  channelMatrix = GetNewChannel(channelParams, table3gpp, aMob, bMob, aAntenna, bAntenna);
1181  channelMatrix->m_antennaPair =
1182  std::make_pair(aAntenna->GetId(),
1183  bAntenna->GetId()); // save antenna pair, with the exact order of s and u
1184  // antennas at the moment of the channel generation
1185 
1186  // store or replace the channel matrix in the channel map
1187  m_channelMatrixMap[channelMatrixKey] = channelMatrix;
1188  }
1189 
1190  return channelMatrix;
1191 }
1192 
1195 {
1196  NS_LOG_FUNCTION(this);
1197 
1198  // Compute the channel key. The key is reciprocal, i.e., key (a, b) = key (b, a)
1199  uint64_t channelParamsKey =
1200  GetKey(aMob->GetObject<Node>()->GetId(), bMob->GetObject<Node>()->GetId());
1201 
1202  if (m_channelParamsMap.find(channelParamsKey) != m_channelParamsMap.end())
1203  {
1204  return m_channelParamsMap.find(channelParamsKey)->second;
1205  }
1206  else
1207  {
1208  NS_LOG_WARN("Channel params map not found. Returning a nullptr.");
1209  return nullptr;
1210  }
1211 }
1212 
1215  const Ptr<const ParamsTable> table3gpp,
1216  const Ptr<const MobilityModel> aMob,
1217  const Ptr<const MobilityModel> bMob) const
1218 {
1219  NS_LOG_FUNCTION(this);
1220  // create a channel matrix instance
1221  Ptr<ThreeGppChannelParams> channelParams = Create<ThreeGppChannelParams>();
1222  channelParams->m_generatedTime = Simulator::Now();
1223  channelParams->m_nodeIds =
1224  std::make_pair(aMob->GetObject<Node>()->GetId(), bMob->GetObject<Node>()->GetId());
1225  channelParams->m_losCondition = channelCondition->GetLosCondition();
1226  channelParams->m_o2iCondition = channelCondition->GetO2iCondition();
1227 
1228  // Step 4: Generate large scale parameters. All LSPS are uncorrelated.
1229  DoubleVector LSPsIndep;
1230  DoubleVector LSPs;
1231  uint8_t paramNum = 6;
1232  if (channelParams->m_losCondition == ChannelCondition::LOS)
1233  {
1234  paramNum = 7;
1235  }
1236 
1237  // Generate paramNum independent LSPs.
1238  for (uint8_t iter = 0; iter < paramNum; iter++)
1239  {
1240  LSPsIndep.push_back(m_normalRv->GetValue());
1241  }
1242  for (uint8_t row = 0; row < paramNum; row++)
1243  {
1244  double temp = 0;
1245  for (uint8_t column = 0; column < paramNum; column++)
1246  {
1247  temp += table3gpp->m_sqrtC[row][column] * LSPsIndep[column];
1248  }
1249  LSPs.push_back(temp);
1250  }
1251 
1252  // NOTE the shadowing is generated in the propagation loss model
1253  double DS;
1254  double ASD;
1255  double ASA;
1256  double ZSA;
1257  double ZSD;
1258  double kFactor = 0;
1259  if (channelParams->m_losCondition == ChannelCondition::LOS)
1260  {
1261  kFactor = LSPs[1] * table3gpp->m_sigK + table3gpp->m_uK;
1262  DS = pow(10, LSPs[2] * table3gpp->m_sigLgDS + table3gpp->m_uLgDS);
1263  ASD = pow(10, LSPs[3] * table3gpp->m_sigLgASD + table3gpp->m_uLgASD);
1264  ASA = pow(10, LSPs[4] * table3gpp->m_sigLgASA + table3gpp->m_uLgASA);
1265  ZSD = pow(10, LSPs[5] * table3gpp->m_sigLgZSD + table3gpp->m_uLgZSD);
1266  ZSA = pow(10, LSPs[6] * table3gpp->m_sigLgZSA + table3gpp->m_uLgZSA);
1267  }
1268  else
1269  {
1270  DS = pow(10, LSPs[1] * table3gpp->m_sigLgDS + table3gpp->m_uLgDS);
1271  ASD = pow(10, LSPs[2] * table3gpp->m_sigLgASD + table3gpp->m_uLgASD);
1272  ASA = pow(10, LSPs[3] * table3gpp->m_sigLgASA + table3gpp->m_uLgASA);
1273  ZSD = pow(10, LSPs[4] * table3gpp->m_sigLgZSD + table3gpp->m_uLgZSD);
1274  ZSA = pow(10, LSPs[5] * table3gpp->m_sigLgZSA + table3gpp->m_uLgZSA);
1275  }
1276  ASD = std::min(ASD, 104.0);
1277  ASA = std::min(ASA, 104.0);
1278  ZSD = std::min(ZSD, 52.0);
1279  ZSA = std::min(ZSA, 52.0);
1280 
1281  // save DS and K_factor parameters in the structure
1282  channelParams->m_DS = DS;
1283  channelParams->m_K_factor = kFactor;
1284 
1285  NS_LOG_INFO("K-factor=" << kFactor << ", DS=" << DS << ", ASD=" << ASD << ", ASA=" << ASA
1286  << ", ZSD=" << ZSD << ", ZSA=" << ZSA);
1287 
1288  // Step 5: Generate Delays.
1289  DoubleVector clusterDelay;
1290  double minTau = 100.0;
1291  for (uint8_t cIndex = 0; cIndex < table3gpp->m_numOfCluster; cIndex++)
1292  {
1293  double tau = -1 * table3gpp->m_rTau * DS * log(m_uniformRv->GetValue(0, 1)); //(7.5-1)
1294  if (minTau > tau)
1295  {
1296  minTau = tau;
1297  }
1298  clusterDelay.push_back(tau);
1299  }
1300 
1301  for (uint8_t cIndex = 0; cIndex < table3gpp->m_numOfCluster; cIndex++)
1302  {
1303  clusterDelay[cIndex] -= minTau;
1304  }
1305  std::sort(clusterDelay.begin(), clusterDelay.end()); //(7.5-2)
1306 
1307  /* since the scaled Los delays are not to be used in cluster power generation,
1308  * we will generate cluster power first and resume to compute Los cluster delay later.*/
1309 
1310  // Step 6: Generate cluster powers.
1311  DoubleVector clusterPower;
1312  double powerSum = 0;
1313  for (uint8_t cIndex = 0; cIndex < table3gpp->m_numOfCluster; cIndex++)
1314  {
1315  double power =
1316  exp(-1 * clusterDelay[cIndex] * (table3gpp->m_rTau - 1) / table3gpp->m_rTau / DS) *
1317  pow(10,
1318  -1 * m_normalRv->GetValue() * table3gpp->m_perClusterShadowingStd / 10.0); //(7.5-5)
1319  powerSum += power;
1320  clusterPower.push_back(power);
1321  }
1322  channelParams->m_clusterPower = clusterPower;
1323 
1324  double powerMax = 0;
1325 
1326  for (uint8_t cIndex = 0; cIndex < table3gpp->m_numOfCluster; cIndex++)
1327  {
1328  channelParams->m_clusterPower[cIndex] =
1329  channelParams->m_clusterPower[cIndex] / powerSum; //(7.5-6)
1330  }
1331 
1332  DoubleVector clusterPowerForAngles; // this power is only for equation (7.5-9) and (7.5-14), not
1333  // for (7.5-22)
1334  if (channelParams->m_losCondition == ChannelCondition::LOS)
1335  {
1336  double kLinear = pow(10, kFactor / 10.0);
1337 
1338  for (uint8_t cIndex = 0; cIndex < table3gpp->m_numOfCluster; cIndex++)
1339  {
1340  if (cIndex == 0)
1341  {
1342  clusterPowerForAngles.push_back(channelParams->m_clusterPower[cIndex] /
1343  (1 + kLinear) +
1344  kLinear / (1 + kLinear)); //(7.5-8)
1345  }
1346  else
1347  {
1348  clusterPowerForAngles.push_back(channelParams->m_clusterPower[cIndex] /
1349  (1 + kLinear)); //(7.5-8)
1350  }
1351  if (powerMax < clusterPowerForAngles[cIndex])
1352  {
1353  powerMax = clusterPowerForAngles[cIndex];
1354  }
1355  }
1356  }
1357  else
1358  {
1359  for (uint8_t cIndex = 0; cIndex < table3gpp->m_numOfCluster; cIndex++)
1360  {
1361  clusterPowerForAngles.push_back(channelParams->m_clusterPower[cIndex]); //(7.5-6)
1362  if (powerMax < clusterPowerForAngles[cIndex])
1363  {
1364  powerMax = clusterPowerForAngles[cIndex];
1365  }
1366  }
1367  }
1368 
1369  // remove clusters with less than -25 dB power compared to the maxim cluster power;
1370  // double thresh = pow(10, -2.5);
1371  double thresh = 0.0032;
1372  for (uint8_t cIndex = table3gpp->m_numOfCluster; cIndex > 0; cIndex--)
1373  {
1374  if (clusterPowerForAngles[cIndex - 1] < thresh * powerMax)
1375  {
1376  clusterPowerForAngles.erase(clusterPowerForAngles.begin() + cIndex - 1);
1377  channelParams->m_clusterPower.erase(channelParams->m_clusterPower.begin() + cIndex - 1);
1378  clusterDelay.erase(clusterDelay.begin() + cIndex - 1);
1379  }
1380  }
1381 
1382  NS_ASSERT(channelParams->m_clusterPower.size() < UINT8_MAX);
1383  channelParams->m_reducedClusterNumber = channelParams->m_clusterPower.size();
1384  // Resume step 5 to compute the delay for LoS condition.
1385  if (channelParams->m_losCondition == ChannelCondition::LOS)
1386  {
1387  double cTau =
1388  0.7705 - 0.0433 * kFactor + 2e-4 * pow(kFactor, 2) + 17e-6 * pow(kFactor, 3); //(7.5-3)
1389  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1390  {
1391  clusterDelay[cIndex] = clusterDelay[cIndex] / cTau; //(7.5-4)
1392  }
1393  }
1394 
1395  // Step 7: Generate arrival and departure angles for both azimuth and elevation.
1396 
1397  double cNlos;
1398  // According to table 7.5-6, only cluster number equals to 8, 10, 11, 12, 19 and 20 is valid.
1399  // Not sure why the other cases are in Table 7.5-2.
1400  switch (table3gpp->m_numOfCluster) // Table 7.5-2
1401  {
1402  case 4:
1403  cNlos = 0.779;
1404  break;
1405  case 5:
1406  cNlos = 0.860;
1407  break;
1408  case 8:
1409  cNlos = 1.018;
1410  break;
1411  case 10:
1412  cNlos = 1.090;
1413  break;
1414  case 11:
1415  cNlos = 1.123;
1416  break;
1417  case 12:
1418  cNlos = 1.146;
1419  break;
1420  case 14:
1421  cNlos = 1.190;
1422  break;
1423  case 15:
1424  cNlos = 1.221;
1425  break;
1426  case 16:
1427  cNlos = 1.226;
1428  break;
1429  case 19:
1430  cNlos = 1.273;
1431  break;
1432  case 20:
1433  cNlos = 1.289;
1434  break;
1435  default:
1436  NS_FATAL_ERROR("Invalid cluster number");
1437  }
1438 
1439  double cPhi = cNlos;
1440 
1441  if (channelParams->m_losCondition == ChannelCondition::LOS)
1442  {
1443  cPhi *= (1.1035 - 0.028 * kFactor - 2e-3 * pow(kFactor, 2) +
1444  1e-4 * pow(kFactor, 3)); //(7.5-10))
1445  }
1446 
1447  switch (table3gpp->m_numOfCluster) // Table 7.5-4
1448  {
1449  case 8:
1450  cNlos = 0.889;
1451  break;
1452  case 10:
1453  cNlos = 0.957;
1454  break;
1455  case 11:
1456  cNlos = 1.031;
1457  break;
1458  case 12:
1459  cNlos = 1.104;
1460  break;
1461  case 15:
1462  cNlos = 1.1088;
1463  break;
1464  case 19:
1465  cNlos = 1.184;
1466  break;
1467  case 20:
1468  cNlos = 1.178;
1469  break;
1470  default:
1471  NS_FATAL_ERROR("Invalid cluster number");
1472  }
1473 
1474  double cTheta = cNlos;
1475  if (channelCondition->IsLos())
1476  {
1477  cTheta *= (1.3086 + 0.0339 * kFactor - 0.0077 * pow(kFactor, 2) +
1478  2e-4 * pow(kFactor, 3)); //(7.5-15)
1479  }
1480 
1481  DoubleVector clusterAoa;
1482  DoubleVector clusterAod;
1483  DoubleVector clusterZoa;
1484  DoubleVector clusterZod;
1485  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1486  {
1487  double logCalc = -1 * log(clusterPowerForAngles[cIndex] / powerMax);
1488  double angle = 2 * sqrt(logCalc) / 1.4 / cPhi; //(7.5-9)
1489  clusterAoa.push_back(ASA * angle);
1490  clusterAod.push_back(ASD * angle);
1491  angle = logCalc / cTheta; //(7.5-14)
1492  clusterZoa.push_back(ZSA * angle);
1493  clusterZod.push_back(ZSD * angle);
1494  }
1495 
1496  Angles sAngle(bMob->GetPosition(), aMob->GetPosition());
1497  Angles uAngle(aMob->GetPosition(), bMob->GetPosition());
1498 
1499  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1500  {
1501  int Xn = 1;
1502  if (m_uniformRv->GetValue(0, 1) < 0.5)
1503  {
1504  Xn = -1;
1505  }
1506  clusterAoa[cIndex] = clusterAoa[cIndex] * Xn + (m_normalRv->GetValue() * ASA / 7.0) +
1507  RadiansToDegrees(uAngle.GetAzimuth()); //(7.5-11)
1508  clusterAod[cIndex] = clusterAod[cIndex] * Xn + (m_normalRv->GetValue() * ASD / 7.0) +
1509  RadiansToDegrees(sAngle.GetAzimuth());
1510  if (channelCondition->IsO2i())
1511  {
1512  clusterZoa[cIndex] =
1513  clusterZoa[cIndex] * Xn + (m_normalRv->GetValue() * ZSA / 7.0) + 90; //(7.5-16)
1514  }
1515  else
1516  {
1517  clusterZoa[cIndex] = clusterZoa[cIndex] * Xn + (m_normalRv->GetValue() * ZSA / 7.0) +
1518  RadiansToDegrees(uAngle.GetInclination()); //(7.5-16)
1519  }
1520  clusterZod[cIndex] = clusterZod[cIndex] * Xn + (m_normalRv->GetValue() * ZSD / 7.0) +
1521  RadiansToDegrees(sAngle.GetInclination()) +
1522  table3gpp->m_offsetZOD; //(7.5-19)
1523  }
1524 
1525  if (channelParams->m_losCondition == ChannelCondition::LOS)
1526  {
1527  // The 7.5-12 can be rewrite as Theta_n,ZOA = Theta_n,ZOA - (Theta_1,ZOA - Theta_LOS,ZOA) =
1528  // Theta_n,ZOA - diffZOA, Similar as AOD, ZSA and ZSD.
1529  double diffAoa = clusterAoa[0] - RadiansToDegrees(uAngle.GetAzimuth());
1530  double diffAod = clusterAod[0] - RadiansToDegrees(sAngle.GetAzimuth());
1531  double diffZsa = clusterZoa[0] - RadiansToDegrees(uAngle.GetInclination());
1532  double diffZsd = clusterZod[0] - RadiansToDegrees(sAngle.GetInclination());
1533 
1534  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1535  {
1536  clusterAoa[cIndex] -= diffAoa; //(7.5-12)
1537  clusterAod[cIndex] -= diffAod;
1538  clusterZoa[cIndex] -= diffZsa; //(7.5-17)
1539  clusterZod[cIndex] -= diffZsd;
1540  }
1541  }
1542 
1543  double sizeTemp = clusterZoa.size();
1544  for (uint8_t ind = 0; ind < 4; ind++)
1545  {
1546  DoubleVector angleDegree;
1547  switch (ind)
1548  {
1549  case 0:
1550  angleDegree = clusterAoa;
1551  break;
1552  case 1:
1553  angleDegree = clusterZoa;
1554  break;
1555  case 2:
1556  angleDegree = clusterAod;
1557  break;
1558  case 3:
1559  angleDegree = clusterZod;
1560  break;
1561  default:
1562  NS_FATAL_ERROR("Programming Error");
1563  }
1564  for (uint8_t nIndex = 0; nIndex < sizeTemp; nIndex++)
1565  {
1566  while (angleDegree[nIndex] > 360)
1567  {
1568  angleDegree[nIndex] -= 360;
1569  }
1570 
1571  while (angleDegree[nIndex] < 0)
1572  {
1573  angleDegree[nIndex] += 360;
1574  }
1575 
1576  if (ind == 1 || ind == 3)
1577  {
1578  if (angleDegree[nIndex] > 180)
1579  {
1580  angleDegree[nIndex] = 360 - angleDegree[nIndex];
1581  }
1582  }
1583  }
1584  switch (ind)
1585  {
1586  case 0:
1587  clusterAoa = angleDegree;
1588  break;
1589  case 1:
1590  clusterZoa = angleDegree;
1591  break;
1592  case 2:
1593  clusterAod = angleDegree;
1594  break;
1595  case 3:
1596  clusterZod = angleDegree;
1597  break;
1598  default:
1599  NS_FATAL_ERROR("Programming Error");
1600  }
1601  }
1602 
1603  DoubleVector attenuationDb;
1604  if (m_blockage)
1605  {
1606  attenuationDb = CalcAttenuationOfBlockage(channelParams, clusterAoa, clusterZoa);
1607  for (uint8_t cInd = 0; cInd < channelParams->m_reducedClusterNumber; cInd++)
1608  {
1609  channelParams->m_clusterPower[cInd] =
1610  channelParams->m_clusterPower[cInd] / pow(10, attenuationDb[cInd] / 10.0);
1611  }
1612  }
1613  else
1614  {
1615  attenuationDb.push_back(0);
1616  }
1617 
1618  // store attenuation
1619  channelParams->m_attenuation_dB = attenuationDb;
1620 
1621  // Step 8: Coupling of rays within a cluster for both azimuth and elevation
1622  // shuffle all the arrays to perform random coupling
1624  channelParams->m_reducedClusterNumber,
1625  DoubleVector(table3gpp->m_raysPerCluster,
1626  0)); // rayAoaRadian[n][m], where n is cluster index, m is ray index
1628  channelParams->m_reducedClusterNumber,
1629  DoubleVector(table3gpp->m_raysPerCluster,
1630  0)); // rayAodRadian[n][m], where n is cluster index, m is ray index
1632  channelParams->m_reducedClusterNumber,
1633  DoubleVector(table3gpp->m_raysPerCluster,
1634  0)); // rayZoaRadian[n][m], where n is cluster index, m is ray index
1636  channelParams->m_reducedClusterNumber,
1637  DoubleVector(table3gpp->m_raysPerCluster,
1638  0)); // rayZodRadian[n][m], where n is cluster index, m is ray index
1639 
1640  for (uint8_t nInd = 0; nInd < channelParams->m_reducedClusterNumber; nInd++)
1641  {
1642  for (uint8_t mInd = 0; mInd < table3gpp->m_raysPerCluster; mInd++)
1643  {
1644  double tempAoa = clusterAoa[nInd] + table3gpp->m_cASA * offSetAlpha[mInd]; //(7.5-13)
1645  double tempZoa = clusterZoa[nInd] + table3gpp->m_cZSA * offSetAlpha[mInd]; //(7.5-18)
1646  std::tie(rayAoaRadian[nInd][mInd], rayZoaRadian[nInd][mInd]) =
1647  WrapAngles(DegreesToRadians(tempAoa), DegreesToRadians(tempZoa));
1648 
1649  double tempAod = clusterAod[nInd] + table3gpp->m_cASD * offSetAlpha[mInd]; //(7.5-13)
1650  double tempZod = clusterZod[nInd] +
1651  0.375 * pow(10, table3gpp->m_uLgZSD) * offSetAlpha[mInd]; //(7.5-20)
1652  std::tie(rayAodRadian[nInd][mInd], rayZodRadian[nInd][mInd]) =
1653  WrapAngles(DegreesToRadians(tempAod), DegreesToRadians(tempZod));
1654  }
1655  }
1656 
1657  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1658  {
1659  Shuffle(&rayAodRadian[cIndex][0], &rayAodRadian[cIndex][table3gpp->m_raysPerCluster]);
1660  Shuffle(&rayAoaRadian[cIndex][0], &rayAoaRadian[cIndex][table3gpp->m_raysPerCluster]);
1661  Shuffle(&rayZodRadian[cIndex][0], &rayZodRadian[cIndex][table3gpp->m_raysPerCluster]);
1662  Shuffle(&rayZoaRadian[cIndex][0], &rayZoaRadian[cIndex][table3gpp->m_raysPerCluster]);
1663  }
1664 
1665  // store values
1666  channelParams->m_rayAodRadian = rayAodRadian;
1667  channelParams->m_rayAoaRadian = rayAoaRadian;
1668  channelParams->m_rayZodRadian = rayZodRadian;
1669  channelParams->m_rayZoaRadian = rayZoaRadian;
1670 
1671  // Step 9: Generate the cross polarization power ratios
1672  // Step 10: Draw initial phases
1673  Double2DVector crossPolarizationPowerRatios; // vector containing the cross polarization power
1674  // ratios, as defined by 7.5-21
1675  Double3DVector clusterPhase; // rayAoaRadian[n][m], where n is cluster index, m is ray index
1676  for (uint8_t nInd = 0; nInd < channelParams->m_reducedClusterNumber; nInd++)
1677  {
1678  DoubleVector temp; // used to store the XPR values
1680  temp2; // used to store the PHI values for all the possible combination of polarization
1681  for (uint8_t mInd = 0; mInd < table3gpp->m_raysPerCluster; mInd++)
1682  {
1683  double uXprLinear = pow(10, table3gpp->m_uXpr / 10.0); // convert to linear
1684  double sigXprLinear = pow(10, table3gpp->m_sigXpr / 10.0); // convert to linear
1685 
1686  temp.push_back(
1687  std::pow(10, (m_normalRv->GetValue() * sigXprLinear + uXprLinear) / 10.0));
1688  DoubleVector temp3; // used to store the PHI values
1689  for (uint8_t pInd = 0; pInd < 4; pInd++)
1690  {
1691  temp3.push_back(m_uniformRv->GetValue(-1 * M_PI, M_PI));
1692  }
1693  temp2.push_back(temp3);
1694  }
1695  crossPolarizationPowerRatios.push_back(temp);
1696  clusterPhase.push_back(temp2);
1697  }
1698  // store the cluster phase
1699  channelParams->m_clusterPhase = clusterPhase;
1700  channelParams->m_crossPolarizationPowerRatios = crossPolarizationPowerRatios;
1701 
1702  uint8_t cluster1st = 0;
1703  uint8_t cluster2nd = 0; // first and second strongest cluster;
1704  double maxPower = 0;
1705  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1706  {
1707  if (maxPower < channelParams->m_clusterPower[cIndex])
1708  {
1709  maxPower = channelParams->m_clusterPower[cIndex];
1710  cluster1st = cIndex;
1711  }
1712  }
1713  channelParams->m_cluster1st = cluster1st;
1714  maxPower = 0;
1715  for (uint8_t cIndex = 0; cIndex < channelParams->m_reducedClusterNumber; cIndex++)
1716  {
1717  if (maxPower < channelParams->m_clusterPower[cIndex] && cluster1st != cIndex)
1718  {
1719  maxPower = channelParams->m_clusterPower[cIndex];
1720  cluster2nd = cIndex;
1721  }
1722  }
1723  channelParams->m_cluster2nd = cluster2nd;
1724 
1725  NS_LOG_INFO("1st strongest cluster:" << +cluster1st
1726  << ", 2nd strongest cluster:" << +cluster2nd);
1727 
1728  // store the delays and the angles for the subclusters
1729  if (cluster1st == cluster2nd)
1730  {
1731  clusterDelay.push_back(clusterDelay[cluster1st] + 1.28 * table3gpp->m_cDS);
1732  clusterDelay.push_back(clusterDelay[cluster1st] + 2.56 * table3gpp->m_cDS);
1733 
1734  clusterAoa.push_back(clusterAoa[cluster1st]);
1735  clusterAoa.push_back(clusterAoa[cluster1st]);
1736 
1737  clusterZoa.push_back(clusterZoa[cluster1st]);
1738  clusterZoa.push_back(clusterZoa[cluster1st]);
1739 
1740  clusterAod.push_back(clusterAod[cluster1st]);
1741  clusterAod.push_back(clusterAod[cluster1st]);
1742 
1743  clusterZod.push_back(clusterZod[cluster1st]);
1744  clusterZod.push_back(clusterZod[cluster1st]);
1745  }
1746  else
1747  {
1748  double min;
1749  double max;
1750  if (cluster1st < cluster2nd)
1751  {
1752  min = cluster1st;
1753  max = cluster2nd;
1754  }
1755  else
1756  {
1757  min = cluster2nd;
1758  max = cluster1st;
1759  }
1760  clusterDelay.push_back(clusterDelay[min] + 1.28 * table3gpp->m_cDS);
1761  clusterDelay.push_back(clusterDelay[min] + 2.56 * table3gpp->m_cDS);
1762  clusterDelay.push_back(clusterDelay[max] + 1.28 * table3gpp->m_cDS);
1763  clusterDelay.push_back(clusterDelay[max] + 2.56 * table3gpp->m_cDS);
1764 
1765  clusterAoa.push_back(clusterAoa[min]);
1766  clusterAoa.push_back(clusterAoa[min]);
1767  clusterAoa.push_back(clusterAoa[max]);
1768  clusterAoa.push_back(clusterAoa[max]);
1769 
1770  clusterZoa.push_back(clusterZoa[min]);
1771  clusterZoa.push_back(clusterZoa[min]);
1772  clusterZoa.push_back(clusterZoa[max]);
1773  clusterZoa.push_back(clusterZoa[max]);
1774 
1775  clusterAod.push_back(clusterAod[min]);
1776  clusterAod.push_back(clusterAod[min]);
1777  clusterAod.push_back(clusterAod[max]);
1778  clusterAod.push_back(clusterAod[max]);
1779 
1780  clusterZod.push_back(clusterZod[min]);
1781  clusterZod.push_back(clusterZod[min]);
1782  clusterZod.push_back(clusterZod[max]);
1783  clusterZod.push_back(clusterZod[max]);
1784  }
1785 
1786  channelParams->m_delay = clusterDelay;
1787  channelParams->m_angle.clear();
1788  channelParams->m_angle.push_back(clusterAoa);
1789  channelParams->m_angle.push_back(clusterZoa);
1790  channelParams->m_angle.push_back(clusterAod);
1791  channelParams->m_angle.push_back(clusterZod);
1792 
1793  // Compute alpha and D as described in 3GPP TR 37.885 v15.3.0, Sec. 6.2.3
1794  // These terms account for an additional Doppler contribution due to the
1795  // presence of moving objects in the surrounding environment, such as in
1796  // vehicular scenarios.
1797  // This contribution is applied only to the delayed (reflected) paths and
1798  // must be properly configured by setting the value of
1799  // m_vScatt, which is defined as "maximum speed of the vehicle in the
1800  // layout".
1801  // By default, m_vScatt is set to 0, so there is no additional Doppler
1802  // contribution.
1803 
1804  DoubleVector dopplerTermAlpha;
1805  DoubleVector dopplerTermD;
1806 
1807  // 2 or 4 is added to account for additional subrays for the 1st and 2nd clusters, if there is
1808  // only one cluster then would be added 2 more subrays (see creation of Husn channel matrix)
1809  uint8_t updatedClusterNumber = (channelParams->m_reducedClusterNumber == 1)
1810  ? channelParams->m_reducedClusterNumber + 2
1811  : channelParams->m_reducedClusterNumber + 4;
1812 
1813  for (uint8_t cIndex = 0; cIndex < updatedClusterNumber; cIndex++)
1814  {
1815  double alpha = 0;
1816  double D = 0;
1817  if (cIndex != 0)
1818  {
1819  alpha = m_uniformRvDoppler->GetValue(-1, 1);
1821  }
1822  dopplerTermAlpha.push_back(alpha);
1823  dopplerTermD.push_back(D);
1824  }
1825  channelParams->m_alpha = dopplerTermAlpha;
1826  channelParams->m_D = dopplerTermD;
1827 
1828  return channelParams;
1829 }
1830 
1833  Ptr<const ParamsTable> table3gpp,
1834  const Ptr<const MobilityModel> sMob,
1835  const Ptr<const MobilityModel> uMob,
1836  Ptr<const PhasedArrayModel> sAntenna,
1837  Ptr<const PhasedArrayModel> uAntenna) const
1838 {
1839  NS_LOG_FUNCTION(this);
1840 
1841  NS_ASSERT_MSG(m_frequency > 0.0, "Set the operating frequency first!");
1842 
1843  // create a channel matrix instance
1844  Ptr<ChannelMatrix> channelMatrix = Create<ChannelMatrix>();
1845  channelMatrix->m_generatedTime = Simulator::Now();
1846  // save in which order is generated this matrix
1847  channelMatrix->m_nodeIds =
1848  std::make_pair(sMob->GetObject<Node>()->GetId(), uMob->GetObject<Node>()->GetId());
1849  // check if channelParams structure is generated in direction s-to-u or u-to-s
1850  bool isSameDirection = (channelParams->m_nodeIds == channelMatrix->m_nodeIds);
1851 
1856 
1857  // if channel params is generated in the same direction in which we
1858  // generate the channel matrix, angles and zenith od departure and arrival are ok,
1859  // just set them to corresponding variable that will be used for the generation
1860  // of channel matrix, otherwise we need to flip angles and zeniths of departure and arrival
1861  if (isSameDirection)
1862  {
1863  rayAodRadian = channelParams->m_rayAodRadian;
1864  rayAoaRadian = channelParams->m_rayAoaRadian;
1865  rayZodRadian = channelParams->m_rayZodRadian;
1866  rayZoaRadian = channelParams->m_rayZoaRadian;
1867  }
1868  else
1869  {
1870  rayAodRadian = channelParams->m_rayAoaRadian;
1871  rayAoaRadian = channelParams->m_rayAodRadian;
1872  rayZodRadian = channelParams->m_rayZoaRadian;
1873  rayZoaRadian = channelParams->m_rayZodRadian;
1874  }
1875 
1876  // Step 11: Generate channel coefficients for each cluster n and each receiver
1877  // and transmitter element pair u,s.
1878  // where n is cluster index, u and s are receive and transmit antenna element.
1879  size_t uSize = uAntenna->GetNumberOfElements();
1880  size_t sSize = sAntenna->GetNumberOfElements();
1881 
1882  // NOTE: Since each of the strongest 2 clusters are divided into 3 sub-clusters,
1883  // the total cluster will generally be numReducedCLuster + 4.
1884  // However, it might be that m_cluster1st = m_cluster2nd. In this case the
1885  // total number of clusters will be numReducedCLuster + 2.
1886  uint16_t numOverallCluster = (channelParams->m_cluster1st != channelParams->m_cluster2nd)
1887  ? channelParams->m_reducedClusterNumber + 4
1888  : channelParams->m_reducedClusterNumber + 2;
1889  Complex3DVector hUsn(uSize, sSize, numOverallCluster); // channel coefficient hUsn (u, s, n);
1890  NS_ASSERT(channelParams->m_reducedClusterNumber <= channelParams->m_clusterPhase.size());
1891  NS_ASSERT(channelParams->m_reducedClusterNumber <= channelParams->m_clusterPower.size());
1892  NS_ASSERT(channelParams->m_reducedClusterNumber <=
1893  channelParams->m_crossPolarizationPowerRatios.size());
1894  NS_ASSERT(channelParams->m_reducedClusterNumber <= rayZoaRadian.size());
1895  NS_ASSERT(channelParams->m_reducedClusterNumber <= rayZodRadian.size());
1896  NS_ASSERT(channelParams->m_reducedClusterNumber <= rayAoaRadian.size());
1897  NS_ASSERT(channelParams->m_reducedClusterNumber <= rayAodRadian.size());
1898  NS_ASSERT(table3gpp->m_raysPerCluster <= channelParams->m_clusterPhase[0].size());
1899  NS_ASSERT(table3gpp->m_raysPerCluster <=
1900  channelParams->m_crossPolarizationPowerRatios[0].size());
1901  NS_ASSERT(table3gpp->m_raysPerCluster <= rayZoaRadian[0].size());
1902  NS_ASSERT(table3gpp->m_raysPerCluster <= rayZodRadian[0].size());
1903  NS_ASSERT(table3gpp->m_raysPerCluster <= rayAoaRadian[0].size());
1904  NS_ASSERT(table3gpp->m_raysPerCluster <= rayAodRadian[0].size());
1905 
1906  double x = sMob->GetPosition().x - uMob->GetPosition().x;
1907  double y = sMob->GetPosition().y - uMob->GetPosition().y;
1908  double distance2D = sqrt(x * x + y * y);
1909  // NOTE we assume hUT = min (height(a), height(b)) and
1910  // hBS = max (height (a), height (b))
1911  double hUt = std::min(sMob->GetPosition().z, uMob->GetPosition().z);
1912  double hBs = std::max(sMob->GetPosition().z, uMob->GetPosition().z);
1913  // compute the 3D distance using eq. 7.4-1
1914  double distance3D = std::sqrt(distance2D * distance2D + (hBs - hUt) * (hBs - hUt));
1915 
1916  Angles sAngle(uMob->GetPosition(), sMob->GetPosition());
1917  Angles uAngle(sMob->GetPosition(), uMob->GetPosition());
1918 
1919  Complex2DVector raysPreComp(channelParams->m_reducedClusterNumber,
1920  table3gpp->m_raysPerCluster); // stores part of the ray expression,
1921  // cached as independent from the u- and s-indexes
1922  Double2DVector sinCosA; // cached multiplications of sin and cos of the ZoA and AoA angles
1923  Double2DVector sinSinA; // cached multiplications of sines of the ZoA and AoA angles
1924  Double2DVector cosZoA; // cached cos of the ZoA angle
1925  Double2DVector sinCosD; // cached multiplications of sin and cos of the ZoD and AoD angles
1926  Double2DVector sinSinD; // cached multiplications of the cosines of the ZoA and AoA angles
1927  Double2DVector cosZoD; // cached cos of the ZoD angle
1928 
1929  // resize to appropriate dimensions
1930  sinCosA.resize(channelParams->m_reducedClusterNumber);
1931  sinSinA.resize(channelParams->m_reducedClusterNumber);
1932  cosZoA.resize(channelParams->m_reducedClusterNumber);
1933  sinCosD.resize(channelParams->m_reducedClusterNumber);
1934  sinSinD.resize(channelParams->m_reducedClusterNumber);
1935  cosZoD.resize(channelParams->m_reducedClusterNumber);
1936  for (uint8_t nIndex = 0; nIndex < channelParams->m_reducedClusterNumber; nIndex++)
1937  {
1938  sinCosA[nIndex].resize(table3gpp->m_raysPerCluster);
1939  sinSinA[nIndex].resize(table3gpp->m_raysPerCluster);
1940  cosZoA[nIndex].resize(table3gpp->m_raysPerCluster);
1941  sinCosD[nIndex].resize(table3gpp->m_raysPerCluster);
1942  sinSinD[nIndex].resize(table3gpp->m_raysPerCluster);
1943  cosZoD[nIndex].resize(table3gpp->m_raysPerCluster);
1944  }
1945  // pre-compute the terms which are independent from uIndex and sIndex
1946  for (uint8_t nIndex = 0; nIndex < channelParams->m_reducedClusterNumber; nIndex++)
1947  {
1948  for (uint8_t mIndex = 0; mIndex < table3gpp->m_raysPerCluster; mIndex++)
1949  {
1950  DoubleVector initialPhase = channelParams->m_clusterPhase[nIndex][mIndex];
1951  NS_ASSERT(4 <= initialPhase.size());
1952  double k = channelParams->m_crossPolarizationPowerRatios[nIndex][mIndex];
1953 
1954  // cache the component of the "rays" terms which depend on the random angle of arrivals
1955  // and departures and initial phases only
1956  auto [rxFieldPatternPhi, rxFieldPatternTheta] = uAntenna->GetElementFieldPattern(
1957  Angles(channelParams->m_rayAoaRadian[nIndex][mIndex],
1958  channelParams->m_rayZoaRadian[nIndex][mIndex]));
1959  auto [txFieldPatternPhi, txFieldPatternTheta] = sAntenna->GetElementFieldPattern(
1960  Angles(channelParams->m_rayAodRadian[nIndex][mIndex],
1961  channelParams->m_rayZodRadian[nIndex][mIndex]));
1962  raysPreComp(nIndex, mIndex) =
1963  std::complex<double>(cos(initialPhase[0]), sin(initialPhase[0])) *
1964  rxFieldPatternTheta * txFieldPatternTheta +
1965  std::complex<double>(cos(initialPhase[1]), sin(initialPhase[1])) *
1966  std::sqrt(1.0 / k) * rxFieldPatternTheta * txFieldPatternPhi +
1967  std::complex<double>(cos(initialPhase[2]), sin(initialPhase[2])) *
1968  std::sqrt(1.0 / k) * rxFieldPatternPhi * txFieldPatternTheta +
1969  std::complex<double>(cos(initialPhase[3]), sin(initialPhase[3])) *
1970  rxFieldPatternPhi * txFieldPatternPhi;
1971 
1972  // cache the component of the "rxPhaseDiff" terms which depend on the random angle of
1973  // arrivals only
1974  double sinRayZoa = sin(rayZoaRadian[nIndex][mIndex]);
1975  double sinRayAoa = cos(rayAoaRadian[nIndex][mIndex]);
1976  double cosRayAoa = cos(rayAoaRadian[nIndex][mIndex]);
1977  sinCosA[nIndex][mIndex] = sinRayZoa * cosRayAoa;
1978  sinSinA[nIndex][mIndex] = sinRayZoa * sinRayAoa;
1979  cosZoA[nIndex][mIndex] = cos(rayZoaRadian[nIndex][mIndex]);
1980 
1981  // cache the component of the "txPhaseDiff" terms which depend on the random angle of
1982  // departure only
1983  double sinRayZod = sin(rayZodRadian[nIndex][mIndex]);
1984  double sinRayAod = cos(rayAodRadian[nIndex][mIndex]);
1985  double cosRayAod = cos(rayAodRadian[nIndex][mIndex]);
1986  sinCosD[nIndex][mIndex] = sinRayZod * cosRayAod;
1987  sinSinD[nIndex][mIndex] = sinRayZod * sinRayAod;
1988  cosZoD[nIndex][mIndex] = cos(rayZodRadian[nIndex][mIndex]);
1989  }
1990  }
1991 
1992  // The following for loops computes the channel coefficients
1993  // Keeps track of how many sub-clusters have been added up to now
1994  uint8_t numSubClustersAdded = 0;
1995  for (uint8_t nIndex = 0; nIndex < channelParams->m_reducedClusterNumber; nIndex++)
1996  {
1997  for (size_t uIndex = 0; uIndex < uSize; uIndex++)
1998  {
1999  Vector uLoc = uAntenna->GetElementLocation(uIndex);
2000 
2001  for (size_t sIndex = 0; sIndex < sSize; sIndex++)
2002  {
2003  Vector sLoc = sAntenna->GetElementLocation(sIndex);
2004  // Compute the N-2 weakest cluster, assuming 0 slant angle and a
2005  // polarization slant angle configured in the array (7.5-22)
2006  if (nIndex != channelParams->m_cluster1st && nIndex != channelParams->m_cluster2nd)
2007  {
2008  std::complex<double> rays(0, 0);
2009  for (uint8_t mIndex = 0; mIndex < table3gpp->m_raysPerCluster; mIndex++)
2010  {
2011  // lambda_0 is accounted in the antenna spacing uLoc and sLoc.
2012  double rxPhaseDiff =
2013  2 * M_PI *
2014  (sinCosA[nIndex][mIndex] * uLoc.x + sinSinA[nIndex][mIndex] * uLoc.y +
2015  cosZoA[nIndex][mIndex] * uLoc.z);
2016 
2017  double txPhaseDiff =
2018  2 * M_PI *
2019  (sinCosD[nIndex][mIndex] * sLoc.x + sinSinD[nIndex][mIndex] * sLoc.y +
2020  cosZoD[nIndex][mIndex] * sLoc.z);
2021 
2022  // NOTE Doppler is computed in the CalcBeamformingGain function and is
2023  // simplified to only account for the center angle of each cluster.
2024  rays += raysPreComp(nIndex, mIndex) *
2025  std::complex<double>(cos(rxPhaseDiff), sin(rxPhaseDiff)) *
2026  std::complex<double>(cos(txPhaseDiff), sin(txPhaseDiff));
2027  }
2028  rays *=
2029  sqrt(channelParams->m_clusterPower[nIndex] / table3gpp->m_raysPerCluster);
2030  hUsn(uIndex, sIndex, nIndex) = rays;
2031  }
2032  else //(7.5-28)
2033  {
2034  std::complex<double> raysSub1(0, 0);
2035  std::complex<double> raysSub2(0, 0);
2036  std::complex<double> raysSub3(0, 0);
2037 
2038  for (uint8_t mIndex = 0; mIndex < table3gpp->m_raysPerCluster; mIndex++)
2039  {
2040  // ZML:Just remind me that the angle offsets for the 3 subclusters were not
2041  // generated correctly.
2042  double rxPhaseDiff =
2043  2 * M_PI *
2044  (sinCosA[nIndex][mIndex] * uLoc.x + sinSinA[nIndex][mIndex] * uLoc.y +
2045  cosZoA[nIndex][mIndex] * uLoc.z);
2046 
2047  double txPhaseDiff =
2048  2 * M_PI *
2049  (sinCosD[nIndex][mIndex] * sLoc.x + sinSinD[nIndex][mIndex] * sLoc.y +
2050  cosZoD[nIndex][mIndex] * sLoc.z);
2051 
2052  std::complex<double> raySub =
2053  raysPreComp(nIndex, mIndex) *
2054  std::complex<double>(cos(rxPhaseDiff), sin(rxPhaseDiff)) *
2055  std::complex<double>(cos(txPhaseDiff), sin(txPhaseDiff));
2056 
2057  switch (mIndex)
2058  {
2059  case 9:
2060  case 10:
2061  case 11:
2062  case 12:
2063  case 17:
2064  case 18:
2065  raysSub2 += raySub;
2066  break;
2067  case 13:
2068  case 14:
2069  case 15:
2070  case 16:
2071  raysSub3 += raySub;
2072  break;
2073  default: // case 1,2,3,4,5,6,7,8,19,20
2074  raysSub1 += raySub;
2075  break;
2076  }
2077  }
2078  raysSub1 *=
2079  sqrt(channelParams->m_clusterPower[nIndex] / table3gpp->m_raysPerCluster);
2080  raysSub2 *=
2081  sqrt(channelParams->m_clusterPower[nIndex] / table3gpp->m_raysPerCluster);
2082  raysSub3 *=
2083  sqrt(channelParams->m_clusterPower[nIndex] / table3gpp->m_raysPerCluster);
2084  hUsn(uIndex, sIndex, nIndex) = raysSub1;
2085  hUsn(uIndex,
2086  sIndex,
2087  channelParams->m_reducedClusterNumber + numSubClustersAdded) = raysSub2;
2088  hUsn(uIndex,
2089  sIndex,
2090  channelParams->m_reducedClusterNumber + numSubClustersAdded + 1) =
2091  raysSub3;
2092  }
2093  }
2094  }
2095  if (nIndex == channelParams->m_cluster1st || nIndex == channelParams->m_cluster2nd)
2096  {
2097  numSubClustersAdded += 2;
2098  }
2099  }
2100 
2101  if (channelParams->m_losCondition == ChannelCondition::LOS) //(7.5-29) && (7.5-30)
2102  {
2103  double lambda = 3.0e8 / m_frequency; // the wavelength of the carrier frequency
2104  std::complex<double> phaseDiffDueToDistance(cos(-2 * M_PI * distance3D / lambda),
2105  sin(-2 * M_PI * distance3D / lambda));
2106 
2107  const double sinUAngleIncl = sin(uAngle.GetInclination());
2108  const double cosUAngleIncl = cos(uAngle.GetInclination());
2109  const double sinUAngleAz = sin(uAngle.GetAzimuth());
2110  const double cosUAngleAz = cos(uAngle.GetAzimuth());
2111  const double sinSAngleIncl = sin(sAngle.GetInclination());
2112  const double cosSAngleIncl = cos(sAngle.GetInclination());
2113  const double sinSAngleAz = sin(sAngle.GetAzimuth());
2114  const double cosSAngleAz = cos(sAngle.GetAzimuth());
2115 
2116  for (size_t uIndex = 0; uIndex < uSize; uIndex++)
2117  {
2118  Vector uLoc = uAntenna->GetElementLocation(uIndex);
2119  double rxPhaseDiff = 2 * M_PI *
2120  (sinUAngleIncl * cosUAngleAz * uLoc.x +
2121  sinUAngleIncl * sinUAngleAz * uLoc.y + cosUAngleIncl * uLoc.z);
2122 
2123  for (size_t sIndex = 0; sIndex < sSize; sIndex++)
2124  {
2125  Vector sLoc = sAntenna->GetElementLocation(sIndex);
2126  std::complex<double> ray(0, 0);
2127  double txPhaseDiff =
2128  2 * M_PI *
2129  (sinSAngleIncl * cosSAngleAz * sLoc.x + sinSAngleIncl * sinSAngleAz * sLoc.y +
2130  cosSAngleIncl * sLoc.z);
2131 
2132  auto [rxFieldPatternPhi, rxFieldPatternTheta] = uAntenna->GetElementFieldPattern(
2133  Angles(uAngle.GetAzimuth(), uAngle.GetInclination()));
2134  auto [txFieldPatternPhi, txFieldPatternTheta] = sAntenna->GetElementFieldPattern(
2135  Angles(sAngle.GetAzimuth(), sAngle.GetInclination()));
2136 
2137  ray = (rxFieldPatternTheta * txFieldPatternTheta -
2138  rxFieldPatternPhi * txFieldPatternPhi) *
2139  phaseDiffDueToDistance *
2140  std::complex<double>(cos(rxPhaseDiff), sin(rxPhaseDiff)) *
2141  std::complex<double>(cos(txPhaseDiff), sin(txPhaseDiff));
2142 
2143  double kLinear = pow(10, channelParams->m_K_factor / 10.0);
2144  // the LOS path should be attenuated if blockage is enabled.
2145  hUsn(uIndex, sIndex, 0) =
2146  sqrt(1.0 / (kLinear + 1)) * hUsn(uIndex, sIndex, 0) +
2147  sqrt(kLinear / (1 + kLinear)) * ray /
2148  pow(10,
2149  channelParams->m_attenuation_dB[0] / 10.0); //(7.5-30) for tau = tau1
2150  for (uint16_t nIndex = 1; nIndex < hUsn.GetNumPages(); nIndex++)
2151  {
2152  hUsn(uIndex, sIndex, nIndex) *=
2153  sqrt(1.0 / (kLinear + 1)); //(7.5-30) for tau = tau2...tauN
2154  }
2155  }
2156  }
2157  }
2158 
2159  NS_LOG_DEBUG("Husn (sAntenna, uAntenna):" << sAntenna->GetId() << ", " << uAntenna->GetId());
2160  for (uint16_t cIndex = 0; cIndex < hUsn.GetNumPages(); cIndex++)
2161  {
2162  for (uint16_t rowIdx = 0; rowIdx < hUsn.GetNumRows(); rowIdx++)
2163  {
2164  for (uint16_t colIdx = 0; colIdx < hUsn.GetNumCols(); colIdx++)
2165  {
2166  NS_LOG_DEBUG(" " << hUsn(rowIdx, colIdx, cIndex) << ",");
2167  }
2168  }
2169  }
2170 
2171  NS_LOG_INFO("size of coefficient matrix (rows, columns, clusters) = ("
2172  << hUsn.GetNumRows() << ", " << hUsn.GetNumCols() << ", " << hUsn.GetNumPages()
2173  << ")");
2174  channelMatrix->m_channel = hUsn;
2175  return channelMatrix;
2176 }
2177 
2178 std::pair<double, double>
2179 ThreeGppChannelModel::WrapAngles(double azimuthRad, double inclinationRad)
2180 {
2181  inclinationRad = WrapTo2Pi(inclinationRad);
2182  if (inclinationRad > M_PI)
2183  {
2184  // inclination must be in [0, M_PI]
2185  inclinationRad -= M_PI;
2186  azimuthRad += M_PI;
2187  }
2188 
2189  azimuthRad = WrapTo2Pi(azimuthRad);
2190 
2191  NS_ASSERT_MSG(0 <= inclinationRad && inclinationRad <= M_PI,
2192  "inclinationRad=" << inclinationRad << " not valid, should be in [0, pi]");
2193  NS_ASSERT_MSG(0 <= azimuthRad && azimuthRad <= 2 * M_PI,
2194  "azimuthRad=" << azimuthRad << " not valid, should be in [0, 2*pi]");
2195 
2196  return std::make_pair(azimuthRad, inclinationRad);
2197 }
2198 
2202  const DoubleVector& clusterAOA,
2203  const DoubleVector& clusterZOA) const
2204 {
2205  NS_LOG_FUNCTION(this);
2206 
2207  DoubleVector powerAttenuation;
2208  uint8_t clusterNum = clusterAOA.size();
2209  for (uint8_t cInd = 0; cInd < clusterNum; cInd++)
2210  {
2211  powerAttenuation.push_back(0); // Initial power attenuation for all clusters to be 0 dB;
2212  }
2213  // step a: the number of non-self blocking blockers is stored in m_numNonSelfBlocking.
2214 
2215  // step b:Generate the size and location of each blocker
2216  // generate self blocking (i.e., for blockage from the human body)
2217  // table 7.6.4.1-1 Self-blocking region parameters.
2218  // Defaults: landscape mode
2219  double phiSb = 40;
2220  double xSb = 160;
2221  double thetaSb = 110;
2222  double ySb = 75;
2223  if (m_portraitMode)
2224  {
2225  phiSb = 260;
2226  xSb = 120;
2227  thetaSb = 100;
2228  ySb = 80;
2229  }
2230 
2231  // generate or update non-self blocking
2232  if (channelParams->m_nonSelfBlocking.empty()) // generate new blocking regions
2233  {
2234  for (uint16_t blockInd = 0; blockInd < m_numNonSelfBlocking; blockInd++)
2235  {
2236  // draw value from table 7.6.4.1-2 Blocking region parameters
2237  DoubleVector table;
2238  table.push_back(m_normalRv->GetValue()); // phi_k: store the normal RV that will be
2239  // mapped to uniform (0,360) later.
2240  if (m_scenario == "InH-OfficeMixed" || m_scenario == "InH-OfficeOpen")
2241  {
2242  table.push_back(m_uniformRv->GetValue(15, 45)); // x_k
2243  table.push_back(90); // Theta_k
2244  table.push_back(m_uniformRv->GetValue(5, 15)); // y_k
2245  table.push_back(2); // r
2246  }
2247  else
2248  {
2249  table.push_back(m_uniformRv->GetValue(5, 15)); // x_k
2250  table.push_back(90); // Theta_k
2251  table.push_back(5); // y_k
2252  table.push_back(10); // r
2253  }
2254  channelParams->m_nonSelfBlocking.push_back(table);
2255  }
2256  }
2257  else
2258  {
2259  double deltaX = sqrt(pow(channelParams->m_preLocUT.x - channelParams->m_locUT.x, 2) +
2260  pow(channelParams->m_preLocUT.y - channelParams->m_locUT.y, 2));
2261  // if deltaX and speed are both 0, the autocorrelation is 1, skip updating
2262  if (deltaX > 1e-6 || m_blockerSpeed > 1e-6)
2263  {
2264  double corrDis;
2265  // draw value from table 7.6.4.1-4: Spatial correlation distance for different
2266  // m_scenarios.
2267  if (m_scenario == "InH-OfficeMixed" || m_scenario == "InH-OfficeOpen")
2268  {
2269  // InH, correlation distance = 5;
2270  corrDis = 5;
2271  }
2272  else
2273  {
2274  if (channelParams->m_o2iCondition == ChannelCondition::O2I) // outdoor to indoor
2275  {
2276  corrDis = 5;
2277  }
2278  else // LOS or NLOS
2279  {
2280  corrDis = 10;
2281  }
2282  }
2283  double R;
2284  if (m_blockerSpeed > 1e-6) // speed not equal to 0
2285  {
2286  double corrT = corrDis / m_blockerSpeed;
2287  R = exp(-1 * (deltaX / corrDis +
2288  (Now().GetSeconds() - channelParams->m_generatedTime.GetSeconds()) /
2289  corrT));
2290  }
2291  else
2292  {
2293  R = exp(-1 * (deltaX / corrDis));
2294  }
2295 
2296  NS_LOG_INFO("Distance change:"
2297  << deltaX << " Speed:" << m_blockerSpeed << " Time difference:"
2298  << Now().GetSeconds() - channelParams->m_generatedTime.GetSeconds()
2299  << " correlation:" << R);
2300 
2301  // In order to generate correlated uniform random variables, we first generate
2302  // correlated normal random variables and map the normal RV to uniform RV. Notice the
2303  // correlation will change if the RV is transformed from normal to uniform. To
2304  // compensate the distortion, the correlation of the normal RV is computed such that the
2305  // uniform RV would have the desired correlation when transformed from normal RV.
2306 
2307  // The following formula was obtained from MATLAB numerical simulation.
2308 
2309  if (R * R * (-0.069) + R * 1.074 - 0.002 <
2310  1) // transform only when the correlation of normal RV is smaller than 1
2311  {
2312  R = R * R * (-0.069) + R * 1.074 - 0.002;
2313  }
2314  for (uint16_t blockInd = 0; blockInd < m_numNonSelfBlocking; blockInd++)
2315  {
2316  // Generate a new correlated normal RV with the following formula
2317  channelParams->m_nonSelfBlocking[blockInd][PHI_INDEX] =
2318  R * channelParams->m_nonSelfBlocking[blockInd][PHI_INDEX] +
2319  sqrt(1 - R * R) * m_normalRv->GetValue();
2320  }
2321  }
2322  }
2323 
2324  // step c: Determine the attenuation of each blocker due to blockers
2325  for (uint8_t cInd = 0; cInd < clusterNum; cInd++)
2326  {
2327  NS_ASSERT_MSG(clusterAOA[cInd] >= 0 && clusterAOA[cInd] <= 360,
2328  "the AOA should be the range of [0,360]");
2329  NS_ASSERT_MSG(clusterZOA[cInd] >= 0 && clusterZOA[cInd] <= 180,
2330  "the ZOA should be the range of [0,180]");
2331 
2332  // check self blocking
2333  NS_LOG_INFO("AOA=" << clusterAOA[cInd] << " Block Region[" << phiSb - xSb / 2.0 << ","
2334  << phiSb + xSb / 2.0 << "]");
2335  NS_LOG_INFO("ZOA=" << clusterZOA[cInd] << " Block Region[" << thetaSb - ySb / 2.0 << ","
2336  << thetaSb + ySb / 2.0 << "]");
2337  if (std::abs(clusterAOA[cInd] - phiSb) < (xSb / 2.0) &&
2338  std::abs(clusterZOA[cInd] - thetaSb) < (ySb / 2.0))
2339  {
2340  powerAttenuation[cInd] += 30; // attenuate by 30 dB.
2341  NS_LOG_INFO("Cluster[" << +cInd
2342  << "] is blocked by self blocking region and reduce 30 dB power,"
2343  "the attenuation is ["
2344  << powerAttenuation[cInd] << " dB]");
2345  }
2346 
2347  // check non-self blocking
2348  for (uint16_t blockInd = 0; blockInd < m_numNonSelfBlocking; blockInd++)
2349  {
2350  // The normal RV is transformed to uniform RV with the desired correlation.
2351  double phiK =
2352  (0.5 * erfc(-1 * channelParams->m_nonSelfBlocking[blockInd][PHI_INDEX] / sqrt(2))) *
2353  360;
2354  while (phiK > 360)
2355  {
2356  phiK -= 360;
2357  }
2358 
2359  while (phiK < 0)
2360  {
2361  phiK += 360;
2362  }
2363 
2364  double xK = channelParams->m_nonSelfBlocking[blockInd][X_INDEX];
2365  double thetaK = channelParams->m_nonSelfBlocking[blockInd][THETA_INDEX];
2366  double yK = channelParams->m_nonSelfBlocking[blockInd][Y_INDEX];
2367 
2368  NS_LOG_INFO("AOA=" << clusterAOA[cInd] << " Block Region[" << phiK - xK << ","
2369  << phiK + xK << "]");
2370  NS_LOG_INFO("ZOA=" << clusterZOA[cInd] << " Block Region[" << thetaK - yK << ","
2371  << thetaK + yK << "]");
2372 
2373  if (std::abs(clusterAOA[cInd] - phiK) < (xK) &&
2374  std::abs(clusterZOA[cInd] - thetaK) < (yK))
2375  {
2376  double A1 = clusterAOA[cInd] - (phiK + xK / 2.0); //(7.6-24)
2377  double A2 = clusterAOA[cInd] - (phiK - xK / 2.0); //(7.6-25)
2378  double Z1 = clusterZOA[cInd] - (thetaK + yK / 2.0); //(7.6-26)
2379  double Z2 = clusterZOA[cInd] - (thetaK - yK / 2.0); //(7.6-27)
2380  int signA1;
2381  int signA2;
2382  int signZ1;
2383  int signZ2;
2384  // draw sign for the above parameters according to table 7.6.4.1-3 Description of
2385  // signs
2386  if (xK / 2.0 < clusterAOA[cInd] - phiK && clusterAOA[cInd] - phiK <= xK)
2387  {
2388  signA1 = -1;
2389  }
2390  else
2391  {
2392  signA1 = 1;
2393  }
2394  if (-1 * xK < clusterAOA[cInd] - phiK && clusterAOA[cInd] - phiK <= -1 * xK / 2.0)
2395  {
2396  signA2 = -1;
2397  }
2398  else
2399  {
2400  signA2 = 1;
2401  }
2402 
2403  if (yK / 2.0 < clusterZOA[cInd] - thetaK && clusterZOA[cInd] - thetaK <= yK)
2404  {
2405  signZ1 = -1;
2406  }
2407  else
2408  {
2409  signZ1 = 1;
2410  }
2411  if (-1 * yK < clusterZOA[cInd] - thetaK &&
2412  clusterZOA[cInd] - thetaK <= -1 * yK / 2.0)
2413  {
2414  signZ2 = -1;
2415  }
2416  else
2417  {
2418  signZ2 = 1;
2419  }
2420  double lambda = 3e8 / m_frequency;
2421  double fA1 =
2422  atan(signA1 * M_PI / 2.0 *
2423  sqrt(M_PI / lambda * channelParams->m_nonSelfBlocking[blockInd][R_INDEX] *
2424  (1.0 / cos(DegreesToRadians(A1)) - 1))) /
2425  M_PI; //(7.6-23)
2426  double fA2 =
2427  atan(signA2 * M_PI / 2.0 *
2428  sqrt(M_PI / lambda * channelParams->m_nonSelfBlocking[blockInd][R_INDEX] *
2429  (1.0 / cos(DegreesToRadians(A2)) - 1))) /
2430  M_PI;
2431  double fZ1 =
2432  atan(signZ1 * M_PI / 2.0 *
2433  sqrt(M_PI / lambda * channelParams->m_nonSelfBlocking[blockInd][R_INDEX] *
2434  (1.0 / cos(DegreesToRadians(Z1)) - 1))) /
2435  M_PI;
2436  double fZ2 =
2437  atan(signZ2 * M_PI / 2.0 *
2438  sqrt(M_PI / lambda * channelParams->m_nonSelfBlocking[blockInd][R_INDEX] *
2439  (1.0 / cos(DegreesToRadians(Z2)) - 1))) /
2440  M_PI;
2441  double lDb = -20 * log10(1 - (fA1 + fA2) * (fZ1 + fZ2)); //(7.6-22)
2442  powerAttenuation[cInd] += lDb;
2443  NS_LOG_INFO("Cluster[" << +cInd << "] is blocked by no-self blocking, the loss is ["
2444  << lDb << "] dB");
2445  }
2446  }
2447  }
2448  return powerAttenuation;
2449 }
2450 
2451 void
2452 ThreeGppChannelModel::Shuffle(double* first, double* last) const
2453 {
2454  for (auto i = (last - first) - 1; i > 0; --i)
2455  {
2456  std::swap(first[i], first[m_uniformRvShuffle->GetInteger(0, i)]);
2457  }
2458 }
2459 
2460 int64_t
2462 {
2463  NS_LOG_FUNCTION(this << stream);
2464  m_normalRv->SetStream(stream);
2465  m_uniformRv->SetStream(stream + 1);
2466  m_uniformRvShuffle->SetStream(stream + 2);
2467  m_uniformRvDoppler->SetStream(stream + 3);
2468  return 4;
2469 }
2470 
2471 } // namespace ns3
#define min(a, b)
Definition: 80211b.c:42
double f(double x, void *params)
Definition: 80211b.c:71
#define max(a, b)
Definition: 80211b.c:43
Class holding the azimuth and inclination angles of spherical coordinates.
Definition: angles.h:118
double GetInclination() const
Getter for inclination angle.
Definition: angles.cc:216
double GetAzimuth() const
Getter for azimuth angle.
Definition: angles.cc:210
AttributeValue implementation for Boolean.
Definition: boolean.h:37
This class can be used to hold variables of floating point type such as 'double' or 'float'.
Definition: double.h:42
Hold a signed integer type.
Definition: integer.h:45
MatrixArray class inherits ValArray class and provides additional interfaces to ValArray which enable...
Definition: matrix-array.h:83
This is an interface for a channel model that can be described by a channel matrix,...
std::vector< double > DoubleVector
Type definition for vectors of doubles.
std::vector< Double2DVector > Double3DVector
Type definition for 3D matrices of doubles.
std::vector< DoubleVector > Double2DVector
Type definition for matrices of doubles.
static uint64_t GetKey(uint32_t a, uint32_t b)
Generate a unique value for the pair of unsigned integer of 32 bits, where the order does not matter,...
A network Node.
Definition: node.h:56
uint32_t GetId() const
Definition: node.cc:117
Hold objects of type Ptr<T>.
Definition: pointer.h:37
Smart pointer class similar to boost::intrusive_ptr.
Definition: ptr.h:78
void SetStream(int64_t stream)
Specifies the stream number for the RngStream.
static Time Now()
Return the current simulation virtual time.
Definition: simulator.cc:199
Hold variables of type string.
Definition: string.h:56
DoubleVector CalcAttenuationOfBlockage(const Ptr< ThreeGppChannelModel::ThreeGppChannelParams > channelParams, const DoubleVector &clusterAOA, const DoubleVector &clusterZOA) const
Applies the blockage model A described in 3GPP TR 38.901.
int64_t AssignStreams(int64_t stream)
Assign a fixed random variable stream number to the random variables used by this model.
bool m_portraitMode
true if portrait mode, false if landscape
bool ChannelParamsNeedsUpdate(Ptr< const ThreeGppChannelParams > channelParams, Ptr< const ChannelCondition > channelCondition) const
Check if the channel params has to be updated.
virtual Ptr< const ParamsTable > GetThreeGppTable(Ptr< const ChannelCondition > channelCondition, double hBS, double hUT, double distance2D) const
Get the parameters needed to apply the channel generation procedure.
Ptr< NormalRandomVariable > m_normalRv
normal random variable
static const uint8_t Y_INDEX
index of the Y value in the m_nonSelfBlocking array
bool m_blockage
enables the blockage model A
Ptr< const ChannelParams > GetParams(Ptr< const MobilityModel > aMob, Ptr< const MobilityModel > bMob) const override
Looks for the channel params associated to the aMob and bMob pair in m_channelParamsMap.
~ThreeGppChannelModel() override
Destructor.
bool ChannelMatrixNeedsUpdate(Ptr< const ThreeGppChannelParams > channelParams, Ptr< const ChannelMatrix > channelMatrix)
Check if the channel matrix has to be updated (it needs update when the channel params generation tim...
static const uint8_t THETA_INDEX
index of the THETA value in the m_nonSelfBlocking array
std::unordered_map< uint64_t, Ptr< ThreeGppChannelParams > > m_channelParamsMap
map containing the common channel parameters per pair of nodes, the key of this map is reciprocal and...
static std::pair< double, double > WrapAngles(double azimuthRad, double inclinationRad)
Wrap an (azimuth, inclination) angle pair in a valid range.
double m_blockerSpeed
the blocker speed
Ptr< const ChannelMatrix > GetChannel(Ptr< const MobilityModel > aMob, Ptr< const MobilityModel > bMob, Ptr< const PhasedArrayModel > aAntenna, Ptr< const PhasedArrayModel > bAntenna) override
Looks for the channel matrix associated to the aMob and bMob pair in m_channelMatrixMap.
void SetFrequency(double f)
Sets the center frequency of the model.
std::unordered_map< uint64_t, Ptr< ChannelMatrix > > m_channelMatrixMap
map containing the channel realizations per pair of PhasedAntennaArray instances, the key of this map...
Ptr< UniformRandomVariable > m_uniformRv
uniform random variable
void DoDispose() override
Destructor implementation.
void SetScenario(const std::string &scenario)
Sets the propagation scenario.
void SetChannelConditionModel(Ptr< ChannelConditionModel > model)
Set the channel condition model.
Ptr< UniformRandomVariable > m_uniformRvDoppler
uniform random variable, used to compute the additional Doppler contribution
uint16_t m_numNonSelfBlocking
number of non-self-blocking regions
std::string GetScenario() const
Returns the propagation scenario.
virtual Ptr< ChannelMatrix > GetNewChannel(Ptr< const ThreeGppChannelParams > channelParams, Ptr< const ParamsTable > table3gpp, const Ptr< const MobilityModel > sMob, const Ptr< const MobilityModel > uMob, Ptr< const PhasedArrayModel > sAntenna, Ptr< const PhasedArrayModel > uAntenna) const
Compute the channel matrix between two nodes a and b, and their antenna arrays aAntenna and bAntenna ...
static const uint8_t PHI_INDEX
index of the PHI value in the m_nonSelfBlocking array
double m_frequency
the operating frequency
double m_vScatt
value used to compute the additional Doppler contribution for the delayed paths
Ptr< ChannelConditionModel > GetChannelConditionModel() const
Get the associated channel condition model.
Ptr< ChannelConditionModel > m_channelConditionModel
the channel condition model
std::string m_scenario
the 3GPP scenario
static const uint8_t R_INDEX
index of the R value in the m_nonSelfBlocking array
static TypeId GetTypeId()
Get the type ID.
void Shuffle(double *first, double *last) const
Shuffle the elements of a simple sequence container of type double.
Ptr< ThreeGppChannelParams > GenerateChannelParameters(const Ptr< const ChannelCondition > channelCondition, const Ptr< const ParamsTable > table3gpp, const Ptr< const MobilityModel > aMob, const Ptr< const MobilityModel > bMob) const
Prepare 3gpp channel parameters among the nodes a and b.
double GetFrequency() const
Returns the center frequency.
Time m_updatePeriod
the channel update period
static const uint8_t X_INDEX
index of the X value in the m_nonSelfBlocking array
Ptr< UniformRandomVariable > m_uniformRvShuffle
uniform random variable used to shuffle array in GetNewChannel
@ NS
nanosecond
Definition: nstime.h:119
bool IsZero() const
Exactly equivalent to t == 0.
Definition: nstime.h:314
AttributeValue implementation for Time.
Definition: nstime.h:1423
a unique identifier for an interface.
Definition: type-id.h:60
TypeId SetGroupName(std::string groupName)
Set the group name.
Definition: type-id.cc:943
TypeId SetParent(TypeId tid)
Set the parent TypeId.
Definition: type-id.cc:935
double GetValue(double min, double max)
Get the next random value drawn from the distribution.
uint32_t GetInteger(uint32_t min, uint32_t max)
Get the next random value drawn from the distribution.
#define NS_ASSERT(condition)
At runtime, in debugging builds, if this condition is not true, the program prints the source file,...
Definition: assert.h:66
#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:86
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:86
Ptr< const AttributeChecker > MakeBooleanChecker()
Definition: boolean.cc:124
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:43
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:46
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:231
Ptr< const AttributeChecker > MakeStringChecker()
Definition: string.cc:30
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:57
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:1424
#define NS_FATAL_ERROR(msg)
Report a fatal error with a message and terminate.
Definition: fatal-error.h:179
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition: log.h:202
#define NS_LOG_DEBUG(msg)
Use NS_LOG to output a message of level LOG_DEBUG.
Definition: log.h:268
#define NS_LOG_FUNCTION(parameters)
If log level LOG_FUNCTION is enabled, this macro will output all input parameters separated by ",...
#define NS_LOG_WARN(msg)
Use NS_LOG to output a message of level LOG_WARN.
Definition: log.h:261
#define NS_LOG_INFO(msg)
Use NS_LOG to output a message of level LOG_INFO.
Definition: log.h:275
#define NS_OBJECT_ENSURE_REGISTERED(type)
Register an Object subclass with the TypeId system.
Definition: object-base.h:46
Time Now()
create an ns3::Time instance which contains the current simulation time.
Definition: simulator.cc:296
Time MilliSeconds(uint64_t value)
Construct a Time in the indicated unit.
Definition: nstime.h:1348
Definition: first.py:1
Every class exported by the ns3 library is enclosed in the ns3 namespace.
static const double offSetAlpha[20]
The ray offset angles within a cluster, given for rms angle spread normalized to 1.
static const double sqrtC_RMa_O2I[6][6]
The square root matrix for RMa O2I, which is generated using the Cholesky decomposition according to ...
static const double sqrtC_UMi_LOS[7][7]
The square root matrix for UMi LOS, which is generated using the Cholesky decomposition according to ...
static const double sqrtC_office_LOS[7][7]
The square root matrix for Indoor-Office LOS, which is generated using the Cholesky decomposition acc...
static const double sqrtC_UMa_O2I[6][6]
The square root matrix for UMa O2I, which is generated using the Cholesky decomposition according to ...
static const double sqrtC_RMa_NLOS[6][6]
The square root matrix for RMa NLOS, which is generated using the Cholesky decomposition according to...
static const double sqrtC_UMa_LOS[7][7]
The square root matrix for UMa LOS, which is generated using the Cholesky decomposition according to ...
static const double sqrtC_UMi_NLOS[6][6]
The square root matrix for UMi NLOS, which is generated using the Cholesky decomposition according to...
Ptr< const AttributeChecker > MakeTimeChecker(const Time min, const Time max)
Helper to make a Time checker with bounded range.
Definition: time.cc:535
static const double sqrtC_RMa_LOS[7][7]
The square root matrix for RMa LOS, which is generated using the Cholesky decomposition according to ...
double DegreesToRadians(double degrees)
converts degrees to radians
Definition: angles.cc:39
static const double sqrtC_UMi_O2I[6][6]
The square root matrix for UMi O2I, which is generated using the Cholesky decomposition according to ...
static const double sqrtC_office_NLOS[6][6]
The square root matrix for Indoor-Office NLOS, which is generated using the Cholesky decomposition ac...
static const double sqrtC_UMa_NLOS[6][6]
The square root matrix for UMa NLOS, which is generated using the Cholesky decomposition according to...
double WrapTo2Pi(double a)
Wrap angle in [0, 2*M_PI)
Definition: angles.cc:102
double RadiansToDegrees(double radians)
converts radians to degrees
Definition: angles.cc:45
Complex3DVector m_channel
Channel matrix H[u][s][n].
std::pair< uint32_t, uint32_t > m_antennaPair
The first element is the ID of the antenna of the s-node (the antenna of the transmitter when the cha...
std::pair< uint32_t, uint32_t > m_nodeIds
The first element is the s-node ID (the transmitter when the channel was generated),...