LibrePilot/flight/libraries/insgps_helper.c

#include "ins.h"
#include "pios.h"
#include "ahrs_spi_comm.h"
#include "insgps.h"
#include "CoordinateConversions.h"

#define INSGPS_GPS_TIMEOUT 2   /* 2 seconds triggers reinit of position */
#define INSGPS_GPS_MINSAT  6   /* 2 seconds triggers reinit of position */
#define INSGPS_GPS_MINPDOP 3.5 /* minimum PDOP for postition updates    */
 /* If GPS is more than this distance on any dimension then wait a few updates */
 /* and reinitialize there */
#define INSGPS_GPS_FAR     10 

//! Contains the data from the mag sensor chip
extern struct mag_sensor mag_data;

//! Contains the data from the accelerometer
extern struct accel_sensor  accel_data;

//! Contains the data from the gyro
extern struct gyro_sensor gyro_data;

//! Conains the current estimate of the attitude
extern struct attitude_solution attitude_data;

//! Contains data from the altitude sensor
extern struct altitude_sensor altitude_data;

//! Contains data from the GPS (via the SPI link)
extern struct gps_sensor gps_data;

//! Offset correction of barometric alt, to match gps data
static float baro_offset = 0;

extern void send_calibration(void);
extern void send_attitude(void);
extern void send_velocity(void);
extern void send_position(void);
extern volatile int8_t ahrs_algorithm;
extern void get_accel_gyro_data();
extern void get_mag_data();

/* INS functions */
/**
 * @brief Update the EKF when in outdoor mode.  The primary difference is using the GPS values.
 */
uint32_t total_far_count = 0;
uint32_t relocated = 0;
void ins_outdoor_update()
{
	static uint32_t ins_last_time;
	float gyro[3], accel[3], vel[3];
	float dT;
	uint16_t sensors;
	static uint32_t gps_far_count = 0;
	
	dT = PIOS_DELAY_DiffuS(ins_last_time) / 1e6;
	ins_last_time = PIOS_DELAY_GetRaw();
	
	// This should only happen at start up or at mode switches
	if(dT > 0.01)
		dT = 0.01;
		
	// format data for INS algo
	gyro[0] = gyro_data.filtered.x;
	gyro[1] = gyro_data.filtered.y;
	gyro[2] = gyro_data.filtered.z;
	accel[0] = accel_data.filtered.x,
	accel[1] = accel_data.filtered.y,
	accel[2] = accel_data.filtered.z,
	
	INSStatePrediction(gyro, accel, dT);
	attitude_data.quaternion.q1 = Nav.q[0];
	attitude_data.quaternion.q2 = Nav.q[1];
	attitude_data.quaternion.q3 = Nav.q[2];
	attitude_data.quaternion.q4 = Nav.q[3];
	send_attitude();  // get message out quickly
	INSCovariancePrediction(dT);

	PositionActualData positionActual;
	PositionActualGet(&positionActual);
	positionActual.North = Nav.Pos[0];
	positionActual.East = Nav.Pos[1];
	positionActual.Down = Nav.Pos[2];
	PositionActualSet(&positionActual);
	
	VelocityActualData velocityActual;
	VelocityActualGet(&velocityActual);
	velocityActual.North = Nav.Vel[0];
	velocityActual.East = Nav.Vel[1];
	velocityActual.Down = Nav.Vel[2];
	VelocityActualSet(&velocityActual);

	sensors = 0;
	
	if (gps_data.updated)
	{
		vel[0] = gps_data.groundspeed * cos(gps_data.heading * DEG_TO_RAD);
		vel[1] = gps_data.groundspeed * sin(gps_data.heading * DEG_TO_RAD);
		vel[2] = 0;
		
		if (abs(gps_data.NED[0] - Nav.Pos[0]) > INSGPS_GPS_FAR ||
		    abs(gps_data.NED[1] - Nav.Pos[1]) > INSGPS_GPS_FAR ||
		    abs(gps_data.NED[2] - Nav.Pos[2]) > INSGPS_GPS_FAR ||
		    abs(vel[0] - Nav.Vel[0]) > INSGPS_GPS_FAR ||
		    abs(vel[1] - Nav.Vel[1]) > INSGPS_GPS_FAR) {
			gps_far_count++;
			total_far_count++;
			gps_data.updated = false;
			
			if(gps_far_count > 30) {
				INSPosVelReset(gps_data.NED,vel);
				relocated++;
			}
		}
		else {			
			sensors |= HORIZ_SENSORS | POS_SENSORS;
			
			/*
			 * When using gps need to make sure that barometer is brought into NED frame
			 * we should try and see if the altitude from the home location is good enough
			 * to use for the offset but for now starting with this conservative filter
			 */
			if(fabs(gps_data.NED[2] + (altitude_data.altitude - baro_offset)) > 10) {
				baro_offset = gps_data.NED[2] + altitude_data.altitude;
			} else {
				/* IIR filter with 100 second or so tau to keep them crudely in the same frame */
				baro_offset = baro_offset * 0.999 + (gps_data.NED[2] + altitude_data.altitude) * 0.001;
			}
			gps_data.updated = false;
		}
	}
	
	if(mag_data.updated) {
		sensors |= MAG_SENSORS;
		mag_data.updated = false;
	}
	
	if(altitude_data.updated) {
		sensors |= BARO_SENSOR;
		altitude_data.updated = false;
	}
	
	/*
	 * TODO: Need to add a general sanity check for all the inputs to make sure their kosher
	 * although probably should occur within INS itself
	 */
	INSCorrection(mag_data.scaled.axis, gps_data.NED, vel, altitude_data.altitude - baro_offset, sensors);
	
	if(fabs(Nav.gyro_bias[0]) > 0.1 || fabs(Nav.gyro_bias[1]) > 0.1 || fabs(Nav.gyro_bias[2]) > 0.1) {
		float zeros[3] = {0,0,0};
		INSSetGyroBias(zeros);
	}
}

