Precise re-entry and landing of propellantless spacecraft

Spacecraft returning scientiﬁc samples or humans from space must be capable of surviving re-entry and landing in a desired location. Traditionally, this has been accomplished via a propulsive de-orbit burn. However, it is not always possible to mount a propulsion sys-tem on board a satellite or a capsule. In the case of small satellites deployed from the International Space Station, for example, on-board propulsion systems are forbidden for safety reasons. Our work proposes a new technological solution for re-entering and landing a spacecraft in a desired location from a low Earth orbit using exclusively aerodynamic drag and eliminating the need for chemical propulsion. First, an iterative procedure is utilized to compute the desired state at the re-entry interface (100 km) such that a propagation of the vehicle dynamics in the nominal re-entry drag conﬁg-uration from this initial state leads to a landing at a desired latitude and longitude on the surface of the Earth. Next, a re-entry point targeting algorithm is utilized to determine the on-orbit ballistic coeﬃcient proﬁle necessary to target the desired re-entry point. Finally, the ballistic coeﬃcient proﬁle during the ﬁnal hours of the trajectory before the re-entry interface is iteratively modiﬁed to correct any remaining along-track error in the landing location. The proposed solution is applied to a small satellite system that is jettisoned from the ISS and is equipped with a deployable heat shield that also serves as a drag device.


Introduction
Controlled spacecraft re-entry has been essential since the beginning of the space program (AviDyne, 1968) to ensure that crewed space capsules land in a desired location for a safe and timely recovery.All space vehicles that have performed such a controlled re-entry and landing have initiated the re-entry with a de-orbit burn.After the burn, to increase landing precision, a control technique called bank angle modulation has been utilized to rotate the vehicle (using thrusters) around the velocity vector to point the aerodynamic lift vector in a desired location (Morth, 1966;Lu, 2012;Tigges et al., 2011;Putnam et al., 2009;Hankey, 1988;Gallais, 2007;Regan, 1993;Bogner, 1966;Wingrove, 1963;De Zaiacomo et al., 2009;Roenneke and Cornwell, 1993;Schierman et al., 2004;Mooij, 2013;D'Souza and Sarigul-Klijn, 2008).To date, there has not been a spacecraft that has successfully re-entered the atmosphere and achieved a desired landing location without the use of thrusters.
However, some other control mechanisms such as drag modulation (DM) can be used instead of the de-orbit burn technique.In DM the effective drag force on the vehicle is varied to achieve the necessary reduction in orbital energy normally associated with a de-orbit burn during the orbital sector before the re-entry interface.Multiple methods exist for achieving the drag variation including deploying and retracting booms behind the vehicle (Omar and Bevilacqua, 2019) or having a mechanism that can continuously change the aerodynamic shape of the vehicle to generate the required drag force (Fortezza et al., 2013).
In this work, we build upon prior work in drag modulation to develop a set of guidance, navigation, and control algorithms that utilize the modulation of a drag device to lead a low-Earth-orbit (LEO) satellite to a desired landing location using exclusively aerodynamic drag.Prior work on drag modulation explored algorithms that guide a satellite only to a desired re-entry interface point at 100-120 km altitude and then have an uncontrolled re-entry with dispersions on the order of thousands of kilometers on the ground (Omar and Bevilacqua, 2019;Carna ´et al., 2019).This level of accuracy is insufficient for applications such as spacecraft sample return or the recovery of crewed capsules in the event of a propulsion system failure.To facilitate these important potential applications, the original reentry point targeting algorithm developed by the authors is used as one stage of a new ground targeting algorithm in the same way that the re-entry burn is but one component of a more elaborate entry, descent, and landing algorithm on traditional space capsules.This ground targeting algorithm is the first of its kind to the authors' knowledge to provide for the guidance of a satellite from its initial orbit to a ground landing location using exclusively aerodynamic forces with highly accurate guidance trajectories on the order of a few tens of kilometers.
In this paper the methodology developed is applied to a micro satellite (50 kg class) released from an ISS orbit with the intent of a re-entry and landing in Kazakhstan for payload recovery.The micro satellite is equipped with a heat shield that can be repeatedly deployed and folded onorbit (Fig. 1)) to provide drag modulation.The heat shield is based on the design of the Mini Irene Flight demonstrator (Bassano et al., 2011Gardi et al., 2017Vernillo et al., 2017Zamprotta et al., 2019Vernillo et al., 2019).The Mini Irene Flight demonstrator has already been manufactured, tested in a plasma wind tunnel facility, and is now ready to be tested in an upcoming ESA suborbital mission (Fedele and Mungiguerra, 2018Fedele et al., 2019Fedele et al., 2020Fedele, 2020).The same heat shield has also been analyzed for use in a Mars aerocapture mission for a small satellite (Isoletta et al., 2021).
The work here presented is organized as follows: Section 3 explains the guidance trajectory generation algorithm that computes the trajectory and drag profile the satellite must follow to land in a desired location.Section 4 presents the tracking algorithms that enable the satellite to follow the guidance trajectory in the face of uncertainties, particularly in the drag force.Finally, Section 5 shows the results of Monte Carlo simulations and case-studies conducted to verify the performances of the aforementioned algorithms applied to the reference capsule.
Ultimately, the ability to re-enter and recover a propellant-free capsule in a desired location using exclusively aerodynamic drag is a game-changing technology that has significant implications.Services such as ISS sample return could be made possible using small capsules with deployable drag devices stored on the station.Satellites containing scientific payloads could return these payloads to Earth without needing propulsion systems or traditional large heat shields.These techniques could also facilitate the controlled landing of a large number of small, propellantless probes deployed from a mother ship on other planets with atmospheres like Mars or could enable aerocapture for science missions at Uranus and Neptune.Finally, the GNC techniques discussed in this paper could be used as backup methods for the controlled re-entry of crewed capsules in the event of a propulsion system malfunction.

