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