src/Radio/radio.cxx

   1 // radio.cxx -- implementation of FGRadio
   2 // Class to manage radio propagation using the ITM model
   3 // Written by Adrian Musceac YO8RZZ, started August 2011.
   4 //
   5 // This program is free software; you can redistribute it and/or
   6 // modify it under the terms of the GNU General Public License as
   7 // published by the Free Software Foundation; either version 2 of the
   8 // License, or (at your option) any later version.
   9 //
  10 // This program is distributed in the hope that it will be useful, but
  11 // WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 // General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU General Public License
  16 // along with this program; if not, write to the Free Software
  17 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  18
  19
  20
  21 #ifdef HAVE_CONFIG_H
  22 #  include <config.h>
  23 #endif
  24
  25 #include <math.h>
  26
  27 #include <stdlib.h>
  28 #include <deque>
  29 #include "radio.hxx"
  30 #include <simgear/scene/material/mat.hxx>
  31 #include <Scenery/scenery.hxx>
  32 #include <boost/scoped_array.hpp>
  33
  34 #define WITH_POINT_TO_POINT 1
  35 #include "itm.cpp"
  36
  37
  38 FGRadioTransmission::FGRadioTransmission() {
  39
  40
  41         _receiver_sensitivity = -105.0; // typical AM receiver sensitivity seems to be 0.8 microVolt at 12dB SINAD or less
  42
  43         /** AM transmitter power in dBm.
  44         *       Typical output powers for ATC ground equipment, VHF-UHF:
  45         *       40 dBm - 10 W (ground, clearance)
  46         *       44 dBm - 20 W (tower)
  47         *       47 dBm - 50 W (center, sectors)
  48         *       50 dBm - 100 W (center, sectors)
  49         *       53 dBm - 200 W (sectors, on directional arrays)
  50         **/
  51         _transmitter_power = 43.0;
  52
  53         _tx_antenna_height = 2.0; // TX antenna height above ground level
  54
  55         _rx_antenna_height = 2.0; // RX antenna height above ground level
  56
  57
  58         _rx_antenna_gain = 1.0; // maximum antenna gain expressed in dBi
  59         _tx_antenna_gain = 1.0;
  60
  61         _rx_line_losses = 2.0;  // to be configured for each station
  62         _tx_line_losses = 2.0;
  63
  64         _polarization = 1; // default vertical
  65
  66         _propagation_model = 2;
  67
  68         _root_node = fgGetNode("sim/radio", true);
  69         _terrain_sampling_distance = _root_node->getDoubleValue("sampling-distance", 90.0); // regular SRTM is 90 meters
  70
  71
  72 }
  73
  74 FGRadioTransmission::~FGRadioTransmission()
  75 {
  76 }
  77
  78
  79 double FGRadioTransmission::getFrequency(int radio) {
  80         double freq = 118.0;
  81         switch (radio) {
  82                 case 1:
  83                         freq = fgGetDouble("/instrumentation/comm[0]/frequencies/selected-mhz");
  84                         break;
  85                 case 2:
  86                         freq = fgGetDouble("/instrumentation/comm[1]/frequencies/selected-mhz");
  87                         break;
  88                 default:
  89                         freq = fgGetDouble("/instrumentation/comm[0]/frequencies/selected-mhz");
  90
  91         }
  92         return freq;
  93 }
  94
  95
  96 void FGRadioTransmission::receiveChat(SGGeod tx_pos, double freq, string text, int ground_to_air) {
  97
  98 }
  99
 100
 101 double FGRadioTransmission::receiveNav(SGGeod tx_pos, double freq, int transmission_type) {
 102
 103         // typical VOR/LOC transmitter power appears to be 100 - 200 Watt i.e 50 - 53 dBm
 104         // vor/loc typical sensitivity between -107 and -101 dBm
 105         // glideslope sensitivity between -85 and -81 dBm
 106         if ( _propagation_model == 1) {
 107                 return LOS_calculate_attenuation(tx_pos, freq, 1);
 108         }
 109         else if ( _propagation_model == 2) {
 110                 return ITM_calculate_attenuation(tx_pos, freq, 1);
 111         }
 112
 113         return -1;
 114
 115 }
 116
 117
 118 double FGRadioTransmission::receiveBeacon(SGGeod &tx_pos, double heading, double pitch) {
 119
 120         // these properties should be set by an instrument
 121         _receiver_sensitivity = _root_node->getDoubleValue("station[0]/rx-sensitivity", _receiver_sensitivity);
 122         _transmitter_power = watt_to_dbm(_root_node->getDoubleValue("station[0]/tx-power-watt", _transmitter_power));
 123         _polarization = _root_node->getIntValue("station[0]/polarization", 1);
 124         _tx_antenna_height += _root_node->getDoubleValue("station[0]/tx-antenna-height", 0);
 125         _rx_antenna_height += _root_node->getDoubleValue("station[0]/rx-antenna-height", 0);
 126         _tx_antenna_gain += _root_node->getDoubleValue("station[0]/tx-antenna-gain", 0);
 127         _rx_antenna_gain += _root_node->getDoubleValue("station[0]/rx-antenna-gain", 0);
 128
 129         double freq = _root_node->getDoubleValue("station[0]/frequency", 144.8);        // by default stay in the ham 2 meter band
 130
 131         double comm1 = getFrequency(1);
 132         double comm2 = getFrequency(2);
 133         if ( !(fabs(freq - comm1) <= 0.0001) &&  !(fabs(freq - comm2) <= 0.0001) ) {
 134                 return -1;
 135         }
 136
 137         double signal = ITM_calculate_attenuation(tx_pos, freq, 1);
 138
 139         return signal;
 140 }
 141
 142
 143
 144 void FGRadioTransmission::receiveATC(SGGeod tx_pos, double freq, string text, int ground_to_air) {
 145
 146         // adjust some default parameters in case the ATC code does not set them
 147         if(ground_to_air == 1) {
 148                 _transmitter_power += 4.0;
 149                 _tx_antenna_height += 30.0;
 150                 _tx_antenna_gain += 2.0;
 151         }
 152
 153         double comm1 = getFrequency(1);
 154         double comm2 = getFrequency(2);
 155         if ( !(fabs(freq - comm1) <= 0.0001) &&  !(fabs(freq - comm2) <= 0.0001) ) {
 156                 return;
 157         }
 158         else {
 159
 160                 if ( _propagation_model == 0) {         // skip propagation routines entirely
 161                         fgSetString("/sim/messages/atc", text.c_str());
 162                 }
 163                 else if ( _propagation_model == 1 ) {           // Use free-space, round earth
 164
 165                         double signal = LOS_calculate_attenuation(tx_pos, freq, ground_to_air);
 166                         if (signal <= 0.0) {
 167                                 return;
 168                         }
 169                         else {
 170                                 fgSetString("/sim/messages/atc", text.c_str());
 171                         }
 172                 }
 173                 else if ( _propagation_model == 2 ) {   // Use ITM propagation model
 174
 175                         double signal = ITM_calculate_attenuation(tx_pos, freq, ground_to_air);
 176                         if (signal <= 0.0) {
 177                                 return;
 178                         }
 179                         if ((signal > 0.0) && (signal < 12.0)) {
 180                                 /** for low SNR values need a way to make the conversation
 181                                 *       hard to understand but audible
 182                                 *       in the real world, the receiver AGC fails to capture the slope
 183                                 *       and the signal, due to being amplitude modulated, decreases volume after demodulation
 184                                 *       the workaround below is more akin to what would happen on a FM transmission
 185                                 *       therefore the correct way would be to work on the volume
 186                                 **/
 187                                 /*
 188                                 string hash_noise = " ";
 189                                 int reps = (int) (fabs(floor(signal - 11.0)) * 2);
 190                                 int t_size = text.size();
 191                                 for (int n = 1; n <= reps; ++n) {
 192                                         int pos = rand() % (t_size -1);
 193                                         text.replace(pos,1, hash_noise);
 194                                 }
 195                                 */
 196                                 //double volume = (fabs(signal - 12.0) / 12);
 197                                 //double old_volume = fgGetDouble("/sim/sound/voices/voice/volume");
 198
 199                                 //fgSetDouble("/sim/sound/voices/voice/volume", volume);
 200                                 fgSetString("/sim/messages/atc", text.c_str());
 201                                 //fgSetDouble("/sim/sound/voices/voice/volume", old_volume);
 202                         }
 203                         else {
 204                                 fgSetString("/sim/messages/atc", text.c_str());
 205                         }
 206                 }
 207         }
 208 }
 209
 210
 211 double FGRadioTransmission::ITM_calculate_attenuation(SGGeod pos, double freq, int transmission_type) {
 212
 213
 214         if((freq < 40.0) || (freq > 20000.0))   // frequency out of recommended range
 215                 return -1;
 216         /** ITM default parameters
 217                 TODO: take them from tile materials (especially for sea)?
 218         **/
 219         double eps_dielect=15.0;
 220         double sgm_conductivity = 0.005;
 221         double eno = 301.0;
 222         double frq_mhz = freq;
 223
 224         int radio_climate = 5;          // continental temperate
 225         int pol= _polarization;
 226         double conf = 0.90;     // 90% of situations and time, take into account speed
 227         double rel = 0.90;
 228         double dbloss;
 229         char strmode[150];
 230         int p_mode = 0; // propgation mode selector: 0 LOS, 1 diffraction dominant, 2 troposcatter
 231         double horizons[2];
 232         int errnum;
 233
 234         double clutter_loss = 0.0;      // loss due to vegetation and urban
 235         double tx_pow = _transmitter_power;
 236         double ant_gain = _rx_antenna_gain + _tx_antenna_gain;
 237         double signal = 0.0;
 238
 239
 240         double link_budget = tx_pow - _receiver_sensitivity - _rx_line_losses - _tx_line_losses + ant_gain;
 241         double signal_strength = tx_pow - _rx_line_losses - _tx_line_losses + ant_gain;
 242         double tx_erp = dbm_to_watt(tx_pow + _tx_antenna_gain - _tx_line_losses);
 243
 244
 245         FGScenery * scenery = globals->get_scenery();
 246
 247         double own_lat = fgGetDouble("/position/latitude-deg");
 248         double own_lon = fgGetDouble("/position/longitude-deg");
 249         double own_alt_ft = fgGetDouble("/position/altitude-ft");
 250         double own_heading = fgGetDouble("/orientation/heading-deg");
 251         double own_alt= own_alt_ft * SG_FEET_TO_METER;
 252
 253
 254
 255
 256         SGGeod own_pos = SGGeod::fromDegM( own_lon, own_lat, own_alt );
 257         SGGeod max_own_pos = SGGeod::fromDegM( own_lon, own_lat, SG_MAX_ELEVATION_M );
 258         SGGeoc center = SGGeoc::fromGeod( max_own_pos );
 259         SGGeoc own_pos_c = SGGeoc::fromGeod( own_pos );
 260
 261
 262         double sender_alt_ft,sender_alt;
 263         double transmitter_height=0.0;
 264         double receiver_height=0.0;
 265         SGGeod sender_pos = pos;
 266
 267         sender_alt_ft = sender_pos.getElevationFt();
 268         sender_alt = sender_alt_ft * SG_FEET_TO_METER;
 269         SGGeod max_sender_pos = SGGeod::fromGeodM( pos, SG_MAX_ELEVATION_M );
 270         SGGeoc sender_pos_c = SGGeoc::fromGeod( sender_pos );
 271
 272
 273         double point_distance= _terrain_sampling_distance;
 274         double course = SGGeodesy::courseRad(own_pos_c, sender_pos_c);
 275         double reverse_course = SGGeodesy::courseRad(sender_pos_c, own_pos_c);
 276         double distance_m = SGGeodesy::distanceM(own_pos, sender_pos);
 277         double probe_distance = 0.0;
 278         /** If distance larger than this value (300 km), assume reception imposssible to spare CPU cycles */
 279         if (distance_m > 300000)
 280                 return -1.0;
 281         /** If above 8000 meters, consider LOS mode and calculate free-space att to spare CPU cycles */
 282         if (own_alt > 8000) {
 283                 dbloss = 20 * log10(distance_m) +20 * log10(frq_mhz) -27.55;
 284                 SG_LOG(SG_GENERAL, SG_BULK,
 285                         "ITM Free-space mode:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm, free-space attenuation");
 286                 //cerr << "ITM Free-space mode:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm, free-space attenuation" << endl;
 287                 signal = link_budget - dbloss;
 288                 return signal;
 289         }
 290
 291
 292         int max_points = (int)floor(distance_m / point_distance);
 293         //double delta_last = fmod(distance_m, point_distance);
 294
 295         deque<double> elevations;
 296         deque<string*> materials;
 297
 298
 299         double elevation_under_pilot = 0.0;
 300         if (scenery->get_elevation_m( max_own_pos, elevation_under_pilot, NULL )) {
 301                 receiver_height = own_alt - elevation_under_pilot;
 302         }
 303
 304         double elevation_under_sender = 0.0;
 305         if (scenery->get_elevation_m( max_sender_pos, elevation_under_sender, NULL )) {
 306                 transmitter_height = sender_alt - elevation_under_sender;
 307         }
 308         else {
 309                 transmitter_height = sender_alt;
 310         }
 311
 312
 313         transmitter_height += _tx_antenna_height;
 314         receiver_height += _rx_antenna_height;
 315
 316         //cerr << "ITM:: RX-height: " << receiver_height << " meters, TX-height: " << transmitter_height << " meters, Distance: " << distance_m << " meters" << endl;
 317         _root_node->setDoubleValue("station[0]/rx-height", receiver_height);
 318         _root_node->setDoubleValue("station[0]/tx-height", transmitter_height);
 319         _root_node->setDoubleValue("station[0]/distance", distance_m / 1000);
 320
 321         unsigned int e_size = (deque<unsigned>::size_type)max_points;
 322
 323         while (elevations.size() <= e_size) {
 324                 probe_distance += point_distance;
 325                 SGGeod probe = SGGeod::fromGeoc(center.advanceRadM( course, probe_distance ));
 326                 const SGMaterial *mat = 0;
 327                 double elevation_m = 0.0;
 328
 329                 if (scenery->get_elevation_m( probe, elevation_m, &mat )) {
 330                         if((transmission_type == 3) || (transmission_type == 4)) {
 331                                 elevations.push_back(elevation_m);
 332                                 if(mat) {
 333                                         const std::vector<string> mat_names = mat->get_names();
 334                                         string* name = new string(mat_names[0]);
 335                                         materials.push_back(name);
 336                                 }
 337                                 else {
 338                                         string* no_material = new string("None");
 339                                         materials.push_back(no_material);
 340                                 }
 341                         }
 342                         else {
 343                                  elevations.push_front(elevation_m);
 344                                  if(mat) {
 345                                          const std::vector<string> mat_names = mat->get_names();
 346                                          string* name = new string(mat_names[0]);
 347                                          materials.push_front(name);
 348                                 }
 349                                 else {
 350                                         string* no_material = new string("None");
 351                                         materials.push_front(no_material);
 352                                 }
 353                         }
 354                 }
 355                 else {
 356                         if((transmission_type == 3) || (transmission_type == 4)) {
 357                                 elevations.push_back(0.0);
 358                                 string* no_material = new string("None");
 359                                 materials.push_back(no_material);
 360                         }
 361                         else {
 362                                 string* no_material = new string("None");
 363                                 elevations.push_front(0.0);
 364                                 materials.push_front(no_material);
 365                         }
 366                 }
 367         }
 368         if((transmission_type == 3) || (transmission_type == 4)) {
 369                 elevations.push_front(elevation_under_pilot);
 370                 //if (delta_last > (point_distance / 2) )                       // only add last point if it's farther than half point_distance
 371                         elevations.push_back(elevation_under_sender);
 372         }
 373         else {
 374                 elevations.push_back(elevation_under_pilot);
 375                 //if (delta_last > (point_distance / 2) )
 376                         elevations.push_front(elevation_under_sender);
 377         }
 378
 379
 380         double num_points= (double)elevations.size();
 381
 382
 383         elevations.push_front(point_distance);
 384         elevations.push_front(num_points -1);
 385
 386         int size = elevations.size();
 387         boost::scoped_array<double> itm_elev( new double[size] );
 388
 389         for(int i=0;i<size;i++) {
 390                 itm_elev[i]=elevations[i];
 391         }
 392
 393         if((transmission_type == 3) || (transmission_type == 4)) {
 394                 // the sender and receiver roles are switched
 395                 ITM::point_to_point(itm_elev.get(), receiver_height, transmitter_height,
 396                         eps_dielect, sgm_conductivity, eno, frq_mhz, radio_climate,
 397                         pol, conf, rel, dbloss, strmode, p_mode, horizons, errnum);
 398                 if( _root_node->getBoolValue( "use-clutter-attenuation", false ) )
 399                         calculate_clutter_loss(frq_mhz, itm_elev.get(), materials, receiver_height, transmitter_height, p_mode, horizons, clutter_loss);
 400         }
 401         else {
 402                 ITM::point_to_point(itm_elev.get(), transmitter_height, receiver_height,
 403                         eps_dielect, sgm_conductivity, eno, frq_mhz, radio_climate,
 404                         pol, conf, rel, dbloss, strmode, p_mode, horizons, errnum);
 405                 if( _root_node->getBoolValue( "use-clutter-attenuation", false ) )
 406                         calculate_clutter_loss(frq_mhz, itm_elev.get(), materials, transmitter_height, receiver_height, p_mode, horizons, clutter_loss);
 407         }
 408
 409         double pol_loss = 0.0;
 410         // TODO: remove this check after we check a bit the axis calculations in this function
 411         if (_polarization == 1) {
 412                 pol_loss = polarization_loss();
 413         }
 414         //SG_LOG(SG_GENERAL, SG_BULK,
 415         //              "ITM:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm, " << strmode << ", Error: " << errnum);
 416         //cerr << "ITM:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm, " << strmode << ", Error: " << errnum << endl;
 417         _root_node->setDoubleValue("station[0]/link-budget", link_budget);
 418         _root_node->setDoubleValue("station[0]/terrain-attenuation", dbloss);
 419         _root_node->setStringValue("station[0]/prop-mode", strmode);
 420         _root_node->setDoubleValue("station[0]/clutter-attenuation", clutter_loss);
 421         _root_node->setDoubleValue("station[0]/polarization-attenuation", pol_loss);
 422         //if (errnum == 4)      // if parameters are outside sane values for lrprop, bail out fast
 423         //      return -1;
 424
 425         // temporary, keep this antenna radiation pattern code here
 426         double tx_pattern_gain = 0.0;
 427         double rx_pattern_gain = 0.0;
 428         double sender_heading = 270.0; // due West
 429         double tx_antenna_bearing = sender_heading - reverse_course * SGD_RADIANS_TO_DEGREES;
 430         double rx_antenna_bearing = own_heading - course * SGD_RADIANS_TO_DEGREES;
 431         double rx_elev_angle = atan((itm_elev[2] + transmitter_height - itm_elev[(int)itm_elev[0] + 2] + receiver_height) / distance_m) * SGD_RADIANS_TO_DEGREES;
 432         double tx_elev_angle = 0.0 - rx_elev_angle;
 433         if (_root_node->getBoolValue("use-tx-antenna-pattern", false)) {
 434                 FGRadioAntenna* TX_antenna;
 435                 TX_antenna = new FGRadioAntenna("Plot2");
 436                 TX_antenna->set_heading(sender_heading);
 437                 TX_antenna->set_elevation_angle(0);
 438                 tx_pattern_gain = TX_antenna->calculate_gain(tx_antenna_bearing, tx_elev_angle);
 439                 delete TX_antenna;
 440         }
 441         if (_root_node->getBoolValue("use-rx-antenna-pattern", false)) {
 442                 FGRadioAntenna* RX_antenna;
 443                 RX_antenna = new FGRadioAntenna("Plot2");
 444                 RX_antenna->set_heading(own_heading);
 445                 RX_antenna->set_elevation_angle(fgGetDouble("/orientation/pitch-deg"));
 446                 rx_pattern_gain = RX_antenna->calculate_gain(rx_antenna_bearing, rx_elev_angle);
 447                 delete RX_antenna;
 448         }
 449
 450         signal = link_budget - dbloss - clutter_loss + pol_loss + rx_pattern_gain + tx_pattern_gain;
 451         double signal_strength_dbm = signal_strength - dbloss - clutter_loss + pol_loss + rx_pattern_gain + tx_pattern_gain;
 452         double field_strength_uV = dbm_to_microvolt(signal_strength_dbm);
 453         _root_node->setDoubleValue("station[0]/signal-dbm", signal_strength_dbm);
 454         _root_node->setDoubleValue("station[0]/field-strength-uV", field_strength_uV);
 455         _root_node->setDoubleValue("station[0]/signal", signal);
 456         _root_node->setDoubleValue("station[0]/tx-erp", tx_erp);
 457
 458         //_root_node->setDoubleValue("station[0]/tx-pattern-gain", tx_pattern_gain);
 459         //_root_node->setDoubleValue("station[0]/rx-pattern-gain", rx_pattern_gain);
 460
 461         for (unsigned i =0; i < materials.size(); i++) {
 462                 delete materials[i];
 463         }
 464
 465         return signal;
 466
 467 }
 468
 469
 470 void FGRadioTransmission::calculate_clutter_loss(double freq, double itm_elev[], deque<string*> &materials,
 471         double transmitter_height, double receiver_height, int p_mode,
 472         double horizons[], double &clutter_loss) {
 473
 474         double distance_m = itm_elev[0] * itm_elev[1]; // only consider elevation points
 475         unsigned mat_size = materials.size();
 476         if (p_mode == 0) {      // LOS: take each point and see how clutter height affects first Fresnel zone
 477                 int mat = 0;
 478                 int j=1;
 479                 for (int k=3;k < (int)(itm_elev[0]) + 2;k++) {
 480
 481                         double clutter_height = 0.0;    // mean clutter height for a certain terrain type
 482                         double clutter_density = 0.0;   // percent of reflected wave
 483                         if((unsigned)mat >= mat_size) { //this tends to happen when the model interferes with the antenna (obstructs)
 484                                 //cerr << "Array index out of bounds 0-0: " << mat << " size: " << mat_size << endl;
 485                                 break;
 486                         }
 487                         get_material_properties(materials[mat], clutter_height, clutter_density);
 488
 489                         double grad = fabs(itm_elev[2] + transmitter_height - itm_elev[(int)itm_elev[0] + 2] + receiver_height) / distance_m;
 490                         // First Fresnel radius
 491                         double frs_rad = 548 * sqrt( (j * itm_elev[1] * (itm_elev[0] - j) * itm_elev[1] / 1000000) / (  distance_m * freq / 1000) );
 492                         if (frs_rad <= 0.0) {   //this tends to happen when the model interferes with the antenna (obstructs)
 493                                 //cerr << "Frs rad 0-0: " << frs_rad << endl;
 494                                 continue;
 495                         }
 496                         //double earth_h = distance_m * (distance_m - j * itm_elev[1]) / ( 1000000 * 12.75 * 1.33 );    // K=4/3
 497
 498                         double min_elev = SGMiscd::min(itm_elev[2] + transmitter_height, itm_elev[(int)itm_elev[0] + 2] + receiver_height);
 499                         double d1 = j * itm_elev[1];
 500                         if ((itm_elev[2] + transmitter_height) > ( itm_elev[(int)itm_elev[0] + 2] + receiver_height) ) {
 501                                 d1 = (itm_elev[0] - j) * itm_elev[1];
 502                         }
 503                         double ray_height = (grad * d1) + min_elev;
 504
 505                         double clearance = ray_height - (itm_elev[k] + clutter_height) - frs_rad * 8/10;
 506                         double intrusion = fabs(clearance);
 507
 508                         if (clearance >= 0) {
 509                                 // no losses
 510                         }
 511                         else if (clearance < 0 && (intrusion < clutter_height)) {
 512
 513                                 clutter_loss += clutter_density * (intrusion / (frs_rad * 2) ) * (freq/100) * (itm_elev[1]/100);
 514                         }
 515                         else if (clearance < 0 && (intrusion > clutter_height)) {
 516                                 clutter_loss += clutter_density * (clutter_height / (frs_rad * 2 ) ) * (freq/100) * (itm_elev[1]/100);
 517                         }
 518                         else {
 519                                 // no losses
 520                         }
 521                         j++;
 522                         mat++;
 523                 }
 524
 525         }
 526         else if (p_mode == 1) {         // diffraction
 527
 528                 if (horizons[1] == 0.0) {       //      single horizon: same as above, except pass twice using the highest point
 529                         int num_points_1st = (int)floor( horizons[0] * itm_elev[0]/ distance_m );
 530                         int num_points_2nd = (int)ceil( (distance_m - horizons[0]) * itm_elev[0] / distance_m );
 531                         //cerr << "Diffraction 1 horizon:: points1: " << num_points_1st << " points2: " << num_points_2nd << endl;
 532                         int last = 1;
 533                         /** perform the first pass */
 534                         int mat = 0;
 535                         int j=1;
 536                         for (int k=3;k < num_points_1st + 2;k++) {
 537                                 if (num_points_1st < 1)
 538                                         break;
 539                                 double clutter_height = 0.0;    // mean clutter height for a certain terrain type
 540                                 double clutter_density = 0.0;   // percent of reflected wave
 541
 542                                 if((unsigned)mat >= mat_size) {
 543                                         //cerr << "Array index out of bounds 1-1: " << mat << " size: " << mat_size << endl;
 544                                         break;
 545                                 }
 546                                 get_material_properties(materials[mat], clutter_height, clutter_density);
 547
 548                                 double grad = fabs(itm_elev[2] + transmitter_height - itm_elev[num_points_1st + 2] + clutter_height) / distance_m;
 549                                 // First Fresnel radius
 550                                 double frs_rad = 548 * sqrt( (j * itm_elev[1] * (num_points_1st - j) * itm_elev[1] / 1000000) / ( num_points_1st * itm_elev[1] * freq / 1000) );
 551                                 if (frs_rad <= 0.0) {
 552                                         //cerr << "Frs rad 1-1: " << frs_rad << endl;
 553                                         continue;
 554                                 }
 555                                 //double earth_h = distance_m * (distance_m - j * itm_elev[1]) / ( 1000000 * 12.75 * 1.33 );    // K=4/3
 556
 557                                 double min_elev = SGMiscd::min(itm_elev[2] + transmitter_height, itm_elev[num_points_1st + 2] + clutter_height);
 558                                 double d1 = j * itm_elev[1];
 559                                 if ( (itm_elev[2] + transmitter_height) > (itm_elev[num_points_1st + 2] + clutter_height) ) {
 560                                         d1 = (num_points_1st - j) * itm_elev[1];
 561                                 }
 562                                 double ray_height = (grad * d1) + min_elev;
 563
 564                                 double clearance = ray_height - (itm_elev[k] + clutter_height) - frs_rad * 8/10;
 565                                 double intrusion = fabs(clearance);
 566
 567                                 if (clearance >= 0) {
 568                                         // no losses
 569                                 }
 570                                 else if (clearance < 0 && (intrusion < clutter_height)) {
 571
 572                                         clutter_loss += clutter_density * (intrusion / (frs_rad * 2) ) * (freq/100) * (itm_elev[1]/100);
 573                                 }
 574                                 else if (clearance < 0 && (intrusion > clutter_height)) {
 575                                         clutter_loss += clutter_density * (clutter_height / (frs_rad * 2 ) ) * (freq/100) * (itm_elev[1]/100);
 576                                 }
 577                                 else {
 578                                         // no losses
 579                                 }
 580                                 j++;
 581                                 mat++;
 582                                 last = k;
 583                         }
 584
 585                         /** and the second pass */
 586                         mat +=1;
 587                         j =1; // first point is diffraction edge, 2nd the RX elevation
 588                         for (int k=last+2;k < (int)(itm_elev[0]) + 2;k++) {
 589                                 if (num_points_2nd < 1)
 590                                         break;
 591                                 double clutter_height = 0.0;    // mean clutter height for a certain terrain type
 592                                 double clutter_density = 0.0;   // percent of reflected wave
 593
 594                                 if((unsigned)mat >= mat_size) {
 595                                         //cerr << "Array index out of bounds 1-2: " << mat << " size: " << mat_size << endl;
 596                                         break;
 597                                 }
 598                                 get_material_properties(materials[mat], clutter_height, clutter_density);
 599
 600                                 double grad = fabs(itm_elev[last+1] + clutter_height - itm_elev[(int)itm_elev[0] + 2] + receiver_height) / distance_m;
 601                                 // First Fresnel radius
 602                                 double frs_rad = 548 * sqrt( (j * itm_elev[1] * (num_points_2nd - j) * itm_elev[1] / 1000000) / (  num_points_2nd * itm_elev[1] * freq / 1000) );
 603                                 if (frs_rad <= 0.0) {
 604                                         //cerr << "Frs rad 1-2: " << frs_rad << " numpoints2 " << num_points_2nd << " j: " << j << endl;
 605                                         continue;
 606                                 }
 607                                 //double earth_h = distance_m * (distance_m - j * itm_elev[1]) / ( 1000000 * 12.75 * 1.33 );    // K=4/3
 608
 609                                 double min_elev = SGMiscd::min(itm_elev[last+1] + clutter_height, itm_elev[(int)itm_elev[0] + 2] + receiver_height);
 610                                 double d1 = j * itm_elev[1];
 611                                 if ( (itm_elev[last+1] + clutter_height) > (itm_elev[(int)itm_elev[0] + 2] + receiver_height) ) {
 612                                         d1 = (num_points_2nd - j) * itm_elev[1];
 613                                 }
 614                                 double ray_height = (grad * d1) + min_elev;
 615
 616                                 double clearance = ray_height - (itm_elev[k] + clutter_height) - frs_rad * 8/10;
 617                                 double intrusion = fabs(clearance);
 618
 619                                 if (clearance >= 0) {
 620                                         // no losses
 621                                 }
 622                                 else if (clearance < 0 && (intrusion < clutter_height)) {
 623
 624                                         clutter_loss += clutter_density * (intrusion / (frs_rad * 2) ) * (freq/100) * (itm_elev[1]/100);
 625                                 }
 626                                 else if (clearance < 0 && (intrusion > clutter_height)) {
 627                                         clutter_loss += clutter_density * (clutter_height / (frs_rad * 2 ) ) * (freq/100) * (itm_elev[1]/100);
 628                                 }
 629                                 else {
 630                                         // no losses
 631                                 }
 632                                 j++;
 633                                 mat++;
 634                         }
 635
 636                 }
 637                 else {  // double horizon: same as single horizon, except there are 3 segments
 638
 639                         int num_points_1st = (int)floor( horizons[0] * itm_elev[0] / distance_m );
 640                         int num_points_2nd = (int)floor(horizons[1] * itm_elev[0] / distance_m );
 641                         int num_points_3rd = (int)itm_elev[0] - num_points_1st - num_points_2nd;
 642                         //cerr << "Double horizon:: horizon1: " << horizons[0] << " horizon2: " << horizons[1] << " distance: " << distance_m << endl;
 643                         //cerr << "Double horizon:: points1: " << num_points_1st << " points2: " << num_points_2nd << " points3: " << num_points_3rd << endl;
 644                         int last = 1;
 645                         /** perform the first pass */
 646                         int mat = 0;
 647                         int j=1; // first point is TX elevation, 2nd is obstruction elevation
 648                         for (int k=3;k < num_points_1st +2;k++) {
 649                                 if (num_points_1st < 1)
 650                                         break;
 651                                 double clutter_height = 0.0;    // mean clutter height for a certain terrain type
 652                                 double clutter_density = 0.0;   // percent of reflected wave
 653                                 if((unsigned)mat >= mat_size) {
 654                                         //cerr << "Array index out of bounds 2-1: " << mat << " size: " << mat_size << endl;
 655                                         break;
 656                                 }
 657                                 get_material_properties(materials[mat], clutter_height, clutter_density);
 658
 659                                 double grad = fabs(itm_elev[2] + transmitter_height - itm_elev[num_points_1st + 2] + clutter_height) / distance_m;
 660                                 // First Fresnel radius
 661                                 double frs_rad = 548 * sqrt( (j * itm_elev[1] * (num_points_1st - j) * itm_elev[1] / 1000000) / (  num_points_1st * itm_elev[1] * freq / 1000) );
 662                                 if (frs_rad <= 0.0) {
 663                                         //cerr << "Frs rad 2-1: " << frs_rad << " numpoints1 " << num_points_1st << " j: " << j << endl;
 664                                         continue;
 665                                 }
 666                                 //double earth_h = distance_m * (distance_m - j * itm_elev[1]) / ( 1000000 * 12.75 * 1.33 );    // K=4/3
 667
 668                                 double min_elev = SGMiscd::min(itm_elev[2] + transmitter_height, itm_elev[num_points_1st + 2] + clutter_height);
 669                                 double d1 = j * itm_elev[1];
 670                                 if ( (itm_elev[2] + transmitter_height) > (itm_elev[num_points_1st + 2] + clutter_height) ) {
 671                                         d1 = (num_points_1st - j) * itm_elev[1];
 672                                 }
 673                                 double ray_height = (grad * d1) + min_elev;
 674
 675                                 double clearance = ray_height - (itm_elev[k] + clutter_height) - frs_rad * 8/10;
 676                                 double intrusion = fabs(clearance);
 677
 678                                 if (clearance >= 0) {
 679                                         // no losses
 680                                 }
 681                                 else if (clearance < 0 && (intrusion < clutter_height)) {
 682
 683                                         clutter_loss += clutter_density * (intrusion / (frs_rad * 2) ) * (freq/100) * (itm_elev[1]/100);
 684                                 }
 685                                 else if (clearance < 0 && (intrusion > clutter_height)) {
 686                                         clutter_loss += clutter_density * (clutter_height / (frs_rad * 2 ) ) * (freq/100) * (itm_elev[1]/100);
 687                                 }
 688                                 else {
 689                                         // no losses
 690                                 }
 691                                 j++;
 692                                 mat++;
 693                                 last = k;
 694                         }
 695                         mat +=1;
 696                         /** and the second pass */
 697                         int last2=1;
 698                         j =1; // first point is 1st obstruction elevation, 2nd is 2nd obstruction elevation
 699                         for (int k=last+2;k < num_points_1st + num_points_2nd +2;k++) {
 700                                 if (num_points_2nd < 1)
 701                                         break;
 702                                 double clutter_height = 0.0;    // mean clutter height for a certain terrain type
 703                                 double clutter_density = 0.0;   // percent of reflected wave
 704                                 if((unsigned)mat >= mat_size) {
 705                                         //cerr << "Array index out of bounds 2-2: " << mat << " size: " << mat_size << endl;
 706                                         break;
 707                                 }
 708                                 get_material_properties(materials[mat], clutter_height, clutter_density);
 709
 710                                 double grad = fabs(itm_elev[last+1] + clutter_height - itm_elev[num_points_1st + num_points_2nd + 2] + clutter_height) / distance_m;
 711                                 // First Fresnel radius
 712                                 double frs_rad = 548 * sqrt( (j * itm_elev[1] * (num_points_2nd - j) * itm_elev[1] / 1000000) / (  num_points_2nd * itm_elev[1] * freq / 1000) );
 713                                 if (frs_rad <= 0.0) {
 714                                         //cerr << "Frs rad 2-2: " << frs_rad << " numpoints2 " << num_points_2nd << " j: " << j << endl;
 715                                         continue;
 716                                 }
 717                                 //double earth_h = distance_m * (distance_m - j * itm_elev[1]) / ( 1000000 * 12.75 * 1.33 );    // K=4/3
 718
 719                                 double min_elev = SGMiscd::min(itm_elev[last+1] + clutter_height, itm_elev[num_points_1st + num_points_2nd +2] + clutter_height);
 720                                 double d1 = j * itm_elev[1];
 721                                 if ( (itm_elev[last+1] + clutter_height) > (itm_elev[num_points_1st + num_points_2nd + 2] + clutter_height) ) {
 722                                         d1 = (num_points_2nd - j) * itm_elev[1];
 723                                 }
 724                                 double ray_height = (grad * d1) + min_elev;
 725
 726                                 double clearance = ray_height - (itm_elev[k] + clutter_height) - frs_rad * 8/10;
 727                                 double intrusion = fabs(clearance);
 728
 729                                 if (clearance >= 0) {
 730                                         // no losses
 731                                 }
 732                                 else if (clearance < 0 && (intrusion < clutter_height)) {
 733
 734                                         clutter_loss += clutter_density * (intrusion / (frs_rad * 2) ) * (freq/100) * (itm_elev[1]/100);
 735                                 }
 736                                 else if (clearance < 0 && (intrusion > clutter_height)) {
 737                                         clutter_loss += clutter_density * (clutter_height / (frs_rad * 2 ) ) * (freq/100) * (itm_elev[1]/100);
 738                                 }
 739                                 else {
 740                                         // no losses
 741                                 }
 742                                 j++;
 743                                 mat++;
 744                                 last2 = k;
 745                         }
 746
 747                         /** third and final pass */
 748                         mat +=1;
 749                         j =1; // first point is 2nd obstruction elevation, 3rd is RX elevation
 750                         for (int k=last2+2;k < (int)itm_elev[0] + 2;k++) {
 751                                 if (num_points_3rd < 1)
 752                                         break;
 753                                 double clutter_height = 0.0;    // mean clutter height for a certain terrain type
 754                                 double clutter_density = 0.0;   // percent of reflected wave
 755                                 if((unsigned)mat >= mat_size) {
 756                                         //cerr << "Array index out of bounds 2-3: " << mat << " size: " << mat_size << endl;
 757                                         break;
 758                                 }
 759                                 get_material_properties(materials[mat], clutter_height, clutter_density);
 760
 761                                 double grad = fabs(itm_elev[last2+1] + clutter_height - itm_elev[(int)itm_elev[0] + 2] + receiver_height) / distance_m;
 762                                 // First Fresnel radius
 763                                 double frs_rad = 548 * sqrt( (j * itm_elev[1] * (num_points_3rd - j) * itm_elev[1] / 1000000) / (  num_points_3rd * itm_elev[1] * freq / 1000) );
 764                                 if (frs_rad <= 0.0) {
 765                                         //cerr << "Frs rad 2-3: " << frs_rad << " numpoints3 " << num_points_3rd << " j: " << j << endl;
 766                                         continue;
 767                                 }
 768
 769                                 //double earth_h = distance_m * (distance_m - j * itm_elev[1]) / ( 1000000 * 12.75 * 1.33 );    // K=4/3
 770
 771                                 double min_elev = SGMiscd::min(itm_elev[last2+1] + clutter_height, itm_elev[(int)itm_elev[0] + 2] + receiver_height);
 772                                 double d1 = j * itm_elev[1];
 773                                 if ( (itm_elev[last2+1] + clutter_height) > (itm_elev[(int)itm_elev[0] + 2] + receiver_height) ) {
 774                                         d1 = (num_points_3rd - j) * itm_elev[1];
 775                                 }
 776                                 double ray_height = (grad * d1) + min_elev;
 777
 778                                 double clearance = ray_height - (itm_elev[k] + clutter_height) - frs_rad * 8/10;
 779                                 double intrusion = fabs(clearance);
 780
 781                                 if (clearance >= 0) {
 782                                         // no losses
 783                                 }
 784                                 else if (clearance < 0 && (intrusion < clutter_height)) {
 785
 786                                         clutter_loss += clutter_density * (intrusion / (frs_rad * 2) ) * (freq/100) * (itm_elev[1]/100);
 787                                 }
 788                                 else if (clearance < 0 && (intrusion > clutter_height)) {
 789                                         clutter_loss += clutter_density * (clutter_height / (frs_rad * 2 ) ) * (freq/100) * (itm_elev[1]/100);
 790                                 }
 791                                 else {
 792                                         // no losses
 793                                 }
 794                                 j++;
 795                                 mat++;
 796
 797                         }
 798
 799                 }
 800         }
 801         else if (p_mode == 2) {         //      troposcatter: ignore ground clutter for now... maybe do something with weather
 802                 clutter_loss = 0.0;
 803         }
 804
 805 }
 806
 807
 808 void FGRadioTransmission::get_material_properties(string* mat_name, double &height, double &density) {
 809
 810         if(!mat_name)
 811                 return;
 812
 813         if(*mat_name == "Landmass") {
 814                 height = 15.0;
 815                 density = 0.2;
 816         }
 817
 818         else if(*mat_name == "SomeSort") {
 819                 height = 15.0;
 820                 density = 0.2;
 821         }
 822
 823         else if(*mat_name == "Island") {
 824                 height = 15.0;
 825                 density = 0.2;
 826         }
 827         else if(*mat_name == "Default") {
 828                 height = 15.0;
 829                 density = 0.2;
 830         }
 831         else if(*mat_name == "EvergreenBroadCover") {
 832                 height = 20.0;
 833                 density = 0.2;
 834         }
 835         else if(*mat_name == "EvergreenForest") {
 836                 height = 20.0;
 837                 density = 0.2;
 838         }
 839         else if(*mat_name == "DeciduousBroadCover") {
 840                 height = 15.0;
 841                 density = 0.3;
 842         }
 843         else if(*mat_name == "DeciduousForest") {
 844                 height = 15.0;
 845                 density = 0.3;
 846         }
 847         else if(*mat_name == "MixedForestCover") {
 848                 height = 20.0;
 849                 density = 0.25;
 850         }
 851         else if(*mat_name == "MixedForest") {
 852                 height = 15.0;
 853                 density = 0.25;
 854         }
 855         else if(*mat_name == "RainForest") {
 856                 height = 25.0;
 857                 density = 0.55;
 858         }
 859         else if(*mat_name == "EvergreenNeedleCover") {
 860                 height = 15.0;
 861                 density = 0.2;
 862         }
 863         else if(*mat_name == "WoodedTundraCover") {
 864                 height = 5.0;
 865                 density = 0.15;
 866         }
 867         else if(*mat_name == "DeciduousNeedleCover") {
 868                 height = 5.0;
 869                 density = 0.2;
 870         }
 871         else if(*mat_name == "ScrubCover") {
 872                 height = 3.0;
 873                 density = 0.15;
 874         }
 875         else if(*mat_name == "BuiltUpCover") {
 876                 height = 30.0;
 877                 density = 0.7;
 878         }
 879         else if(*mat_name == "Urban") {
 880                 height = 30.0;
 881                 density = 0.7;
 882         }
 883         else if(*mat_name == "Construction") {
 884                 height = 30.0;
 885                 density = 0.7;
 886         }
 887         else if(*mat_name == "Industrial") {
 888                 height = 30.0;
 889                 density = 0.7;
 890         }
 891         else if(*mat_name == "Port") {
 892                 height = 30.0;
 893                 density = 0.7;
 894         }
 895         else if(*mat_name == "Town") {
 896                 height = 10.0;
 897                 density = 0.5;
 898         }
 899         else if(*mat_name == "SubUrban") {
 900                 height = 10.0;
 901                 density = 0.5;
 902         }
 903         else if(*mat_name == "CropWoodCover") {
 904                 height = 10.0;
 905                 density = 0.1;
 906         }
 907         else if(*mat_name == "CropWood") {
 908                 height = 10.0;
 909                 density = 0.1;
 910         }
 911         else if(*mat_name == "AgroForest") {
 912                 height = 10.0;
 913                 density = 0.1;
 914         }
 915         else {
 916                 height = 0.0;
 917                 density = 0.0;
 918         }
 919
 920 }
 921
 922
 923 double FGRadioTransmission::LOS_calculate_attenuation(SGGeod pos, double freq, int transmission_type) {
 924
 925         double frq_mhz = freq;
 926         double dbloss;
 927         double tx_pow = _transmitter_power;
 928         double ant_gain = _rx_antenna_gain + _tx_antenna_gain;
 929         double signal = 0.0;
 930
 931         double sender_alt_ft,sender_alt;
 932         double transmitter_height=0.0;
 933         double receiver_height=0.0;
 934         double own_lat = fgGetDouble("/position/latitude-deg");
 935         double own_lon = fgGetDouble("/position/longitude-deg");
 936         double own_alt_ft = fgGetDouble("/position/altitude-ft");
 937         double own_alt= own_alt_ft * SG_FEET_TO_METER;
 938
 939
 940         double link_budget = tx_pow - _receiver_sensitivity - _rx_line_losses - _tx_line_losses + ant_gain;
 941
 942         //cerr << "ITM:: pilot Lat: " << own_lat << ", Lon: " << own_lon << ", Alt: " << own_alt << endl;
 943
 944         SGGeod own_pos = SGGeod::fromDegM( own_lon, own_lat, own_alt );
 945
 946         SGGeod sender_pos = pos;
 947
 948         sender_alt_ft = sender_pos.getElevationFt();
 949         sender_alt = sender_alt_ft * SG_FEET_TO_METER;
 950
 951         receiver_height = own_alt;
 952         transmitter_height = sender_alt;
 953
 954         double distance_m = SGGeodesy::distanceM(own_pos, sender_pos);
 955
 956
 957         transmitter_height += _tx_antenna_height;
 958         receiver_height += _rx_antenna_height;
 959
 960
 961         /** radio horizon calculation with wave bending k=4/3 */
 962         double receiver_horizon = 4.12 * sqrt(receiver_height);
 963         double transmitter_horizon = 4.12 * sqrt(transmitter_height);
 964         double total_horizon = receiver_horizon + transmitter_horizon;
 965
 966         if (distance_m > total_horizon) {
 967                 return -1;
 968         }
 969         double pol_loss = 0.0;
 970         if (_polarization == 1) {
 971                 pol_loss = polarization_loss();
 972         }
 973         // free-space loss (distance calculation should be changed)
 974         dbloss = 20 * log10(distance_m) +20 * log10(frq_mhz) -27.55;
 975         signal = link_budget - dbloss + pol_loss;
 976
 977         //cerr << "LOS:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm " << endl;
 978         return signal;
 979
 980 }
 981
 982 /*** calculate loss due to polarization mismatch
 983 *       this function is only reliable for vertical polarization
 984 *       due to the V-shape of horizontally polarized antennas
 985 ***/
 986 double FGRadioTransmission::polarization_loss() {
 987
 988         double theta_deg;
 989         double roll = fgGetDouble("/orientation/roll-deg");
 990         if (fabs(roll) > 85.0)
 991                 roll = 85.0;
 992         double pitch = fgGetDouble("/orientation/pitch-deg");
 993         if (fabs(pitch) > 85.0)
 994                 pitch = 85.0;
 995         double theta = fabs( atan( sqrt(
 996                 pow(tan(roll * SGD_DEGREES_TO_RADIANS), 2) +
 997                 pow(tan(pitch * SGD_DEGREES_TO_RADIANS), 2) )) * SGD_RADIANS_TO_DEGREES);
 998
 999         if (_polarization == 0)
1000                 theta_deg = 90.0 - theta;
1001         else
1002                 theta_deg = theta;
1003         if (theta_deg > 85.0)   // we don't want to converge into infinity
1004                 theta_deg = 85.0;
1005
1006         double loss = 10 * log10( pow(cos(theta_deg * SGD_DEGREES_TO_RADIANS), 2) );
1007         //cerr << "Polarization loss: " << loss << " dBm " << endl;
1008         return loss;
1009 }
1010
1011
1012 double FGRadioTransmission::watt_to_dbm(double power_watt) {
1013         return 10 * log10(1000 * power_watt);   // returns dbm
1014 }
1015
1016 double FGRadioTransmission::dbm_to_watt(double dbm) {
1017         return exp( (dbm-30) * log(10.0) / 10.0);       // returns Watts
1018 }
1019
1020 double FGRadioTransmission::dbm_to_microvolt(double dbm) {
1021         return sqrt(dbm_to_watt(dbm) * 50) * 1000000;   // returns microvolts
1022 }
1023
1024