ardupilot/libraries/AP_NavEKF3/derivation/main.py

# Copied from https://github.com/PX4/ecl/commit/264c8c4e8681704e4719d0a03b848df8617c0863
# and modified for ArduPilot
from sympy import *
from code_gen import *
import numpy as np

# q: quaternion describing rotation from frame 1 to frame 2
# returns a rotation matrix derived form q which describes the same
# rotation
def quat2Rot(q):
    q0 = q[0]
    q1 = q[1]
    q2 = q[2]
    q3 = q[3]

    # This form is the one normally used in flight dynamics and inertial navigation texts, eg
    # Aircraft Control and Simulation, Stevens,B.L, Lewis,F.L, Johnson,E.N, Third Edition, eqn 1.8-18
    # It does produce second order terms in the covariance prediction that can be problematic
    # with single precision processing.
    # It requires the quternion to be unit length.
    # Rot = Matrix([[q0**2 + q1**2 - q2**2 - q3**2, 2*(q1*q2 - q0*q3), 2*(q1*q3 + q0*q2)],
    #               [2*(q1*q2 + q0*q3), q0**2 - q1**2 + q2**2 - q3**2, 2*(q2*q3 - q0*q1)],
    #                [2*(q1*q3-q0*q2), 2*(q2*q3 + q0*q1), q0**2 - q1**2 - q2**2 + q3**2]])

    # This form removes q1 from the 0,0, q2 from the 1,1 and q3 from the 2,2 entry and results
    # in a covariance prediction that is better conditioned.
    # It requires the quaternion to be unit length and is mathematically identical
    # to the alternate form when q0**2 + q1**2 + q2**2 + q3**2 = 1
    # See https://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm
    Rot = Matrix([[1 - 2*(q2**2 + q3**2), 2*(q1*q2 - q0*q3)    , 2*(q1*q3 + q0*q2)    ],
                 [2*(q1*q2 + q0*q3)     , 1 - 2*(q1**2 + q3**2), 2*(q2*q3 - q0*q1)    ],
                 [2*(q1*q3-q0*q2)       , 2*(q2*q3 + q0*q1)    , 1 - 2*(q1**2 + q2**2)]])

    return Rot

def create_cov_matrix(i, j):
    if j >= i:
        # return Symbol("P(" + str(i) + "," + str(j) + ")", real=True)
        # legacy array format
        return Symbol("P[" + str(i) + "][" + str(j) + "]", real=True)
    else:
        return 0

def create_yaw_estimator_cov_matrix():
    # define a symbolic covariance matrix
    P = Matrix(3,3,create_cov_matrix)

    for index in range(3):
        for j in range(3):
            if index > j:
                P[index,j] = P[j,index]

    return P

def create_Tbs_matrix(i, j):
    # return Symbol("Tbs(" + str(i) + "," + str(j) + ")", real=True)
    # legacy array format
    return Symbol("Tbs[" + str(i) + "][" + str(j) + "]", real=True)

def quat_mult(p,q):
    r = Matrix([p[0] * q[0] - p[1] * q[1] - p[2] * q[2] - p[3] * q[3],
                p[0] * q[1] + p[1] * q[0] + p[2] * q[3] - p[3] * q[2],
                p[0] * q[2] - p[1] * q[3] + p[2] * q[0] + p[3] * q[1],
                p[0] * q[3] + p[1] * q[2] - p[2] * q[1] + p[3] * q[0]])

    return r

def create_symmetric_cov_matrix(n):
    # define a symbolic covariance matrix
    P = Matrix(n,n,create_cov_matrix)

    for index in range(n):
        for j in range(n):
            if index > j:
                P[index,j] = P[j,index]

    return P

# generate equations for observation vector innovation variances
def generate_observation_vector_innovation_variances(P,state,observation,variance,n_obs):
    H = observation.jacobian(state)
    innovation_variance = zeros(n_obs,1)
    for index in range(n_obs):
        H[index,:] = Matrix([observation[index]]).jacobian(state)
        innovation_variance[index] = H[index,:] * P * H[index,:].T + Matrix([variance])

    IV_simple = cse(innovation_variance, symbols("IV0:1000"), optimizations='basic')

    return IV_simple

