# Patent

## Title
Orbit insertion device for artificial satellite aimed to explore a planet of unknown characteristics

## Abstract

An orbit insertion device for artificial satellite aiming to explore a planet with unknown characteristics. A plurality of sensors send a plurality of measured states of the artificial satellite. A planetary attraction constant estimator calculates a planetary attraction constant based on the plurality of measured states. A drag force coefficient estimator calculates a drag force coefficient based on the plurality of measured states. A rate of convergence calculation unit calculates a rate of convergence based on the drag force coefficient. A fictitious control calculation unit calculates a fictitious control based on the planetary attraction constant, the drag force coefficient and the rate of convergence. A normalized thrust calculation unit calculates a normalized thrust based on the fictitious control. A plurality of thrusters generate a thrust based on the normalized thrust.

## Background
BACKGROUND OF THE INVENTION
The present invention relates to an orbit insertion device for artificial satellite aimed to explore a planet of unknown characteristics such as planetary mass density, planetary radius and atmospheric composition.
U.S. Pat. No. 5,595,360 discloses an apparatus for translating a spacecraft from one orbit to another by generating a thrust of a force in timing sequences, both of which are predetermined based on known characteristics of a planet. U.S. Pat. No. 6,213,432 discloses a method for inserting a space vehicle into a target orbit of a planet through a spiral trajectory pre-calculated based on known characteristics of a planet. U.S. Pat. No. 6,341,749 B1 discloses an apparatus for translating a spacecraft from an orbit to another using a nonlinearly varying thrust minimizing transfer time and fuel calculated based on known characteristics of a planet.
The prior arts in U.S. Pat. Nos. 5,595,360, 6,213,432, and 6,341,749 B1, if applied to orbit insertion for artificial satellite aimed to explore a planet of unknown characteristics, fail to insert the artificial satellite into a target orbit since the prior arts are not able to calculate a thrust required by an orbit insertion without a pre-knowledge on the characteristics of the planet.
SUMMARY OF THE INVENTION
The present invention provides an orbit insertion device for artificial satellite aimed to explore a planet of unknown characteristics.
The orbit insertion device includes:

    
    
        a planetary attraction constant estimator receiving the plurality of measured states, and calculating the planetary attraction constant, which is a product of gravitational constant and planetary mass, as a polynomial of the radial coordinate, a radial acceleration and an angular speed, where the radial acceleration is a second order time derivative of the radial coordinate, and the angular speed is a first order time derivative of the angular coordinate;
        a drag force coefficient estimator receiving the plurality of measured states, and calculating the drag force coefficient, which is a drag force acting on the artificial satellite in radial direction with respect to a planetary center divided by a product of planetary mass, an opposite sign of radial speed and squared radial speed, as a polynomial of the radial speed and a function of the radial coordinate monotonically decreasing with respect to increasing radial coordinate;
        a rate of convergence calculation unit receiving the drag force coefficient, and calculating the rate of convergence, with which the radial coordinate of the artificial satellite approaches to radius of the target orbit, as a linear combination of the drag force coefficient divided by the radial speed and a small quantity;
        a fictitious control calculation unit receiving the planetary attraction constant, the drag force coefficient and the rate of convergence, and calculating the fictitious control, which is a centripetal acceleration required to propel the artificial satellite to the target orbit, as a linear combination of the planetary attraction constant divided by the radial coordinate squared, the drag force coefficient, a product of the rate of convergence and the radial speed, and a small quantity; and
        a normalized thrust calculation unit receiving the fictitious control, and calculating a normalized thrust, the thrust normalized with respect to a satellite mass, as a product of the radial coordinate and a first order time derivative of a square root of the fictitious control divided by the radial coordinate.

This construction utilizing the estimated planetary attraction constant, the estimated drag force coefficient and the fictitious control make possible an orbit insertion without a prior knowledge about the planet explored with more certainty than the prior arts. With this construction, the time required to complete an orbit insertion is minimized, since this construction allows a direct orbit insertion without passing through intermediate orbits before arriving at the target orbit.

    
    





BRIEF DESCRIPTION OF THE DRAWINGS
The novel features of the invention will be readily understood and appreciated by reference to the following detailed description taken in conjunction with the drawings, in which:
FIG. 1 is a block diagram of an orbit insertion device for artificial satellite in accordance with an exemplary embodiment of the invention;
FIG. 2 is a flowchart showing a processing flow of the orbit insertion device for artificial satellite in accordance with an exemplary embodiment of the invention;
FIG. 3 is a simulation result of an exemplary embodiment of the invention showing radial coordinate due to orbit insertion control with the present invention; and
FIG. 4 is a simulation result of an exemplary embodiment of the invention showing hyper plane variable.