Simulation environment
The high fidelity simulation environment utilized to test the algorithms discussed in this papers consists of two com- ponents.The orbital simulator is utilized to compute the spacecraft's trajectory while in orbit (above 100 km) and contains high fidelity gravitational and density models.The re-entry simulator models the trajectory below 100 km and more precisely characterizes the re-entry aero-thermodynamic environment experienced by the satellite.

Orbital Environment
During the orbital phase of a satellite's trajectory, aerothermodynamic heating is not significant, the dynamic pressure experienced by the satellite is low, and the aerodynamic environment is one of free molecular flow where the mean free path between particles is large and particles do not interact with each other.While the drag coefficient for a given satellite geometry is not expected to change significantly during this regime, large atmospheric density variations can occur due to changes in solar and geomagnetic activity (Vallado and Finkleman, 2014).In addition, because a satellite may spend weeks or months on orbit during a period of interest, the perturbations due to the non-uniform gravitational field of the Earth can have a significant long-term effect and must be taken into account.In the low Earth orbits considered in this work, solar pressure, solar gravity, lunar gravity, relativity, atmospheric winds, tidal effects, and Earth precession and nutation were not found not be significant and were not considered.Only aerodynamic drag and Earth's gravity were precisely modeled in this orbital simulator, but all algorithms in this work were designed to be usable in a higher fidelity simulation environment such as STK if a higher accuracy is desired.

Gravitational Model
The simplest gravitational model assumes that the Earth is a point mass and that gravity always acts toward the center of earth.With this two-body model, the acceleration due to gravity is based on the ECI position of the satellite (r) and is given by Eq. 2.1 (pp.15) in Ref. Montenbruck and Gill, 2005 as where l is Earth's gravitational parameter.A more realistic gravitational model provides a g as the cumulative effect of gravitational perturbations.These perturbations are divided into zonal harmonics which capture variations in Earth's gravity at different latitudes, sectorial harmonics which capture longitude-dependent gravitational effects, and tesseral harmonics which capture gravitational effects that are dependent on both longitude and latitude Madden, 2006.The coefficients from the EGM2008 gravitational model Pavlis et al., 2008 and the technique in  are utilized in this work to accurately compute the acceleration due to the Earth's non-uniform gravity field at each point in time by calculat-ing the gradient of the potential function (Montenbruck and Gill, 2005) 2.1.2.Aerodynamic Model Aerodynamic drag force is discussed in  and is calculated by (Vallado, 2013) where C d is the drag coefficient, A is a reference area, q is the ambient density, and v rel is the velocity of the spacecraft relative to the atmosphere.Aerodynamic lift is not significant during the orbital phase of the trajectory and is not considered in the orbit simulator.Substituting the ballistic coefficient defined as into Eq.( 3), the acceleration due to drag can be written as Drag is by far the most difficult force to accurately predict due to uncertainties in C d and q.Because the atmosphere tends to rotate with Earth due to viscous forces, the velocity vector of the satellite relative to the rotating atmosphere can be approximated as (Vallado, 2013) where r is the spacecraft position vector measured in the ECI frame.
The NRLMSISE-00 atmospheric density model was utilized in this work because it is modern, high performing (Vallado, 2013), and an implementation is readily available in the MATLAB aerospace toolbox.Based on 69,932 density measurements on satellite between 200 and 620 km altitudes, the NRLMSISE-00 model exhibited a mean ratio of measured to actual density of.9949with a standard deviation of.1717(Marcos et al., 2006).In addition to latitude, longitude, altitude, and time, the NRLMSISE-00 model takes as inputs the F10.7 solar indices and the Ap geomagnetic indices.Details about the inputs and implementation of MATLAB's atmosnrlmsise00 function are provided on the MathWorks Website (MathWorks, 2016).45 day forecasts of F10.7 and Ap are available online from the National Oceanic and Atmospheric Administration, 2017.

Numerical Simulation Technique
To simulate the orbit of the spacecraft, the equations of motion are first written in state space form and numerically integrated.The spacecraft state vector x during the orbital phase consists of the ECI (Earth-Centered Inertial) position and velocity. where The ECI frame is defined as aligned with the Earth-Centered-Earth-Fixed (ECEF) frame at the simulation epoch, and the ECEF frame is assumed to rotate about the ECI zaxis (through the North pole) at a constant rate of x È ¼ 7:292 Â 10 À5 rad/s.The derivative of the state vector in the ECI frame can be written as where a is the summation of the accelerations induced by all forces acting on the spacecraft.a can be computed by Newton's second law (assuming spacecraft mass is not changing) as Eq. ( 8) along with an initial state value are numerically integrated using MATLAB's ode113 function (MathWorks, 2017).The numerical integration process provides the evolution of the state vector over time in the ECI frame.

Re-Entry Environment
A three degree of freedom simulator was developed in SIMULINK to simulate the re-entry environment.The NRMLSISE-00 atmosphere model has been used to evaluate air density q and sound speed a values depending on the current altitude, position, and epoch.An aerodynamic database based on CIRA heritage projects is also used to help evaluate aerodynamic forces.The core of the reentry simulator is composed of an "Equations of Motion" block which takes aerodynamic forces as an input and integrates the differential equations of motion given by Eq. 14a and Eq.15a in order to compute the spacecraft's state over time.