# generate equations for observation Jacobian and Kalman gain
def generate_observation_equations(P,state,observation,variance,varname="HK"):
    H = Matrix([observation]).jacobian(state)
    innov_var = H * P * H.T + Matrix([variance])
    assert(innov_var.shape[0] == 1)
    assert(innov_var.shape[1] == 1)
    K = P * H.T / innov_var[0,0]
    extension="0:1000"
    var_string = varname+extension
    HK_simple = cse(Matrix([H.transpose(), K]), symbols(var_string), optimizations='basic')

    return HK_simple

# generate equations for observation vector Jacobian and Kalman gain
# n_obs is the vector dimension and must be >= 2
def generate_observation_vector_equations(P,state,observation,variance,n_obs):
    K = zeros(24,n_obs)
    H = observation.jacobian(state)
    HK = zeros(n_obs*48,1)
    for index in range(n_obs):
        H[index,:] = Matrix([observation[index]]).jacobian(state)
        innov_var = H[index,:] * P * H[index,:].T + Matrix([variance])
        assert(innov_var.shape[0] == 1)
        assert(innov_var.shape[1] == 1)
        K[:,index] = P * H[index,:].T / innov_var[0,0]
        HK[index*48:(index+1)*48,0] = Matrix([H[index,:].transpose(), K[:,index]])

    HK_simple = cse(HK, symbols("HK0:1000"), optimizations='basic')

    return HK_simple

# write single observation equations to file
def write_equations_to_file(equations,code_generator_id,n_obs):
    if (n_obs < 1):
        return

    if (n_obs == 1):
        code_generator_id.print_string("Sub Expressions")
        code_generator_id.write_subexpressions(equations[0])
        code_generator_id.print_string("Observation Jacobians")
        code_generator_id.write_matrix(Matrix(equations[1][0][0:24]), "Hfusion", False)
        code_generator_id.print_string("Kalman gains")
        code_generator_id.write_matrix(Matrix(equations[1][0][24:]), "Kfusion", False)
    else:
        code_generator_id.print_string("Sub Expressions")
        code_generator_id.write_subexpressions(equations[0])
        for axis_index in range(n_obs):
            start_index = axis_index*48
            code_generator_id.print_string("Observation Jacobians - axis %i" % axis_index)
            code_generator_id.write_matrix(Matrix(equations[1][0][start_index:start_index+24]), "Hfusion", False)
            code_generator_id.print_string("Kalman gains - axis %i" % axis_index)
            code_generator_id.write_matrix(Matrix(equations[1][0][start_index+24:start_index+48]), "Kfusion", False)

    return

# derive equations for sequential fusion of optical flow measurements
def optical_flow_observation(P,state,R_to_body,vx,vy,vz):
    flow_code_generator = CodeGenerator("./generated/flow_generated.cpp")
    range = symbols("range", real=True) # range from camera focal point to ground along sensor Z axis
    obs_var = symbols("R_LOS", real=True) # optical flow line of sight rate measurement noise variance

    # Define rotation matrix from body to sensor frame
    Tbs = Matrix(3,3,create_Tbs_matrix)

    # Calculate earth relative velocity in a non-rotating sensor frame
    relVelSensor = Tbs * R_to_body * Matrix([vx,vy,vz])

    # Divide by range to get predicted angular LOS rates relative to X and Y
    # axes. Note these are rates in a non-rotating sensor frame
    losRateSensorX = +relVelSensor[1]/range
    losRateSensorY = -relVelSensor[0]/range

    # calculate the observation Jacobian and Kalman gains for the X axis
    equations = generate_observation_equations(P,state,losRateSensorX,obs_var)

    flow_code_generator.print_string("X Axis Equations")
    write_equations_to_file(equations,flow_code_generator,1)

    # calculate the observation Jacobian and Kalman gains for the Y axis
    equations = generate_observation_equations(P,state,losRateSensorY,obs_var)

    flow_code_generator.print_string("Y Axis Equations")
    write_equations_to_file(equations,flow_code_generator,1)

    flow_code_generator.close()

    # calculate a combined result for a possible reduction in operations, but will use more stack
    observation = Matrix([relVelSensor[1]/range,-relVelSensor[0]/range])
    equations = generate_observation_vector_equations(P,state,observation,obs_var,2)
    flow_code_generator_alt = CodeGenerator("./generated/flow_generated_alt.cpp")
    write_equations_to_file(equations,flow_code_generator_alt,2)
    flow_code_generator_alt.close()

    return