/**
 * @brief Update the EKF when in indoor mode
 */
void ins_indoor_update()
{
	static uint32_t updated_without_gps = 0;
	
	float gyro[3], accel[3];
	float zeros[3] = {0, 0, 0};
	static uint32_t ins_last_time = 0;
	uint16_t sensors = 0;
	float dT;

	dT = PIOS_DELAY_DiffuS(ins_last_time) / 1e6;
	ins_last_time = PIOS_DELAY_GetRaw();
	
	// This should only happen at start up or at mode switches
	if(dT > 0.01)
		dT = 0.01;

	// format data for INS algo
	gyro[0] = gyro_data.filtered.x;
	gyro[1] = gyro_data.filtered.y;
	gyro[2] = gyro_data.filtered.z;
	accel[0] = accel_data.filtered.x,
	accel[1] = accel_data.filtered.y,
	accel[2] = accel_data.filtered.z,
	
	INSStatePrediction(gyro, accel, dT);
	attitude_data.quaternion.q1 = Nav.q[0];
	attitude_data.quaternion.q2 = Nav.q[1];
	attitude_data.quaternion.q3 = Nav.q[2];
	attitude_data.quaternion.q4 = Nav.q[3];
	send_attitude();  // get message out quickly
	INSCovariancePrediction(dT);
	
	/* Indoors, update with zero position and velocity and high covariance */
	sensors = HORIZ_SENSORS | VERT_SENSORS;
	
	if(mag_data.updated && (ahrs_algorithm == INSSETTINGS_ALGORITHM_INSGPS_INDOOR)) {
		sensors |= MAG_SENSORS;
		mag_data.updated = false;
	}
	
	if(altitude_data.updated) {
		sensors |= BARO_SENSOR;
		altitude_data.updated = false;
	}
	
	if(gps_data.updated) {
		PositionActualData positionActual;
		PositionActualGet(&positionActual);
		positionActual.North = gps_data.NED[0];
		positionActual.East = gps_data.NED[1];
		positionActual.Down = Nav.Pos[2];
		PositionActualSet(&positionActual);
		
		VelocityActualData velocityActual;
		VelocityActualGet(&velocityActual);
		velocityActual.North = gps_data.groundspeed * cos(gps_data.heading * DEG_TO_RAD);
		velocityActual.East = gps_data.groundspeed * sin(gps_data.heading * DEG_TO_RAD);
		velocityActual.Down = Nav.Vel[2];
		VelocityActualSet(&velocityActual);
		
		updated_without_gps = 0;
		gps_data.updated = false;
	} else {
		PositionActualData positionActual;
		PositionActualGet(&positionActual);

		VelocityActualData velocityActual;
		VelocityActualGet(&velocityActual);

		positionActual.Down = Nav.Pos[2];
		velocityActual.Down = Nav.Vel[2];

		if(updated_without_gps > 500) {
			// After 2-3 seconds without a GPS update set velocity estimate to NAN
			positionActual.North = NAN;
			positionActual.East = NAN;
			velocityActual.North = NAN;
			velocityActual.East = NAN;
		} else
			updated_without_gps++;

		PositionActualSet(&positionActual);
		VelocityActualSet(&velocityActual);
	}
	
	/*
	 * TODO: Need to add a general sanity check for all the inputs to make sure their kosher
	 * although probably should occur within INS itself
	 */
	INSCorrection(mag_data.scaled.axis, zeros, zeros, altitude_data.altitude, sensors);
	
	if(fabs(Nav.gyro_bias[0]) > 0.1 || fabs(Nav.gyro_bias[1]) > 0.1 || fabs(Nav.gyro_bias[2]) > 0.1) {
		float zeros[3] = {0,0,0};
		INSSetGyroBias(zeros);
	}

}

