]> git.mxchange.org Git - flightgear.git/blob - src/Instrumentation/newnavradio.cxx
Improve gyro heading indicator failure simulation
[flightgear.git] / src / Instrumentation / newnavradio.cxx
1 // navradio.cxx -- class to manage a nav radio instance
2 //
3 // Written by Curtis Olson, started April 2000.
4 // Rewritten by Torsten Dreyer, August 2011
5 //
6 // Copyright (C) 2000 - 2011  Curtis L. Olson - http://www.flightgear.org/~curt
7 //
8 // This program is free software; you can redistribute it and/or
9 // modify it under the terms of the GNU General Public License as
10 // published by the Free Software Foundation; either version 2 of the
11 // License, or (at your option) any later version.
12 //
13 // This program is distributed in the hope that it will be useful, but
14 // WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16 // General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with this program; if not, write to the Free Software
20 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
21 //
22
23 #ifdef HAVE_CONFIG_H
24 #  include <config.h>
25 #endif
26
27 #include "newnavradio.hxx"
28
29 #include <assert.h>
30 #include <boost/foreach.hpp>
31
32 #include <simgear/math/interpolater.hxx>
33 #include <simgear/sg_inlines.h>
34 #include <simgear/props/propertyObject.hxx>
35 #include <simgear/misc/strutils.hxx>
36 #include <simgear/sound/sample_group.hxx>
37
38 #include <Main/fg_props.hxx>
39 #include <Navaids/navlist.hxx>
40 #include <Sound/audioident.hxx>
41
42 #include "navradio.hxx"
43
44
45 namespace Instrumentation {
46
47 using simgear::PropertyObject;
48
49 /* --------------The Navigation Indicator ----------------------------- */
50
51 class NavIndicator {
52 public:
53     NavIndicator( SGPropertyNode * rootNode ) :
54       _cdi( rootNode->getNode("heading-needle-deflection", true ) ),
55       _cdiNorm( rootNode->getNode("heading-needle-deflection-norm", true ) ),
56       _course( rootNode->getNode("radials/selected-deg", true ) ),
57       _toFlag( rootNode->getNode("to-flag", true ) ),
58       _fromFlag( rootNode->getNode("from-flag", true ) ),
59       _signalQuality( rootNode->getNode("signal-quality-norm", true ) ),
60       _hasGS( rootNode->getNode("has-gs", true ) ),
61       _gsDeflection(rootNode->getNode("gs-needle-deflection", true )),
62       _gsDeflectionDeg(rootNode->getNode("gs-needle-deflection-deg", true )),
63       _gsDeflectionNorm(rootNode->getNode("gs-needle-deflection-norm", true ))
64   {
65   }
66
67   virtual ~NavIndicator() {}
68
69   /**
70    * set the normalized CDI deflection
71    * @param norm the cdi deflection normalized [-1..1]
72    */
73   void setCDI( double norm )
74   {
75       _cdi = norm * 10.0;
76       _cdiNorm = norm;
77   }
78
79   /**
80    * set the normalized GS deflection
81    * @param norm the gs deflection normalized to [-1..1]
82    */
83   void setGS( double norm )
84   {
85       _gsDeflectionNorm = norm;
86       _gsDeflectionDeg = norm * 0.7;
87       _gsDeflection = norm * 3.5;
88   }
89
90   void setGS( bool enabled )
91   {
92       _hasGS = enabled;
93       if( !enabled ) {
94         setGS( 0.0 );
95       }
96   }
97
98   void showFrom( bool on )
99   {
100       _fromFlag = on;
101   }
102
103   void showTo( bool on )
104   {
105       _toFlag = on;
106   }
107       
108   void setSelectedCourse( double course )
109   {
110       _course = course;
111   }
112
113   double getSelectedCourse() const
114   {
115       return SGMiscd::normalizePeriodic(0.0, 360.0, _course );
116   }
117
118   void setSignalQuality( double signalQuality )
119   {
120       _signalQuality = signalQuality;
121   }
122
123 private:
124   PropertyObject<double> _cdi;
125   PropertyObject<double> _cdiNorm;
126   PropertyObject<double> _course;
127   PropertyObject<double> _toFlag;
128   PropertyObject<double> _fromFlag;
129   PropertyObject<double> _signalQuality;
130   PropertyObject<double> _hasGS;
131   PropertyObject<double> _gsDeflection;
132   PropertyObject<double> _gsDeflectionDeg;
133   PropertyObject<double> _gsDeflectionNorm;
134 };
135
136 /* ---------------------------------------------------------------- */
137
138 class NavRadioComponent {
139 public:
140   NavRadioComponent( const std::string & name, SGPropertyNode_ptr rootNode );
141   virtual ~NavRadioComponent();
142
143   virtual void   update( double dt, const SGGeod & aircraftPosition );
144   virtual void   search( double frequency, const SGGeod & aircraftPosition );
145   virtual double getRange_nm( const SGGeod & aircraftPosition );
146   virtual void   display( NavIndicator & navIndicator ) = 0;
147   virtual bool   valid() const { return NULL != _navRecord && true == _serviceable; }
148   virtual const std::string getIdent() const { return _ident; }
149
150 protected:
151   virtual double computeSignalQuality_norm( const SGGeod & aircraftPosition );
152   virtual FGNavList::TypeFilter* getNavaidFilter() = 0;
153
154   // General-purpose sawtooth function.  Graph looks like this:
155   //         /\                                    .
156   //       \/
157   // Odd symmetry, inversion symmetry about the origin.
158   // Unit slope at the origin.
159   // Max 1, min -1, period 4.
160   // Two zero-crossings per period, one with + slope, one with - slope.
161   // Useful for false localizer courses.
162   static double sawtooth(double xx)
163   {
164     return 4.0 * fabs(xx/4.0 + 0.25 - floor(xx/4.0 + 0.75)) - 1.0;
165   }
166
167   SGPropertyNode_ptr _rootNode;
168   const std::string _name;
169   FGNavRecord * _navRecord;
170   PropertyObject<bool>   _serviceable;
171   PropertyObject<double> _signalQuality_norm;
172   PropertyObject<double> _trueBearingTo_deg;
173   PropertyObject<double> _trueBearingFrom_deg;
174   PropertyObject<double> _trackDistance_m;
175   PropertyObject<double> _slantDistance_m;
176   PropertyObject<double> _heightAboveStation_ft;
177   PropertyObject<string> _ident;
178   PropertyObject<bool>   _inRange;
179   PropertyObject<double> _range_nm;
180 };
181
182 class NavRadioComponentWithIdent : public NavRadioComponent {
183 public:
184   NavRadioComponentWithIdent( const std::string & name, SGPropertyNode_ptr rootNode, AudioIdent * audioIdent );
185   virtual ~NavRadioComponentWithIdent();
186   void update( double dt, const SGGeod & aircraftPosition );
187 protected:
188   static std::string getIdentString( const std::string & name, int index );
189 private:
190   AudioIdent * _audioIdent;
191   PropertyObject<double> _identVolume;
192   PropertyObject<bool>   _identEnabled;
193 };
194
195 std::string NavRadioComponentWithIdent::getIdentString( const std::string & name, int index )
196 {
197   std::ostringstream temp;
198   temp << name << "-ident-" << index;
199   return temp.str();
200 }
201
202 NavRadioComponentWithIdent::NavRadioComponentWithIdent( const std::string & name, SGPropertyNode_ptr rootNode, AudioIdent * audioIdent ) :
203   NavRadioComponent( name, rootNode ),
204   _audioIdent( audioIdent ),
205   _identVolume( rootNode->getNode(name,true)->getNode("ident-volume",true) ),
206   _identEnabled( rootNode->getNode(name,true)->getNode("ident-enabled",true) )
207 {
208   _audioIdent->init();
209
210 }
211 NavRadioComponentWithIdent::~NavRadioComponentWithIdent()
212 {
213   delete _audioIdent;
214 }
215
216 void NavRadioComponentWithIdent::update( double dt, const SGGeod & aircraftPosition )
217 {
218   NavRadioComponent::update( dt, aircraftPosition );
219   _audioIdent->update( dt );
220
221   if( false == ( valid() && _identEnabled && _signalQuality_norm > 0.1 ) ) {
222       _audioIdent->setIdent("", 0.0 );
223       return;
224   }
225   _audioIdent->setIdent( _ident, SGMiscd::clip(_identVolume, 0.0, 1.0) );
226 }
227
228 NavRadioComponent::NavRadioComponent( const std::string & name, SGPropertyNode_ptr rootNode ) :
229   _rootNode(rootNode),
230   _name(name),
231   _navRecord(NULL),
232   _serviceable( rootNode->getNode(name,true)->getNode("serviceable",true) ),
233   _signalQuality_norm( rootNode->getNode(name,true)->getNode("signal-quality-norm",true) ),
234   _trueBearingTo_deg( rootNode->getNode(name,true)->getNode("true-bearing-to-deg",true) ),
235   _trueBearingFrom_deg( rootNode->getNode(name,true)->getNode("true-bearing-from-deg",true) ),
236   _trackDistance_m( rootNode->getNode(name,true)->getNode("track-distance-m",true) ),
237   _slantDistance_m( rootNode->getNode(name,true)->getNode("slant-distance-m",true) ),
238   _heightAboveStation_ft( rootNode->getNode(name,true)->getNode("height-above-station-ft",true) ),
239   _ident( rootNode->getNode(name,true)->getNode("ident",true) ),
240   _inRange( rootNode->getNode(name,true)->getNode("in-range",true) ),
241   _range_nm( rootNode->getNode(_name,true)->getNode("range-nm",true) )
242 {
243   simgear::props::Type typ = _serviceable.node()->getType();
244   if ((typ == simgear::props::NONE) || (typ == simgear::props::UNSPECIFIED))
245     _serviceable = true;
246 }
247
248 NavRadioComponent::~NavRadioComponent()
249 {
250 }
251
252 double NavRadioComponent::getRange_nm( const SGGeod & aircraftPosition )
253
254   if( _navRecord == NULL ) return 0.0; // no station: no range
255   double d = _navRecord->get_range();
256   if( d <= SGLimitsd::min() ) return 25.0; // no configured range: arbitrary number
257   return d; // configured range
258 }
259
260 void NavRadioComponent::search( double frequency, const SGGeod & aircraftPosition )
261 {
262   _navRecord = FGNavList::findByFreq(frequency, aircraftPosition, getNavaidFilter() );
263   if( NULL == _navRecord ) {
264     SG_LOG(SG_INSTR,SG_ALERT, "No " << _name << " available at " << frequency );
265     _ident = "";
266     return;
267   }
268   SG_LOG(SG_INSTR,SG_ALERT, "Using " << _name << "'" << _navRecord->get_ident() << "' for " << frequency );
269   _ident = _navRecord->ident();
270 }
271
272 double NavRadioComponent::computeSignalQuality_norm( const SGGeod & aircraftPosition )
273 {
274   if( false == valid() ) return 0.0;
275
276   double distance_nm = _slantDistance_m * SG_METER_TO_NM;
277   double range_nm = _range_nm;
278
279   // assume signal quality is 100% up to the published range and 
280   // decay with the distance squared further out
281   if ( distance_nm <= range_nm ) return 1.0;
282   return range_nm*range_nm/(distance_nm*distance_nm);
283 }
284
285 void NavRadioComponent::update( double dt, const SGGeod & aircraftPosition )
286 {
287     if( false == valid() ) {
288       _signalQuality_norm = 0.0;
289       _trueBearingTo_deg = 0.0;
290       _trueBearingFrom_deg = 0.0;
291       _trackDistance_m = 0.0;
292       _slantDistance_m = 0.0;
293       return;
294     } 
295
296     _slantDistance_m = dist(_navRecord->cart(), SGVec3d::fromGeod(aircraftPosition));
297
298     double az1 = 0.0, az2 = 0.0, dist = 0.0;
299     SGGeodesy::inverse(aircraftPosition, _navRecord->geod(), az1, az2, dist );
300     _trueBearingTo_deg = az1; _trueBearingFrom_deg = az2; _trackDistance_m = dist;
301     _heightAboveStation_ft = SGMiscd::max(0.0, aircraftPosition.getElevationFt() - _navRecord->get_elev_ft());
302
303     _range_nm = getRange_nm(aircraftPosition);
304     _signalQuality_norm = computeSignalQuality_norm( aircraftPosition );
305     _inRange = _signalQuality_norm > 0.2;
306 }
307
308 /* ---------------------------------------------------------------- */
309
310 static std::string VORTablePath( const char * name )
311 {
312     SGPath path( globals->get_fg_root() );
313     path.append( "Navaids" );
314     path.append(name);
315     return path.str();
316 }
317
318 class VOR : public NavRadioComponentWithIdent {
319 public:
320   VOR( SGPropertyNode_ptr rootNode);
321   virtual ~VOR();
322   virtual void update( double dt, const SGGeod & aircraftPosition );
323   virtual void display( NavIndicator & navIndicator );
324   virtual double getRange_nm(const SGGeod & aircraftPosition);
325 protected:
326   virtual double computeSignalQuality_norm( const SGGeod & aircraftPosition );
327   virtual FGNavList::TypeFilter* getNavaidFilter();
328
329 private:
330   double _totalTime;
331   class ServiceVolume {
332   public:
333     ServiceVolume() :
334       term_tbl(VORTablePath("range.term")),
335       low_tbl(VORTablePath("range.low")),
336       high_tbl(VORTablePath("range.high")) {
337     }
338     double adjustRange( double height_ft, double nominalRange_nm );
339
340   private:
341     SGInterpTable term_tbl;
342     SGInterpTable low_tbl;
343     SGInterpTable high_tbl;
344   } _serviceVolume;
345
346   PropertyObject<double> _radial;
347   PropertyObject<double> _radialInbound;
348 };
349
350 // model standard VOR/DME/TACAN service volumes as per AIM 1-1-8
351 double VOR::ServiceVolume::adjustRange( double height_ft, double nominalRange_nm )
352 {
353     if (nominalRange_nm < SGLimitsd::min() )
354       nominalRange_nm = FG_NAV_DEFAULT_RANGE;
355     
356     // extend out actual usable range to be 1.3x the published safe range
357     const double usability_factor = 1.3;
358
359     // assumptions we model the standard service volume, plus
360     // ... rather than specifying a cylinder, we model a cone that
361     // contains the cylinder.  Then we put an upside down cone on top
362     // to model diminishing returns at too-high altitudes.
363
364     if ( nominalRange_nm < 25.0 + SG_EPSILON ) {
365         // Standard Terminal Service Volume
366         return term_tbl.interpolate( height_ft ) * usability_factor;
367     } else if ( nominalRange_nm < 50.0 + SG_EPSILON ) {
368         // Standard Low Altitude Service Volume
369         // table is based on range of 40, scale to actual range
370         return low_tbl.interpolate( height_ft ) * nominalRange_nm / 40.0
371             * usability_factor;
372     } else {
373         // Standard High Altitude Service Volume
374         // table is based on range of 130, scale to actual range
375         return high_tbl.interpolate( height_ft ) * nominalRange_nm / 130.0
376             * usability_factor;
377     }
378 }
379
380 VOR::VOR( SGPropertyNode_ptr rootNode) :
381   NavRadioComponentWithIdent("vor", rootNode, new VORAudioIdent(getIdentString(string("vor"), rootNode->getIndex()))),
382   _totalTime(0.0),
383   _radial( rootNode->getNode(_name,true)->getNode("radial",true) ),
384   _radialInbound( rootNode->getNode(_name,true)->getNode("radial-inbound",true) )
385 {
386 }
387
388 VOR::~VOR()
389 {
390 }
391
392 double VOR::getRange_nm( const SGGeod & aircraftPosition )
393 {
394   return _serviceVolume.adjustRange( _heightAboveStation_ft, _navRecord->get_range() );
395 }
396
397 FGNavList::TypeFilter* VOR::getNavaidFilter()
398 {
399   static FGNavList::TypeFilter filter(FGPositioned::VOR);
400   return &filter;
401 }
402
403 double VOR::computeSignalQuality_norm( const SGGeod & aircraftPosition )
404 {
405   // apply cone of confusion. Some sources say it's opening angle is 53deg, others estimate
406   // a diameter of 1NM per 6000ft (approx. 45deg). ICAO Annex 10 says minimum 40deg.
407   // We use 1NM@6000ft and a distance-squared
408   // function to make signal-quality=100% 0.5NM@6000ft from the center and zero overhead
409   double cone_of_confusion_width = 0.5 * _heightAboveStation_ft / 6000.0 * SG_NM_TO_METER;
410   if( _trackDistance_m < cone_of_confusion_width ) {
411     double d = cone_of_confusion_width <= SGLimitsd::min() ? 1 : 
412               (1 - _trackDistance_m/cone_of_confusion_width);
413     return 1-d*d;
414   } 
415
416   // use default decay function outside the cone of confusion
417   return NavRadioComponentWithIdent::computeSignalQuality_norm( aircraftPosition );
418 }
419
420 void VOR::update( double dt, const SGGeod & aircraftPosition )
421 {
422   _totalTime += dt;
423   NavRadioComponentWithIdent::update( dt, aircraftPosition );
424
425   if( false == valid() ) {
426       _radial = 0.0;
427       return;
428   }
429
430   // an arbitrary error function
431   double error = 0.5*(sin(_totalTime/11.0) + sin(_totalTime/23.0));
432
433   // add 1% error at 100% signal-quality
434   // add 50% error at  0% signal-quality
435   // of full deflection (+/-10deg)
436   double e = 10.0 * ( 0.01 + (1-_signalQuality_norm) * 0.49 ) * error;
437
438   // compute magnetic bearing from the station (aka current radial)
439   double r = SGMiscd::normalizePeriodic(0.0, 360.0, _trueBearingFrom_deg - _navRecord->get_multiuse() + e );
440
441   _radial = r;
442   _radialInbound = SGMiscd::normalizePeriodic(0.0,360.0, 180.0 + _radial);
443 }
444
445 void VOR::display( NavIndicator & navIndicator )
446 {
447   if( false == valid() ) return;
448
449   double offset = SGMiscd::normalizePeriodic(-180.0,180.0,_radial - navIndicator.getSelectedCourse());
450   bool to = fabs(offset) >= 90.0;
451
452   if( to ) offset = -offset + copysign(180.0,offset);
453
454   navIndicator.showTo( to );
455   navIndicator.showFrom( !to );
456   // normalize to +/- 1.0 for +/- 10deg, decrease deflection with decreasing signal
457   navIndicator.setCDI( SGMiscd::clip( -offset/10.0, -1.0, 1.0 ) * _signalQuality_norm );
458   navIndicator.setSignalQuality( _signalQuality_norm );
459 }
460
461 /* ---------------------------------------------------------------- */
462 class LOC : public NavRadioComponentWithIdent {
463 public:
464   LOC( SGPropertyNode_ptr rootNode );
465   virtual ~LOC();
466   virtual void update( double dt, const SGGeod & aircraftPosition );
467   virtual void search( double frequency, const SGGeod & aircraftPosition );
468   virtual void display( NavIndicator & navIndicator );
469   virtual double getRange_nm(const SGGeod & aircraftPosition);
470
471 protected:
472   virtual double computeSignalQuality_norm( const SGGeod & aircraftPosition );
473   virtual FGNavList::TypeFilter* getNavaidFilter();
474
475 private:
476   class ServiceVolume {
477   public:
478       ServiceVolume();
479       double adjustRange( double azimuthAngle_deg, double elevationAngle_deg );
480   private:
481       SGInterpTable _azimuthTable;
482       SGInterpTable _elevationTable;
483   } _serviceVolume;
484   PropertyObject<double> _localizerOffset_norm;
485   PropertyObject<double> _localizerWidth_deg;
486 };
487
488 LOC::ServiceVolume::ServiceVolume()
489 {
490 // maybe this: http://www.tpub.com/content/aviation2/P-1244/P-12440125.htm
491   // ICAO Annex 10 - 3.1.3.2.2: The emission from the localizer
492   // shall be horizontally polarized
493   // very rough abstraction of a 5-element yagi antenna's
494   // E-plane radiation diagram
495   _azimuthTable.addEntry(   0.0, 1.0 );
496   _azimuthTable.addEntry(  10.0, 1.0 );
497   _azimuthTable.addEntry(  30.0, 0.75 );
498   _azimuthTable.addEntry(  40.0, 0.50 );
499   _azimuthTable.addEntry(  50.0, 0.20 );
500   _azimuthTable.addEntry(  60.0, 0.10 );
501   _azimuthTable.addEntry(  70.0, 0.20 );
502   _azimuthTable.addEntry(  80.0, 0.10 );
503   _azimuthTable.addEntry(  90.0, 0.05 );
504   _azimuthTable.addEntry( 105.0, 0.10 );
505   _azimuthTable.addEntry( 130.0, 0.05 );
506   _azimuthTable.addEntry( 150.0, 0.30 );
507   _azimuthTable.addEntry( 160.0, 0.40 );
508   _azimuthTable.addEntry( 170.0, 0.50 );
509   _azimuthTable.addEntry( 180.0, 0.50 );
510
511   _elevationTable.addEntry(   0.0, 0.1 );
512   _elevationTable.addEntry(  1.05, 1.0 );
513   _elevationTable.addEntry(  7.00, 1.0 );
514   _elevationTable.addEntry(  45.0, 0.3 );
515   _elevationTable.addEntry(  90.0, 0.1 );
516   _elevationTable.addEntry( 180.0, 0.01 );
517 }
518
519 double LOC::ServiceVolume::adjustRange( double azimuthAngle_deg, double elevationAngle_deg )
520 {
521     return _azimuthTable.interpolate( fabs(azimuthAngle_deg) ) * 
522         _elevationTable.interpolate( fabs(elevationAngle_deg) );
523 }
524
525 LOC::LOC( SGPropertyNode_ptr rootNode) :
526   NavRadioComponentWithIdent("loc", rootNode, new LOCAudioIdent(getIdentString(string("loc"), rootNode->getIndex()))),
527   _serviceVolume(),
528   _localizerOffset_norm( rootNode->getNode(_name,true)->getNode("offset-norm",true) ),
529   _localizerWidth_deg( rootNode->getNode(_name,true)->getNode("width-deg",true) )
530 {
531 }
532
533 LOC::~LOC()
534 {
535 }
536
537 FGNavList::TypeFilter* LOC::getNavaidFilter()
538 {
539   return FGNavList::locFilter();
540 }
541
542 void LOC::search( double frequency, const SGGeod & aircraftPosition )
543 {
544   NavRadioComponentWithIdent::search( frequency, aircraftPosition );
545   if( false == valid() ) {
546       _localizerWidth_deg = 0.0;
547       return;
548   }
549
550   // cache slightly expensive value, 
551   // sanitized in FGNavRecord::localizerWidth() to  never become zero
552   _localizerWidth_deg = _navRecord->localizerWidth();
553 }
554
555 /* Localizer coverage (ICAO Annex 10 Volume I 3.1.3.3 
556   25NM within +/-10 deg from the front course line
557   17NM between 10 and 35deg from the front course line
558   10NM outside of +/- 35deg  if coverage is provided
559   at and above a height of 2000ft above threshold or
560   1000ft above the highest point within intermediate
561   and final approach areas. Upper limit is a surface
562   extending outward from the localizer and inclined at
563   7 degrees above the horizontal
564  */
565 double LOC::getRange_nm(const SGGeod & aircraftPosition)
566 {
567   double elevationAngle = ::atan2(_heightAboveStation_ft*SG_FEET_TO_METER, _trackDistance_m)*SG_RADIANS_TO_DEGREES;
568   double azimuthAngle = SGMiscd::normalizePeriodic( -180.0, 180.0, _trueBearingFrom_deg + 180.0 - _navRecord->get_multiuse() );
569
570   // looks like our navrecord declared range is based on 10NM?
571   return  _navRecord->get_range() * _serviceVolume.adjustRange( azimuthAngle, elevationAngle );
572 }
573
574 double LOC::computeSignalQuality_norm( const SGGeod & aircraftPosition )
575 {
576   return NavRadioComponentWithIdent::computeSignalQuality_norm( aircraftPosition );
577 }
578
579 void LOC::update( double dt, const SGGeod & aircraftPosition )
580 {
581   NavRadioComponentWithIdent::update( dt, aircraftPosition );
582
583   if( false == valid() ) {
584     _localizerOffset_norm = 0.0;
585     return;
586   }
587
588   double offset = SGMiscd::normalizePeriodic( -180.0, 180.0, _trueBearingFrom_deg + 180.0 - _navRecord->get_multiuse() );
589
590   // The factor of 30.0 gives a period of 120 which gives us 3 cycles and six 
591   // zeros i.e. six courses: one front course, one back course, and four 
592   // false courses. Three of the six are reverse sensing.
593   offset = 30.0 * sawtooth(offset / 30.0);
594
595   // normalize offset to the localizer width, scale and clip to [-1..1]
596   offset = SGMiscd::clip( 2.0 * offset / _localizerWidth_deg, -1.0, 1.0 );
597   
598   _localizerOffset_norm = offset;
599 }
600
601 void LOC::display( NavIndicator & navIndicator )
602 {
603   if( false == valid() ) 
604     return;
605
606   navIndicator.showTo( true );
607   navIndicator.showFrom( false );
608
609   navIndicator.setCDI( _localizerOffset_norm * _signalQuality_norm );
610   navIndicator.setSignalQuality( _signalQuality_norm );
611 }
612
613 class GS : public NavRadioComponent {
614 public:
615   GS( SGPropertyNode_ptr rootNode);
616   virtual ~GS();
617   virtual void update( double dt, const SGGeod & aircraftPosition );
618   virtual void search( double frequency, const SGGeod & aircraftPosition );
619   virtual void display( NavIndicator & navIndicator );
620
621   virtual double getRange_nm(const SGGeod & aircraftPosition);
622 protected:
623   virtual FGNavList::TypeFilter* getNavaidFilter();
624 private:
625   class ServiceVolume {
626   public:
627       ServiceVolume();
628       double adjustRange( double azimuthAngle_deg, double elevationAngle_deg );
629   private:
630       SGInterpTable _azimuthTable;
631       SGInterpTable _elevationTable;
632   } _serviceVolume;
633   static SGVec3d tangentVector(const SGGeod& midpoint, const double heading);
634
635   PropertyObject<double>  _targetGlideslope_deg;
636   PropertyObject<double>  _glideslopeOffset_norm;
637   SGVec3d _gsAxis;
638   SGVec3d _gsVertical;
639 };
640
641 GS::ServiceVolume::ServiceVolume()
642 {
643 // maybe this: http://www.tpub.com/content/aviation2/P-1244/P-12440125.htm
644   // ICAO Annex 10 - 3.1.5.2.2: The emission from the glide path equipment
645   // shall be horizontally polarized
646   // very rough abstraction of a 5-element yagi antenna's
647   // E-plane radiation diagram
648   _azimuthTable.addEntry(   0.0, 1.0 );
649   _azimuthTable.addEntry(  10.0, 1.0 );
650   _azimuthTable.addEntry(  30.0, 0.75 );
651   _azimuthTable.addEntry(  40.0, 0.50 );
652   _azimuthTable.addEntry(  50.0, 0.20 );
653   _azimuthTable.addEntry(  60.0, 0.10 );
654   _azimuthTable.addEntry(  70.0, 0.20 );
655   _azimuthTable.addEntry(  80.0, 0.10 );
656   _azimuthTable.addEntry(  90.0, 0.05 );
657   _azimuthTable.addEntry( 105.0, 0.10 );
658   _azimuthTable.addEntry( 130.0, 0.05 );
659   _azimuthTable.addEntry( 150.0, 0.30 );
660   _azimuthTable.addEntry( 160.0, 0.40 );
661   _azimuthTable.addEntry( 170.0, 0.50 );
662   _azimuthTable.addEntry( 180.0, 0.50 );
663
664   _elevationTable.addEntry(   0.0, 0.1 );
665   _elevationTable.addEntry(  1.05, 1.0 );
666   _elevationTable.addEntry(  7.00, 1.0 );
667   _elevationTable.addEntry(  45.0, 0.3 );
668   _elevationTable.addEntry(  90.0, 0.1 );
669   _elevationTable.addEntry( 180.0, 0.01 );
670 }
671
672 double GS::ServiceVolume::adjustRange( double azimuthAngle_deg, double elevationAngle_deg )
673 {
674     return _azimuthTable.interpolate( fabs(azimuthAngle_deg) ) * 
675         _elevationTable.interpolate( fabs(elevationAngle_deg) );
676 }
677
678 GS::GS( SGPropertyNode_ptr rootNode) :
679   NavRadioComponent("gs", rootNode ),
680   _targetGlideslope_deg( rootNode->getNode(_name,true)->getNode("slope",true) ),
681   _glideslopeOffset_norm( rootNode->getNode(_name,true)->getNode("offset-norm",true) ),
682   _gsAxis(SGVec3d::zeros()),
683   _gsVertical(SGVec3d::zeros())
684 {
685 }
686
687 GS::~GS()
688 {
689 }
690
691 FGNavList::TypeFilter* GS::getNavaidFilter()
692 {
693   static FGNavList::TypeFilter filter(FGPositioned::GS);
694   return &filter;
695 }
696
697 double GS::getRange_nm(const SGGeod & aircraftPosition)
698 {
699   double elevationAngle = ::atan2(_heightAboveStation_ft*SG_FEET_TO_METER, _trackDistance_m)*SG_RADIANS_TO_DEGREES;
700   double azimuthAngle = SGMiscd::normalizePeriodic( -180.0, 180.0, _trueBearingFrom_deg + 180.0 - fmod(_navRecord->get_multiuse(), 1000.0) );
701   return  _navRecord->get_range() * _serviceVolume.adjustRange( azimuthAngle, elevationAngle );
702 }
703
704 // Calculate a Cartesian unit vector in the
705 // local horizontal plane, i.e. tangent to the 
706 // surface of the earth at the local ground zero.
707 // The tangent vector passes through the given  <midpoint> 
708 // and points forward along the given <heading>.
709 // The <heading> is given in degrees.
710 SGVec3d GS::tangentVector(const SGGeod& midpoint, const double heading)
711 {
712   // move 100m away from the midpoint - arbitrary number
713   const double delta(100.0);
714   SGGeod head, tail;
715   double az2;                   // ignored
716   SGGeodesy::direct(midpoint, heading,     delta, head, az2);
717   SGGeodesy::direct(midpoint, 180+heading, delta, tail, az2);
718   head.setElevationM(midpoint.getElevationM());
719   tail.setElevationM(midpoint.getElevationM());
720   SGVec3d head_xyz = SGVec3d::fromGeod(head);
721   SGVec3d tail_xyz = SGVec3d::fromGeod(tail);
722 // Awkward formula here, needed because vector-by-scalar
723 // multiplication is defined, but not vector-by-scalar division.
724   return (head_xyz - tail_xyz) * (0.5/delta);
725 }
726
727 void GS::search( double frequency, const SGGeod & aircraftPosition )
728 {
729   NavRadioComponent::search( frequency, aircraftPosition );
730   if( false == valid() ) {
731       _gsAxis = SGVec3d::zeros();
732       _gsVertical = SGVec3d::zeros();
733       _targetGlideslope_deg = 3.0;
734       return;
735   }
736   
737   double gs_radial = SGMiscd::normalizePeriodic(0.0, 360.0, fmod(_navRecord->get_multiuse(), 1000.0) );
738
739   _gsAxis = tangentVector(_navRecord->geod(), gs_radial);
740   SGVec3d gsBaseline = tangentVector(_navRecord->geod(), gs_radial + 90.0);
741   _gsVertical = cross(gsBaseline, _gsAxis);
742
743   int tmp = (int)(_navRecord->get_multiuse() / 1000.0);
744   // catch unconfigured glideslopes here, they will cause nan later
745   _targetGlideslope_deg = SGMiscd::max( 1.0, (double)tmp / 100.0 );
746 }
747
748 void GS::update( double dt, const SGGeod & aircraftPosition )
749 {
750   NavRadioComponent::update( dt, aircraftPosition );
751   if( false == valid() ) {
752       _glideslopeOffset_norm = 0.0;
753       return;
754   }
755   
756   SGVec3d pos = SGVec3d::fromGeod(aircraftPosition) - _navRecord->cart(); // relative vector from gs antenna to aircraft
757   // The positive GS axis points along the runway in the landing direction,
758   // toward the far end, not toward the approach area, so we need a - sign here:
759   double comp_h = -dot(pos, _gsAxis);      // component in horiz direction
760   double comp_v = dot(pos, _gsVertical);   // component in vertical direction
761   //double comp_b = dot(pos, _gsBaseline);   // component in baseline direction
762   //if (comp_b) {}                           // ... (useful for debugging)
763
764 // _gsDirect represents the angle of elevation of the aircraft
765 // as seen by the GS transmitter.
766   double gsDirect = atan2(comp_v, comp_h) * SGD_RADIANS_TO_DEGREES;
767 // At this point, if the aircraft is centered on the glide slope,
768 // _gsDirect will be a small positive number, e.g. 3.0 degrees
769
770 // Aim the branch cut straight down 
771 // into the ground below the GS transmitter:
772   if (gsDirect < -90.0) gsDirect += 360.0;
773
774   double offset = _targetGlideslope_deg - gsDirect;
775   if( offset < 0.0 )
776     offset = _targetGlideslope_deg/2 * sawtooth(2.0*offset/_targetGlideslope_deg);
777   assert( false == isnan(offset) );
778 // GS is documented to be 1.4 degrees thick, 
779 // i.e. plus or minus 0.7 degrees from the midline:
780   _glideslopeOffset_norm = SGMiscd::clip(offset/0.7, -1.0, 1.0);
781 }
782
783 void GS::display( NavIndicator & navIndicator )
784 {
785   if( false == valid() ) {
786     navIndicator.setGS( false );
787     return;
788   }
789   navIndicator.setGS( true );
790   navIndicator.setGS( _glideslopeOffset_norm );
791 }
792
793 /* ------------- A NAV/COMM Frequency formatter ---------------------- */
794
795 class FrequencyFormatter : public SGPropertyChangeListener {
796 public:
797   FrequencyFormatter( SGPropertyNode_ptr freqNode, SGPropertyNode_ptr fmtFreqNode, double channelSpacing ) :
798     _freqNode( freqNode ),
799     _fmtFreqNode( fmtFreqNode ),
800     _channelSpacing(channelSpacing)
801   {
802     _freqNode->addChangeListener( this );
803     valueChanged(_freqNode);
804   }
805   ~FrequencyFormatter()
806   {
807     _freqNode->removeChangeListener( this );
808   }
809
810   void valueChanged (SGPropertyNode * prop)
811   {
812     // format as fixed decimal "nnn.nn"
813     std::ostringstream buf;
814     buf << std::fixed 
815         << std::setw(5) 
816         << std::setfill('0') 
817         << std::setprecision(2)
818         << getFrequency();
819     _fmtFreqNode->setStringValue( buf.str() );
820   }
821
822   double getFrequency() const 
823   {
824     double d = SGMiscd::roundToInt(_freqNode->getDoubleValue() / _channelSpacing) * _channelSpacing;
825     // strip last digit, do not round
826     return ((int)(d*100))/100.0;
827   }
828
829 private:
830   SGPropertyNode_ptr _freqNode;
831   SGPropertyNode_ptr _fmtFreqNode;
832   double _channelSpacing;
833 };
834
835
836 /* ------------- The NavRadio implementation ---------------------- */
837
838 class NavRadioImpl : public NavRadio {
839 public:
840   NavRadioImpl( SGPropertyNode_ptr node );
841   virtual ~NavRadioImpl();
842
843   virtual void update( double dt );
844   virtual void init();
845 private:
846   void search();
847
848   class Legacy {
849   public:
850       Legacy( NavRadioImpl * navRadioImpl ) : _navRadioImpl( navRadioImpl ) {}
851
852       void init();
853       void update( double dt );
854   private:
855       NavRadioImpl * _navRadioImpl;
856       SGPropertyNode_ptr is_valid_node;
857       SGPropertyNode_ptr nav_serviceable_node;
858       SGPropertyNode_ptr nav_id_node;
859       SGPropertyNode_ptr id_c1_node;
860       SGPropertyNode_ptr id_c2_node;
861       SGPropertyNode_ptr id_c3_node;
862       SGPropertyNode_ptr id_c4_node;
863   } _legacy;
864
865   const static int VOR_COMPONENT = 0;
866   const static int LOC_COMPONENT = 1;
867   const static int GS_COMPONENT  = 2;
868
869   std::string _name;
870   int         _num;
871   SGPropertyNode_ptr _rootNode;
872   FrequencyFormatter _useFrequencyFormatter;
873   FrequencyFormatter _stbyFrequencyFormatter;
874   std::vector<NavRadioComponent*> _components;
875   NavIndicator _navIndicator;
876   double _stationTTL;
877   double _frequency;
878   PropertyObject<bool> _cdiDisconnected;
879 };
880
881 NavRadioImpl::NavRadioImpl( SGPropertyNode_ptr node ) :
882   _legacy( this ),
883   _name(node->getStringValue("name", "nav")),
884   _num(node->getIntValue("number", 0)),
885   _rootNode(fgGetNode( string("/instrumentation/") + _name, _num, true)),
886   _useFrequencyFormatter( _rootNode->getNode("frequencies/selected-mhz",true), _rootNode->getNode("frequencies/selected-mhz-fmt",true), 0.05 ),
887   _stbyFrequencyFormatter( _rootNode->getNode("frequencies/standby-mhz",true), _rootNode->getNode("frequencies/standby-mhz-fmt",true), 0.05 ),
888   _navIndicator(_rootNode),
889   _stationTTL(0.0),
890   _frequency(-1.0),
891   _cdiDisconnected(_rootNode->getNode("cdi-disconnected",true))
892 {
893 }
894
895 NavRadioImpl::~NavRadioImpl()
896 {
897   BOOST_FOREACH( NavRadioComponent * p, _components ) {
898     delete p;
899   }
900 }
901
902 void NavRadioImpl::init()
903 {
904   if( 0 < _components.size() )
905     return;
906
907   _components.push_back( new VOR(_rootNode) );
908   _components.push_back( new LOC(_rootNode) );
909   _components.push_back( new GS(_rootNode) );
910
911   _legacy.init();
912 }
913
914 void NavRadioImpl::search()
915 {
916 }
917
918 void NavRadioImpl::update( double dt )
919 {
920   if( dt < SGLimitsd::min() ) return;
921
922   SGGeod position;
923
924   try {
925     position = globals->get_aircraft_position();
926   }
927   catch( std::exception & ) {
928     return;
929   }
930
931   _stationTTL -= dt;
932   if( _frequency != _useFrequencyFormatter.getFrequency() ) {
933       _frequency = _useFrequencyFormatter.getFrequency();
934       _stationTTL = 0.0;
935   }
936
937   BOOST_FOREACH( NavRadioComponent * p, _components ) {
938       if( _stationTTL <= 0.0 )
939           p->search( _frequency, position );
940       p->update( dt, position );
941
942       if( false == _cdiDisconnected )
943           p->display( _navIndicator );
944   }
945
946   if( _stationTTL <= 0.0 )
947       _stationTTL = 30.0;
948
949   _legacy.update( dt );
950 }
951
952 void NavRadioImpl::Legacy::init()
953 {
954     is_valid_node = _navRadioImpl->_rootNode->getChild("data-is-valid", 0, true);
955     nav_serviceable_node = _navRadioImpl->_rootNode->getChild("serviceable", 0, true);
956
957     nav_id_node = _navRadioImpl->_rootNode->getChild("nav-id", 0, true );
958     id_c1_node = _navRadioImpl->_rootNode->getChild("nav-id_asc1", 0, true );
959     id_c2_node = _navRadioImpl->_rootNode->getChild("nav-id_asc2", 0, true );
960     id_c3_node = _navRadioImpl->_rootNode->getChild("nav-id_asc3", 0, true );
961     id_c4_node = _navRadioImpl->_rootNode->getChild("nav-id_asc4", 0, true );
962
963 }
964
965 void NavRadioImpl::Legacy::update( double dt )
966 {
967     is_valid_node->setBoolValue( 
968         _navRadioImpl->_components[VOR_COMPONENT]->valid() || _navRadioImpl->_components[LOC_COMPONENT]->valid()  
969         );
970
971     string ident = _navRadioImpl->_components[VOR_COMPONENT]->getIdent();
972     if( ident.empty() )
973         ident = _navRadioImpl->_components[LOC_COMPONENT]->getIdent();
974
975     nav_id_node->setStringValue( ident );
976
977     ident = simgear::strutils::rpad( ident, 4, ' ' );
978     id_c1_node->setIntValue( (int)ident[0] );
979     id_c2_node->setIntValue( (int)ident[1] );
980     id_c3_node->setIntValue( (int)ident[2] );
981     id_c4_node->setIntValue( (int)ident[3] );
982 }
983
984
985 SGSubsystem * NavRadio::createInstance( SGPropertyNode_ptr rootNode )
986 {
987     // use old navradio code by default
988     if( fgGetBool( "/instrumentation/use-new-navradio", false ) )
989         return new NavRadioImpl( rootNode );
990
991     return new FGNavRadio( rootNode );
992 }
993
994 } // namespace Instrumentation
995