# Derive equations for sequential fusion of body frame velocity measurements
def body_frame_velocity_observation(P,state,R_to_body,vx,vy,vz):
    obs_var = symbols("R_VEL", real=True) # measurement noise variance

    # Calculate earth relative velocity in a non-rotating sensor frame
    vel_bf = R_to_body * Matrix([vx,vy,vz])

    vel_bf_code_generator = CodeGenerator("./generated/vel_bf_generated.cpp")
    axes = [0,1,2]
    H_obs = vel_bf.jacobian(state) # observation Jacobians
    K_gain = zeros(24,3)
    for index in axes:
        equations = generate_observation_equations(P,state,vel_bf[index],obs_var)

        vel_bf_code_generator.print_string("axis %i" % index)
        vel_bf_code_generator.write_subexpressions(equations[0])
        vel_bf_code_generator.write_matrix(Matrix(equations[1][0][0:24]), "H_VEL", False)
        vel_bf_code_generator.write_matrix(Matrix(equations[1][0][24:]), "Kfusion", False)

    vel_bf_code_generator.close()

    # calculate a combined result for a possible reduction in operations, but will use more stack
    equations = generate_observation_vector_equations(P,state,vel_bf,obs_var,3)

    vel_bf_code_generator_alt = CodeGenerator("./generated/vel_bf_generated_alt.cpp")
    write_equations_to_file(equations,vel_bf_code_generator_alt,3)
    vel_bf_code_generator_alt.close()

# derive equations for fusion of dual antenna yaw measurement
def gps_yaw_observation(P,state,R_to_body):
    obs_var = symbols("R_YAW", real=True) # measurement noise variance
    ant_yaw = symbols("ant_yaw", real=True) # yaw angle of antenna array axis wrt X body axis

    # define antenna vector in body frame
    ant_vec_bf = Matrix([cos(ant_yaw),sin(ant_yaw),0])

    # rotate into earth frame
    ant_vec_ef = R_to_body.T * ant_vec_bf

    # Calculate the yaw angle from the projection
    observation = atan(ant_vec_ef[1]/ant_vec_ef[0])

    equations = generate_observation_equations(P,state,observation,obs_var)

    gps_yaw_code_generator = CodeGenerator("./generated/gps_yaw_generated.cpp")
    write_equations_to_file(equations,gps_yaw_code_generator,1)
    gps_yaw_code_generator.close()

    return

# derive equations for fusion of declination
def declination_observation(P,state,ix,iy):
    obs_var = symbols("R_DECL", real=True) # measurement noise variance

    # the predicted measurement is the angle wrt magnetic north of the horizontal
    # component of the measured field
    observation = atan(iy/ix)

    equations = generate_observation_equations(P,state,observation,obs_var)

    mag_decl_code_generator = CodeGenerator("./generated/mag_decl_generated.cpp")
    write_equations_to_file(equations,mag_decl_code_generator,1)
    mag_decl_code_generator.close()

    return

