- double accelN = _north_accel_node->getDoubleValue();
- double accelE = _east_accel_node->getDoubleValue();
- double accelU = _down_accel_node->getDoubleValue() - 32.0; // why?
-
- // force vector towards magnetic north pole
- double var = _variation_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS;
- double dip = _dip_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS;
- double cosdip = cos(dip);
- double forceN = cosdip * cos(var);
- double forceE = cosdip * sin(var);
- double forceU = sin(dip);
-
- // rotation is around acceleration axis
- // (magnitude doesn't matter)
- double accel = accelN * accelN + accelE * accelE + accelU * accelU;
- if (accel > 1.0)
- accel = sqrt(accel);
- else
- accel = 1.0;
-
- // North marking on compass card
- double edgeN = cos(_error_deg * SGD_DEGREES_TO_RADIANS);
- double edgeE = sin(_error_deg * SGD_DEGREES_TO_RADIANS);
- double edgeU = 0.0;
-
- // apply the force to that edge to get torques
- double torqueN = edgeE * forceU - edgeU * forceE;
- double torqueE = edgeU * forceN - edgeN * forceU;
- double torqueU = edgeN * forceE - edgeE * forceN;
-
- // get the component parallel to the axis
- double torque = (torqueN * accelN +
- torqueE * accelE +
- torqueU * accelU) * 5.0 / accel;
-
- // the compass has angular momentum,
- // so apply a torque and wait
- if (delta_time_sec < 1.0) {
- _rate_degps = _rate_degps * (1.0 - delta_time_sec) - torque;
- _error_deg += delta_time_sec * _rate_degps;
+ /*
+ Formula for northernly turning error from
+ http://williams.best.vwh.net/compass/node4.html:
+
+ Hc: compass heading
+ psi: magnetic heading
+ theta: bank angle (right positive; should be phi here)
+ mu: dip angle (down positive)
+
+ Hc = atan2(sin(Hm)cos(theta)-tan(mu)sin(theta), cos(Hm))
+
+ This function changes the variable names to the more common psi
+ for the heading, theta for the pitch, and phi for the roll (and
+ target_deg for Hc). It also modifies the equation to
+ incorporate pitch as well as roll, as suggested by Chris
+ Metzler.
+ */
+
+ // bank angle (radians)
+ double phi = _roll_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS;
+
+ // pitch angle (radians)
+ double theta = _pitch_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS
+ + _pitch_offset_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS;
+
+ // magnetic heading (radians)
+ double psi = _heading_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS;
+
+ // magnetic dip (radians)
+ double mu = _dip_node->getDoubleValue() * SGD_DEGREES_TO_RADIANS;
+
+
+ /*
+ Tilt adjustments for accelerations.
+
+ The magnitudes of these are totally made up, but in real life,
+ they would depend on the fluid level, the amount of friction,
+ etc. anyway. Basically, the compass float tilts forward for
+ acceleration and backward for deceleration. Tilt about 4
+ degrees (0.07 radians) for every G (32 fps/sec) of
+ acceleration.
+
+ TODO: do something with the vertical acceleration.
+ */
+ double x_accel_g = _x_accel_node->getDoubleValue() / 32;
+ double y_accel_g = _y_accel_node->getDoubleValue() / 32;
+ //double z_accel_g = _z_accel_node->getDoubleValue() / 32;
+
+ theta -= 0.07 * x_accel_g;
+ phi -= 0.07 * y_accel_g;
+
+ ////////////////////////////////////////////////////////////////////
+ // calculate target compass heading degrees
+ ////////////////////////////////////////////////////////////////////
+
+ // these are expensive: don't repeat
+ double sin_phi = sin(phi);
+ double sin_theta = sin(theta);
+ double sin_mu = sin(mu);
+ double cos_theta = cos(theta);
+ double cos_psi = cos(psi);
+ double cos_mu = cos(mu);
+
+ double a = cos(phi) * sin(psi) * cos_mu
+ - sin_phi * cos_theta * sin_mu
+ - sin_phi* sin_theta * cos_mu * cos_psi;
+
+ double b = cos_theta * cos_psi * cos(mu)
+ - sin_theta * sin_mu;
+
+ // This is the value that the compass
+ // is *trying* to display.
+ double target_deg = atan2(a, b) * SGD_RADIANS_TO_DEGREES;
+
+ if( _deviation_node ) {
+ target_deg -= _deviation_node->getDoubleValue();
+ } else if( _deviation_table ) {
+ target_deg -= _deviation_table->interpolate( SGMiscd::normalizePeriodic( 0.0, 360.0, target_deg ) );