Reference Frames
Since the objective of the simulation is to reproduce the motion of the re-entry capsule, the equations of motion must be integrated with respected to an inertial (nonmoving) reference frame.The geocentric-equatorial Earth-centered inertial (ECI) reference frame ðO; X ; Y ; ZÞ with its origin at the center of the planet, the x-axis pointing in the direction of the vernal equinox, and the z-axis pointing through the North Pole is used for this.The yaxis completes the right-handed coordinate system and lies in the planet's equatorial plane.
Assuming the planet rotates with a constant velocity x È around the z-axis, X ¼ x È Dt represents the rotation rate between the inertial and the Earth-centered-Earth-fixed (ECEF) ðX R ; Y R ; Z R Þ reference frames.The two frames are aligned when Dt ¼ 0.
In the ECEF reference frame, h and / represent the longitude and the geocentric latitude of the vehicle's position, respectively, while r represents the vehicle's geocentric altitude (see Fig. 2).
Note that the geocentric latitude / is different from the geodetic latitude / G . while the first is defined as the angle between the radius r and the equatorial plane, the latter is defined as the angle between the local surface normal vector and the equatorial plane.This definition is related to distinction between geocentric and geodetic altitude (r and h respectively) as illustrated in Fig. 3.
If a local horizontal plane is introduced as the plane perpendicular to the vector from the surface of the Earth to the satellite at any given instant, the flight-path angle c can be defined as the angle between the local horizontal plane and the capsule velocity vector Ṽ .Similarly, the heading w is the angle between the local parallel of latitude and the projection of Ṽ on the horizontal plane.
By convention, c is positive when Ṽ is above the local horizontal plane while w is increased when turns are made toward the left and is zero when facing east.ðCM; X B ; Y B ; Z B Þ represents the body reference frame for the considered spacecraft (see Fig. 4).

Atmosphere and Gravitational Model
Different analytical and empirical atmospheric models exist for the computation of capsule re-entry trajectories.Their major objective is to provide air density q and temperature T values since these are essential for the computation of aerodynamic forces, Mach number, and other coefficients.
As in the orbital simulator, an empirical global reference atmospheric model called NRMLSISE-00 has been implemented in the 3DOF simulator in order to model density and airspeed values (see Fig. 5).
As discussed in Section 2.1.1,the gravitational acceleration vector is given by the gradient of Eq. 2. Since the J 2 term is significantly larger than the higher order coefficients, a first or second order approximation is sufficient for re-entry modeling.When the J 2 effect is considered, the gravity vector remains in the plane containing the satellite position vector and the north pole but is given by two components.The first component is directed opposite the position vector ẽr , while the second is perpendicular to the position vector along the Àẽ / direction: The two components are given by: Comparing these two components, g r (or z g considering the vehicle-pointing frame) is considerably larger than g / (or x g ).The Earth-related constants used inside the simulation are listed in Table 1.

Aerodynamic Forces
Aerodynamic forces can be decomposed into three main components: drag (D), lift (L) and side force (Y).Drag acts opposite the direction of the velocity vector.Lift is perpendicular to the velocity vector, in the symmetry plane of the vehicle.Side force acts in the horizontal plane.While only drag is considered during the orbital simulation, all three aerodynamic force components are significant during the re-entry and are therefore considered.
The aerodynamic forces can be defined as:    of attack and sideslip angles respectively.S represents the spacecraft reference surface area.

Equations of Motion
A set of first order nonlinear ordinary differential equations is needed to simulate vehicle's 3DOF dynamics Hankey, 1988.
A set of 3 dynamic equations is needed to take into account forces acting on the spacecraft.Velocity, flight path angle, and heading angle are the variables which describe vehicle's dynamic as: where m is the spacecraft's mass.
A set of 3 kinematic equations is needed to compute vehicle's position.Geocentric altitude, latitude, and longitude describe the spacecraft's position over time as: This system of equations is numerically integrated over time using a 4th order Runge-Kutta solver.

Algorithm Overview
The goal of this algorithm is to compute a trajectory and corresponding ballistic coefficient profile, subject to constraints, that if followed leads a spacecraft to impact the ground in a desired location.For any given spacecraft geometry and drag modulation capability, the minimum and maximum ballistic coefficient can be calculated a priori based on drag coefficient, mass, and cross-sectional area and provided as inputs to the targeting algorithm.The drag coefficient is the most challenging to compute and uncertain of the parameters in the ballistic coefficient formula, and a detailed analysis of drag coefficient computation for each regime (free molecular, transitional, continuum) used in this work can be found in Fedele and Mungiguerra, 2018.The algorithm will compute a guidance trajectory with ballistic coefficients within the specified range.This makes the algorithm agnostic to the specific spacecraft design and facilitates the propellantfree landing and recovery of an arbitrary spacecraft in a precise location provided that there is some known ballistic coefficient modulation capability.The methods discussed in this section build on prior work by the authors on space-   J 2 1:082629 Á 10 À3 craft de-orbit point targeting (Omar and Bevilacqua, 2019;Omar et al., 2017;Omar and Bevilacqua, 2019).The previous work only provided a guidance trajectory to a specified de-orbit latitude and longitude at 100-120 km altitude but did not provide means of generating a trajectory from the initial conditions to a specified ground impact point.It also did not consider requirements on the velocity and flight path angle at the re-entry interface.The key contribution of this work is the ability to generate and track a trajectory all the way to a ground impact point.
The first novel component of the ground targeting algorithm involves the selection of the de-orbit point.The desired de-orbit state is defined such that a numerical orbit propagation of the spacecraft from this state results in a ground impact in a desired location.The de-orbit point is the latitude, longitude, and altitude at the deorbit state which in this work is at 100 km geocentric altitude.The next phase of the algorithm involves computing a trajectory that guides the satellite to the de-orbit point from the initial conditions utilizing the range of ballistic coefficient modulation achievable by the satellite.Finally, while a guidance trajectory to the de-orbit point may be achieved, the final velocity, flight path angle, and atmospheric conditions may vary slightly from those in the originally selected de-orbit state.This may result in some along-track error in the landing location.For this reason, the last novel component of the ground targeting algorithm performs a correction to the portion of the generated guidance trajectory between 140 km and the deorbit point at 100 km to correct for this error and ensure that the desired landing point is achieved.After this correction, the complete guidance trajectory consists of a ballistic coefficient, position, and velocity time-series from the spacecraft initial conditions to the ground impact point, though the tracking of this trajectory can be terminated at any point prior to the ground impact point if desired.