# derive equations for fusion of lateral body acceleration (multirotors only)
def body_frame_accel_observation(P,state,R_to_body,vx,vy,vz,wx,wy):
    obs_var = symbols("R_ACC", real=True) # measurement noise variance
    Kaccx = symbols("Kaccx", real=True) # measurement noise variance
    Kaccy = symbols("Kaccy", real=True) # measurement noise variance

    # use relationship between airspeed along the X and Y body axis and the
    # drag to predict the lateral acceleration for a multirotor vehicle type
    # where propulsion forces are generated primarily along the Z body axis

    vrel = R_to_body*Matrix([vx-wx,vy-wy,vz]) # predicted wind relative velocity

    # Use this nonlinear model for the prediction in the implementation only
    # It uses a ballistic coefficient for each axis
    # accXpred = -0.5*rho*vrel[0]*vrel[0]*BCXinv # predicted acceleration measured along X body axis
    # accYpred = -0.5*rho*vrel[1]*vrel[1]*BCYinv # predicted acceleration measured along Y body axis

    # Use a simple viscous drag model for the linear estimator equations
    # Use the the derivative from speed to acceleration averaged across the
    # speed range. This avoids the generation of a dirac function in the derivation
    # The nonlinear equation will be used to calculate the predicted measurement in implementation
    observation = Matrix([-Kaccx*vrel[0],-Kaccy*vrel[1]])

    acc_bf_code_generator  = CodeGenerator("./generated/acc_bf_generated.cpp")
    H = observation.jacobian(state)
    K = zeros(24,2)
    axes = [0,1]
    for index in axes:
        equations = generate_observation_equations(P,state,observation[index],obs_var)
        acc_bf_code_generator.print_string("Axis %i equations" % index)
        write_equations_to_file(equations,acc_bf_code_generator,1)

    acc_bf_code_generator.close()

    # calculate a combined result for a possible reduction in operations, but will use more stack
    equations = generate_observation_vector_equations(P,state,observation,obs_var,2)

    acc_bf_code_generator_alt  = CodeGenerator("./generated/acc_bf_generated_alt.cpp")
    write_equations_to_file(equations,acc_bf_code_generator_alt,3)
    acc_bf_code_generator_alt.close()

    return

# yaw fusion
def yaw_observation(P,state,R_to_earth):
    yaw_code_generator = CodeGenerator("./generated/yaw_generated.cpp")

    # Derive observation Jacobian for fusion of 321 sequence yaw measurement
    # Calculate the yaw (first rotation) angle from the 321 rotation sequence
    # Provide alternative angle that avoids singularity at +-pi/2 yaw
    angMeasA = atan(R_to_earth[1,0]/R_to_earth[0,0])
    H_YAW321_A = Matrix([angMeasA]).jacobian(state)
    H_YAW321_A_simple = cse(H_YAW321_A, symbols('SA0:200'))

    angMeasB = pi/2 - atan(R_to_earth[0,0]/R_to_earth[1,0])
    H_YAW321_B = Matrix([angMeasB]).jacobian(state)
    H_YAW321_B_simple = cse(H_YAW321_B, symbols('SB0:200'))

    yaw_code_generator.print_string("calculate 321 yaw observation matrix - option A")
    yaw_code_generator.write_subexpressions(H_YAW321_A_simple[0])
    yaw_code_generator.write_matrix(Matrix(H_YAW321_A_simple[1]).T, "H_YAW", False)

    yaw_code_generator.print_string("calculate 321 yaw observation matrix - option B")
    yaw_code_generator.write_subexpressions(H_YAW321_B_simple[0])
    yaw_code_generator.write_matrix(Matrix(H_YAW321_B_simple[1]).T, "H_YAW", False)

    # Derive observation Jacobian for fusion of 312 sequence yaw measurement
    # Calculate the yaw (first rotation) angle from an Euler 312 sequence
    # Provide alternative angle that avoids singularity at +-pi/2 yaw
    angMeasA = atan(-R_to_earth[0,1]/R_to_earth[1,1])
    H_YAW312_A = Matrix([angMeasA]).jacobian(state)
    H_YAW312_A_simple = cse(H_YAW312_A, symbols('SA0:200'))

    angMeasB = pi/2 - atan(-R_to_earth[1,1]/R_to_earth[0,1])
    H_YAW312_B = Matrix([angMeasB]).jacobian(state)
    H_YAW312_B_simple = cse(H_YAW312_B, symbols('SB0:200'))

    yaw_code_generator.print_string("calculate 312 yaw observation matrix - option A")
    yaw_code_generator.write_subexpressions(H_YAW312_A_simple[0])
    yaw_code_generator.write_matrix(Matrix(H_YAW312_A_simple[1]).T, "H_YAW", False)

    yaw_code_generator.print_string("calculate 312 yaw observation matrix - option B")
    yaw_code_generator.write_subexpressions(H_YAW312_B_simple[0])
    yaw_code_generator.write_matrix(Matrix(H_YAW312_B_simple[1]).T, "H_YAW", False)

    yaw_code_generator.close()

    return

