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