/** * Aircraft Navigation Functions * @author: Robert John Morton UK-YE572246C * @version: 1997.01 These 'C' functions do not constitute a running program. They are merely snippets of code used in the prototyping of the Moving Map Navigator I wrote to navigate an air vehicle completely automatically, along a pre- scribed route of waypoints, from any origin to any destination on the planet. An applet illustrating this Moving Map Navigator was written in Java and appears within this section of the web site. Because the major web browsers no longer support Java applets, the applet has been re- written as a Java Web Start application, which can be launched from the same web page. Details of the aircraft's position is maintained by the flight program in a structure: */ struct AirData { double Lat, // aircraft's latitude Lng, // aircraft's longitude Ht, // height above sea-level (Earth Radians) GndHt, // estimated height of ground beneath aircraft Hdg, // ground-track heading (compass + wind − drift angle) Des, // aircraft's angle of descent Spd, // aircraft's true horizontal ground speed Pitch, // aircraft's pitch angle Roll, // aircraft's roll angle Rot, // rate-of-turn in radians per second TD; // diameter of aircraft's minimum turning circle int Sw; // Id of NAV receiver switched to navigation computer } Air; /* The record for each station within GS.DAT has the following structure: */ struct StnData { char Type; // NDB,RRG,VOR,ILS,DME,AMK,BMK,OMK,MMK,IMK double Lat, // Latitude Lng, // Longitude Ht, // Ht of stn above sea level (Earth radians) Rng, // Effective Range Hdg; // Runway or airway heading int RadPat; // 0=omni,1=RRG,2=Marker lobe,3=ILS lobe long Frq; // Frequency char CallSign[16]; // call sign coding string time_t TimeOut, // Time before next dist & bearing check OldTime; // Time when dist & bearing last checked BOOL Receivable; // Flag saying if station is in range } Stn; // pointer to an instance of the above structure /* The first 9 items are constants for a given station, the last 3 vary with aircraft position. GS.DAT is generated from a ground station data source file GSS.DAT into which details of each ground station are typed initially. */ /* The data structure for each currently-active station has the following contents. */ struct TxData { // DATA STRUCTURE FOR AN ACTIVE GROUND STN int StnNum, // Stn's record No within GS.DAT file StnType; // see below double Lat, // Latitude Lng, // Longitude Ht, // Ht of stn > sea level in Earth radians Rng, // Effective Range Hdg, // Runway or airway heading Len, // Runway length (for ILS only) GS; // Glide Slope Angle (for ILS only) int RadPat; // 0=omni, 1=RRG, 2=Marker lobe, 3=ILS lobe long Frq; // Frequency char CallSign[16]; // call sign coding string BOOL Single; // morse version of callsign in Morse[] double Dst, PrevDst, // Current & previous dist from aircraft Brg, // Current true bearing from aircraft Sig; // Signal strength of tx at aircraft char *pCallSign, // ptr to current posn in callsign string Morse[48], // Morse sequence in timer durations *pMorse; // points to current Morse element int Repeat, // Nº of morse seq repeats still to do Pitch, // audio freq of morse tone (see below) MorseTimer; // duration of current Morse Code element BOOL MorseKey; // current state of Morse Key } *Tx[NumTx], // ptrs to NumTx active stn structures *(*TxHi); // ptr to Tx[highest occupied element+1] /* StnType: 1 = NDB, 2 = RRG, 3 = VOR, 4 = VOR/DME, 5 = ILS, 6 = ILS/DME, 7 = DME, 8 = AMK, 9 = BMK, 10 = OMK, 11 = MMK, 12 = IMK Pitch 0 = 1020Hz, 1 = 1300Hz, 2 = 3000Hz The items above the blank line are fixed items copied from the GS.DAT file. The items below the blank line are updated by the real-time NAV software. When a station comes within receivable range, this function copies its details into a newly-allocated TxData structure. Apart from simply copying the fixed data items from the appropriate GS.DAT disk record, SetUpTxData() also possibly sets up the station's callsign for keying. If a station has a simple callsign sequence (ie only one semi-colon appears in it - see later in the section on callsign keying), SetUpTxData() calls SetUpCallSign() which translates the callsign sequence into Morse code once for the whole time the station is active. If the station has a multiple callsign (ie more than one semi-colon appears in it) then each segment of the callsign is set up dynamically during the on-going callsign keying cycle. */ void SetUpTxData(struct TxData *t) { char *u = Stn.CallSign, // ptr to callsign data in disk stream *v = t->CallSign; // ptr to callsign in TxData structure t->Type = Stn.Type; // Type of Station t->Lat = Stn.Lat; // Its Latitude t->Lng = Stn.Lng; // Its Longitude t->Ht = Stn.Ht; // Its height > sea level (Earth rads) t->Rng = Stn.Rng - RNG; // True transmission range t->Hdg = Stn.Hdg; // runway or airway heading t->RadPat = Stn.RadPat; // antenna radiation pattern or gain t->Frq = Stn.Frq; // Frequency if(*u++ = 'S') // S = single, M = multiple callsign t->Single = TRUE; // single callsign else t->Single = FALSE; // multiple callsign for(;*u > '\0';u++,v++) // copy all but S/M byte of callsign in *v = *u; // to TxData callsign buffer CallSign[] *v = '\0'; // and terminate it with a null char t->pCallSign = Stn.CallSign; // set ptr to start of callsign t->MorseKey = TRUE; // for first Morse element t->pMorse = t->Morse; // set ptr to start of Morse Code if(t->Single == TRUE) // if it is a single-part callsign SetUpCallSign(t); // translate it to Morse once-only here } /* Stations are added to and deleted from the active stations list by the function: [StnNum = Stn's Record Nº in GS.DAT file, flag = TRUE if new stn just come in range] */ BOOL SetStn(int StnNum, // Stn's Record Nº in GS.DAT file BOOL flag) { // TRUE if new stn just come in range #define RNG 0.0157079632675 // 100km in Great Circle radians BOOL result = FALSE; static struct TxData *(*TX) = Tx + NumTx, // ptr to last element in ptr array *t, // ptr to subject Tx Data structure **tx; // ptr to ptr to Tx Data structure if(flag == TRUE) { // IF A STATION HAS JUST COME INTO RANGE for(tx = Tx;tx < TX;tx++) { // find spare slot in ptr array if(*tx == NULL) { // then allocate memory for station data *tx = (struct TxData *)malloc(size_t sizeof(struct TxData)); if(*tx != NULL) result = TRUE; break; } } if(result == TRUE) { // if enough memory could be allocated if(TxHi < tx) // keep track of the TxHi = tx; // highest used element SetUpTxData(*tx);} // copy stn data to TxData structure } else { // ELSE A STATION HAS JUST GONE OUT OF RANGE for (tx = Tx;tx < TX;tx++) { // find its Tx Data if((*tx)->StnNum == StnNum) { free(*tx); // free the memory used by this station *tx = NULL; // clear ptr in TxData ptr array if(TxHi == tx) // keep track of the TxHi--; // highest used element result = TRUE; break; } } } return(result); } /* When a station's Stn.TimeOut expires we call the function GetDst() to compute its true accurate great-circle distance as follows: */ double GetDst() { // COMPUTE THE DISTANCE OF AN #define Pi 3.1415926535 // OBSERVED FROM THE OBSERVER double dlong = Stn.Lng - Air.Lng, // Longitudinal difference SinAirLat = sin(Air.Lat), // Sine of latitude of observer CosAirLat = cos(Air.Lat), // Cosine latitude of observer SinStnLat = sin(Stn.Lat), // Sine of latitude of observed CosStnLat = cos(Stn.Lat); // Cosine latitude of observed if(dlong > Pi) dlong -= Pi; // change the range of dlong if(dlong < -Pi) dlong += Pi; // from 0 - 2p to -p - +p CosDst = CosAirLat * CosStnLat * cos(dlong) + SinAirLat * SinStnLat; if(CosDst > +1) CosDst = +1; // avoid acos() domain errors if(CosDst < -1) CosDst = -1; return(acos(CosDst)); // return distance in radians } /* The following 'C' function determines whether or not a ground station is currently receivable at the aircraft. It is called about every 200 milliseconds. */ void StnPoll() { #define VMAX 0.4363323129861111E-4 // 1000km/hr in radians/sec static int StnNum = 0; // Station No preserved from prev pass double d; // distance of station from aircraft if (StnNum == NumStn) // Station's record number in the StnNum = 0; // Ground Station Data file GS.DAT GetStn(++StnNum); // get data of next stn in GS.DAT file if (Stn.TimeOut > 0) { // if not yet time to recalculate distance /* Decrement the station's distance re-calculation timeout, and advance the station's ghost towards the aircraft. */ Stn.TimeOut -= SysTime - Stn.OldTime; Stn.OldTime = SysTime; } else { // Else station's ghost is now in range of aircraft, so... /* If aircraft is within station's signal range, and station is not currently on the 'receivable stations' list then.. */ if(((d = GetDst()) < Stn.Rng) && (Stn.Receivable == FALSE)) { /*Attempt to put the station on the 'receivable' list and if successful store its updated status to the GS.DAT file*/ if((Stn.Receivable = SetStn(StnNum, TRUE)) == TRUE) PutStn(StnNum); } else { // Else the station is out of receivable range, so... /* If station currently on 'receivable' list, delete it from 'active stations' list and set it as unreceivable. */ if (Stn.Receivable == TRUE) { SetStn(StnNum, FALSE); Stn.Receivable = FALSE; } /* Set the timeout to next distance check and then store the updated station data to the GS.DAT file. */ Stn.TimeOut = (d - Stn.Rng) / VMAX; PutStn(StnNum); } } } /* The function returns the directional attenuation factor of the station's signal computed according to a radio range station's 4-leaf clover radiation pattern given by the formula: Signal = MaxSignal * sin( 2 * Bearing ). */ double RadioRange( struct TxData *t ) { #define Pi 3.1415926535 #define HalfPi 1.5707963267 char s; // contains 'A' or 'N' according to quadrant double b = t->Brg; // bearing of aircraft from station if(b > Pi) // if aircraft in 3rd or 4th quadrants, b -= Pi; // consider it in the first or second if(b > HalfPi) // if aircraft in the second or 4th quadrants s = 'N'; // key the letter 'N' else // if aircraft in the first or 3rd quadrants s = 'A'; // key the letter 'A' // Put letter in 4th character position of first callsign segment *( t->CallSign + 3 ) = s; // Compute directional reduction in station's effective range double a = t->Rng * sin( b + b ); // return the current signal strength of station at aircraft return( ( a - t->Dst ) / a ); } /* The following function computes the signal strength within the field of a high-gain antenna as used by a marker beacon or an ILS station. The square root expression is the distance 'R' along the central axis of the station's radiation lobe at which the signal is the same strength as it is at the aircraft. */ double RadiationLobe ( double rsq, //square of distance between aircraft and station double R, //maximum effective range of station double theta; //relative bearing between aircraft and lobe axis ) { return((R - sqrt(rsq /(cos(3.5 * theta) * cos(theta))))/ R); } /* The function below finds a marker's signal strength at the aircraft. The range of a marker beacon is given as a vertical range or height. This is because a marker's radiation lobe points upwards rather than horizontally like a Radio Range or ILS. The function finds the signal strength at the aircraft by finding the height on the central vertical axis of the marker's radiation lobe at which the signal is the same strength as at the aircraft. It then compares this with the maximum range of the marker (vertically) to find the relative signal strength at the aircraft. */ double MarkerBeacon( struct TxData *t ) { double // relative bearing of aircraft from marker axis x = sin(t->Brg - t->Hdg), d = t->Dst, // horizontal distance of aircraft from marker h = Air.Ht, // aircraft height /* Compute the square of the distance along the airway at the height of the aircraft at which marker's signal strength is the same as it is at the aircraft. */ asq = (d * d) / (1 + 3 * x * x), /* Compute the angle from the vertical axis of the marker's radiation lobe to the point along the airway at which the marker's signal strength is the same as it is at the aircraft. */ theta = atan(sqrt(asq) / h ); return(RadiationLobe(asq + h * h, t->Rng, theta)); } /* The following function uses the above to compute the signal strength at the aircraft which is somewhere within range of an ILS localiser station. */ double ILSlocaliser( struct TxData *t ) { double r = t->Dst; // distance between aircraft and station return ( RadiationLobe( r * r, t->Rng, t->Brg - t->Hdg ) ); } /* The function below uses spherical geometry to compute the distance between the aircraft and the station whose TxData is pointed to by 'q'. It then computes the bearing as follows: if FromAircraft = FALSE it computes bearing of aircraft as observed from station if FromAircraft = TRUE it computes bearing of station as observed from aircraft. */ BOOL DandB ( struct TxData *t, // pointer to relevant station's details double MaxDst, // maximum relevant distance BOOL FromAircraft // if TRUE swap observer and observed ) { #define Pi 3.1415926535 double dlong = t->Lng - Air.Lng, // Longitudinal difference SinAirLat = sin(Air.Lat), // Sine of latitude of observer CosAirLat = cos(Air.Lat), // Cosine of latitude of observer SinStnLat = sin(t->Lat), // Sine of latitude of observed CosStnLat = cos(t->Lat); // Cosine of latitude of observed if(dlong > Pi) dlong -= Pi; // change dlong's range if(dlong < -Pi) dlong += Pi; // from 0 - 2p to -p - +p double CosDst = CosAirLat * CosStnLat * cos(dlong) + SinAirLat * SinStnLat; if(CosDst > +1) CosDst = +1; // avoid possible acos() if(CosDst < -1) CosDst = -1; // domain errors if((double Dst = acos(CosDst)) > MaxDst) return( FALSE ); double Brg = 0; // avoid possibility of division by zero if((double SinDst = sin(Dst)) > 0) { if(FromAircraft == TRUE) { // if doing bearing of the station double a = SinStnLat; // from the aircraft, then swap SinStnLat = SinAirLat; // the aircraft and station SinAirLat = a; // co-ordinates over. } double CosBrg = (SinStnLat - SinAirLat * CosDst) / (CosAirLat * SinDst); if(CosBrg < +1 && CosBrg > -1) { // avoid possible acos() Brg = acos(CosBrg); // domain errors if(dlong < 0) Brg += Pi; } else if(CosBrg < 1) Brg = Pi; } t->Dst = Dst; t->Brg = Brg; // put dist & brg in station's return(TRUE); // TxData array } /* Please note: the full implementation of the Moving Map Global Navigator does not use the above method of computing distance and bearing. Instead, it uses the "Solution to the geodetic inverse problem after T. Vincenty, which is a modified Rainsford's Method with Helmert's elliptical terms." This gives far more accurate results and offers a choice of elliptical geoids. This more accurate method is implemented in the applet on this web site. */ /* A function to compute the signal strength a transmitter delivers at the aircraft's current position according to distance and line-of-sight criteria is show below: */ double TxSig(struct TxData *t) { #define LOS .005811946408975 // 37km Line-Of-Sight // (expressed in Earth radians) double (*pf[ ])(struct TxData *) = // array of ptrs to functions { OmniSignal, // omni-directional radiation patterns RadioRange, // 4-leaf clover radiation pattern MarkerBeacon, // vertical lobe radiation pattern ILSlocaliser, // horizontal lobe radiation pattern } Los, // line-of-sight distance limit Sig; // tx's signal strength at aircraft if(t->Frq > 30000) { // If it is a VHF station (freq > 30MHz) /* Compute line-of-sight range limit. If less than horizon scatter range (20 nautical miles or 37 km) set line-of-sight range = scatter range. */ if((Los = sqrt(1.5 * (Air.Ht - t->Ht))) < LOS) Los = LOS; } else Los = Pi; // else set it to half Earth circumference /* If station is effectively within 'line-of-sight', compute signal strength of station at aircraft. If signal strength returned is negative, make it zero. */ if((DandB(t, Los, FALSE) == TRUE) && ((Sig = (*(pf + t->RadPat))(t)) < 0)) Sig = 0; else Sig = 0; // zero signal strength if beyond line-of-sight return(Sig); // return signal strength at aircraft } /* FUNCTION TO SET UP THE CALL SIGN OF A GIVEN STATION READY FOR KEYING @author: Robert John Morton UK-YE572246C @version: 17 January 1997. The parameter 't' is a pointer to the Tx Data structure of the station currently being keyed. */ void SetUpCallSign(struct TxData *t) { /* Array letters of the Morse Code alphabet where the code for each letter is expressed as timer durations: dot=1, dash = 3. */ char *MorseCode[] = { "\3\1", "\3\1\1\3", "\3\1\3\1", "\3\1\1", "\1", "\1\1\3\1", "\3\3\1", "\1\1\1\1", "\1\1", "\1\3\3\3", "\3\1\3", "\1\3\1\1", "\3\3", "\3\1", "\3\3\3", "\1\3\3\1", "\3\3\1\3", "\1\3\1", "\1\1\1", "\3", "\1\1\3", "\1\1\1\3", "\1\3\3", "\3\1\1\3", "\3\1\3\3", "\3\3\1\1" }; char MC, // Morse Code element (dot/dash/gap) *pMC, // ptr to Morse element in MorseCode[ ] *pC = t->pCallSign, // ptr to start of next callsign segment *pM = t->pMorse; // ptr to start of stn's Morse string if(*pC == '\0') // if we've reached end of callsign data pC = t->CallSign; // reset ptr to start of callsign data t->Pitch = (int)*pC++; // audio frequency of keying tone t->Repeat = (int)*pC++; // Nº of times it is to be repeated if((int delay = (int)*pC++) < 0) { // if the delay is negative t->Key = TRUE; // set morse key to the on state delay = -delay; // reverse its sign (make it positive) } else // else delay is zero or positive t->Key = FALSE; // so set morse key to the off state /* While the next call sign character is not a semicolon get pointer to next character's Morse sequence string. */ while((char Letter = *pC) != ';') { pMC = MorseCode + Letter - 65; // For each Morse element (dot or dash) of the character while((MC = *pMC++) > '\0') { *pM++ = MC; // store duration of this Morse element *pM++ = '\1'; // store duration of inter-element pause } *--pM = '\3'; // overwrite with an inter-letter pause pC++; // advance to next letter to be translated } *pM++ = delay; // overwrite with delay at end of callsign *pM = '\0'; // terminating null for Morse Buffer t->pCallSign = ++pC; // start of next segment of callsign } /* A 'C' function which returns the signal attenuation divisor, A, given station frequency and receiver frequency is given below. It calculates the attenuation of a received signal when the receiver is (ω0 − ω) radians per second off tune as effected by 3 cascaded tuned circuits of Q = 100 as would be produced by the receiver's IF amplifier chain. */ double DeTune(long sf, long rf) { // stn freq, rx freq #define L .0000003422 // tuned circuit's inductance in Henries #define R .01 // circuit resistance to give a Q = 100 double w = (double)sf, // Compute the reactance X of single w0 = (double)rf, // tuned circuit. X = L * (w - w0 * w0 / w), Z = X / R, // impedance of single tuned circuit A = sqrt(1 + Z * Z); // attenuation of single tuned circuit return(A * A * A); // attenuation of 3 cascaded circuits } /* The frequency to which the ADF receiver is currently tuned can therefore be calculated as shown in the following 'C' function: */ long ADFfrq( struct RxData *s ) { int BandStart[] = {189, 396, 830}, // kHz BandExtent[] = {221, 474, 990}, // kHz syn = r->Syn, // current position of tuning synchro wb = r->WB, // currently selected waveband hi = r->SynHi[wb], /* synchro position for high frequency limit of the selected band */ lo = r->SynLo[wb]; /* synchro position for low frequency limit of the selected band */ if(syn < lo) syn = lo; if(syn > hi) syn = hi; return((long)(BandStart + ((syn - lo) / (hi - lo)) / BandExtent[wb])); } struct RxData { // DATA STRUCTURE FOR A RECEIVER int Syn; // current position of tuning synchro int SynHi[3]; // synchro position for highest frequency in each band int SynLo[3]; // synchro position for lowest frequency in each band int WB; // currently selected waveband int Sig; // signal strength delivered by receiver's IF amp int Tone[3]; // tone levels delivered to receiver's AF stage input int Valid; // NAV data validation bits long Frq; // frequency to which receiver is currently tuned double ISR; // selected inbound VOR radial or ILS runway heading double OSR; // selected outbound VOR radial double Brg; // bearing of tuned station from receiver double RadDev; // deviation of aircraft from selected radial double StrOff; // required steering offset from station double ReqHdg; // heading the aircraft must fly double DistTD; // ILS distance to touch-down double PrvDst; // dist to/from stn at previous program pass double HdgErr; // horizontal (VOR/ILS) deviation correction command double GsErr; // aircraft's ILS glideslope deviation double DesErr; // difference between required and actual angle of descent double Atten; // signal attenuation divisor due to receiver being off-tune BOOL dFrq; // set TRUE by Receivers() when receiver is re-tuned BOOL capture; // indicates a selected radial has been captured BOOL OnOff; // receiver's on/off switch BOOL CctBrk; // receiver's d/c power supply circuit breaker BOOL SynBrk; // receiver's synchro a/c power supply circuit breaker struct TxData *t; // pointer to the TxData structure of the } // dominant currently on-tune station. *Rx[10]; // Array of pointers to NumRx receiver data structures /* The following function determines the signal strength for all receivers of the type which receives the type of station whose TxData is pointed to by 't'. It also calculates the possible separate audio tone levels for each of the three possible pitches of 1020Hz, 1300Hz and 3000Hz. */ int RxSig(t) { int RxType[] = { // Type of Stn received by Type of Rx 0, // 0 NDB 0 ADF 0, // 1 RRG 0 ADF 1, // 2 VOR 1 NAV 1, // 3 VOR/DME 1 NAV 1, // 4 ILS 1 NAV 1, // 5 ILS/DME 1 NAV 1, // 6 DME 1 NAV 2, // 7 AMK 2 MKR 2, // 8 BMK 2 MKR 2, // 9 OMK 2 MKR 2, // 10 MMK 2 MKR 2 // 11 IMK 2 MKR }; // An array of pointers to pointers to RxData structures struct RxData **pRx[] = { Rx, // points to Rx[] element ptg to RxData for ADF1 Rx + 4, // points to Rx[] element ptg to RxData for NAV1 Rx + 8, // points to Rx[] element ptg to RxData for MKR1 Rx + 10 // points to the element beyond end of Rx array }, *r, // points to RxData struct of rx being dealt with **rx, // points to element of Rx[] containing 'r' **RX, /* This pointer points to the element of Rx[] after the one containing the pointer to last RxData structure for the type of receiver concerned. */ /* The final pointer points to the element within pRx[] (above) which points to the element within Rx[ ] which points to the RxData structure of the first receiver of the type which receives the type of station whose TxData structure is pointed to by 't'. */ ***prx = pRx + *( RxType + t->StnType ); /* FOR EACH RECEIVER OF THIS TYPE, IF Rx[rx] contains a pointer to an exist- ing RxData structure AND if this receiver has been re-tuned since the last station scan, THEN cancel the 'frequency changed' flag and get its detuning attenuation divisor. */ for(rx = *prx, RX = *(prx + 1); rx < RX; rx++) { if(((r = *rx) > NULL)) && ((r->dFrq == TRUE)) { r->dFrq = FALSE; r->Atten = DeTune(t->Frq, r-Frq); } /* Get the signal strength of this station at the aircraft and divide it by the de-tuning attenuation of receiver's passband. Rescale it from a floating point double to an interger ranging from 0 to 255. If the res- ulting signal is the strongest so far encountered during this station scan THEN set it as receiver's current signal strength, and set this station as dominant on-tune station. */ if((Sig = (int)(255 * TxSig(t)/r->Atten)) > r->Sig) { r->Sig = Sig; r->t = t; /* If the station's Morse key is down, and this is the strongest station keying at this pitch, make this the sound level for this pitch of tone. */ if((t->MorseKey == TRUE) && (r->Tone[t->pitch] < Sig)) r->Tone[t->pitch] = Sig; } } } /* This function checks each receiver's on/off switch and dc power circuit breaker and then determines the frequency to which it is currently tuned. If the receiver has been retuned since the last pass, the new frequency is stored and a 'frequency changed' flag is set in the receiver's RxData structure. If the receiver is switched off or its circuit breaker is out, its frequency is set to zero. The signal and tone levels of all receivers are reset to zero ready for the next station scan. */ void RxInput() { long(*GetFrq[])(struct RxData *) = { ADFfrq, ADFfrq, // function returning ADF rx frequency ADFfrq, ADFfrq, NAVfrq, NAVfrq, // function returning NAV rx frequency NAVfrq, NAVfrq, MKRfrq,MKRfrq // function returning Marker frequency }; struct RxData *r; // points to RxData of current rx **rx; // ptr to Rx[] element holding 'r' **RX = Rx + 10; static int r = 0; // element number of Rx[] range 0 to 9 long frq; // current receiver frequency register n = 0; // element number of Rx[] for(rx = Rx; rx <= RX; rx++, n++) { // for all possible receivers if((r = *rx) > NULL) { // if an rx is installed in this 'position' if(r->CctBrk == TRUE // if rx's dc circuit breaker closed && r->OnOff == TRUE) { // and rx is switched on if((frq = (*(GetFrq + rx))(r)) != r->Frq) { // if rx has been retuned r->Frq = frq; // store new frequency r->dFrq = TRUE; // indicate that it has changed } } else r->Frq = 0; // frequency = 0 means rx inoperative } r->Sig = 0; // clear the dominant signal strength for(register x = 0; x < 3; x++) r->Tone[x] = 0; // clear the Morse tone levels } } /* The function below simulates the activity of all active stations by controlling the keying of their callsigns and computing their signal strengths and audio tone levels presented to each receiver. */ void TxScan() { // NAVIGATION STATION SCANNER static double Sigma, InvDst; // For computing height of static int N; // ground under aircraft. static struct TxData **tx; // Pointer to an element in // Tx[ ] pointer array. struct TxData *t; // Pointer to a TxData (ie a // station data) structure. for(tx = Tx; tx <= TxHi; tx++) { // for each element in Tx[] pointer array if(((t = *tx) > NULL) // If it contains a ptr to a TxData struct && (--(t->MorseTimer) < 0)) { // & it's the end of current Morse element t->MorseKey = !(t->MorseKey); // reverse the state of the Morse key. if(*(t->pMorse) == '\0') { // if end of current dots/dashes sequence t->pMorse = t->Morse; // reset ptr to start of Morse Code buffer. if(t->Single == FALSE) // if station has more than 1 callsig && ((--(t->Repeat) < 0)) // avoid possible division by zero SetUpCallSign(t); // get the next part of the call sign } t->MorseTimer = *(t->pMorse)++; // Set timer to duration of next dash/dot } RxSig(t); // compute stn's signal strength for each rx if((t->StnType == 4 // if the station is an ILS || t->StnType == 5) // or an ILS/DME && (t->Dst > 0)) { // and its distance is > zero InvDst += 1 / t->Dst; // add inverse dist Sigma += t->Ht / t->Dst; // add dist/height } } Air.GndHt = Sigma / (N * InvDst); // Compute estimated height of ground } // beneath aircraft /* The following 'C' function simulates the navigation computer's function which computes the required heading. The required heading r->ReqHdg has nothing to do with the aircraft's actual current heading. */ double GetReqHdg(struct TxData *t, struct RxData *r) { #define R 40 // radius of aircraft's minimum turning circle #define RR 1600 // square of this radius double a, // used for intermediate angular quantities b, c, e, // sides of construction triangles see diagram d = t->Dst; // distance between aircraft and waypoint BOOL f = FALSE; // to indicate whether to use + or − root /*If the aircraft has not so far reached the station and it is now closer than (within) the radius of the minimum turning circle and it is now receding from the station then it is deemed now to have passed station. */ if(r->capture == FALSE && d < R && d > t->PrevDst) r->capture = TRUE; /*If the aircraft has not yet reached the station then get its deviation from the selected radial, determine which side of the radial it is, and finally find b and c (see Pre-capture diagram).*/ if(r->capture == FALSE) { r->RadDev = (a = BipAng(t->Brg - r->ISR)); if (a < 0) { a += Pi; f = TRUE; } b = R * cos(a); c = d - R * sin(a); /*If the aircraft is outside the minimum turning circle compute the square root of the appropriate sign, then compute the rationalised required heading it must fly.*/ if((e = b * b + c * c - RR) > 0) { e = sqrt(e); if(f == TRUE) e = -e; r->ReqHdg = RatAng(r->Brg + (r->StrOff = atan(R/e)- atan(b/c))); } } else { // else the aircraft has passed the station // compute a and b (see the Post-capture diagram). r->RadDev = (a = BipAng(t->Brg - r->OSR)); a = r->Brg - a - a - Pi; r->StrOff = (b = BipAng(r->OSR - a)); // Get steering offset if (b > +1) a = r->OSR + 1; // and limit it to +/- if (b < -1) a = r->OSR - 1; // one radian. r->ReqHdg = RatAng(a); } return(r->ReqHdg); // return the aircraft's required heading } /* When adding and subtracting a number of angles, some of which could be compass bearings and some of which could be deviations from reference lines, we can end up with very large angular quantities which can be positive or negative. Many mathematical functions can only accept angles from 0 to 2π or from -π/2 to +π/2. So we need to be able to rationalise and polarise angles. These two tasks are done by RatAng() and BipAng() below which should appear before (or at least prototyped before) GetReqHdg(): Put a wild angle in range 0 to 2π: */ double RatAng(double theta) { while(theta < 0) theta += TwoPi; while(theta > TwoPi) theta -= TwoPi; return(theta); } // Make a wild angle bipolar: double BipAng(double theta) { while(theta < -Pi) theta += TwoPi; while(theta > Pi) theta -= TwoPi; return(theta); } /* Although not used in GetReqHdg(), later on we will also need to be able to make sure that an angle such as a corrective deviation cannot exceed an operating limit. This is done by LimAng(), which for completeness is shown now below: Impose an upper and lower limit on a bipolar angle: */ double LimAng(double theta, double limit) { if(theta < 0) { if(theta < -limit) theta = -limit; } else if(theta > limit) theta = limit; return(theta); } double ElapsedTime() { // elapsed time since previous pass static struct timeb T; // System time at last program pass static double PrevTime = 0; // absolute system time at last pass ElapsedTime = 0; // elapsed time since previous pass double CurrentTime; // absolute system time this pass ftime(&T); // get current system time CurrentTime = T.time + ( (double)T.millitm ) / 1000; if(PrevTime > 0) ElapsedTime = CurrentTime - PrevTime; PrevTime = CurrentTime; // previous time ready for next pass return(ElapsedTime); // return number of seconds +- 1ms } /* The following function stands in place of a full autopilot, airframe, engines and flight dynamics simulation for this type of aircraft. We pass to it the current heading error (HdgErr) and the elapsed time (et) since it was last called, and it returns the aircraft's average rate-of -turn since it was last called. */ double GetRot(double et, double HdgErr) { // get rate of turn #define RollHdgErr 0.03 // Required roll per unit hdg error #define RollMax 0.03 // Max allowed roll angle (radians) #define RollDamp 0.02 // Rate-of-change of roll damping factor // Rate-of-turn per unit roll (radians per second per radian) #define RotRoll 1.00 double dRoll = RollDamp * et * (LimAng(RollHdgErr * HdgErr, RollMax) − Air.Roll); Air.Roll += dRoll; // update aircraft's roll angle // Return average rate of turn since last program pass return(RotRoll * BipAng(Air.Roll - dRoll / 2)); } /* The following function uses this formula to return the appropriate instrument needle deviation (on a -100 - 0 - +100 scale) for a given angular error ranging from -π to +π: Get Deviation Angle, imposing a scaling factor of -100 to +100. */ int GetDev(double theta) { return( (int)( 100 * sin(theta / 2 ) ); } /* The function below, in effect, simulates the action of a basic GPS receiver. Given the aircraft's current position, speed and heading, it computes the aircraft's new position, heading and rate-of-turn based on the time that has elapsed since the previous program pass. Compute the aircraft's heading error and store it as a non-linear instrument deviation: */ void HorzSteer(struct TxData *t, struct RxData *r, double et) { r->HdgErr = GetDev(double HdgErr = BipAng(ReqHdg(t, r) - Air.Hdg)); /*Compute & store aircraft's average rate-of-turn and from that compute its incremental change in heading since last time this function was called.*/ double dHdg = (Air.Rot = GetRot(et, HdgErr)) * et; double AveHdg = Air.Hdg + dHdg / 2; /* Compute aircraft's average heading since last program pass. */ Air.Hdg = RatAng(Air.Hdg + dHdg); // Update the aircraft's actual heading double dDst = Air.Spd * et; // Compute distance aircraft has moved since last program pass. */ double dLat = dDst * cos(AveHdg); /* Compute change in aircraft's latitude since last pass. */ Air.Lat += dLat; // Update aircraft's latitude and longitude Air.Lng += dDst * sin(AveHdg) / cos(Air.Lat + dLat / 2); } /* The following ground height function is illustrative only. The ground height computation is in fact built in as part of the function TxScan(). */ void GndHt( struct TxData *t, // pointer to this station's data structure int n // total number of in-range stations ) { static double Sigma, InvDst; // summation variables static int N; // ILS station count N++; // increment Nº of ILS stns considered if(t->Dst == 0) // if aircraft directly over an ILS stn Air.GndHt = t->Ht; // height of ground = height of station else { // else continue the summations InvDst += 1 / t->Dst; Sigma += t->Ht / t->Dst; Air.GndHt = Sigma / (N * InvDst); } } /* The Gyrocompass is simulated by the following 'C' function (excluding circuit breaker logic): */ void Gyro(double et) { #define T ???? // annunciator's tick-tock constant #define FSD ???? // annunciator's full-scale deflection #define HALFDEG .008725 // ½° in radians #define ERRDEC .002908333 /* 10°/min card error decrement constant in radians/second */ static int AnnDecr = T; /* Use Annunciator tick-tock constant as iterative decrement */ /* If the compass card is within half a degree of indicating the correct magnetic heading, make the annunciator tick-tock, otherwise set it to full scale deflection. */ if((double CardErr = BipAng(Air.CardHdg - Air.Hdg - Air.MagVar)) < HALFDEG) { if(Air.Ann > +FSD) // if annunciator > full scale +ve AnnDecr = -T; // make deflection decrementer -ve else if(Air.Ann < -FSD) // if annunciator < full scale -ve AnnDecr = +T; // make deflection decrementer +ve Air.Ann += AnnDecr * et; // Decrement annunciator by the } // amount of the tick-tock constant. else Ann = FSD; // else set annunciator to full scale deflection double ErrDec = ERRDEC * et; /* The amount by which to decrement the card error this program pass */ if(CardErr > 0) // a +ve error requires a -ve ErrDec = -ErrDec; // decrement and vice versa Air.CardHdg += ErrDec; // advance card heading towards magnetic heading } /* The following data has to be maintained independently for each gyro system. The data structures for the auxiliary gyro (system 0), the captain's main gyro (system 1) and the first officer's main gyro (system 2) are accessed through the array of pointers *Gx[ ]. Data relating to the horizon instruments is in the captain's (NAV1) or the first officer's (NAV2) RxData structure as appropriate. */ struct GyroData { int Status; // stage in the gyro precess cycle (see above) double Pre[2]; // current gyro pitch & roll precess angles BOOL GyroCB; // gyro's 115volt a.c. circuit breaker state BOOL Erect; // TRUE = gyro erected OK, FALSE = gyro failed } *Gx[3]; /* A 'C' function which simulates the gyro's precessing cycle is shown below. This function is entered separately for pitch and for roll for each gyro system. */ double GyroPrecess( struct GyroData *g, // pointer to data struct of Gyro 0, 1, 2 int i, // pitch/roll index: 0 = pitch, 1 = roll double et // elapsed time since last program pass ) { #define INIT .0001 // gyro angle at start of precess cycle #define K 1.5 // precess iterative multiplying factor double pre = *(g->Pre + i); // gyro's current precess angle if(g->Status == 0) { // If this gyro has just been switched on or reset pre = INIT; // set its precess angle to its initial 'seed' value g->Status = 1; // advance the gyro to its 'precessing' state } else // Otherwise if((g->Status == 1) // if this gyro is currently in its precessing cycle && (pre < HalfPi)) { // and it has not yet precessed 90 degrees pre *= (K * et); // multiply its current precess angle by 90 degrees if(pre > HalfPi) { // If, as a result, the precess angle has now over- pre = HalfPi; // shot its 90° limit, set it to its 90° limit and g->Status = 2; // advance the gyro to its 'precess completed' state } } *(g->Pre + i) = pre; // Store and return the updated precess return(pre); // angle ready for next pass } /* The circuit breaker and main/auxiliary gyro switching logic is dealt with by the following function. It returns an index number 'j' which indicates which gyro platform is supplying the attitude information. If we are dealing with the captain's instruments, j can be 0 or 1; if we are deal- ing with the first officer's instruments, j can be 0 or 2 according to whether the pilot concerned is switched to the auxiliary or his main gyro. However, if the gyro is inoperative for some reason, it returns a negative value of 'j' according to the nature of the problem. */ int GyroLogic( int i, // index: 0 for captain's or 1 for first officer's instruments double et // elapsed time since last program pass ) { struct RxData *r = *(Rx + 4 + i); /* Declare pointer to the captain's or first officer's RxData */ if(r->TranCB == TRUE) // If 'transfer switch' circuit breaker closed && (r->TranSw == True)) { // & transfer switch is set to 'auxiliary gyro' r->TranLmp = TRUE; // light the 'transferred to auxiliary' lamp j = 0; // set gyro index = 0 [for auxiliary gyro] } else { // Else r->TranLmp = FALSE; // kill the 'transferred to auxiliary' lamp j = i; // set the gyro's index the same as } // the captain's or first officer's struct GyroData *g = *(Gx + j); /* Declare a pointer to the data structure for the gyro in use. */ if(g->GyroCB == FALSE) // If relevant gyro's circuit breaker is out, set j = -1; // gyro index to show gyro platform inoperative. if(r->SphCB == FALSE // if horizon circuit breaker is out || g->Erect == FALSE) // or gyro platform failed to erect j = -2; // set gyro index to show attitude system failed return(j); // return valid gyro platform index or failure code. } /* The main attitude system function below uses the previous two functions to provide the attitude (pitch and roll) indications on the captain's or first officer's horizon display, and also displays the appropriate instr- ument warning flags and warning lamps. */ void AttitudeSystem ( int i, // index: 0 = captain's or 1 = first officer's instruments double et // elapsed time since last program pass ) { if((int j = GyroLogic(i, et)) < -1) { // IF CURRENT GYRO IS INOPERABLE r->GyroFlag = TRUE; // display its warning flag if(j < -1) // If gyro platform failed r->AttFailLmp = TRUE; // light attitude system fail lamp } else { // ELSE THE GYRO CONCERNED IS WORKING r->AttFailLmp = FALSE; // kill attitude system fail lamp r->Pitch = Air.Pitch; // set sphere pitch to aircraft pitch r->Roll = Air.Roll; // set sphere roll to aircraft roll // If the pilot concerned is holding down 'precess test' button, then if(r->PreReq == TRUE) { r->GyroFlag = TRUE; // show the gyro warning flag r->Pitch += GyroPrecess(g,0,et); // add gyro precess pitch increment r->Roll += GyroPrecess(g,1,et); // add gyro precess roll increment } else r->GyroFlag = FALSE; // else hide warning flag } } double ReqDes(struct TxData *t, struct RxData *r) { #define FlairHeight .0000015 // ~9.55 metres in Earth radians double Hdev, // aircraft's angular displacement from runway centre line ReqDes, // Required angle of descent gs, // runway's glideslope angle GsErr, // linear version of glideslope error RunLen = t->Len; // runway length if(Air.Ht < FlairHeight) { // if aircraft is now below flairout height RunLen *= .75; // advance touchdown point 1/4 way down runway gs = 0; // and set glideslope angle to 0 } else gs = t->GS; // else, use runway's actual glideslope angle Hdev = BipAng(t->Brg - r->ISR); /* Compute aircraft's angular displace- ment from runway centre line. */ /*Provided aircraft's distance along the extended runway centre line from current point of touchdown is non-zero: */ if((double x = (double D = t->Dst) * cos(Hdev)) - RunLen) != 0) { /* Compute the glide slope error as a linear angular quantity and store it in RxData as an appropriately scaled instrument needle deviation: */ r->GsErr = GetDev(GsErr = BipAng(double a = atan(Air.Ht / x)) - gs)); double y = D * sin(Hdev); // compute perpendicular dist. to runway r->DstTD = sqrt(x * x + y * y); // line & distance to touchdown point ReqDes = BipAng(a + GsErr); // required angle of descent } else { // ELSE, DISTANCE TO TOUCHDOWN IS ZERO r->GsErr = 0; // glideslope error indication r->DstTD = 0; // distance to touchdown point ReqDes = 0; // required angle of descent } return(ReqDes); // return required angle of descent } double GetdDes(double et, double DesErr) { #define PitchDesErr .03 // change in pitch per unit angle of descent error #define RopLim .01 // rate of change of pitch limit in radians/second #define MaxPitch .5 // maximum permitted pitch in radians #define DesPitch 1.00 // change in Air.Des per unit change in Air.Pitch dPitch = LimAng(K1 * BipAng(ReqDes(t, r) - Air.Des), RopLim * et); Air.Pitch = LimAng(BipAng(Air.Pitch + dPitch), MaxPitch); return(DesPitch * dPitch); } /* The whole vertical steering task is managed by the following function which, from the required angle of descent and the actual angle of des- cent, updates the aircraft's current height and angle of descent: */ void VertSteer(struct TxData *t, struct RxData, *r, double et) { /*Compute angle of descent error and store it also as a non-linear instrument deviation: */ r->DesErr = GetDev(double DesErr = BipAng(ReqDes(t, r) - Air.Des)); /*Compute the required incremental change in the angle of descent since the previous pass, then update the aircraft's angle of descent and hence its height since previous pass. Note that Air.Spd is the aircraft's horizontal ground track speed. Hence tan(). */ Air.Ht -= Air.Spd * et * tan(Air.Des = RatAng(Air.Des + GetdDes(et, DesErr))); } void GndTrk(HWND hWnd, HDC hDC) { #define X 400 // in logical pixels #define Y 152 // in logical pixels #define TD 140 // distance of touchdown beyond station #define R 40 // Radius of turning circle extern int th; extern HPEN hBlackPen, hWhitePen, hYellowPen; extern HBRUSH hWhiteBrush, hBlueBrush; struct TxData *t = &Tx; // can't pass externally-declared ptrs struct RxData *r = &Rx; // internally-declared ptrs are auto! int SelRad, x, y; // SR in degrees, co-ords of aircraft double a, b, isr, xx, yy, x1, y1, x2, y2; char s[24], ToFrom; TEXTMETRIC tm; GetTextMetrics(hDC, &tm); th = tm.tmExternalLeading + tm.tmHeight; y = 10; strcpy(s,"AIRCRAFT TRACK FOR"); TextOut(hDC, 10, y, s, strlen(s)); y += th; strcpy(s,"AUTOMATIC CAPTURE "); TextOut(hDC, 10, y, s, strlen(s)); y += th; strcpy(s,"OF WAYPOINT RADIAL"); TextOut(hDC, 10, y, s, strlen(s)); y += th; y += th; strcpy(s,"Selected Radial" ); TextOut(hDC, 10, y, s, strlen(s)); y += th; strcpy(s,"Aircraft Bearing"); TextOut(hDC, 10, y, s, strlen(s)); y += th; strcpy(s,"Radial Deviation"); TextOut(hDC, 10, y, s, strlen(s)); y += th; y += th; strcpy(s,"Station Bearing" ); TextOut(hDC, 10, y, s, strlen(s)); y += th; strcpy(s,"Steering Offset" ); TextOut(hDC, 10, y, s, strlen(s)); y += th; strcpy(s,"Required Heading"); TextOut(hDC, 10, y, s, strlen(s)); y += th; y += th; strcpy(s,"Distance"); TextOut(hDC, 10, y, s, strlen(s)); for(SelRad = 0; SelRad <= 360; SelRad++) { // for each of the 360 degrees r->ISR = (isr = SelRad / DegRad); r->OSR = RatAng(isr + Pi); // FORM AND CLEAR THE GROUND TRACK DISPLAY RECTANGLE SelectObject(hDC, hBlackPen); SelectObject(hDC, hBlueBrush); Rectangle(hDC, X - 220, Y - 150, X + 150, Y + 150); // DRAW THE RED & GREEN CIRCLES a = isr + PiBy2; x1 = R * sin(a); y1 = -R * cos(a); a = isr - PiBy2; x2 = R * sin(a); y2 = -R * cos(a); for (b = 0; b <= TwoPi; b += .017453292) { xx = X + R * sin(b); yy = Y + R * cos(b); SetPixel(hDC, x = x1 + xx, y = y1 + yy, RGB(255,0,0)); SetPixel(hDC, x = x2 + xx, y = y2 + yy, RGB(0,255,0)); } // DRAW THE SELECTED RADIAL LINE IN YELLOW SelectObject(hDC, hYellowPen); x = TD * sin(isr); y = -TD * cos(isr); MoveTo(hDC, X - x, Y - y); LineTo(hDC, X + x / 2, Y + y / 2); // DRAW A WHITE BLOB TO REPRESENT THE WAYPOINT SelectObject(hDC, hBlackPen); SelectObject(hDC, hWhiteBrush); Ellipse(hDC, X - 5, Y - 5, X + 5, Y + 5); // PLOT THE AIRCRAFT TRACK xx = -210; yy = 0; // starting position of aircraft t->PrevDst = (t->Dst = sqrt(xx * xx + yy * yy) + 1) + 1; SelectObject(hDC, hWhitePen); MoveTo(hDC, x = X + xx, y = Y + yy);// move to starting posn r->capture = FALSE; // not yet reached the station while(((a = t->Dst) < t->PrevDst && a > R) || a < TD) { t->PrevDst = a; t->Dst = sqrt(xx * xx + yy * yy); //initial dist from stn r->Brg = RatAng((t->Brg = RatAng(atan2(xx, yy))) + Pi); xx += sin(a = GetReqHdg(t, r)); // update aircraft's posn yy += cos(a); LineTo(hDC, x = X + xx, y = Y - yy); // plot its position // on the screen ShowDeg(hDC, 4, r->ISR); ShowDeg(hDC, 5, t->Brg); ShowDeg(hDC, 6, r->RadDev); // Radial Deviation ShowDeg(hDC, 8, r->Brg); ShowDeg(hDC, 9, r->StrOff); ShowDeg(hDC, 10, r->ReqHdg); a = t->Dst; if (r->capture == FALSE) a = -a; ShowDeg(hDC, 12, a); } for(b = 0; b < 3000; b++); // delay } }