# 3D magnetometer fusion
def mag_observation_variance(P,state,R_to_body,i,ib):
    obs_var = symbols("R_MAG", real=True)  # magnetometer measurement noise variance

    m_mag = R_to_body * i + ib

    # separate calculation of innovation variance equations for the y and z axes
    m_mag[0]=0
    innov_var_equations = generate_observation_vector_innovation_variances(P,state,m_mag,obs_var,3)
    mag_innov_var_code_generator = CodeGenerator("./generated/3Dmag_innov_var_generated.cpp")
    write_equations_to_file(innov_var_equations,mag_innov_var_code_generator,3)
    mag_innov_var_code_generator.close()

    return

# 3D magnetometer fusion
def mag_observation(P,state,R_to_body,i,ib):
    obs_var = symbols("R_MAG", real=True)  # magnetometer measurement noise variance

    m_mag = R_to_body * i + ib

    # calculate a separate set of equations for each axis
    mag_code_generator = CodeGenerator("./generated/3Dmag_generated.cpp")

    axes = [0,1,2]
    label="HK"
    for index in axes:
        if (index==0):
            label="HKX"
        elif (index==1):
            label="HKY"
        elif (index==2):
            label="HKZ"
        else:
            return
        equations = generate_observation_equations(P,state,m_mag[index],obs_var,varname=label)
        mag_code_generator.print_string("Axis %i equations" % index)
        write_equations_to_file(equations,mag_code_generator,1)

    mag_code_generator.close()

    # calculate a combined set of equations for a possible reduction in operations, but will use slighlty more stack
    equations = generate_observation_vector_equations(P,state,m_mag,obs_var,3)

    mag_code_generator_alt  = CodeGenerator("./generated/3Dmag_generated_alt.cpp")
    write_equations_to_file(equations,mag_code_generator_alt,3)
    mag_code_generator_alt.close()

    return

# airspeed fusion
def tas_observation(P,state,vx,vy,vz,wx,wy):
    obs_var = symbols("R_TAS", real=True) # true airspeed measurement noise variance

    observation = sqrt((vx-wx)*(vx-wx)+(vy-wy)*(vy-wy)+vz*vz)

    equations = generate_observation_equations(P,state,observation,obs_var)

    tas_code_generator = CodeGenerator("./generated/tas_generated.cpp")
    write_equations_to_file(equations,tas_code_generator,1)
    tas_code_generator.close()

    return

# sideslip fusion
def beta_observation(P,state,R_to_body,vx,vy,vz,wx,wy):
    obs_var = symbols("R_BETA", real=True) # sideslip measurement noise variance

    v_rel_ef = Matrix([vx-wx,vy-wy,vz])
    v_rel_bf = R_to_body * v_rel_ef
    observation = v_rel_bf[1]/v_rel_bf[0]

    equations = generate_observation_equations(P,state,observation,obs_var)

    beta_code_generator = CodeGenerator("./generated/beta_generated.cpp")
    write_equations_to_file(equations,beta_code_generator,1)
    beta_code_generator.close()

    return