De-Orbit Point Selection
The first novel component of the ground targeting algorithm involves computing the de-orbit location at 100 km altitude that the satellite must pass through to achieve a desired ground impact location.This is accomplished by initializing a simulation in a 140 km circular orbit defined by initial mean orbital elements x 140 = ða; e; i; X; x; hÞ = (6518 km, 0, i init , 0, 0) where a is semi-major axis, e is eccentricity, i is inclination, X is right ascension of the ascending node, x is argument of the perigee, and h is true anomaly.The time is set equal to the approximate expected time of de-orbit.A satellite trajectory is simulated from this initial state using an orbit propagator until a geocentric altitude of 100 km and a re-entry simulator from 100 km until ground impact.Initiating the trajectory at an altitude of 140 km instead of at 100 km directly helps ensure that the velocity and flight path angle as the satellite passes through 100 km accurately reflect the appropriate values for a re-entry trajectory.This is important as velocity and flight path angle at the re-entry interface are the primary parameters that affect landing location.After trajectory simulation, the increases in X and h at x 140 necessary to achieve a desired impact location are computed.This is done as follows.

Ground Impact Latitude Targeting
First, the argument of latitude (sum of argument of perigee and true anomaly) at both the desired and actual (from numerical simulation) ground impact locations can be computed based on the orbital inclination and the latitude of each location.During the ascending (northbound) phase of the orbit, the argument of latitude u at a given latitude l is computed by The argument of latitude during the descending portion of the orbit is given by To compute the desired argument of latitude u des , the user must specify whether they want the re-entry trajectory to occur during the ascending or descending portion of the orbit.To determine the actual argument of latitude u act after a re-entry trajectory propagation, it is necessary to utilize the impact latitude in Eq. ( 16) and determine whether the trajectory is ascending or descending by examining trajectory points before the impact point to determine whether the latitude of each point is increasing or decreasing.The true anomaly increase necessary in x 140 is to arrive at an impact point of the correct latitude is

Ground Impact Longitude Targeting
The change in right ascension (DX) necessary to correct for the longitude error can be computed by first calculating the longitude increase Dk u that occurs when a satellite travels from an argument of latitude u act to u des .This is most easily done by setting h ¼ u des in x 140 and computing the longitude associated with the orbital element set.h ¼ u act can then be set in x 140 and the longitude can be once again computed.The longitude at u act can be subtracted from the longitude at u des to get the longitude increase that occurs as the satellite traverses a true anomaly increase of Du.The right ascension change in x 140 necessary to correct the longitude error is then where k des is the longitude of the desired ground impact location and k act is the impact-point longitude after the trajectory simulation.

Initial Conditions Update
The desired orbital element increases DX and Dh ¼ Du can be applied to the initial state at x 140 and the trajectory re-propagated.The process of computing the change in x 140 necessary for a desired landing location and re-propagating the trajectory can be continued until a x 140 is found that results in a trajectory that impacts that ground in the desired location.When the proper x 140 is found, the latitude, longitude, and altitude where the trajectory crosses 100 km geocentric altitude is returned as the target deorbit point.The negative flight path angle that is seen for a passive de-orbit trajectory at 100 km is caused by aerodynamic drag.If trajectory propagation is initiated from a high enough altitude so that aerodynamic forces have sufficient time to shape the trajectory, the drag-induced flight path angle at 100 km will be nearly the same regardless of the precise initial conditions.For this reason, trajectory propagation is initiated at 140 km geocentric altitude with zero flight path angle to ensure that the trajectory has approximately the correct flight path angle at 100 km.
Note that in some cases, directly applying the desired Du to x 140 and propagating the trajectory results in a convergence failure.On one trajectory simulation, an overshoot of the desired landing location may occur, but on the next an undershoot occurs.This overshoot-undershoot cycle continues indefinitely in some cases.To remedy this, a change in argument of latitude of kDu instead of Du can applied to x 140 where k is a scaling factor initially set to 1.Each time the sign of the desired Du changes between iterations of trajectory propagation, k is reduced by a factor of two.This reduction in k breaks the undershootovershoot cycle and helps ensure convergence.

