CubeSat Vehicle Sun Tracking — Robotics/Spacecraft Dynamics

CubeSat Vehicle (Sun Tracking)

Robotics / Spacecraft Dynamics

Flies a CubeSat: its orbit and its attitude together, with the attitude controlled to point a body axis at the Sun, read from JPL's Development Ephemeris DE421: t1 = C_bi·1000·S(JD), S the Sun relative to the Earth in km, at JD = t/86400 + Start Date. Each sample is one RK4 step of the Simulink block's own equations, with the Earth's orientation re-evaluated at every stage:

r′ = v,  v′ = R_woPM^T(ω_E × (ω_E × r_p) + 2ω_E × v_p + R_pm^Ta_ecef)
q′ = ½Ω(w − C_bfω_E) q/|q|,  w′ = M

r and v are the ICRF state, q the ECEF-to-body quaternion and w the inertial body rate; R_woPM is the IAU-2000/2006 ECI-to-ECEF rotation without polar motion and R_pm the polar motion, both from the IERS table at the date; ω_E = [0 0 7.2921158553×10⁻⁵] rad/s. The torque M is a discrete PID (P 10⁻⁵, I 10⁻¹⁰, D 0.01, no derivative filter) on the XYZ Euler angles of the alignment quaternion, which points the body primary alignment a1 along the target t1 and turns the secondary alignment a2 as close as it can to the secondary constraint t2 = C_c·(Secondary Constraint). The body inertia is the identity, as in the Simulink block, so the torque is the angular acceleration.

Ports

• a_ecef – the vehicle's acceleration relative to the rotating Earth, in ECEF axes, [3,1], m/s². Not a force and not gravity: see Notes.
• X_ecef – position in ECEF (ITRF), [3,1], metres.
• V_ecef – velocity relative to ECEF, [3,1], m/s.
• q_eci2b – the ICRF-to-body quaternion, [4,1], scalar first.
• q_ecef2b – the ECEF-to-body quaternion, [4,1], scalar first.

Parameters

• Start Date – the UTC Julian date of t = 0 (default 2458488, which is juliandate(2019, 1, 4, 12, 0, 0)). It dates the Earth's orientation every sample, and the Sun's position.
• Input Method – how the initial orbit is given:
  • Keplerian Orbital Elements – Semi-major Axis (m), Eccentricity, Inclination, RAAN, Argument of Periapsis and True Anomaly (degrees), with μ = 3.986004418×10¹⁴ m³/s². A circular orbit (eccentricity below 1e-12) reads Argument of Latitude, a circular equatorial one True Longitude and an elliptical equatorial one Longitude of Periapsis in place of the angles they replace. Then rotated by the IAU-1976 precession from the ECI Epoch.
  • ECI Position and Velocity – ECI Position (m) and ECI Velocity (m/s) at the ECI Epoch, rotated by the same precession.
  • ECEF Position and Velocity – ECEF Position (m) and ECEF Velocity (m/s) at the Start Date.
  • Geodetic LatLonAlt and Velocity in NED – Geodetic Position (latitude and longitude in degrees, WGS84 altitude in metres) and NED Velocity (m/s) at the Start Date.
• ECI Epoch – the Julian date of the ECI frame the first two methods are given in (default 2451545, J2000).
• Semi-major Axis (m), Eccentricity, Inclination (deg), RAAN (deg), Argument of Periapsis (deg), True Anomaly (deg), True Longitude (deg), Argument of Latitude (deg), Longitude of Periapsis (deg) – the Keplerian elements above, one number each. The orbit must be elliptical (eccentricity at least 0 and below 1).
• ECI Position (m), ECI Velocity (m/s), ECEF Position (m), ECEF Velocity (m/s), Geodetic Position, NED Velocity (m/s) – three numbers each, read by the method that names them.
• Euler Angles (deg) – roll, pitch and yaw of the body relative to the NED frame of the starting point, at t = 0.
• Body Angular Rates (deg/s) – the body rate relative to that NED frame at t = 0, in body axes (default [0 0 -0.05168]).
• Primary Alignment – the body vector a1 the pointing law points at the target (default [0 0 -1]), any nonzero length.
• Secondary Alignment – the body vector a2 that is turned as close as it can come to the secondary constraint (default [0 1 0]).
• Secondary Constraint – the direction t2 for a2, in the constraint frame (default [0 0 1]).
• Constraint Coordinates – the frame the constraints are given in: ECI Axes, ECEF Axes, NED Axes (the local north-east-down frame under the spacecraft) or Body-Fixed Axes.
• Ephemeris Folder – the folder holding JPL's published ASCII DE421 files (header.421 and ascp*.421), as MATLAB's support package holds its own; relative to the project. The Simulink block reads DE421 and nothing else, so neither does this one. Empty refuses the run.
• Sampling Time (s) – zero or less inherits the solver's rate; a positive value runs the block at that period.

