src/Radio/radio.cxx

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