# yaw estimator prediction and observation code
def yaw_estimator():
    dt = symbols("dt", real=True)  # dt (sec)
    psi = symbols("psi", real=True)  # yaw angle of body frame wrt earth frame
    vn, ve = symbols("vn ve", real=True)  # velocity in world frame (north/east) - m/sec
    daz = symbols("daz", real=True)  # IMU z axis delta angle measurement in body axes - rad
    dazVar = symbols("dazVar", real=True) # IMU Z axis delta angle measurement variance (rad^2)
    dvx, dvy = symbols("dvx dvy", real=True)  # IMU x and y axis delta velocity measurement in body axes - m/sec
    dvxVar, dvyVar = symbols("dvxVar dvyVar", real=True)   # IMU x and y axis delta velocity measurement variance (m/s)^2

    # derive the body to nav direction transformation matrix
    Tbn = Matrix([[cos(psi) , -sin(psi)],
                [sin(psi) ,  cos(psi)]])

    # attitude update equation
    psiNew = psi + daz

    # velocity update equations
    velNew = Matrix([vn,ve]) + Tbn*Matrix([dvx,dvy])

    # Define the state vectors
    stateVector = Matrix([vn,ve,psi])

    # Define vector of process equations
    newStateVector = Matrix([velNew,psiNew])

    # Calculate state transition matrix
    F = newStateVector.jacobian(stateVector)

    # Derive the covariance prediction equations
    # Error growth in the inertial solution is assumed to be driven by 'noise' in the delta angles and
    # velocities, after bias effects have been removed.

    # derive the control(disturbance) influence matrix from IMU noise to state noise
    G = newStateVector.jacobian(Matrix([dvx,dvy,daz]))

    # derive the state error matrix
    distMatrix = Matrix([[dvxVar , 0 , 0],
                        [0 , dvyVar , 0],
                        [0 , 0 , dazVar]])

    Q = G * distMatrix * G.T

    # propagate covariance matrix
    P = create_yaw_estimator_cov_matrix()

    P_new = F * P * F.T + Q

    P_new_simple = cse(P_new, symbols("S0:1000"), optimizations='basic')

    yaw_estimator_covariance_generator = CodeGenerator("./generated/yaw_estimator_covariance_prediction_generated.cpp")
    yaw_estimator_covariance_generator.print_string("Equations for covariance matrix prediction")
    yaw_estimator_covariance_generator.write_subexpressions(P_new_simple[0])
    yaw_estimator_covariance_generator.write_matrix(Matrix(P_new_simple[1]), "_ekf_gsf[model_index].P", True)
    yaw_estimator_covariance_generator.close()

    # derive the covariance update equation for a NE velocity observation
    velObsVar = symbols("velObsVar", real=True) # velocity observation variance (m/s)^2
    H = Matrix([[1,0,0],
                [0,1,0]])

    R = Matrix([[velObsVar , 0],
                [0 , velObsVar]])

    S = H * P * H.T + R
    S_det_inv = 1 / S.det()
    S_inv = S.inv()
    K = (P * H.T) * S_inv
    P_new = P - K * S * K.T

    # optimize code
    t, [S_det_inv_s, S_inv_s, K_s, P_new_s] = cse([S_det_inv, S_inv, K, P_new], symbols("t0:1000"), optimizations='basic')

    yaw_estimator_observation_generator = CodeGenerator("./generated/yaw_estimator_measurement_update_generated.cpp")
    yaw_estimator_observation_generator.print_string("Intermediate variables")
    yaw_estimator_observation_generator.write_subexpressions(t)
    yaw_estimator_observation_generator.print_string("Equations for NE velocity innovation variance's determinante inverse")
    yaw_estimator_observation_generator.write_matrix(Matrix([[S_det_inv_s]]), "_ekf_gsf[model_index].S_det_inverse", False)
    yaw_estimator_observation_generator.print_string("Equations for NE velocity innovation variance inverse")
    yaw_estimator_observation_generator.write_matrix(Matrix(S_inv_s), "_ekf_gsf[model_index].S_inverse", True)
    yaw_estimator_observation_generator.print_string("Equations for NE velocity Kalman gain")
    yaw_estimator_observation_generator.write_matrix(Matrix(K_s), "K", False)
    yaw_estimator_observation_generator.print_string("Equations for covariance matrix update")
    yaw_estimator_observation_generator.write_matrix(Matrix(P_new_s), "_ekf_gsf[model_index].P", True)
    yaw_estimator_observation_generator.close()

