C: Directional Gyro drift

There is only a single directional gyro class in FSX and P3D (PLANE HEADING DEGREES GYRO), so if you have multiple directional gyros they will all drift at the same rate. The following code implements a directional gyro class that can be applied to as many directional gyros as the aircraft needs. All will drift at marginally different rates which match real world standards. The code is copy'n'paste; the only section you need to change is below the 'Implementation' header. This has been tested on an aircraft with three directional gyros (captain, first officer and flight engineer) on a four hour flight. It has also been tested using figures from Ed William's Aviation Formulary ([1])

// -----------------------------------------------------------------------
// Helper functions
// -----------------------------------------------------------------------
#define PI_LARGE                 3.1415926535897932384626433832795
#define RADIANS_TO_DEGREE_FACTOR (180.0/PI_LARGE)
#define RAD_TO_DEG(val)          (val)*RADIANS_TO_DEGREE_FACTOR
#define DEG_TO_RAD(val)          (val)/RADIANS_TO_DEGREE_FACTOR

//General to get the current latitude and longitude of the aircraft
double dblAcftCurrLat=0;
double dblAcftCurrLon=0;

// -------------------------------------------------
// Aircraft longitude in radians
double AcftCurrLonRad()
{
  execute_calculator_code("(A:PLANE LONGITUDE,Radians)",&dblAcftCurrLon,NULL,NULL);
  return dblAcftCurrLon;
}
// -------------------------------------------------
// Aircraft latitude in radians
double AcftCurrLatRad()
{
  execute_calculator_code("(A:PLANE LATITUDE,Radians)",&dblAcftCurrLat,NULL,NULL);
  return dblAcftCurrLat;
}
// -------------------------------------------------
// Aircraft longitude in degrees
double AcftCurrLonDeg()
{
  execute_calculator_code("(A:PLANE LONGITUDE,Radians)", &dblAcftCurrLon, NULL, NULL);
  return RAD_TO_DEG(dblAcftCurrLon);
}
// -------------------------------------------------
// Aircraft latitude in degrres
double AcftCurrLatDeg()
{
  execute_calculator_code("(A:PLANE LATITUDE,Radians)", &dblAcftCurrLat, NULL, NULL);
  return RAD_TO_DEG(dblAcftCurrLat);
}

// -----------------------------------------------------------------------
// Drift class explanation
// -----------------------------------------------------------------------
// Definitions of variables
// Real wander (RW)
// Earth rotation (ER)
// Latitude nut correction (LN)
// Transport wander easterly or westerly (TW)
//
// Quadrant sign Table
//         |N. Hem |S. Hem |
//         |-------|-------|
// ER      |   -   |   +   |
// LN      |   +   |   -   |
// TW east |   -   |   +   |
// TW west |   +   |   -   |
//
// Total Drift (TD) = RW + ER + LN + TW
// BUT: if we assume that we're dealing with a perfect gyro then RW = 0, so
// TD = ER + LN + TW
// Period must be in fractions of an hour e.g. 90 mins = 1.5 hours
// Assumption: the latitude nut is always set in the northern hemisphere (positive sign)
// Note also that in this implementation when calculating TW, (x) is already a sine from calculating ER, so TW = y*sin(x) becomes TW = y*x.
// -----------------------------------------------------------------------
class gyro_drift
{
  double current_drift = 0;

public:

  double calc(double lat_nut)
 {
   static double flight_start = 0, prev_lat = 0, prev_lon = 0, prev_lat_deg = 0, prev_lon_deg = 0, flight_progress = 0;
  double x = 0, y = 0, er = 0, ln = 0, tw = 0, curr_lat = 0, curr_lon = 0, curr_lat_deg = 0, curr_lon_deg = 0, flight_time = 0, curr_ew = 0, curr_ns = 0,  travel_dir = 0;
  // Travelling east = 1
  static double prev_ew = 0; // Previous longitude point was east (prev_ew = 1) or west (prev_ew = 0). See also curr_ew
  static double prev_ns = 0; // Previous longitude point was north (prev_ns = 1) or south (prev_ns = 0). See also curr_ns

  // Only update if we are flying
  if (!acft_on_gnd.var_value.n)
  {
   if (!flight_start)
   {
    flight_start = elapsedsecs.var_value.n + 30;
    flight_progress = flight_start;
   }

    // Calculate drift every thirty seconds
   if (elapsedsecs.var_value.n > flight_progress + 30)
   {
    flight_progress = elapsedsecs.var_value.n;
    // Flight time running total in seconds
    flight_time = flight_progress - flight_start;
    // Delay calculation start until the aircraft has been flying for two minutes
    if (flight_time > 120)
    {
     // ----------
     // Current position in radians
     curr_lat = AcftCurrLatRad();
     curr_lon = AcftCurrLonRad();

     // ----------
     // Get the quadrant we are in
     curr_ns = AcftCurrLatDeg();
     curr_lat_deg = curr_ns;    // Needed in degrees for ER calculation
     if (curr_ns > 0)curr_ns = 1;
     curr_ew = AcftCurrLonDeg();
     curr_lon_deg = curr_ew;    // Needed in degrees for TW calculation
     if (curr_ew > 0)curr_ew = 1;

     // ----------
     // Calculate earth rotation - use the mean of the two known points
     // er = 15 x sin of latitude per hour

     x = int(abs(curr_lat_deg) + abs(prev_lat_deg)) / 2;
     x = sin(DEG_TO_RAD(x));
     // Convert flight time to fractions of an hour
     y = minsToDec((int)flight_time / 60);
     // x is already a sin so don't convert it
     er = (15 * x) * y;

     // Negative if in the northern hemisphere
     if (curr_ns)er = -abs(er);

      // ----------
     // Latitude nut error
     // ln = 15 x Sin (Latitude) in degrees per hour
     ln = (15 * sin(lat_nut)*flight_time);
     // Negative if in the southern hemisphere
     if (!curr_ns)ln = -abs(ln);

     // ----------
     // Transport wander
     // Get direction of travel - default travel_dir = 0 i.e. travelling west
     prev_lon_deg = abs(prev_lon_deg);
     curr_lon_deg = abs(curr_lon_deg);
     // Heading east in the eastern hemisphere
     if ((prev_ew && curr_ew) && (curr_lon_deg > prev_lon_deg))
     {
      travel_dir = 1;
      y = DEG_TO_RAD(curr_lon_deg - prev_lon_deg);
     }
     // Heading east in the western hemisphere
     if ((!prev_ew && !curr_ew) && (curr_lon_deg < prev_lon_deg))
     {
      travel_dir = 1;
      y = DEG_TO_RAD(prev_lon_deg - curr_lon_deg);
     }
     // Heading east from the western hemisphere to the eastern hemisphere
     if (!prev_ew && curr_ew)
     {
      travel_dir = 1;
      y = DEG_TO_RAD(prev_lon_deg + curr_lon_deg);
     }
     // Heading west in the eastern hemisphere
     if ((prev_ew && curr_ew) && (curr_lon_deg > prev_lon_deg))
     {
      y = DEG_TO_RAD(prev_lon_deg - curr_lon_deg);
     }
     // Heading west in the western hemisphere
     if ((!prev_ew && !curr_ew) && (curr_lon_deg < prev_lon_deg))
     {
      y = DEG_TO_RAD(curr_lon_deg - prev_lon_deg);
     }
     // Heading west from the eastern hemisphere to the western hemisphere
     if (prev_ew && !curr_ew)
     {
      y = DEG_TO_RAD(prev_lon_deg + curr_lon_deg);
     }

     // Normal calculation is tw = y*sin(x), but x is already a sine (see calculation for er)
     tw = y*x;
     tw = RAD_TO_DEG(tw);
     // Negative if the direction of travel is west to east
     if (travel_dir)tw = -abs(tw);

     // ----------
     // Calculate the current drift
     current_drift = (er)+(ln)+(tw);

     // Done with the calculations, so set curr into prev
     prev_lat = curr_lat;
     prev_lon = curr_lon;
     prev_ns = curr_ns;
     prev_ew = curr_ew;
     prev_lat_deg = curr_lat_deg;
     prev_lon_deg = curr_lon_deg;
    }
   }
  }
  else
  {
   // Reset on landing / flight reload
   flight_start = 0;
   flight_progress = 0;
  }

  // Return drift
  return current_drift;
 }
};

// -----------------------------------------------------------------------
// Directional gyro.
// Calculate the drift to update the display
// Offset is the value assigned by the adjustment knob
// -----------------------------------------------------------------------
class directional_gyro
{
 double dgi_hdg;

public: 

 double get_hdg(double lat_nut, double offset)
 {
  double curr_hdg = 0;

  // Get the current magnetic heading
  curr_hdg = hdg_mag.var_value.n;
  // Instantiate a temporary drift calculation class
  gyro_drift drift;
  // Calculate the current drift
  dgi_hdg = drift.calc(lat_nut);
  // Calculate the current drift offset from the magnetic heading
  dgi_hdg = dgi_hdg + offset + curr_hdg;

  return dgi_hdg;
 }
};

// *************************************************************************
// Implementation
// *************************************************************************

// Instantiate a class per directional gyro
directional_gyro gyro1;
directional_gyro gyro2;
(etc.)

// Randomise a latitude nut setting between 39 degrees north and 56 degrees north per gyro
// Make sure this is only done once during aircraft load, otherwise you'll never get a stable setting
latitude_nut1 = getRand(39, 56);
latitude_nut2 = getRand(39, 56);
(etc.)

// ----------------------------------------------------------
// Directional gyro indicators - 
// ----------------------------------------------------------

// These variables must be global
double gyro1DGAdjustment = 0;  // Gyro1 adjustment knob in degrees 0 - 359
double gyro1DGCard = 0;        // Gyro1 heading display card 0 - 359

double gyro2DGAdjustment = 0;
double gyro2DGCard = 0;
(etc.)

// ----------------------------------------------------------
// Update the directional gyro drift and display
// Called from the _update section of the relevant DG gauge
// ----------------------------------------------------------
void dgUpdate()
{
 double x = 0, y = 0;

 // Captain's directional gyro
 // Current degrees of the directional gyro adjustment knob
 x = gyro1DGAdjustment;
 // Calculate the offset and drift
 y = gyro1.get_hdg(latitude_nut1, x);
 // Correct any overshoot or undershoot in degrees
 if (y < 0)y += 359;
 if (y > 359)y -= 359;
 // Set the card to the corrected heading
 gyro1DGcard = y;

 // First Officer
 x = gyro2DGAdjustment;
 y = gyro2.get_hdg(latitude_nut2, x);
 if (y < 0)y += 359;
 if (y > 359)y -= 359;
 gyro2DGcard = y;

 return;
}