diff src/compass.c @ 604:ca4556fb60b9

bump to 2.99beta, work on 3.00 stable
author heinrichsweikamp
date Thu, 22 Nov 2018 19:47:26 +0100
parents 7d9edd3b8c86
children cd986267a5ca
line wrap: on
line diff
--- a/src/compass.c	Thu Oct 11 21:06:29 2018 +0200
+++ b/src/compass.c	Thu Nov 22 19:47:26 2018 +0100
@@ -1,6 +1,7 @@
 //////////////////////////////////////////////////////////////////////////////
 /// compass.c
-/// Compute north direction from magnetic/acceleration measures.
+/// Compute north direction from magnetic/acceleration measures
+///
 /// Copyright (c) 2012-2015, JD Gascuel, HeinrichsWeikamp, all right reserved.
 //////////////////////////////////////////////////////////////////////////////
 // HISTORY
@@ -12,189 +13,197 @@
 #include "compass.h"
 
 //////////////////////////////////////////////////////////////////////////////
-// mH: Crude work-around, needs to be made right
+// mH: crude work-around, needs to be fixed up
 #ifndef UNIX
-#   pragma udata overlay bank8=0x800
-    static char	  C_STACK[256];  // Overlay C-code data stack here.
-#       define RESET_C_STACK    \
-        _asm                    \
-            LFSR    1, 0x800    \
-            LFSR    2, 0x800    \
-        _endasm
-#   pragma udata overlay bank9_compass
-#   pragma code compass_run
+#	pragma udata overlay bank8=0x800
+	static char C_STACK[256];			// overlay C-code data stack here
+#	define RESET_C_STACK		\
+		_asm					\
+			LFSR	1, 0x800	\
+			LFSR	2, 0x800	\
+		_endasm
+#	pragma udata overlay bank9_compass
+#	pragma code compass_run
 #else
-#   define RESET_C_STACK
+#	define RESET_C_STACK
 #endif
 
 //////////////////////////////////////////////////////////////////////////////
 // fifth order of polynomial approximation of atan(), giving 0.05 deg max error
 //
-#define K1  (5701)  // Needs K1/2**16
-#define K2  (1645)  // Needs K2/2**48  WAS NEGATIV
-#define K3  ( 446)  // Needs K3/2**80
+#define K1	(5701)			// needs K1/2**16
+#define K2	(1645)			// needs K2/2**48  WAS NEGATIV
+#define K3	( 446)			// needs K3/2**80
+
 
 //////////////////////////////////////////////////////////////////////////////
-// Interface to assembleur multiplies
+// Interface to assembler multiplies
+
 Int16 umul(PARAMETER Int16 a, PARAMETER Int16 b)
 {
-    extern Int16 compass_umul(void);
-    extern Int16 compass_a, compass_b;
-    compass_a = a;
-    compass_b = b;
-    return compass_umul();
+	extern Int16 compass_umul(void);
+	extern Int16 compass_a, compass_b;
+	compass_a = a;
+	compass_b = b;
+	return compass_umul();
 }
 
 Int16 imul(PARAMETER Int16 a, PARAMETER Int16 b)
 {
-    extern Int16 compass_imul(void);
-    extern Int16 compass_a, compass_b;
-    compass_a = a;
-    compass_b = b;
-    return compass_imul();
+	extern Int16 compass_imul(void);
+	extern Int16 compass_a, compass_b;
+	compass_a = a;
+	compass_b = b;
+	return compass_imul();
 }
 
 //////////////////////////////////////////////////////////////////////////////
 /// Returns a / b * 2**16
 ///
-/// A 16/16 -> 16 bits divide, returning a scalled result.
+/// A 16/16 -> 16 bits divide, returning a scaled result.
 /// Used to multiply fractional numbers in the range 0..1,
 /// represented as 0..32767.
