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