def quaternion_error_propagation():
    # define quaternion state vector
    q0, q1, q2, q3 = symbols("q0 q1 q2 q3", real=True)
    q = Matrix([q0, q1, q2, q3])

    # define truth gravity unit vector in body frame
    R_to_earth = quat2Rot(q)
    R_to_body = R_to_earth.T
    gravity_ef = Matrix([0,0,1])
    gravity_bf = R_to_body * gravity_ef

    # define perturbations to quaternion state vector q
    dq0, dq1, dq2, dq3 = symbols("dq0 dq1 dq2 dq3", real=True)
    q_delta = Matrix([dq0, dq1, dq2, dq3])

    # apply perturbations
    q_perturbed = q + q_delta

    # gravity unit vector in body frame after quaternion perturbation
    R_to_earth_perturbed = quat2Rot(q_perturbed)
    R_to_body_perturbed = R_to_earth_perturbed.T
    gravity_bf_perturbed = R_to_body_perturbed * gravity_ef

    # calculate the angular difference between the perturbed and unperturbed body frame gravity unit vectors
    # assuming small angles
    tilt_error_bf = gravity_bf.cross(gravity_bf_perturbed)

    # calculate the derivative of the perturbation rotation vector wrt the quaternion perturbations
    J = tilt_error_bf.jacobian(q_delta)

    # remove second order terms
    # we don't want the error deltas to appear in the final result
    J.subs(dq0,0)
    J.subs(dq1,0)
    J.subs(dq2,0)
    J.subs(dq3,0)

    # define covaraince matrix for quaternion states
    P = create_symmetric_cov_matrix(4)

    # discard off diagonals
    P_diag = diag(P[0,0],P[1,1],P[2,2],P[3,3])

    # rotate quaternion covariances into rotation vector state space
    P_rot_vec = J * P_diag * J.transpose()
    P_rot_vec_simple = cse(P_rot_vec, symbols("PS0:400"), optimizations='basic')

    quat_code_generator = CodeGenerator("./generated/tilt_error_cov_mat_generated.cpp")
    quat_code_generator.write_subexpressions(P_rot_vec_simple[0])
    quat_code_generator.write_matrix(Matrix(P_rot_vec_simple[1]), "tiltErrCovMat", False, "[", "]")
    quat_code_generator.close()