DETAILED DESCRIPTION OF THE INVENTION
FIG. 1 is a block diagram of an orbit insertion device for artificial satellite in accordance with an exemplary embodiment of the invention. The artificial satellite includes: an artificial satellite structure 110, a plurality of sensors 120, an orbit insertion device 130 and a plurality of thrusters 140.
The artificial satellite structure 110 is a mechanical structure of an artificial satellite including a plurality of solar arrays, a plurality of propellant tanks, a plurality of propulsion systems that forms a basic shape of the artificial satellite, and mounts the plurality of sensors 120, the orbit insertion device 130 and the plurality of thrusters 140.
The plurality of sensors 120 measure radial coordinate and angular coordinate of the artificial satellite in a polar coordinate system with respect to planetary center, and send them to the orbit insertion device 130 as a plurality of measured states. Examples of the plurality of sensors 120 include but are not limited to static horizon sensor, scanning horizon sensor, and star sensor where the static horizon sensor and the scanning horizon sensor measure position of the artificial satellite with respect to the planet explored using infrared radiation from the planet, and the star sensor measures the position of the artificial satellite with respect to a star around which non-stellar objects including the planet explored are gravitationally bound in orbit.
The orbit insertion device 130 estimates planetary attraction constant which is a product of gravitational constant and planetary mass as a function of the radial coordinate and the angular coordinate, estimates drag force coefficient as a monotonic function of radial coordinate, calculates rate of convergence of the radial coordinate to target orbit radius as a function of estimated drag force coefficient, calculates fictitious control as a function of the estimated planetary attraction constant and the estimated drag force coefficient, calculates normalized thrust as a function of the fictitious control and sends a control signal to cause the plurality of thrusters 140 generate thrust, magnitude of which is a product of the normalized thrust and the satellite mass, in the direction of travel. Detailed description of the orbit insertion device 130 is given together with its functional configuration in what follows.
The orbit insertion device 130 further includes: a planetary attraction constant estimator 131, a drag force coefficient estimator 132, a rate of convergence calculation unit 133, a fictitious control calculation unit 134 and a normalized thrust calculation unit 135.
The planetary attraction constant estimator 131 receives the plurality of measured states, and calculates the planetary attraction constant, which is a product of gravitational constant and planetary mass, as a polynomial of the radial coordinate, a radial acceleration and an angular speed, where the radial acceleration is a second order time derivative of the radial coordinate, and the angular speed is a first order time derivative of the angular coordinate.
The drag force coefficient estimator 132 receives the plurality of measured states, and calculates the drag force coefficient, which is a drag force acting on the artificial satellite in radial direction with respect to a planetary center divided by a product of planetary mass, an opposite sign of radial speed and squared radial speed, as a polynomial of the radial speed and a function of the radial coordinate monotonically decreasing with respect to increasing radial coordinate.
The rate of convergence calculation unit 133 receives the drag force coefficient, and calculates the rate of convergence, with which the radial coordinate of the artificial satellite approaches to radius of the target orbit, as a linear combination of the drag force coefficient divided by the radial speed and a small quantity.
The fictitious control calculation unit receives the planetary attraction constant, the drag force coefficient and the rate of convergence, and calculates the fictitious control, which is a centripetal acceleration required to propel the artificial satellite to the target orbit, as a linear combination of the planetary attraction constant divided by the radial coordinate squared, the drag force coefficient, a product of the rate of convergence and the radial speed, and a small quantity.
The normalized thrust calculation unit 135 receives the fictitious control, and calculates the normalized thrust as a product of the radial coordinate and a first order time derivative of a square root of the fictitious control divided by the radial coordinate.
The plurality of thrusters 140 generate a thrust which is a product of the normalized thrust and a satellite mass in a direction of travel, and propels the artificial satellite to a target orbit.
Derived in what follows is a principle whereby the orbit insertion device 130 calculates the normalized thrust.
Equations of motion of an artificial satellite are given as (1) and (2).




  
    
      
        
          
            r
             ̈
          
          -
          
            
              
                θ
                .
              
              2
            
            ⁢
            r
          
        
        =
        
          {
          
            
              
                
                  -
                  
                    k
                    
                      r
                      2
                    
                  
                
              
              
                
                  r
                  ⪢
                  
                    r
                    0
                  
                
              
            
            
              
                
                  
                    -
                    
                      k
                      
                        r
                        2
                      
                    
                  
                  +
                  
                    d
                    r
                  
                
              
              
                
                  r
                  ≈
                  
                    r
                    0
                  
                
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  
  
    
      
        
          
            θ
             ̈
          
          ⁢
          r
        
        =
        
          f
          θ
        
      
    
    
      
        (
        2
        )
      
    
  




where r is a radial coordinate of the artificial satellite in a polar coordinate system with respect to planetary center (hereinafter referred to as “radial coordinate”) [m], θ is an angular coordinate of the artificial satellite in the polar coordinate system with respect to planetary center (hereinafter referred to as “angular coordinate”) [rad], k is a product of gravitational constant and planetary mass (hereinafter referred to as “planetary attraction constant”) [N·m2/kg] where the gravitational constant is 6.674×10−11 [N·m2/kg2], dr is a radial drag force due to atmosphere of a planet normalized with respect to satellite mass (hereinafter referred to as “normalized radial drag force” [m/s2], fθ is thrust of the artificial satellite in the direction of travel normalized with respect to satellite mass (hereinafter referred to as “normalized thrust”) [m/s2], and r0 is target orbit radius [m].

For r>>r0 the planetary attraction constant is estimated using (1) as (3).

{circumflex over (k)}=r2({dot over (θ)}2r−{umlaut over (r)})  (3)

where {circumflex over (k)} is estimated planetary attraction constant [N·m2/kg].

The normalized radial drag force dr in (1) is expressed as (4).

dr=−sgn({dot over (r)}){dot over (r)}2c(r), c(r)>0  (4)

where c(r) is a function of the radial coordinate monotonically decreasing with respect to increasing radial coordinate (hereinafter referred to as “drag force coefficient”) [N/m2].

The drag force coefficient is estimated using (1), (3), (4) as (5).




  
    
      
        
          
            c
            ^
          
          ⁡
          
            (
            r
            )
          
        
        =
        
          
            -
            
              sgn
              ⁡
              
                (
                
                  r
                  .
                
                )
              
            
          
          ⁢
          
            1
            
              
                r
                .
              
              2
            
          
          ⁢
          
            (
            
              
                r
                 ̈
              
              -
              
                
                  
                    θ
                    .
                  
                  2
                
                ⁢
                r
              
              +
              
                
                  k
                  ^
                
                
                  r
                  2
                
              
            
            )
          
        
      
    
    
      
        (
        5
        )
      
    
  




where ĉ(r) is estimated drag force coefficient [N/m2].

Hereinafter, the normalized thrust fθ in (2) is designed such that the trajectory of the artificial satellite converges to a circular orbit with the target orbit radius around the planetary center. A variable in (6) is introduced.

s={dot over (r)}+σr−σr0  (6)

where s is a hyper plane variable [m/s], vanishing of which causes the trajectory of the artificial satellite to converge to a circular orbit with the target orbit radius around the planetary center, and σ is a rate of convergence [s−1].

Next a candidate for Lyapunov function is introduced as in (7).

V=1⁄2s2>0, ∀s≠0

where V is the candidate for Lyapunov function [m2/s2].

Taking a first order time derivative of (7) yields (8).




  
    
      
        
          
            
              
                V
                .
              
              =
              
              ⁢
              
                s
                ⁢
                
                  s
                  .
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                s
                ⁡
                
                  (
                  
                    
                      
                        
                          θ
                          .
                        
                        2
                      
                      ⁢
                      r
                    
                    -
                    
                      k
                      
                        r
                        2
                      
                    
                    +
                    
                      d
                      r
                    
                    +
                    
                      σ
                      ⁢
                      
                        
                      
                      ⁢
                      
                        r
                        .
                      
                    
                  
                  )
                
              
            
          
        
      
    
    
      
        (
        8
        )
      
    
  




According to Lyapunov's second theorem the hyper plane variable s converges to zero if (9) holds.




  
    
      
        
          
            sgn
            ⁡
            
              (
              
                
                  
                    
                      θ
                      .
                    
                    2
                  
                  ⁢
                  r
                
                -
                
                  k
                  
                    r
                    2
                  
                
                +
                
                  d
                  r
                
                +
                
                  σ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    r
                    .
                  
                
              
              )
            
          
          =
          
            -
            
              sgn
              ⁡
              
                (
                s
                )
              
            
          
        
        ,
        
          
            ∀
            s
          
          ⁢
          
            ≠
            0
          
        
      
    
    
      
        (
        9
        )
      
    
  




Since (2) can be rewritten as (10), centripetal acceleration {dot over (θ)}2r is treated as a fictitious control to make (9) hold.




  
    
      
        
          f
          θ
        
        =
        
          r
          ⁢
          
            d
            
              d
              t
            
          
          ⁢
          
            
              
                
                  
                    θ
                    .
                  
                  2
                
                ⁢
                r
              
              r
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




Sufficient condition for (9) is derived from (9) as (11).




  
    
      
        
          
            
              
                
                  
                    θ
                    .
                  
                  2
                
                ⁢
                r
              
              <
              
                
                  k
                  
                    r
                    2
                  
                
                +
                
                  d
                  r
                
                +
                
                  σ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      r
                      .
                    
                    ⋀
                    σ
                  
                
              
              ≥
              
                -
                
                  
                    d
                    r
                  
                  
                    r
                    .
                  
                
              
            
          
          
            
              
                for
                ⁢
                
                  
                
                ⁢
                s
              
              >
              
                0
                ⋀
                
                  r
                  .
                
              
              >
              0
            
          
        
        
          
            
              
                
                  
                    θ
                    .
                  
                  2
                
                ⁢
                r
              
              <
              
                
                  k
                  
                    r
                    2
                  
                
                +
                
                  d
                  r
                
                +
                
                  σ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      r
                      .
                    
                    ⋀
                    σ
                  
                
              
              ≤
              
                -
                
                  
                    d
                    r
                  
                  
                    r
                    .
                  
                
              
            
          
          
            
              
                for
                ⁢
                
                  
                
                ⁢
                s
              
              >
              
                0
                ⋀
                
                  r
                  .
                
              
              <
              0
            
          
        
        
          
            
              
                
                  
                    θ
                    .
                  
                  2
                
                ⁢
                r
              
              <
              
                k
                
                  r
                  2
                
              
            
          
          
            
              
                
                  for
                  ⁢
                  
                    
                  
                  ⁢
                  s
                
                >
                
                  0
                  ⋀
                  
                    r
                    .
                  
                
              
              =
              0
            
          
        
        
          
            
              
                
                  
                    θ
                    .
                  
                  2
                
                ⁢
                r
              
              >
              
                
                  k
                  
                    r
                    2
                  
                
                +
                
                  d
                  r
                
                +
                
                  σ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    r
                    .
                  
                
              
            
          
          
            
              
                for
                ⁢
                
                  
                
                ⁢
                s
              
              <
              0
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




From (10) and (11), the normalized thrust fθ and the rate of convergence σ satisfying (9) are obtained as (12).




  
    
      
        
          
            
              
                f
                θ
              
              =
              
                r
                ⁢
                
                  d
                  
                    d
                    t
                  
                
                ⁢
                
                  
                    
                      f
                      r
                    
                    r
                  
                
              
            
          
          
            
              
            
          
        
        
          
            
              
                
                  f
                  r
                
                =
                
                  
                    
                      k
                      ^
                    
                    
                      r
                      2
                    
                  
                  +
                  
                    d
                    r
                  
                  +
                  
                    σ
                    ⁢
                    
                      
                    
                    ⁢
                    
                      r
                      .
                    
                  
                  -
                  
                    ɛ
                    f
                  
                
              
              ,
              
                σ
                =
                
                  
                    -
                    
                      
                        d
                        r
                      
                      
                        r
                        .
                      
                    
                  
                  +
                  
                    ɛ
                    σ
                  
                
              
            
          
          
            
              
                for
                ⁢
                
                  
                
                ⁢
                s
              
              >
              
                0
                ⋀
                
                  r
                  .
                
              
              >
              0
            
          
        
        
          
            
              
                
                  f
                  r
                
                =
                
                  
                    
                      k
                      ^
                    
                    
                      r
                      2
                    
                  
                  +
                  
                    d
                    r
                  
                  +
                  
                    σ
                    ⁢
                    
                      
                    
                    ⁢
                    
                      r
                      .
                    
                  
                  -
                  
                    ɛ
                    f
                  
                
              
              ,
              
                σ
                =
                
                  
                    -
                    
                      
                        d
                        r
                      
                      
                        r
                        .
                      
                    
                  
                  -
                  
                    ɛ
                    σ
                  
                
              
            
          
          
            
              
                for
                ⁢
                
                  
                
                ⁢
                s
              
              >
              
                0
                ⋀
                
                  r
                  .
                
              
              <
              0
            
          
        
        
          
            
              
                f
                r
              
              =
              
                
                  
                    k
                    ^
                  
                  
                    r
                    2
                  
                
                -
                
                  ɛ
                  f
                
              
            
          
          
            
              
                
                  for
                  ⁢
                  
                    
                  
                  ⁢
                  s
                
                >
                
                  0
                  ⋀
                  
                    r
                    .
                  
                
              
              =
              0
            
          
        
        
          
            
              
                f
                r
              
              =
              
                
                  
                    k
                    ^
                  
                  
                    r
                    2
                  
                
                +
                
                  d
                  r
                
                +
                
                  σ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    r
                    .
                  
                
                +
                
                  ɛ
                  f
                
              
            
          
          
            
              
                for
                ⁢
                
                  
                
                ⁢
                s
              
              <
              0
            
          
        
        
          
            
              
                
                  ɛ
                  f
                
                >
                0
              
              ,
              
                
                  ɛ
                  σ
                
                >
                0
              
            
          
          
            
              
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




where fr is the fictitious control [m/s2], εf and εσ are small quantities.

With the normalized thrust fθ and the rate of convergence σ in (12), the hyper plane variable in (6) converges to zero and the radial coordinate r converges to the target orbit radius r0.
Next a processing flow of the orbit insertion device for artificial satellite is described using FIG. 2. FIG. 2 is a flowchart showing a processing flow of the orbit insertion device for artificial satellite in accordance with an exemplary embodiment of the invention. The plurality of sensors 120 measure radial coordinate r and angular coordinate θ of the artificial satellite in a polar coordinate system with respect to planetary center, and send them to the orbit insertion device 130 as a plurality of measured states (S201).
In the orbit insertion device 130, the planetary attraction constant estimator 131 estimates the planetary attraction constant using (3) and send it to the fictitious control calculation unit 134 as the estimated planetary attraction constant {circumflex over (k)} (S202).
The drag force coefficient estimator 132 estimates the drag force coefficient using (5) and sends it to the rate of convergence calculation unit 133 and the fictitious control calculation unit 134 as the estimated drag force coefficient ĉ(r) (S203).
The rate of convergence calculation unit 133 calculates the rate of convergence σ using (12) and sends it to the fictitious control calculation unit 134 (S204).
The fictitious control calculation unit 134 calculates the fictitious control fr using (12) and sends it to the normalized thrust calculation unit 135 (S205).
The normalized thrust calculation unit 135 calculates the normalized thrust fθ using (12) and sends it to the plurality of thrusters 140 (S206).
The plurality of thrusters 140 generate a thrust which is a product of the normalized thrust fθ and a satellite mass in a direction of travel (S207).
The plurality of sensors 120 measure radial coordinate r and angular coordinate θ of the artificial satellite, and steps S203 through S207 are repeated to propel the artificial satellite to the target orbit (S208).
Hereinafter, simulation results of aforementioned control processes of the orbit insertion device for artificial satellite in accordance with an exemplary embodiment of the invention are described in detail. The simulation assumes a space mission wherein an artificial satellite attempts to enter into a stable orbit around a planet. Parameters used in the simulation are given in (13).

mp=5.68×1026[kg], G=6.674×10−11[N·m2/kg2], k=mpG[N·m2/kg], 

rp=60×106[m], r0=1.01rp[m], r(0)=1.02r0[m], T=0.1[s]  (13)

where mp is the planetary mass [kg], rp is planetary radius [m], r(0) is an initial condition of the radial coordinate [m], and T is a sampling time of the orbit insertion device 130 [s].

FIG. 3 is a simulation result of an exemplary embodiment of the invention showing radial coordinate due to orbit insertion control with the present invention. In FIG. 3 a solid line represents altitude of the artificial satellite (hereinafter referred to as “the artificial satellite altitude”) with the present invention, a broken line represents the one with the prior art and a dashed line represents altitude of the target orbit (hereinafter referred to as “the target orbit altitude”) where altitude is the difference between the radial coordinate r and the planetary radius rp. The prior art in this simulation utilizes an orbit transfer control described in U.S. Pat. No. 5,595,360, wherein propulsion thrusters are fired in timing sequences predetermined based on a planetary mass incorrectly estimated as 3.12×1022 [kg] in relation to measured position of the artificial satellite relative to the planet. Such incorrect estimation is due to limited information available about the planetary density, and in this simulation the planetary density is estimated 5.49×10−5 times as small as a true value. As in FIG. 3 the artificial satellite altitude with the present invention approaches and converges to the target orbit altitude, since the present invention is based on the planetary attraction constant correctly estimated by (3), and the normalized radial drag force correctly estimated by (4) and (5). On the other hand, the artificial satellite altitude with the prior art continues to decrease after reaching the target orbit altitude at 492 [s] until the artificial satellite plummets to the ground of the planet at 683 [s], since the prior art generates thrust based on the planetary attraction constant derived from incorrectly pre-calculated planetary mass that is insufficient to lift the artificial satellite back to the target orbit.
FIG. 4 is a simulation result of an exemplary embodiment of the invention showing hyper plane variable. In FIG. 4 a solid line represents the hyper plane variable with the present invention and the dashed line represents a target value of the hyper plane variable with which the candidate of Lyapunov function in (7) monotonically decreases. As in FIG. 4 the hyper plane variable reaches the target value at 492 [s] and stays in the vicinity of the target value.
Although, in the foregoing embodiment, a circular target orbit is considered, the present invention can also be applied to an elliptical orbit with arbitrary eccentricity by replacing the target orbit radius r0 with a target distance from a center of the planet, which is a function of angular coordinate θ.
The orbit insertion device for artificial satellite in the foregoing embodiment can be realized as an electrical/electronic/programmable electronic system. Examples of the electrical/electronic/programmable electronic system include but are not limited to a control system comprising an application specific integrated circuit (ASIC) and a microcontroller.
In the above simulation, the planet, an orbit of which the artificial satellite is deployed on, is a solid planet. The present invention can be applied without any modifications to a gas planet.
While, in the foregoing embodiment, a single stage orbit insertion, wherein the artificial satellite is inserted to the target orbit without passing through intermediate orbits, is discussed, successive application of aforementioned control procedure allows the present invention to propel the artificial satellite to the target orbit by a multiple stage orbit insertion wherein the artificial satellite passes through a plurality of intermediate orbits before arriving at the target orbit.
The described embodiments of the current invention are only exemplary and numerous variations thereof apparent to those skilled in the art are intended to be within the scope of the present invention as defined in the following claims.

## Claims

1. An orbit insertion device for an artificial satellite comprising an artificial satellite structure including a plurality of solar arrays, a plurality of propellant tanks, and a plurality of propulsion systems that forms a basic shape of the artificial satellite; a plurality of sensors for measuring a radial coordinate and an angular coordinate of the artificial satellite in a polar coordinate system with respect to a planetary center and sending them as a plurality of measured states that is a plurality of electrical signals representing a measured radial coordinate and a measured angular coordinate of the artificial satellite; and a plurality of thrusters for generating a thrust for propelling the artificial satellite to a target orbit, the orbit insertion device comprising:
a planetary attraction constant estimator receiving the plurality of measured states, and calculating a planetary attraction constant, which is a product of a gravitational constant and a planetary mass, as a polynomial of the radial coordinate, a radial acceleration and an angular speed, where the radial acceleration is a second order time derivative of the radial coordinate;
a drag force coefficient estimator receiving the plurality of measured states, and calculating a drag force coefficient, which is a drag force acting on the artificial satellite in radial direction with respect to the planetary center divided by a product of the planetary mass, an opposite sign of a radial speed and a squared radial speed, as a polynomial of the radial speed and a function of the radial coordinate monotonically decreasing with respect to increasing radial coordinate;
a rate of convergence calculation unit receiving the drag force coefficient, and calculating a rate of convergence, with which the radial coordinate of the artificial satellite approaches to a radius of the target orbit, as a linear combination of the drag force coefficient divided by the radial speed and a small quantity that is a positive number smaller than an absolute value of the drag force coefficient divided by the radial speed;
a fictitious control calculation unit receiving the planetary attraction constant, the drag force coefficient and the rate of convergence, and calculating a fictitious control, which is a centripetal acceleration required to propel the artificial satellite to the target orbit, as a linear combination of the planetary attraction constant divided by the radial coordinate squared, the drag force coefficient, a product of the rate of convergence and the radial speed, and a small quantity that is a positive number smaller than an absolute value of the planetary constant divided by the radial coordinate squared; and
a normalized thrust calculation unit receiving the fictitious control, and calculating a normalized thrust, the thrust normalized with respect to a satellite mass, as a product of the radial coordinate and a first order time derivative of a square root of the fictitious control divided by the radial coordinate.