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