def generate_code():
    print('Starting code generation:')
    print('Creating symbolic variables ...')

    dt = symbols("dt", real=True)  # dt
    g = symbols("g", real=True) # gravity constant

    r_hor_vel = symbols("R_hor_vel", real=True) # horizontal velocity noise variance
    r_ver_vel = symbols("R_vert_vel", real=True) # vertical velocity noise variance
    r_hor_pos = symbols("R_hor_pos", real=True) # horizontal position noise variance

    # inputs, integrated gyro measurements
    # delta angle x y z
    d_ang_x, d_ang_y, d_ang_z = symbols("dax day daz", real=True)  # delta angle x
    d_ang = Matrix([d_ang_x, d_ang_y, d_ang_z])

    # inputs, integrated accelerometer measurements
    # delta velocity x y z
    d_v_x, d_v_y, d_v_z = symbols("dvx dvy dvz", real=True)
    d_v = Matrix([d_v_x, d_v_y,d_v_z])

    u = Matrix([d_ang, d_v])

    # input noise
    d_ang_x_var, d_ang_y_var, d_ang_z_var = symbols("daxVar dayVar dazVar", real=True)

    d_v_x_var, d_v_y_var, d_v_z_var = symbols("dvxVar dvyVar dvzVar", real=True)

    var_u = Matrix.diag(d_ang_x_var, d_ang_y_var, d_ang_z_var, d_v_x_var, d_v_y_var, d_v_z_var)

    # define state vector

    # attitude quaternion
    qw, qx, qy, qz = symbols("q0 q1 q2 q3", real=True)
    q = Matrix([qw,qx,qy,qz])
    R_to_earth = quat2Rot(q)
    R_to_body = R_to_earth.T

    # velocity in NED local frame (north, east, down)
    vx, vy, vz = symbols("vn ve vd", real=True)
    v = Matrix([vx,vy,vz])

    # position in NED local frame (north, east, down)
    px, py, pz = symbols("pn pe pd", real=True)
    p = Matrix([px,py,pz])

    # delta angle bias x y z
    d_ang_bx, d_ang_by, d_ang_bz = symbols("dax_b day_b daz_b", real=True)
    d_ang_b = Matrix([d_ang_bx, d_ang_by, d_ang_bz])
    d_ang_true = d_ang - d_ang_b

    # delta velocity bias x y z
    d_vel_bx, d_vel_by, d_vel_bz = symbols("dvx_b dvy_b dvz_b", real=True)
    d_vel_b = Matrix([d_vel_bx, d_vel_by, d_vel_bz])
    d_vel_true = d_v - d_vel_b

    # earth magnetic field vector x y z
    ix, iy, iz = symbols("magN magE magD", real=True)
    i = Matrix([ix,iy,iz])

    # earth magnetic field bias in body frame
    ibx, iby, ibz = symbols("ibx iby ibz", real=True)

    ib = Matrix([ibx,iby,ibz])

    # wind in local NE frame (north, east)
    wx, wy = symbols("vwn, vwe", real=True)
    w = Matrix([wx,wy])

    # state vector at arbitrary time t
    state = Matrix([q, v, p, d_ang_b, d_vel_b, i, ib, w])

    print('Defining state propagation ...')
    # kinematic processes driven by IMU 'control inputs'
    q_new = quat_mult(q, Matrix([1, 0.5 * d_ang_true[0],  0.5 * d_ang_true[1],  0.5 * d_ang_true[2]]))
    v_new = v + R_to_earth * d_vel_true + Matrix([0,0,g]) * dt
    p_new = p + v * dt

    # static processes
    d_ang_b_new = d_ang_b
    d_vel_b_new = d_vel_b
    i_new = i
    ib_new = ib
    w_new = w

    # predicted state vector at time t + dt
    state_new = Matrix([q_new, v_new, p_new, d_ang_b_new, d_vel_b_new, i_new, ib_new, w_new])

    print('Computing state propagation jacobian ...')
    A = state_new.jacobian(state)
    G = state_new.jacobian(u)

    P = create_symmetric_cov_matrix(24)

    print('Computing covariance propagation ...')
    P_new = A * P * A.T + G * var_u * G.T

    for index in range(24):
        for j in range(24):
            if index > j:
                P_new[index,j] = 0

    print('Simplifying covariance propagation ...')
    P_new_simple = cse(P_new, symbols("PS0:400"), optimizations='basic')

    print('Writing covariance propagation to file ...')
    cov_code_generator = CodeGenerator("./generated/covariance_generated.cpp")
    cov_code_generator.print_string("Equations for covariance matrix prediction, without process noise!")
    cov_code_generator.write_subexpressions(P_new_simple[0])
    cov_code_generator.write_matrix(Matrix(P_new_simple[1]), "nextP", True, "[", "]")

    cov_code_generator.close()


    # derive autocode for other methods
    print('Computing tilt error covariance matrix ...')
    quaternion_error_propagation()
    print('Generating heading observation code ...')
    yaw_observation(P,state,R_to_earth)
    print('Generating gps heading observation code ...')
    gps_yaw_observation(P,state,R_to_body)
    print('Generating mag observation code ...')
    mag_observation_variance(P,state,R_to_body,i,ib)
    mag_observation(P,state,R_to_body,i,ib)
    print('Generating declination observation code ...')
    declination_observation(P,state,ix,iy)
    print('Generating airspeed observation code ...')
    tas_observation(P,state,vx,vy,vz,wx,wy)
    print('Generating sideslip observation code ...')
    beta_observation(P,state,R_to_body,vx,vy,vz,wx,wy)
    print('Generating optical flow observation code ...')
    optical_flow_observation(P,state,R_to_body,vx,vy,vz)
    print('Generating body frame velocity observation code ...')
    body_frame_velocity_observation(P,state,R_to_body,vx,vy,vz)
    print('Generating body frame acceleration observation code ...')
    body_frame_accel_observation(P,state,R_to_body,vx,vy,vz,wx,wy)
    print('Generating yaw estimator code ...')
    yaw_estimator()
    print('Code generation finished!')


if __name__ == "__main__":
    generate_code()