GNC Unit and Integration Test Documentation - spacecraft-design-lab-2019/documentation GitHub Wiki

Detumble Algorithms - detumble_algorithms.cpp

• detumble_B_cross(Vector3d omega, Vector3d B, double k) -> Vector3d

  • Calculates torque for detumbling a spacecraft based on the B_cross algorithm. Returns a 3x1 vector in the body frame [N-m] opposing the given angular rate vector.

  • Inputs:

    • omega - angular rate; the rotation rate of the spacecraft in body coordinates [radians/second].
    • B - magnetic field; best estimate of the Earth's magnetic field in body coordinates [nanoTesla].
    • k - B_cross gain; gain relating angular rate to the appropriate torque. Units work out to [kg-m^2].

  • Test: verify_B_cross.py

  • References: Avanzini and Giulietti

  • TODO:

    • Implement dipole saturation effects on output torque. (Paul DeTrempe, 12/4/2019).
  • Find the optimal gain, k, for our system. This may end up being a bang-bang control law (i.e. maximally actuate the dipoles in the direction given by the output torque). ( Paul DeTrempe, 12/4/2019).
    • Create and use matrix describing dipole arrangement in torque calculation. (Paul DeTrempe, 12/4/2019).

---

• detumble_B_dot_bang_bang(Vector3d B_dot, Vector3d max_dipoles) -> Vector3d

  • Calculates torque for detumbling a spacecraft based on the B_dot algorithm. Returns the command moment vector in the body frame [Amp-meters^2].

  • Inputs:

    • B_dot - the time rate of change of the magnetic field, in spacecraft body coordinates, [nanoTesla/second]. It should be noted that nanoTesla/second are listed as the units here, but units of the measured magnetic field don't matter, as only the sign of the rate of change of the magnetic field is used in the bang-bang control law.
    • max_dipoles - the maximum possible value of the magnetorquer dipole moments in each body coordinate direction [Amp-meters^2]

  • Test Script/Function: example_orbit_prop_detumble.py

  • References: Wertz, Fundamentals of Spacecraft Attitude Determination and Control, Equation 7.54

---

• get_bias_estimate(MatrixXd B_mat) -> Vector3d

  • Estimates the magnetometer bias using least squares ellipsoid fit over a time history of the magnetic field.

  • Inputs:

    • B_mat - an m x 3 matrix of the time history of the magnetic field. Each row (1x3) should represent a measurement of the magnetic field as measured in the body frame [nanoTesla]. It should be noted that while nanoTesla are the units specified here, any units can be entered and the bias will be returned in the same units.

  • Test: test_detumble_algorithms.py, verify_magnetometer_bias_estimation.py

  • References:

    • Ellipsoid fit, MATLAB stack exchange: https://www.mathworks.com/matlabcentral/fileexchange/24693-ellipsoid-fit
    • LEAST SQUARES FITTING OF ELLIPSOID USING ORTHOGONAL DISTANCES: http://dx.doi.org/10.1590/S1982-21702015000200019 (Note, this function uses the non-iterative method described in section 4 of the listed paper).

  • TODO:

    • Add ability to optionally return orientation and dimensions of ellipsoid for estimation of T matrix (combined scaling and misallignment errors). This was left out as it is assumed that T matrix estimation can be performed before launch and will not change on orbit. (Paul DeTrempe 12/4/2019).

---

• get_B_dot(Vector3d B1, Vector3d B2, double dt) -> Vector3d

  • Estimates the time rate of change of the magnetic field in body coordinates.

  • Inputs:

    • B1 - magnetic field in the spacecraft body frame [nanoTesla]. This is the first measurement sequentially (i.e. comes before B2). Units are noted as nanoTesla here, but the function will return a rate of change with dimensions given by those of the magnetic field measurements divided by the units of the time interval [B]/[dt]
    • B2 - The second measurement sequentially (i.e. comes after B1).
    • dt - scalar giving the time interval between the two measurements [seconds].

---

---

Orbit Propagation - orbit.py

• get_orbit_pos

• get_orbit_state

• get_orbit_magnetic

• get_B_field_at_point

---

---

Attitude Dynamics & Kinematics - euler_cpp.cpp

• get_w_dot(Vector3d w, Vector3d M, MatrixXd I) -> Vector3d

  • Takes in applied torque M [Nm], inertia matrix I [kg*m^2], and the current angular velocity w [rad/s] to find the angular acceleration in body frame [rad/s^2].

  • Test: test_euler_quat.py

---

• get_q_dot(Vector4d q, Vector3d w) -> Vector4d

  • Takes in a quaternion q (scalar first) representing a rotation from body to inertial frame, and the current angular velocity w in the body frame. Returns the derivative of the quaternion.

  • Test: test_euler_quat.py

  • References: Zac's Notes

---

• get_attitude_derivative(VectorXd x, Vector3d M, MatrixXd I) -> VectorXd

  • Uses the current state, applied torques, and inertia matrix to calculate the derivative of the state. Returns [q_dot; w_dot]

  • Inputs:

    • x - [q; w]. q is body to ECI, w in the body frame.
    • M - applied torques about the spacecraft in the body frame (Nm)
    • I - inertia matrix about the rotation point (if the CM, I should be diagonal) (kg*m^2)

  • Test: test_euler_quat.py

  • References: Zac's Notes and Wertz