De-Orbit Point Targeting
Once the desired de-orbit point at 100 km geocentric altitude is selected, a trajectory and corresponding drag profile can be computed to guide the satellite to this deorbit point using the techniques detailed in (Omar and Bevilacqua, 2019).A flow chart of this de-orbit point targeting algorithm is shown in Fig. 6.Note that a geocentric instead of geodetic altitude is used to define the de-orbit point since the targeting algorithm is designed to work with geocentric altitude.The details of the de-orbit point targeting algorithm are summarized for the reader as follows.
Define the spacecraft ballistic coefficient as in Eq. 5.In this work, the reference area is assumed to be the projected area perpendicular to velocity.In the case of calculation with the heat shield deployed, this area is equal to the diameter of the open heat shield.In the case of calculation with the heat shield in a folded configuration, this area is equal to the frontal area of the satellite (see Fig. 1).It can be shown that given the ability to vary an initial C b1 , a second ballistic coefficient C b2 , and the time t swap at which the ballistic coefficient is changed from C b1 to C b2 , it is possible to target any point on the Earth with latitude below the orbit inclination if maneuvering is initiated early enough.We will define (C b1 ; C b2 ; t swap ) as the control parameters.Note that C b2 is maintained until some terminal semi-major axis at which point some pre-set ballistic coefficient ðC b term Þ is maintained until de-orbit.
To determine the control parameters needed to target a desired de-orbit location, a trajectory is first numerically propagated in a high fidelity simulation environment using initial guesses for the control parameters denoted by C b10 ; C b20 , and t s0 .An effective initial guess is to set C b10 and C b20 equal to the average achievable spacecraft ballistic coefficient and to set t s0 equal to half the expected orbit lifetime remaining with this drag configuration.The resulting trajectory can then be characterized by two sections.Section one is from the initial conditions to the C b swap point and section two is from the swap point to the terminal point.This analytical solution can be summarized as follows.Assume a satellite with ballistic coefficient C b1 takes time Dt 1 to achieve some drag-induced change in semi major axis Da and undergoes argument of latitude change Du 1 while achieving this Da.The time and argument of latitude change a satellite with the same initial conditions and some different C b2 will undergo to achieve the same Da (same orbital decay) are given by Since the average rate of change of right ascension ( _ X avg ) is independent of C b , the change in X experienced during the orbital decay can be calculated by As shown in Fig. 7, if the trajectory of a satellite with some initial set of control parameters has been numerically propagated (initial trajectory), the de-orbit location of a new trajectory corresponding to the same initial conditions but a different set of control parameters can be analytically estimated by dividing the trajectories into phases where the C b is not changing in either trajectory.In each trajectory, the phases are demarcated by the ballistic coefficient swap point (t swap ), the point at which the semi major axis (orbit energy level) is the same as at the swap point of the other trajectory (t eq ), and the terminal point (t term ).For the three phases before the terminal point, Eqs. ( 20)-( 22) can be utilized to calculate the changes in time and orbital elements experienced in each phase of the new trajectory.All changes in time and orbital elements can be added to calculate the final time and orbital elements, and hence the latitude and longitude, at the de-orbit point.
Through an algebraic manipulation of these same equations, the control parameters (C b1 ; C b2 , and t swap ) needed to achieve a desired total time (Dt t ) and total change in argument of latitude (Du t ) to the terminal point can be computed based on C b10 ; C b20 ; Dt 10 ; Dt 20 ; Du 10 , and u 20 from the original numerically propagated trajectory via the following equations The total argument of latitude u des required to de-orbit at the desired latitude can be calculated by Eq. ( 16) if the orbit is ascending and Eq. ( 17) if the orbit is descending.The total argument of latitude increase for the trajectory to land at the correct latitude is then where k is an integer.The total time to the terminal point required to achieve the correct longitude can be computed by Where Dt t0 is the time to the terminal point in the original numerically propagated trajectory, k deo is longitude of the de-orbit trajectory at the appropriate latitude, k des is the desired de-orbit point longitude, and X e is the rotation rate of the Earth.If the optimal Dt t cannot be achieved within the drag control capability of the satellite, Dt t can be adjusted to ensure that the minimum targeting error is achieved within the range of feasible satellite ballistic coefficients.
The analytically computed ballistic coefficient profile is then refined via a shrinking horizon and drag-workenforcement method until a C b profile is found whose corresponding trajectory leads to the desired de-orbit location.In the shrinking horizon approach, when the trajectory is propagated to the de-orbit point with the ballistic coefficient profile dictated by the analytical solution, the first t g seconds of this trajectory are stored as a part of the guidance solution and the trajectory after t g is utilized to compute another analytical ballistic coefficient profile that will be numerically propagated and will be approximately t g seconds shorter than the previously propagated trajectory.This process continues until a certain error threshold is reached or a trajectory is propagated that has less than a certain amount of orbit lifetime remaining.With the drag-work-enforcement enhancement, the work done by aerodynamic drag is recorded during the trajectory propagation, and the ballistic coefficient of the satellite during the first t g seconds of propagation is varied so that the total work done by drag at t g is equal to the work that should have been done by this time according to the analytical solution.

Terminal Orbit Adjustment
Once a guidance trajectory that leads from the spacecraft initial conditions to the desired de-orbit point is computed, propagation of this trajectory from the de-orbit point to the ground impact point may not necessarily result in a ground impact at the desired location.This is because the final velocity and flight path angle at the de-orbit point as well as the density profile may not be identical to those at the initially calculated de-orbit state, even if the latitude, longitude, and altitude are the same.Fortunately, the impact-point errors that result from this are almost entirely in the along track direction and can be corrected by modifying the ballistic coefficient during the final portion of the guidance trajectory before the de-orbit point.This ballistic coefficient adjustment allows for control of the along-track impact location without the need to control re-entry velocity or flight path angle directly.
To determine the necessary C b adjustment, the argument of latitude increase Du des needed to achieve the desired ground impact location must first be computed using the same novel procedure introduced in Section 3.2.1.Consider for example the desire to de-orbit at latitude l d in the ascending portion of an orbit with inclination i.Given an actual ground impact latitude after numerical simulation, the effective argument of latitude u imp at ground impact can be computed using Eq.16 or Eq. 17 depending whether the orbit was ascending or descending at the time of impact.Du des can then be computed as To achieve this desired Du des , all ballistic coefficient values between time t adj and the de-orbit point (t deo ) can be modified according to where C b0 is the original ballistic coefficient and Du 0 is the argument of latitude change between t adj and the de-orbit point in the initially propagated trajectory.To compute the t adj that should be used to ensure that the ballistic coefficient change is not excessive, the user can first define the maximum allowable percentage of ballistic coefficient change q.In this work q ¼ :04 was used.If T is the approximate orbital period at the semi-major axis of the de-orbit point, t adj can be computed by where the orbital period T is given by (Curtis, 2009) This ensures that the percentage change in ballistic coefficient associated with Eq. ( 29) will be less than q.A new trajectory from t adj to the ground is simulated with the modified C b profile prescribed by Eq. ( 29).The process of updating the C b profile and numerically propagating to the ground continues until a trajectory is found with an impact latitude error below a specified tolerance.Note that the novel process explained in the second paragraph of Section 3.2.3 can be utilized to reduce Du des by a shrinking scaling factor k to aid in the convergence to a trajectory with the desired impact latitude.If a certain number of terminal C b computation iterations take place without convergence to the desired impact latitude, the ballistic coefficient in the next iteration can be computed by where k < 1.The trajectory and C b profile after t adj are appended to the trajectory and C b profile before t adj to obtain a complete guidance trajectory and C b profile from the spacecraft initial conditions to the desired ground impact point.Fig. 8 contains a flowchart of the terminal orbit adjustment algorithm.