Code export

Python, MATLAB, Java, Rust, C and C++, from one program spelled six ways, so every target runs the simulation's arithmetic in the same order: the full 2941-term IAU-2000/2006 series at each RK4 stage, the IERS rows of the days the run covers, from the Start Date to the model's stop time, and the DE421 records of those days. A sample past the stop time reads no IERS row and answers NaN. The initial state is computed at export and written as twenty constants. The other four are not offered and refuse the export by name: VHDL, Verilog, SystemVerilog and PLC Structured Text – a Julian date (2.4 million) and an orbit (7 million metres) do not fit a Q16.16 port, and the series alone is eighteen thousand constants.

Simulink bridge

Import and export, mapped to Aerospace Blockset's aerolibcubesatveh/CubeSat Vehicle (Sun Tracking), one for one: Start Date → sim_t0, Input Method → method, ECI Epoch → epoch, Semi-major Axis (m) → a, Eccentricity → ecc, Inclination (deg) → incl, RAAN (deg) → omega, Argument of Periapsis (deg) → argp, True Anomaly (deg) → nu, True Longitude (deg) → truelon, Argument of Latitude (deg) → arglat, Longitude of Periapsis (deg) → lonper, ECI Position (m) → r_eci, ECI Velocity (m/s) → v_eci, ECEF Position (m) → r_ecef, ECEF Velocity (m/s) → v_ecef, Geodetic Position → lla, NED Velocity (m/s) → v_ned, Euler Angles (deg) → euler, Body Angular Rates (deg/s) → pqr, Primary Alignment → firstAlign, Secondary Alignment → secondAlign, Secondary Constraint → secondRef, Constraint Coordinates → constraintCoord. Input Method and Constraint Coordinates use Simulink's own option spellings, so they cross losslessly. Always written: pointingMode Sun Tracking and firstAlignExt, secondAlignExt, firstRefExt and secondRefExt Dialog, which fix this block's ports. The mission-analysis parameters (missionRT and the rest) do not cross: they only feed the Simulink block's report button, and change nothing it computes (measured, bit for bit). Ephemeris Folder does not cross: Simulink reads its DE421 from the support package on the MATLAB path. The Simulink block has no SampleTime: its PID runs at the model's base rate, which is this block's Sampling Time.

Notes

• ⚠ There is no gravity in it. a_ecef is the acceleration relative to the rotating Earth, and the block adds the Coriolis and centrifugal terms itself, so with a_ecef = 0 it flies a straight line in ECEF. To fly an orbit, feed back a_ecef = g − 2ω_E×V_ecef − ω_E×(ω_E×X_ecef), with g from a gravity block. That is the Simulink block's own contract.
• The outputs are re-derived from the state every sample: X and V through the frames at that instant, and both quaternions through dcm2quat, so either may come out with a negative scalar part.
• The PID's derivative has no filter, so the first sample carries a derivative kick of D·e/Ts – the Simulink block's own, and the reason its first second is lively.
• The ECI-to-ECEF rotation the dynamics use takes ΔAT = 0, while the ECEF and geodetic initialisation takes ΔAT = 37 s; both are the Simulink block's, kept rather than corrected.
• ⚠ The Sun is read at the UTC date as if it were TDB, about 69 s off, exactly as the Simulink block does; it moves the pointing by at most a few microradians.
• Discrete: twenty states (the orbit, the attitude, the PID and a sample count), advanced once per sample. No direct feedthrough: the outputs come from the state.