---

• Lq(Vector4d q) -> MatrixXd

  • Creates the left multiply matrix for quaternion multiplication.

  • Inputs:

    • q - scalar-first quaternion

  • Test: test_euler_quat.py

  • References: Zac's Notes

---

• Rq(Vector4d q) -> MatrixXd

  • Creates the right multiply matrix for quaternion multiplication.

  • Inputs:

    • q - scalar-first quaternion

  • Test: test_euler_quat.py

  • References: Zac's Notes

---

• hat(Vector3d w) -> MatrixXd

  • Create the skew symmetric matrix used for the cross product.

  • Test: Test

---

• rotate_vec(Vector3d xb, Vector4d q) -> Vector3d

  • This function takes in a quaternion and a vector (3x1) in the BODY frame and rotates it to the INERTIAL frame

  • Inputs:

    • q - scalar first quaternion, representing a rotation from BODY to INERTIAL
    • xb - vector in BODY frame

  • Test: Test

  • References: Zac's Notes

---

• get_inverse_quaternion(Vector4d q) -> Vector4d

  • Returns the inverse of a scalar-first quaternion q. The inverse represents a rotation equal and opposite to q.

  • Test: Test

  • References: Zac's Notes

---

---

Magnetic Field Models - magnetic_field.cpp

• get_magnetic_field(double lat, double lon, double alt, double year, int order) -> VectorXd

  • Return B vector in NED.

  • Inputs:

    • lat - latitude (east)
    • lon - longitude
    • alt - altitude (in km)
    • year - (must be at least 2015 and can use decimals)
    • order (degree of model)

  • Test: test_mag_field.py

  • References:

    • MATLAB inspiration for model: https://www.mathworks.com/matlabcentral/fileexchange/34388-international-geomagnetic-reference-field-igrf-model?s_tid=mwa_osa_a

  • TODO:

    • should return Vector3d

---

• get_P_coefficients(x, order)

  • Inputs:

    • x = cos(co-latitude)
    • order

  • Outputs:

    • legendre polynomials

---

• get_Pd_coefficients(P, x, order)

  • Inputs:

    • P - legengre polynomials
    • x = cos(co-latitude)
    • order

  • Outputs:

    • derivative of legendre polynomials

---

• Other functions:
  • get_g_coefficients, get_h_coefficients, get_g_sv_coefficients, get_h_sv_coefficients

  • Inputs:

    • None

  • Outputs:

    • igrf coefficients up to order 10

    • References:

      • Can download excel spreadsheet here: https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html

---

---

Sun Functions - sun_utils.cpp

• sun_position(double MJD) -> VectorXd

  • From MJD, computes the position of the Sun in ECI.

  • Test: Test

  • References:

    • Montenbruck and Gill

---

• sat_sun_vect(VectorXd r, double MJD) -> VectorXd

  • This function takes in the position of the satellite in ECI and a MJD to return the unit vector from the satellite to the Sun in ECI.

  • Test: Test

  • References:

    • Montenbruck and Gill

---

---

TRIAD - deterministic_ad.cpp

• triad_ad(MatrixXd M, MatrixXd V) -> MatrixXd

  • From 3xn measurements matrix (in BODY frame), M, and 3xn modelled vectors matrix (in INERTIAL frame), V, computes the rotation matrix from INERTIAL to BODY.

  • Vectors in M and V should be unit vectors and columns should correspond.

    • Example: Column 1 in M is the unit vector to the Sun given by the Sun Sensor in the BODY frame and Column 1 in V should be the predicted unit vector to the Sun in the INERTIAL frame.

  • Test: Test

  • References:

    • 279C Notes

---

---

Frame Conversion Functions - frame_conversions.cpp

• eci2ecef(double GMST) -> MatrixXd

  • Rotation matrix from ECI to ECEF coordinates given the GMST (0, 2π)

  • Test: Test

  • References: Rotation matrix about z-axis

---

• ecef2enu(double lat, double lon) -> MatrixXd

  • Rotation matrix from ECEF to ENU coordinates given the lat and lon of the satellite (in radians)

  • Test: Test

---

• ecef2lla(Vector3d r) -> (double, double, double)

  • Gives longitude, latitude, and altitude from ECEF position vector

  • Inputs:

    • r - position vector in ECEF

  • Outputs:

    • lat - geocentric latitude, (rad)
    • lon - longitude, (rad)
    • alt - altitude (km)

  • Test: Test

---

---

Time Functions - time_functions.cpp

• MJD2GMST(double MJD) -> double

  • Converts from MJD to GMST

  • Test: Test

  • References:

    • Monetenbruck and Gill

---

• date2MJD(int M, int D, int Y, int HH, int MM, double SS) -> double

  • Takes in the current date and time and returns the MJD

  • Inputs:

    • M - month
    • D - day
    • Y - year
    • HH - hour
    • MM - minute
    • SS - second (can be fractional)

  • Outputs:

    • MJD - modified Julian date

  • Test: Test

  • References:

    • Monetenbruck and Gill

---

• valid_date(int M, int D, int Y, int HH, int MM, double SS) -> bool

  • Checks if a date/time is valid or not

  • Inputs:

    • (see date2MJD)

  • Test:

    • Not currently tested due to how asserts work in C++. It just ends the program so a different use will need to be implemented