+
 Int16 udiv(PARAMETER Int16 a, PARAMETER Int16 b)
 {
-    OVERLAY Int16 d, r;
+	OVERLAY Int16 d, r;
 
-    //---- Pre-scale both numerator and denominator --------------------------
-    while( (((a>>8) | (b>>8)) & 0xC0) == 0 )
-    {
-        a <<= 1;
-        b <<= 1;
-    }
+	//---- Pre-scale both numerator and denominator --------------------------
+	while( (((a>>8) | (b>>8)) & 0xC0) == 0 )
+	{
+		a <<= 1;
+		b <<= 1;
+	}
 
-    //---- Make division trials ----------------------------------------------
-    d = 0x4000;             // Starts with 0.5, because 1.0 is sign bit.
-    b >>= 1;                // Hence pre-shift b.
-    r = 0;
-    do {
-        if( a >= b ) {      // a is big enough ?
-            a -= b;         // then count d times b out of it.
-            r |= d;         // and accumulate that bit.
-        }
-        b >>= 1;            // then loop trying twice smaller.
-        d >>= 1;
-    } while( b );
-    return r;
+	//---- Make division trials ----------------------------------------------
+	d = 0x4000;				// start with 0.5, because 1.0 is sign bit
+	b >>= 1;				// hence pre-shift b
+	r = 0;
+	do {
+		if( a >= b )		// a is big enough ?
+		{
+			a -= b;			// then count d times b out of it
+			r |= d;			// and accumulate that bit
+		}
+		b >>= 1;			// then loop trying twice smaller
+		d >>= 1;
+	} while( b );
+	return r;
 }
 
 //////////////////////////////////////////////////////////////////////////////
 /// Computes atan(y/x) in Angle, for x, y in range 0..32767
 ///
 /// Results a single quadrant Angle, in the range 0 .. Q_PI/2
+
 Angle utan(PARAMETER Int16 y, PARAMETER Int16 x)
 {
-    OVERLAY Int16 ratio, angle, x2, x3;
+	OVERLAY Int16 ratio, angle, x2, x3;
 
-    //---- Handle zero divisor -----------------------------------------------
-    if( x == 0 )
-        return (y == 0) ? 0 : Q_PIO2;
+	//---- Handle zero divisor -----------------------------------------------
+	if( x == 0 )
+		return (y == 0) ? 0 : Q_PIO2;
 
-    //---- Make it half-quadrant : 0 .. 45 deg -------------------------------
-    ratio = (x > y) ? udiv(y, x) : udiv(x, y);
+	//---- Make it half-quadrant : 0 .. 45 deg -------------------------------
+	ratio = (x > y) ? udiv(y, x) : udiv(x, y);
 
-    //---- Then apply the polynomial approximation ---------------------------
-    angle = umul(K1, ratio);            // r*K1 / 2**16
-    x2 = umul(ratio, ratio);            // r**2 / 2**16
-    x3 = umul(x2, ratio);               // r**3 / 2**32
-    angle -= umul(x3, K2);              // K2*r**3 / 2**48: NEGATIV.
+	//---- Then apply the polynomial approximation ---------------------------
+	angle = umul(K1, ratio);			// r*K1 / 2**16
+	x2 = umul(ratio, ratio);			// r**2 / 2**16
+	x3 = umul(x2, ratio);				// r**3 / 2**32
+	angle -= umul(x3, K2);				// K2*r**3 / 2**48: NEGATIV.
 
-    x3 = umul(x3, x2);                  // r**5 / 2**64
-    angle += umul(x3, K3);              // K3*r**5 / 2**80
+	x3 = umul(x3, x2);					// r**5 / 2**64
+	angle += umul(x3, K3);				// K3*r**5 / 2**80
 
-    //---- Recover the full quadrant -----------------------------------------
-    return (x < y) ? (Angle)(Q_PIO2 - angle)
-                   : (Angle)(angle);
+	//---- Recover the full quadrant -----------------------------------------
+	return (x < y) ? (Angle)(Q_PIO2 - angle)
+				   : (Angle)(angle);
 }
 
 //////////////////////////////////////////////////////////////////////////////
 /// Computes atan2(y/x) in Angle, for x, y in range -32768 to 32767
 ///
 /// Results a four quadrant Angle, in the range -Q_PI .. +Q_PI