LQR Feedback Control
Due to uncertainties in the drag force, the spacecraft will eventually drift from the guidance trajectory if the desired ballistic coefficient profile is applied open loop.While it is possible to re-generate the guidance trajectory once the drift exceeds a given threshold, guidance trajectory generation is computationally expensive and there is no guarantee that a new guidance trajectory with equally low error will be found from the new spacecraft initial conditions.For this reason, feedback control techniques must be utilized to vary the commanded spacecraft C b based on the difference between the actual and desired state to ensure that the computed guidance trajectory is followed.While the spacecraft is tracking an initial guidance trajectory, new guidance trajectories can be periodically generated and tracked to take into account updated density forecasts.
An LQR-based guidance tracking algorithm can be utilized to ensure the spacecraft remains on the guidance trajectory despite uncertainties.In summary, this tracking algorithm revolves around the use of the Schweighart Sedwick (SS) equations of relative motion (Schweighart and Sedwick, 2002) to specify the evolution of the in-plane position and velocity of the spacecraft relative to the guidance trajectory in the classic state-space form The use of relative dynamics is advantageous because it provides an accurate linear approximation of otherwise nonlinear absolute dynamics which facilitates the implementation of feedback control.Considering only the in-plane relative position and velocity (dx; dy; d_ x; d _ y) in the LVLH frame (Fig. 9) because aerodynamic drag cannot be used for out-of-plane control, and considering a relative d€ y due to a difference in the C b between the spacecraft and the guidance trajectory, the SS linearization can be written as (Schweighart and Sedwick, 2002) where With the equations of relative motion in state-space form, it is possible to use the linear quadratic regulator (LQR) (Franklin et al., 2002) feedback control approach  (Curtis, 2009).
to drive the relative position and velocity to zero with K selected to minimize the cost functional where Q and R are square weighting matrices of appropriate dimension.Because the state is four-dimensional and the control is one-dimensional, Q and R will be 4 by 4 and 1 by 1 matrices respectively and K will be a 1 by 4 gain matrix with the control given by Recall that dx and dy are the radial and along track position and velocity of the satellite with respect to the guidance trajectory in the LVLH frame shown in Fig. 9.Because along-track error is far greater than radial error in general, it makes sense to consider radial error only in terms of its contribution to along-track error.For this reason, setting yields superior tracking performance for a given value of R than any other value of Q with a comparable matrix 2norm.With Q fixed, R must be set based on the desired magnitude of the spacecraft response to deviations from the guidance trajectory.This can be done by first defining Dr sat as the desired dy at which the magnitude of the commanded change in ballistic coefficient will be equal to C bmax À C b min .That is, the controller will be guaranteed to saturate at dy ¼ Dr sat .The LQR gain K can first be computed using an arbitrary initial R (R ¼ 10; 000 used in this work) and then recomputed after updating R based on the initially obtained K to enforce controller saturation at dy ¼ Dr sat .The equation to update R is

Kalman Filter for State Estimation
Because the LQR controller needs an estimate of the relative, not absolute, position and velocity, the LVLH position and velocity relative to the guidance can be computed for each GPS measurement and treated as the "observation" for estimation purposes.A standard Extended Kalman Filter (EKF) formulation is used to remove noise from this observation (Ljung, 1979).First, an initial rela-tive state estimate ðx þ iÀ1 Þ is converted to the ECI frame, propagated to the time of the next available measurement and converted back to a relative position and velocity ðx À i Þ.Let this numerical propagation process be denoted by The estimation error covariance (P) is then updated using the state transition matrix (U) derived from the SS dynamics as Such that This is the predict phase and can be described mathematically as where Q is a user-defined process noise covariance matrix.
The actual (noisy) GPS measurement converted to a relative in-plane position and velocity (z) is finally used to update the state and error covariance matrices as shown in Eq. 46.
In Eq. 46, G ¼ I 4x4 is the mapping from the state that is being estimated to the measurement (z i ¼ Gx i ), W is a user defined measurement noise covariance matrix, K i is the Kalman gain, and S is an intermediate term used in the calculations.K is an anti-smugness term set slightly greater than 1 to ensure that the process noise covariance (P) does not become too small and lead to a lack of responsiveness of the filter to new inputs.

Results
In this section, the results of Monte Carlo simulations and case-studies conducted to verify the performances of the algorithms are presented.The reference mission scenario is a micro satellite (50 kg class) released from an ISS orbit with the intent of a re-entry and landing in Kazakhstan for payload recovery.The micro satellite is equipped with a heat shield that can be repeatedly deployed and folded on-orbit to provide drag modulation and orbital maneuvering capabilities.The simulations discussed utilized high fidelity orbit and re-entry propagators with gravity calculated using the EGM2008 model (Montenbruck and Gill, 2005) and atmospheric density calculated using the NRLMSISE-00 model (Vallado and Finkleman, 2014).Gravitational spherical harmonics through degree and order 4 were used in the EGM2008 model, and historic F10.7 and Ap index values were used for the NRLMSISE-00 density model.