/**
 * @brief Initialize the EKF assuming stationary
 */
bool inited = false;
float init_q[4];
void ins_init_algorithm()
{
	inited = true;
	float Rbe[3][3], q[4], accels[3], rpy[3], mag;
	float ge[3]={0,0,-9.81}, zeros[3]={0,0,0}, Pdiag[16]={25,25,25,5,5,5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-4,1e-4,1e-4};
	bool using_mags, using_gps;
	
	INSGPSInit();
	
	HomeLocationData home;
	HomeLocationGet(&home);
	
	accels[0]=accel_data.filtered.x;
	accels[1]=accel_data.filtered.y;
	accels[2]=accel_data.filtered.z;
	
	using_mags = (ahrs_algorithm == INSSETTINGS_ALGORITHM_INSGPS_OUTDOOR) || (ahrs_algorithm == INSSETTINGS_ALGORITHM_INSGPS_INDOOR);
	using_mags &= (home.Be[0] != 0) || (home.Be[1] != 0) || (home.Be[2] != 0);  /* only use mags when valid home location */
	
	using_gps = (ahrs_algorithm == INSSETTINGS_ALGORITHM_INSGPS_OUTDOOR) && (gps_data.quality >= INSGPS_GPS_MINSAT);
	
	/* Block till a data update */
	get_accel_gyro_data();
	
	/* Ensure we get mag data in a timely manner */
	uint16_t fail_count = 50; // 50 at 200 Hz is up to 0.25 sec
	while(using_mags && !mag_data.updated && fail_count--) {
		get_mag_data();
		get_accel_gyro_data();
		AhrsPoll();
		PIOS_DELAY_WaituS(2000);
	}
	using_mags &= mag_data.updated;
	
	if (using_mags) {		
		RotFrom2Vectors(accels, ge, mag_data.scaled.axis, home.Be, Rbe);
		R2Quaternion(Rbe,q);
		if (using_gps)
			INSSetState(gps_data.NED, zeros, q, zeros, zeros);
		else
			INSSetState(zeros, zeros, q, zeros, zeros);
	} else {
		// assume yaw = 0
		mag = VectorMagnitude(accels);
		rpy[1] = asinf(-accels[0]/mag);
		rpy[0] = atan2(accels[1]/mag,accels[2]/mag);
		rpy[2] = 0;
		RPY2Quaternion(rpy,init_q);
		if (using_gps)
			INSSetState(gps_data.NED, zeros, init_q, zeros, zeros);
		else {
			for (uint32_t i = 0; i < 5; i++) {
				INSSetState(zeros, zeros, init_q, zeros, zeros);
				ins_indoor_update();
			}
		}
	}
	
	INSResetP(Pdiag);
	
	// TODO: include initial estimate of gyro bias?
}