+
 Angle itan(PARAMETER Int16 y, PARAMETER Int16 x)
 {
-    // Beware: -32768 is not properly handled (sgn error).
-    if( x == -32768 ) x = -32767;
-    if( y == -32768 ) y = -32767;
+	// Beware: -32768 is not properly handled (sign error)
+	if( x == -32768 ) x = -32767;
+	if( y == -32768 ) y = -32767;
 
-    if( x >= 0 )
-        if( y >= 0 )    // First quadrant: 0..90 deg.
-            return utan(y,x);
-        else            // Fourth quadrant: 0..-90 deg
-            return -utan(-y,x);
-    else
-        if( y >= 0 )    // Second quadrant: 90..180 deg
-            return Q_PI - utan(y, -x);
-        else            // Third quadrant: -90..-180 deg;
-            return -Q_PI + utan(-y, -x);
+	if( x >= 0 )
+		if( y >= 0 )					// first  quadrant: 0..90 deg.
+			return utan(y,x);
+		else							// fourth quadrant: 0..-90 deg
+			return -utan(-y,x);
+	else
+		if( y >= 0 )					// second quadrant: 90..180 deg
+			return Q_PI - utan(y, -x);
+		else							// third  quadrant: -90..-180 deg;
+			return -Q_PI + utan(-y, -x);
 }
 
 //////////////////////////////////////////////////////////////////////////////
 /// Computes cos(theta) = sqrtf(x2/h2),
 /// when theta = atan(y/x) and h2=x*x+y*y
 ///
+
 Int16 cosxh(PARAMETER Int16 x2, PARAMETER Int16 h2)
 {
-    OVERLAY Int16 r = 0;
-    OVERLAY Int16 d = 0x4000;
+	OVERLAY Int16 r = 0;
+	OVERLAY Int16 d = 0x4000;
 
-    do {
-        OVERLAY Int16 a = r + d;
-        a = umul(a, a);
-        a = umul(a, h2);
-        if( a <= x2 ) r += d;
-        d >>= 1;
-    } while( d );
+	do {
+		OVERLAY Int16 a = r + d;
+		a = umul(a, a);
+		a = umul(a, h2);
+		if( a <= x2 ) r += d;
+		d >>= 1;
+	} while( d );
 
-    return r;
+	return r;
 }
 
 //////////////////////////////////////////////////////////////////////////////
 /// Computes both sin and cos of angle y/x,
 /// with h = sqrt(x**2+y**2).
 ///
+
 void sincos(PARAMETER Int16 x, PARAMETER Int16 y, Int16* sin, Int16* cos)
 {
-    OVERLAY Int16 x2, y2, h2;
+	OVERLAY Int16 x2, y2, h2;
 
-    //---- Fold into one quadant ---------------------------------------------
-    OVERLAY char neg = 0;
-    if( x < 0 )
-    {
-        neg |= 1;
-        x = -x;
-    }
-    if( y < 0 )
-    {
-        neg |= 2;
-        y = -y;
-    }
+	//---- Fold into one quadant ---------------------------------------------
+	OVERLAY char neg = 0;
+	if( x < 0 )
+	{
+		neg |= 1;
+		x = -x;
+	}
+	if( y < 0 )
+	{
+		neg |= 2;
+		y = -y;
+	}
 
-    //---- Pre-scale both numerator and denominator ----------------------
-    while( (((x>>8) | (y>>8)) & 0xE0) == 0 )
-    {
-        x <<= 1;
-        y <<= 1;
-    }
+	//---- Pre-scale both numerator and denominator ----------------------
+	while( (((x>>8) | (y>>8)) & 0xE0) == 0 )
+	{
+		x <<= 1;
+		y <<= 1;
+	}
 
-    //---- Uses trig() to do the stuff one on quadrant -------------------
-    x2 = umul(x,x);
-    y2 = umul(y,y);
-    h2 = x2 + y2;
-    x2 = cosxh(x2, h2);
+	//---- Uses trig() to do the stuff one on quadrant -------------------
+	x2 = umul(x,x);
+	y2 = umul(y,y);
+	h2 = x2 + y2;
+	x2 = cosxh(x2, h2);
 
-    //---- Results back in four quadrants --------------------------------
-    *cos = (neg & 1) ? -x2 : x2;
-    y2 = cosxh(y2, h2);
-    *sin = (neg & 2) ? -y2 : y2;
+	//---- Results back in four quadrants --------------------------------
+	*cos = (neg & 1) ? -x2 : x2;
+	y2   = cosxh(y2, h2);
+	*sin = (neg & 2) ? -y2 : y2;
 }
 
 //////////////////////////////////////////////////////////////////////////////