Case Specific Simulation Results
Once a desired de-orbit point at 100 km altitude is selected, the guidance trajectory generation algorithm is utilized to generate a trajectory from the satellite initial conditions to the de-orbit point.The ballistic coefficient during the last part of this trajectory is then adjusted such that the satellite lands as close as possible to a desired location after propagation to the ground.The complete guidance trajectory is then a position, velocity, and ballistic coefficient profile all the way from the spacecraft initial conditions to the desired ground impact location.In this study, the desired capsule landing location was a geodetic latitude and longitude of 47.06 degrees and 58.3 degrees respectively.In one re-entry scenario, the initial mean orbital elements were ða; e; X; x; h; iÞ = (6795.6km,.0007 0, 290.92 deg, 130.06 deg, 80.81 deg, 51.64 deg).The minimum and maximum ballistic coefficients were 0.008715 m 2 /kg and 0.03064 m 2 /kg respectively.The deorbit point selected by the algorithm described in Section 3.2 was at a latitude of 30.54 degrees and a longitude of 29.87 degrees at a 100 km geocentric altitude.This led to a guidance trajectory with an error of under 4 km at the ground impact location.Three sinusoidally varying error factors were included in the simulated tracking of this trajectory and were multiplied by the density to model the effects of uncertainties in the aerodynamic drag force.These error factors had periods of 26 days (solar rotation period), 1 day (Earth rotation period), and 5400 s (approximate orbital period) and amplitudes of 1.25, 1.1, and 1.1 respectively.These density error profiles were selected based on prior literature showing that NRLMSISE-00 density estimation errors have periodicity strongly influenced by the orbital period, Earth period, and Solar period and are on bounded by the aforementioned magnitudes (Lean et al., 2006).The tracking was also simulated with an actuator deadband of 5% and the assumption that the tracker could command a maximum ballistic coefficient equal to 1.5 times the maximum guidance ballistic coefficient and a minimum ballistic coefficient equal to the minimum guidance ballistic coefficient divided by 1.5.This ensures that there is always a sufficient C b margin to correct for the simulated drag uncertainties and any tracking errors are a result of suboptimal controller performance rather than a complete saturation of the actuator.Ballistic coefficient was maintained as a state variable during the numerical integration process with the assumption that 240 s were required for the drag device to move from the minimum to the maximum drag configuration and vice versa _ C b ¼ C bmax ÀC bmin 240 À Á .The deadband prevents the commanded ballistic coefficient from changing unless the desired ballistic coefficient from the feedback controller differs from the latest commanded ballistic coefficient by more than 5%.This minimizes actuator chattering in the face of sensor noise and greatly reduces the total actuator run-time necessary for tracking.Fig. 10 shows the ballistic coefficient and tracking error over time to the de-orbit point during the guidance trajectory tracking simulation.This highlight's the tracker's ability to maintain the satellite on the guidance trajectory despite drag force uncertainties.Beyond the de-orbit point, the satellite can be allowed to re-enter without further control or an alternative set of tracking algorithms can be utilized to follow the guidance trajectory all the way to the desired ground impact point.

Guidance Algorithm Performance
Five hundred Monte Carlo simulations of the guidance trajectory generation algorithm were conducted to verify the ability to calculate an achievable drag profile and corresponding trajectory that if followed, will allow the spacecraft to re-enter and land in a desired location.The initial mean orbital elements were randomly selected to represent potential ISS initial conditions and the epoch for each run was randomly selected from within an eleven year period corresponding to a complete solar activity cycle.The distributions from which the simulation initial conditions and scenario parameters were selected are shown in Table 2.
In Fig. 11, a map with the landing locations of the simulated trajectories is shown.As better shown in Fig. 12, in 100% of the cases the guidance trajectory landing error is below 20 km and in 73.2% of the cases the error is below 10 km.The average error is 7.146 km and the standard deviation is 4.05 km, showcasing the innovative algorithm's ability to reliably generate highly accurate guidance trajectories down to the ground.Previous iterations of the guidance trajectory generation algorithm were only capable of generating trajectories to a de-orbit point latitude and longitude without consideration of the landing location, leading to ground impact errors on the order of hundreds of kilometers (Carna ´et al., 2019).
The overall mission duration varied between fifty and three hundred days (see Fig. 13) depending on the simulation epoch.This is because density can vary by up to two orders of magnitude at a given location within the eleven year solar activity cycle.These time scales are reasonable for many potential Space Station sample return missions and could be further reduced to less then a month using a larger drag device or a bigger initial delta-v at the deployment from the ISS (using for example the tethered technique explained in Brunello et al., 2021).The overall mission duration reduction is achievable because for all the guidance trajectories generated, the drag device is generally only actuated during the last ten to twenty days of the mission (see Fig. 14).Prior to this, the satellite is in its maximum drag configuration (heat shield fully deployed) to decay as fast as possible and no significant operator inputs are required.

Tracker Algorithm Performance
In this section, the ability of the system to track a reference trajectory is analyzed.The tracking of each guidance trajectory presented in Section 5.2 is simulated using sinusoidally varying drag errors during the orbital phase and a fixed drag error during the re-entry phase.As in the example case in Section 5.1, the sinusoidally varying error factors have periods of 26 days, 1 day, and 5400 s and amplitudes of 1.25, 1.1, and 1.1 respectively.Uncertainties of 10% in the density, vehicle mass, and aerodynamic coefficient during the re-entry phase are modeled by a single random error term between [.9 and 1.1] that remains constant during each re-entry simulation and is multiplied by the total aerodynamic drag force.This uncertainty parameter has been selected based on a review of the literature on re-entry modeling Desai et al., 1997;Desai and Cheatwood, 2001;Desai et al., 2008.The 10% uncertainty on the drag force during the re-entry phase is conservative as ambient density fluctuates significantly less and is better known at lower altitudes than at higher altitudes.The re-entry phase is also short enough that a single time-invariant uncertainty term can be used.It is worth underlining that many other factors could be taken into account for a more detailed analysis.A sensitivity analysis could be used to identify the dominants parameters, but the objective of the Monte Carlo analysis reported in this work is to provide an order of magnitude for landing dispersions.All tracking simulations are run assuming that the maximum C b achievable by the spacecraft is a factor of 1.5 greater than the maximum allowable guidance C b and the minimum achievable C b is a factor of 1.5 less than the smallest allowable guidance C b .Active feedback control is terminated at the de-orbit point and the drag device is assumed "locked" to a predetermined C b at this point.However, feedback control beyond the de-orbit point (Fedele et al., 2021) could be employed for actuators that support it as the guidance trajectory is valid all the way to the ground impact point.

Tracker Performance
Five hundred Monte Carlo simulations were conducted to verify the ability to track a drag profile and correspond-ing trajectory that if followed, lead the satellite to a desired de-orbit point.The position error at the de-orbit point is under 1 km in most of the cases (see Fig. 15), highlighting the robustness and performance of the tracking algorithm.
Assuming that 240 s are required to fully deploy or retract the drag device, the actuator run time was below 0.3% of the total mission time in each simulated scenario (see Fig. 16), indicating that power consumption is unlikely to be a concern with the proposed aerodynamically controlled re-entry approach.To give an example, the Drag De-Orbit Device (D3) is a drag actuator for CubeSats that uses 16.4 Watts peak power while deploying or retracting (Omar and Bevilacqua, 2019).At a duty cycle of .3%, this comes to 50 milliwatts of orbit-averaged power draw which is generally well below the power generation capability of even the smallest CubeSats.
In Fig. 17, a map with the landing location dispersion is shown.As better shown in Figs.18 and 19, in 99.6% of the cases the error is below 150 km and in 66.6% of the cases the error is below 50 km.The average error is 42.148 km and the standard deviation is 31.314km.Note that no control is performed beyond the 100 km altitude de-orbit point in these cases.While such performance is sufficient for many missions, greater landing point accuracy may be necessary for the recovery of sensitive payloads.To further increase accuracy at landing, a control system based on aerodynamic flaps can be used beyond the 100 km altitude.This propellantless method involves using the flaps to control the vehicle's angle-of-attack to generate lift and side forces (Fedele et al., 2021).

