src/Radio/radio.cxx

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