@@ -202,47 +211,48 @@
 
 void compass(void)
 {
-    OVERLAY Int16 sin, cos;
-    OVERLAY Int16 iBfx, iBfy, Gz;
-    OVERLAY Int16 iBpx, iBpy, iBpz;
-    RESET_C_STACK;
+	OVERLAY Int16 sin, cos;
+	OVERLAY Int16 iBfx, iBfy, Gz;
+	OVERLAY Int16 iBpx, iBpy, iBpz;
+
+	RESET_C_STACK;
 
-    //---- Make hard iron correction -----------------------------------------
-    // Measured magnetometer orientation, measured ok.
-    // From matthias drawing: (X,Y,Z) --> (X,Y,Z) : no rotation.
-    iBpx = compass_DX_f - compass_CX_f; // X
-    iBpy = compass_DY_f - compass_CY_f; // Y
-    iBpz = compass_DZ_f - compass_CZ_f; // Z
+	//---- Make hard iron correction -----------------------------------------
+	// Measured magnetometer orientation, measured ok
+	// From matthias drawing: (X,Y,Z) --> (X,Y,Z) : no rotation
+	iBpx = compass_DX_f - compass_CX_f; // X
+	iBpy = compass_DY_f - compass_CY_f; // Y
+	iBpz = compass_DZ_f - compass_CZ_f; // Z
 
-    //---- Calculate sine and cosine of roll angle Phi -----------------------
-    sincos(accel_DZ_f, accel_DY_f, &sin, &cos);
+	//---- Calculate sine and cosine of roll angle Phi -----------------------
+	sincos(accel_DZ_f, accel_DY_f, &sin, &cos);
 
-    //---- rotate by roll angle (-Phi) ---------------------------------------
-    iBfy = imul(iBpy, cos) - imul(iBpz, sin);
-    iBpz = imul(iBpy, sin) + imul(iBpz, cos);
-    Gz = imul(accel_DY_f, sin) + imul(accel_DZ_f, cos);
+	//---- rotate by roll angle (-Phi) ---------------------------------------
+	iBfy = imul(iBpy, cos) - imul(iBpz, sin);
+	iBpz = imul(iBpy, sin) + imul(iBpz, cos);
+	Gz = imul(accel_DY_f, sin) + imul(accel_DZ_f, cos);
 
-    //---- calculate sin and cosine of pitch angle Theta ---------------------
-    sincos(Gz, -accel_DX_f, &sin, &cos);     // NOTE: changed sin sign.
+	//---- calculate sin and cosine of pitch angle Theta ---------------------
+	sincos(Gz, -accel_DX_f, &sin, &cos);	 // NOTE: changed sin sign
 
-    /* correct cosine if pitch not in range -90 to 90 degrees */
-    if( cos < 0 ) cos = -cos;
+	/* correct cosine if pitch not in range -90 to 90 degrees */
+	if( cos < 0 ) cos = -cos;
 
-    ///---- de-rotate by pitch angle Theta -----------------------------------
-    iBfx = imul(iBpx, cos) + imul(iBpz, sin);
+	///---- de-rotate by pitch angle Theta -----------------------------------
+	iBfx = imul(iBpx, cos) + imul(iBpz, sin);
 
-    //---- Detect uncalibrated compass ---------------------------------------
-    if( !compass_CX_f && !compass_CY_f && !compass_CZ_f )
-    {
-        compass_heading = -1;
-        return;
-    }
+	//---- Detect uncalibrated compass ---------------------------------------
+	if( !compass_CX_f && !compass_CY_f && !compass_CZ_f )
+	{
+		compass_heading = -1;
+		return;
+	}
 
-    //---- calculate current yaw = e-compass angle Psi -----------------------
-    // Result in degree (no need of 0.01 deg precision...
-    compass_heading = itan(-iBfy, iBfx) / 100;
+	//---- calculate current yaw = e-compass angle Psi -----------------------
+	// Result in degree (no need of 0.01 deg precision...
+	compass_heading = itan(-iBfy, iBfx) / 100;
 
-    // Result in 0..360 range:
-    if( compass_heading < 0 )
-        compass_heading += 360;
+	// Result in 0..360 range:
+	if( compass_heading < 0 )
+		compass_heading += 360;
 }