Tracker Performance with Sensor Noise, Actuator Delays, and Kalman filter
In this section, the ability of the system to track the reference trajectory simulating both sensor noise and actuator delays is analyzed.A Kalman filter as discussed in Sec-tion 4.2 is used for this analysis.A GPS position and velocity measurement is assumed to be available every second and is generated by adding noise terms to the true ECI position and velocity.The noise terms consist of a Gaussian component with standard deviation of 5 m for position where t is the time since simulation epoch.Fig. 20 shows the tracking error at 100 km geocentric altitude (the re-entry interface).In Fig. 21 a map with the relative dispersion on the ground is shown assuming tracking with a Kalman filter to 100 km altitude and an uncontrolled trajectory thereafter.Fig. 22 shows the latitudinal and longitudinal components of the final position error on the ground.As shown in Fig. 23, in 81.8% of the cases the error is below 50 km.The average error is 30.5613km and the standard deviation is 23.6871.In simulations including noise and Kalman filtering, the error at the deorbit point is slightly larger and the error on the ground is slightly lower than for the case with no noise and no filtering.This occurs because the Kalman filter acts as a sort of low-pass filter on the error signal and prevents the feedback controller from "chasing" small perceived errors caused by either sensor noise or minor deviations from the guidance trajectory.Near the de-orbit point where dynamics are changing rapidly, error is slightly larger when the filter is included because the feedback controller does not act as aggressively to correct for observed alongtrack errors.However, this lack of rapid maneuvering minimizes the magnitude of the drag changes induced by the feedback controller.This means that the commanded ballistic coefficient is closer to the guidance value and thus the final flight path angle is closer to what is expected in the guidance trajectory.Because of the sensitivity of the re-entry trajectory to the flight path angle at the de-orbit point, this in turn results in a smaller ground impact error than in the scenario without the Kalman filter.
Fig. 24 shows the drag device actuation motor run times as a percentage of total orbit lifetime during the tracking phase with sensor noise, Kalman filtering, and the assumption that 240 s are required for the drag device to fully deploy or retract.Due to erroneous tracking error signals caused by sensor noise that persists despite the Kalman filter, actuator run time is about 3 times greater than in the noise free case.However, average actuator run time is still less than 1% of the total orbit lifetime and is substantially less than for cases with sensor noise and no Kalman filter.

Conclusions
This paper presents a means of targeting a precise landing location on Earth for a spacecraft by modulating the aerodynamic forces the spacecraft experiences.First, a de-orbit point at 100 km altitude is selected such that a re-entry trajectory propagated from this point leads to a desired landing location.Next, a trajectory and corresponding drag profile are computed that guide the satellite from the initial conditions to this de-orbit point.Finally, the ballistic coefficient during the last part of this trajectory is adjusted to correct for any remaining along-track error in the landing location.During the orbital phase, the total aerodynamic drag is modulated based on a feedback control scheme to ensure the spacecraft remains on the desired guidance trajectory.Monte Carlo simulations were conducted to verify the performance and robustness of the ground targeting algorithms.In each of the 500 tested cases, a guidance trajectory and corresponding drag profile were generated from the initial conditions to the desired landing location.In all cases, the guidance trajectory landing error was less than 20 km.This represents an order of magnitude improvement over prior algorithms that can not generate trajectories with landing errors this low.The tracking of each guidance trajectory was then simulated using realistic sinusoidally varying drag errors during the orbital phase and fixed aerodynamic force errors during the re-entry phase.In all cases, the position error at the de-orbit point was less than 2 km.The average position error on the ground was 42.2 km with 99.6% of errors less than 150 km with no control beyond the de-orbit point (100 km altitude).The average error was 30.56 km with 81.8% of errors below 50 km when trajectory tracking was simulated with sensor noise and a Kalman filter in the loop.Since the guidance trajectory is valid all the way to the ground, landing accuracy can be further improved by performing trajectory tracking beyond the de-orbit point.
Ultimately, the work in this paper demonstrates that the landing of a spacecraft in a desired recovery location is indeed feasible using exclusively aerodynamic forces.The algorithms are sufficiently robust for a flight mission and are capable of maintaining the nominal satellite trajectory in the face of uncertainties in the aerodynamic model.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Example of a Re-entry Satellite in Folded and Deployed Heat Shield Configurations.
13cÞ where C D ; C L and C Y are the non-dimensional drag, lift and side force coefficients while a and b are the spacecraft angle

Fig. 7 .
Fig. 7. Characterizing Behavior of New Trajectory Based on Old Trajectory.

Fig. 20 .
Fig. 20.Tracking Algorithm Performance with Sensor Noise, Actuator Delay and Kalman Filter: Tracking Error at 100 km De-Orbit Point.

Fig. 21 .
Fig. 21.Tracking Algorithm Performance with sensor noise, actuator delay and Kalman filter: Dispersion on Ground.

Fig. 22 .
Fig. 22. Tracking Algorithm Performance with sensor noise, actuator delay and Kalman filter: Latitude and Longitude Errors on Ground.

Fig. 23 .
Fig. 23.Tracking Algorithm Performance with Sensor Noise, Actuator Delay and Kalman Filter: Tracking Errors on Ground for 500 MC Runs.

Table 1
Planet Earth Constants.

Table 2
Monte Carlo simulation parameters.