Fran De Aquino

The Bipolar Linear Momentum transported by the Electromagnetic Waves: Origin of the Gravitational Interaction

Keywords: Gravity, Gravitation, Electromagnetic Waves, Radiation Pressure

Besides energy, the electromagnetic waves transport linear momentum. Then, if this momentum is absorbed by a surface, pressure is exerted on the surface. This is the so-called Radiation Pressure. Here we show that this pressure has a negative component (opposite to the direction of propagation of the radiation) due to the existence of the negative linear momentum transported by the electromagnetic waves. This fact leads to an important theoretical discovery: the velocity of the electromagnetic waves in free space is not a constant. In addition, a generalized equation of the Newton's law of Gravitation, is deduced starting from the concept of negative radiation pressure applied on the Gravitational Interaction.

Introduction

Electromagnetic waves transport energy as well as linear momentum. Then, if this momentum is absorbed by a surface, pressure is exerted on the surface. Maxwell showed that, if the incident energy U is totally absorbed by the surface during a time t , then the total momentum q transferred to the surface is v U q = [START_REF] Henry | Radiation Pressure[END_REF]. Then, the pressure, p (defined as force per unit area F A ), exerted on the surface, is given by ( ) ( )

1 1 1 1 A dt dU v v U dt d A dt dq A A F p = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = = =
whence we recognize the term ( ) A dt dU as the radiation power density, , (in watts/m D 2 ) arriving at the surface 1 . Thus, if c v = the radiation pressure exerted on the surface is ( )

2 c D p =
Here we show that this pressure has a negative component (opposite to the direction of propagation of the photons) due to the existence of the negative linear momentum transported by the photons. This fact leads to an important theoretical discovery: the velocity of the electromagnetic waves in free space is not a constant. In addition, a generalized equation of the Newton's law of Gravitation is deduced starting from the concept of negative radiation pressure applied on the Gravitational Interaction.

.

Theory

The energy of a harmonic oscillator is quantized in multiples of , and given by hf 1 This value is also called of Poynting vector. . This value is called energy of the zero point. Thus, the energy of the harmonic oscillator, at equilibrium with the surrounding medium, does not tend to zero when temperature approaches to absolute zero, but stays equal to . 0 E In the particular case of massless oscillators (photons, for example), does not correspond to 0 E the lowest value of the energy, which the oscillator can have, because, when temperature approaches to absolute zero the oscillator frequency becomes dependent of the temperature T , as show the well-known expression of the thermal De Broglie wavelength ( ) Λ for massless particles [START_REF] Kittel | Thermal Physics[END_REF][START_REF] Yan | [END_REF], which is given by ( ) , in the case of the photon, for example, will be a fraction of the value hf , can express satisfactorily the mentioned fraction. Thus, according to Eq. ( 4), we can write that 2 ( )

5 3 1 2 1 kT f f hf g π = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛
Above , the photon absorbs energy from the surrounding medium c T2 [START_REF] Beiser | Concepts of Modern Physics, Portuguese version[END_REF], and its energy becomes equal to hf . Therefore, the energy absorbed by the photon is ( )

6 2 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = g f f hf hf U ( ) g f f hf 2 1 hf photon T =T c photon T >T c Fig. 1 -Above c
T , the photon absorbs energy from the surrounding medium, and its energy becomes equal to hf .

The absorbed energy is that thrust the photon and gives to it its velocity . Consequently, the momentum transported by the photon with velocity will be expressed by v r q v r ( )

7 2 1 1 2 1 1 2 1 1 2 1 1 2 1 r g g g g g n c hf f f v c c hf f f c c v hf f f v hf f f v f f hf hf v U q r r r r r r r ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = =
Equation above shows the existence of a bipolar linear momentum transported by the electromagnetic waves. For the resultant momentum transported by the photon is positive, i.e., If this momentum is absorbed by a surface, pressure is exerted on the surface,

g f f 2 <
in the same direction of propagation of the photon. These photons are well-known. However, Eq. ( 7) point to a new type of photons when . In this

case g f f 2 = 0 = q
, i.e., this type of photon does not exert pressure when it incides on a surface. What means that it does not interact with matter. Obviously, this corresponds to a special type of photon, which we will call of neutral photon. Finally, if the resultant momentum transported by the photon is negative. If this momentum is absorbed by a surface, pressure is exerted on the surface,

g f f 2 >
in the opposite direction of propagation of the photon. This special type of photon will be denominated of attractive photon.

The quantization of gravity shows that the gravitational mass m g and inertial mass m i are correlated by means of the following factor [START_REF] De Aquino | Mathematical Foundations of the Relativistic Theory of Quantum Gravity[END_REF]:

( ) 8 1 1 2 1 2 0 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Δ + - = = c m p m m i i g χ
where is the rest inertial mass of the particle and 0 i m p Δ is the variation in the particle's kinetic momentum.

Another important equation obtained in the quantization theory of gravity is the new expression for the momentum of a particle with gravitational mass and velocity , which is given by

q g M v ( ) 9 v M q g r r = where 2 2 1 c v m M g g - =
; is given by Eq.( 8), i.e., .

g m i g m m χ =
By comparing Eq. ( 9) with ( 7) we obtain ( )

10 2 1 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = g g f f c M hf c v
Mass-energy equivalence principle states that anything having mass has an equivalent amount of energy and vice versa. In the particular case of photons, the energy of the photons, , has a corresponding equivalent mass, hf E = g M , given by its energy E divided by the speed of light squared c², i.e., ( )

11 2 c E M g ≡ or 3 ( ) 12 2 hf c M g ≡
Considering this expression, Eq. ( 10) can be rewritten as follows ( ) 13) shows that the speed of the electromagnetic waves in free space is not a constant. Thus, also the speed of light is not a constant.

13 2 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = g f f c v Equation (
Theories proposing a varying speed of light have recently been widely proposed under the claim that they offer a solution to cosmological puzzles [7 , 8].

It is known that the interactions are communicated by means of the changing of "virtual" quanta. The maximum velocity of these quanta is a constant called maximum velocity of propagation of the interactions. Currently, it is assumed that this velocity is equal to the velocity of the electromagnetic waves in free space ( ) c . This is the reason of the constant c to appear in the relativistic factor (Eq. ( 10)). However, Eq. ( 13) shows that the velocity of the electromagnetic waves in free space is not a constant. In addition, Eq. ( 13) shows that for the velocity of the photon is imaginary.

g f f 2 >
This means that the attractive photons are virtual photons. Now we will apply the concept of negative radiation pressure, here developed, to the Gravitational Interaction.

According to Eq. ( 7), the resultant momentum transported by the photons with frequency is negative. 6), are respectively given by ( ) 14) and (15) can be rewritten as follows

14 2 1 f f f hf N P g A A ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ( ) 15 2 1 f f f hf N P g B B ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - =
( ) ( ) ( ) 16 1 1 2 1 2 2 2 f c M f E f f f hf N P A g c A V A g g A A ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ - ≅ = = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ≅ ( ) ( ) ( ) 17 1 1 2 1 2 2 2 f c M f E f f f hf N P B g c B V B g g B B ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ - ≅ = = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ≅
where and are respectively the total energies of the particles ( )

A g E ( ) B g E
A and B , given respectively by:

( ) ( ) 2 2 2 2 1 c V c m c M E A A A g A g - = = and ( ) ( ) 2 2 2 2 1 c V c m c M E B B B g B g - = = [6];
is the velocity of the particle

A V
A in respect to the particle B and is the velocity of the particle B V B in respect to the particle A . Obviously,

B A V V = .
Thus, if r is the distance between the mentioned particles, then the power densities of the attractive radiation in A and B are respectively, given by ( ) ( ) ( ) It is known that the electric force on an electric charge A, due to another electric charge B is related to the product of the charge of A by the flux density on A, due to the charge of B,

18 1 1 4 1 1 4 2 0 2 2 2 2 2 2 2 r M r f c M r P D A g c V A g c V A A ξ π π ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = = ( ) ( ) ( ) 19 1 1 4 1 1 4 2 0 2 2 2 2 2 2 2 r M r f c M r P D B g c V B g c V B B ξ π π ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = = where B A V V V = = ; 0 ξ is expressed by ( ) 20
( A q ) ( ) 2 4 r q D B B π = , i.e., A B BA q D F r r ∝
. By analogy, the charge A exerts an opposite electric force on the charge B, which is related to the product of the charge of B ( by the flux density on B, due to the charge of A,

) B q ( ) 2 4 r q D A A π = , i.e., AB AB F F r r - = , B A AB q D F r r ∝
. These proportionalities are usually written by means of the following equations:

( ) In the case of the forces produced by the action of attractive photons emitted reciprocally from the particles A and B, their expressions can be deduced by using the same argument previously shown in order to obtain the expressions of the electric forces and magnetic forces. That is, the force exerted on the particle A (whose gravitational mass is ), by another particle B (whose gravitational mass is

21 4 2 0 0 0 μ πε ε ε r q q q D q D F F B A A B B A BA AB = - = = - = r r r r
) is related to the product of by the flux density (power density) on A, due to the mass , ( )

A g M ( ) B g M (A g M ) ( ) B g M ( ) ( ) 19 . . 4 2 Eq See r P D B B π = , i.e., ( ) A g B BA M D F r r ∝
. By analogy, the particle A exerts an opposite force on the particle B, which is related to the product of the mass

B g A AB M D F r r ∝
. Thus, we can write that 18) and ( 19) into Eq. ( 23) we get

( ) ( ) ( ) 23 0 0 μ μ k M D k M D F F B b B A g A BA AB r r r r - = = - = Substitution of Eqs. (
( ) ( ) ( ) 24 11 2 0 0 2 2 μ ξ r M M k F F B g A g c V BA AB ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = - = r r For 0 = V Eq. (24) reduces to ( ) 25 ˆ2 2 0 0 μ μ ξ r m m G r m m k F F B A B A BA AB = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = - = r r
where is the Universal constant of Gravitation.

27 4 2 0 0 G f c G k π ξ = =
From the above exposed, we can then conclude that the gravitational interaction is caused by the interchange of virtual photons with frequencies (attractive photons). In this way, the called graviton must have spin 1 and not 2. Consequently, the gravitational forces are also gauge forces because they are yielded by the exchange of virtual quanta of spin 1, such as the electromagnetic forces and the weak and strong nuclear forces. , we get ( )

28 1 1 2 1 2 0 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = = cv m U m m i i g χ
Substitution of Eq. ( 6) and (13) into Eq. ( 28) gives ( )

29 1 2 1 1 2 1 2 2 0 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = g i i g f f c m Nhf m m
In case of protons, for example, the number of attractive photons emitted from a proton, , can be expressed by means of the following relation:

p N ( )

0 0 i p p m m N N =
, where is the number of attractive photons emitted from the particle with mass ; is the rest inertial mass of the proton. In this case, the equation (29) will be rewritten as follows:

N 0 i m 0 p m ( ) 30 1 2 1 1 2 1 1 2 1 1 2 1 2 2 0 2 2 0 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = g i g p p p gp f f c m hf N f f c m hf N m m
In the case of, electrons, for example, the number of attractive photons emitted from an electron , can be expressed by means of the following relation: ), where is the rest inertial mass of the electron. In this case, the equation (29) will be rewritten as 0 e m ( )

31 1 2 1 1 2 1 1 2 1 1 2 1 2 2 0 2 2 0 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = g i g e e e ge f f c m hf N f f c m hf N m m
This is exactly the same expression for the proton (Eq. (30).

In the case of, neutrinos, for example, the number of attractive photons, , emitted from a neutrino can be expressed by means of the following relation:

n N ( 0 0 i n n m m N N =
), where is the rest inertial mass of the neutrino. In this case, the equation (29) will be rewritten as follows:

0 n m ( ) 32 1 2 1 1 2 1 1 2 1 1 2 1 2 2 0 2 2 0 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = g i g n n n gn f f c m hf N f f c m hf N m m
Also, in this case, the obtained expression is exactly the same expression for the proton and the electron (the term 0 i m N is a constant). In short, the result is the same for any particle with non-null mass.

By solving equation below ( ) ( )

33 1 2 1 1 2 1 1 2 1 1 2 1 2 tan 2 2 0 2 2 0 χ = ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - g t cons i g i f f f c m h N f f c m hf N
35 1 2 1 1 2 1 max max 2 max tan 2 2 0 χ = ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - g t cons i f f f c m h N 4 3 4 2 1
By solving this equation, we obtain: ( )

36 4 5 6 2 1 max 2 max 2 4 2 0 max 2 max ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - χ χ h c m f f f i g
The energy density at a distance D r from the mentioned particle can be expressed by: ( )( ) ( )

i i ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = χ χ π
By comparing Eq. (38) with equations ( 18) and (19), we can conclude that ( )

39 4 5 6 4 0 max 2 max 4 0 ξ χ χ π = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ h c m i
By comparing this equation with Eq.(20), we obtain ( )

40 4 5 6 2 0 max 2 max max h c m f i ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = χ χ
Now, taking Eq. ( 8), where the term c m p i0 Δ is putted in the following form: should be very close it, because the spectrum of the attractive photons should make limit with the gamma rays spectrum (See Fig. 2). Thus, we can write that It is very unlikely that there are gamma rays in the Nature with frequency much greater than the aforementioned value, but if they exist, they would only show that the value of would be situated above the value indicated by Eq (43).

( )( ) 2 1 2 2 0 1 - - = Δ c v c v c m p i [6], we get ( ) 41 1 1 1 2 1 2 2 2 0 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - = = c v c v

Gamma rays X rays

Gravific photons 2f g 3 The largest air shower detected is from a particle of around but this is thought to be from a cosmic ray particle rather than photon.

  frequency of oscillation, and h is the Planck's constant [

  lowest value of the energy hf 2 1

  absolute zero. The mentioned fraction must be only related to the frequency , and a frequency limit, , whose value must be extremely large. Just a simple algebraic form, the

  photons

ε

  is the so-called permissivity constant for free space () magnetic force that a magnetic pole B exerts on another magnetic pole A is related to the product of the pole intensity of the pole A by the flux density on A, due to the pole intensity of the pole B, pole A exerts an opposite magnetic force on the pole B, which is related to the product of the pole intensity of the pole B by the flux density on B, due to the pole intensity of the pole A,

ξ

  given by Eq. (20) into this expression, we obtain ( )

  becomes equal to . But, the minimum frequency of the gravific photons is very close to , then we can write that f terms of the Eq. (34) become approximately equal to zero. On the other hand, for 1 min = = N N (one gravific photon) the value of in Eq. (34) becomes equal to (the maximum frequency of the gravific photons) and Eq. (33) shows that, in this case,

.

  the maximum frequency of the gravific photons (See Eq. (40)) is very greater than f Thus, we can define the frequency spectrum of the gravific photons by means of the following expression:

  This expression shows then that the frequency spectrum of the gravific photons must be above the spectrum of the gamma rays (neutrino mass:

g f 2 Fig. 2 -

 22 Fig. 2 -The Gravific Photons Spectrum (above Hz f g 28 10 2 ≈

  If this momentumNow consider two particles A and B with

					f 2 >	f	g
	is absorbed by a surface, pressure is exerted on
	the surface,	in the opposite direction of
	propagation of the photon.
	gravitational masses	m	A	and	m	B	, respectively.
	If both particles emit attractive radiation
	(	f 2 >	f	g	), then the powers	A P	and	B P	emitted
	from	A and B , according to Eq. (

In order to provide the equilibrium the harmonic oscillator absorbs energy from the surrounding medium.


Zaaraoui
email: hind.zaaraoui@orange.com

Zwi Altman
email: zwi.altman@orange.com

Sana Ben Jemaa
email: sana.benjemaa@orange.com

Eitan Altman
email: eitan.altman@sophia.inria.fr

Tania Jimenez
email: tania.jimenez@univ-avignon.fr

Heuristic approach for forecast scheduling

Keywords: Forecast scheduler, alpha-fair, high mobility, Radio Environment Maps, geo-localized measurements

. This paper proposes a fast heuristic for the FS problem based on two FS users' scheduling. Furthermore, it is shown that in the case of two users, the FS problem can be solved analytically, making the heuristic computationally very efficient. Numerical results illustrate the throughput gain brought about by the scheduling solution. a

I. INTRODUCTION

The concept of FS has been introduced in [START_REF] Zaaraoui | Forecast scheduling for mobile users[END_REF], where it is shown that the knowledge of the present and future rates along the users' trajectories can be exploited by the scheduler in order to significantly improve the average user throughput. The throughput gain is achieved by exploiting long term time and spatial diversity along the users' trajectory. The FS allocation is posed as a convex optimization problem that can be solved using fast convex optimziation solvers. The randomness of the traffic, i.e. arrival and departure of communications on the one hand and randomenss in the user trajectories on the other hand can be incorporated into the FS solution [START_REF]Forecast scheduling and its extensions to account for random events[END_REF].

The SINR or rate along the mobile users' trajectories can be provided using a REM. The generation of REMs from Geo-Localized Measurements (GLM) has been the topic of recent research activity in both industry and a This work has been partially carried out in the framework of IDEFIX project, funded by the ANR under the contract number ANR-13-INFR-0006. academia and is decribed in more details in Section II. While abundant work on REM construction has been published, little work is available on how to exploit GLM and REMs to better manage the network.

The reasons for the growing interest in exploiting GLM are the following: First, it allows network operators to better assess the actual Quality of Service (QoS) of mobile users at any location; Second, it provides significant levers to optimize the network performance; and third, it allows personalized, user centric type of optimization. GLM can feed new Radio Resource Management (RRM) algorithms, and in general, a management entity to (self-) configure, optimize and troubleshoot the network. The use of GLM corresponds to a general trend in 5G networks, namely the exploitation of data from different sources in order to improve the network operation.

The purpose of this paper is to introduce a fast heuristic for the problem of calculating the FS allocation. Such a heuristic can simplify the implementation of the FS in real equipment. The contributions of this paper are the following:

• Provides an analytical solution for the FS allocation problem for the case of two mobile users; • Develops a computationally fast FS allocation heuristic based on the solution for two mobile users. The paper is organized as follows: Section II provides an overview of REMs and illustrates how the FS can exploit the REM to achieve throughput gain. Section III recalls the basic FS formulation and provides the main results for a closed form allocation rule for the case of two users. Section IV presents a two users based FS heuristic. Numerical results are described in Section V followed by concluding remarks in Section VI.

II. RADIO ENVIRONMENT MAPS

With the availability of Global Positioning System (GPS) information in most mobile terminals, it has become possible to retrieve radio measurements together with their corresponding location. GLMs can be reported to the network using User Equipments (UEs) with chipsets that implement the Minimization of Drive Testing (MDT) feature standardized by the 3rd Generation Partnership Project (3GPP) [START_REF] Hapsari | Minimization of drive tests solution in 3GPP[END_REF]. GLM comprises a rich source of information that allows network operators to optimize the network with finer granularity than the classical cell level, and to go one step further towards personalized optimization of the end-user experience.

As the measurements are performed by UEs, the operator does not have the control on the positions where these measurements are performed. Hence the first challenge has been to complement the information provided by GLM namely predicting the considered quantity in the entire area of interest and to build the complete REM.

The concept of REM was first introduced in cognitive radio and TV whitespace database systems [START_REF] Fette | Cognitive radio technology[END_REF], [START_REF] Zhao | Network support-The radio environment map[END_REF]. In mobile networks, the main research has focused on building reliable maps that could be trusted by radio engineers and that could be integrated in automated optimization processes and tools.

Several works ([6], [START_REF] Konak | A kriging approach to predicting coverage in wireless networks[END_REF], [START_REF] Riihijarvi | Characterization and modelling of spectrum for dynamic spectrum access with spatial statistics and random fields[END_REF]) have considered the geostatistical interpolation technique known as Kriging [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty[END_REF] to build REMs. This technique exploits the spatial correlation of the signal propagation to provide higher prediction accuracy than classical propagation models. The accuracy of the prediction can be enhanced by increasing the density of the measurement samples. However, the computational complexity of the algorithm increases rapidly with the number of measurement points n, namely as O(n 3 ). The reduction of the computational complexity has been among the main research challenges. The Fixed Rank Kriging (FRK) [START_REF] Cressie | Fixed rank kriging for very large spatial data sets[END_REF] solution brings the computational complexity of the REM construction to a linear level with respect to n, while maintaining high quality. The latter can be measured in terms of the root-mean square error (RMSE) of the interpolated GLM with respect to a reference test set of GLM for example. The coverage map is obtained using a ray tracing propagation model.

Figure 1(a) shows GLM of the received signal strength over 10 percent of the meshes in a cell area, and the predicted coverage area obtained using the FRK algorithm in Figure 1(b). In this example, the FRK is able to build the entire map with a root-mean square error not exceeding 3 dB.

Once the received signal strength coverage map is built with a satisfactory accuracy, it can serve as input to construct other radio measurement maps such as the SINR map.

III. FORECAST SCHEDULING MODEL AND ANALYTICAL RESULTS

In this Section we first briefly recall the basic formulation of the FS as presented in [START_REF] Zaaraoui | Forecast scheduling for mobile users[END_REF]. Then, using Karush-Kuhn-Tucker (KKT) conditions, we provide the main results for the closed form solution for n = 2 users. Due to space limitation, the complete formulation is given in a supplementary document [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF].

A. Basic formulation

Consider a macro-cell (Base Station (BS)) surrounded by interfering BSs. We suppose that the SINR of the users in mobility at any location is provided to the BS by a REM. Consider n full buffer users moving at a constant speed during a time interval T -the scheduling period, over which n is considered constant. Suppose that time is in a discrete space: t ∈ {1, 2, .., T } = [|1, T |] and let i denote the user number, i ∈ {1, 2, .., n} = [|1, n|]. We suppose that during the scheduling duration there are no arrivals or departures of users (it is noted that this assumption can be relaxed [START_REF]Forecast scheduling and its extensions to account for random events[END_REF]).

A scheduling period T (typically of the order of seconds) is divided into scheduling time slots e.g. of 1ms, during which the bandwidth is shared among the scheduled users. In practice, the time resolution of the FS is longer than the system time slot (e.g. of 1ms), and is referred to hereafter as time slot (see Section V). Let a i (t) denote the bandwidth proportion allocated to the user i at time t, a i (t) ∈ [0, 1] according to the scheduling strategy, where ∀t, n i=1 a i (t) = 1, and W is the total bandwidth. Using the Shannon equation, the throughput as a function φ of the SINR of user i reads

ai(t)φ(SIN Ri(t)) = ai(t)W log 2 (1 + SIN Ri(t)). (1)
Denote by S t i the predicted SINR (i.e. the one provided by the REM). The FS allocation policy is defined by the following optimization problem, with α = 1:

maximize : f (a) = n i=1 ( T t=1 ai(t)φ(S t i )) 1-α 1 -α (2) subject to : ∀i, ∀t, ai(t) ≥ 0 ∀t, n i=1 ai(t) = 1.
α is a positive fixed parameter of the optimization problem. For α → 1, the optimization problem with the same constraints reads:

maximize : f (a) = n i=1 log( T t=1 ai(t)φ(S t i )). (3)
Both equations ( 2) and (3) have concave functions f for α ≥ 0, and can be solved using convex optimization (e.g. CVX solver [START_REF] Research | CVX: Matlab software for disciplined convex programming, version 2.0[END_REF], see more details in Section V). The size of the optimization problem is defined by the number of unknown variables, namely n × T (with T being the number of time slots).

The interpretation of ( 2) and ( 3) is the following: resources are shared fairly among the users according to the data-rate variation along their trajectories. For example, if a user has a large enough coverage hole in his future trajectory, the FS may allocate to this user as much data as possible before reaching the coverage hole so as to remain fair with respect to the other users.

B. KKT conditions

We first derive the KKT conditions for the problem (2). The optimization problem (2) has a nonlinear objective function f with regular equality and inequality conditions i.e. differentiable constraint functions. If the objective and constraint functions (2) are continuously differentiable at any a = (a 1 (1), ...a 1 (T ), a 2 (1), ..., a n (T )) ∈ R nT , then there exist multipliers λ k,j and ν j , where k ∈ [|1, n|] and j ∈ [|1, T |], called KKT multipliers ( [START_REF] Boyd | Convex optimization[END_REF], Chap.5) with the following Lagrangian function:

L(a, ν, λ) = f (a)+ T j=1 ν j ( n k=1 a k (j)-1)+ k,j λ k,j a k (j), (4)
where λ k,j ≥ 0.

We define the Lagrange dual function as the maximum value of L over a. Let a * maximize the Lagrangian function (4) for the optimal multipliers λ * k,j and ν * j , where k ∈ [|1, n|] and j ∈ [|1, T |], and therefore its gradient is null at this point:

∇L(a * , ν * , λ * ) = 0. (5)
From ( 5) one obtains for all i ∈ [|1, n|] and t ∈ [|1, T |] the corresponding KKT conditions:

φ(S t i )( T j=1 a * i (j)φ(S j i )) -α + ν * t + λ * i,t = 0, (6)
a * i (t) ≥ 0, (7)
n k=1 a * k (t) = 1, (8)
λ * i,t a * i (t) = 0 (9) λ * i,t ≥ 0. ( 10
)
We note from (6) that for any time t, for all users i and w with w = i φ(S t i )(

T j=1 a * i (j)φ(S j i )) -α + λ * i,t = φ(S t w )( T j=1 a * w (j)φ(S j w )) -α + λ * w,t (11)
and also from [START_REF] Ojaniemi | Optimal field measurement design for radio environment mapping[END_REF], for all user i, we have for any time t and u:

λ * i,t + ν * t φ(S t i ) = λ * i,u + ν * u φ(S u i ) (12)
as ν * t does not depend on the users, and T j=1 a * i (j)φ(S j i ) does not depend on the time. Equality [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF] gives explicitly the resource balancing among users at each time relatively to α in the sense of equalizing the two expressions for each user.

C. Analytical solution for n = 2 users

In the calculations that follow, we rewrite a = a * , λ = λ * and ν = ν * . We assume two users in the cell, i.e. n = 2 in problem [START_REF]Forecast scheduling and its extensions to account for random events[END_REF].

KKT analysis shows that there exists at most one time slot K ∈ [1, T ] for which two users can be scheduled simultaneously [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF]. For all other times in [1, T ] one user is scheduled at a time. We formulate hereafter the more frequent case where there is no K ∈ [1, T ] for which two users scheduled at the same time. The complete analysis is available in [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF].

We suppose in this section that at time u = 1 user 1 is scheduled and therefore λ 1,1 = 0. We need to verify this assumption, namely that user 1 is indeed selected at time u = 1 (see Theorem 2). For all time t, the sign of λ 1,tλ 2,t determines the scheduling decision: if it is positive then user 2 is scheduled at time t (λ 1,t ≥ λ 2,t ⇐⇒ a 2 (t) ≥ a 1 (t) by equations (9-10)). One of the lambdas must be null and the other one positive for the case of two users.

Using equation ( 12), as λ 1,1 = 0, we have:

λ 1,t -λ 2,t = ν 1 ψ t/1 1 -ν 1 ψ t/1 2 -λ 2,1 ψ t/1 2 (13)
where

ψ t/u i = φ(S t i ) φ(S u i ) . We have then: λ 1,t -λ 2,t = ν 1 (ψ t/1 1 -ψ t/1 2 ) -λ 2,1 ψ t/1 2 (14)
As λ 2,1 ψ t/1 2 > 0 we divide equation ( 14) by λ 2,1 ψ t/1 2 and hence the sign studied is the same as the sign of

ψ t/1 2 -ψ t/1 1 ψ t/1 2 - λ 2,1 ν 1 -1. Note that from (6), ν 1 is negative. If ψ t/1 2 < ψ t/1 1
then λ 1,tλ 2,t < 0 hence user 1 will be scheduled at time t. We cannot say more using this method.

Theorem 1. If user 1 is supposed scheduled at time u = 1, then:

ψ t/1 2 < ψ t/1 1
=⇒ user 1 is scheduled at time t.

Denote by A u

1 and A u 2 the following sets:

A u 1 = {t = u, ψ t/u 1 > ψ t/u 2 } ∪ {u} A u 2 = {t = u, ψ t/u 1 < ψ t/u 2 }
For u = 1 and using Theorem 1 i.e. if t ∈ A u 1 then we have that user 1 is surely scheduled at time t > 1 with the assumption that the user is scheduled at time 1. We have the following result: Theorem 2. User 1 is scheduled at time u if and only if:

( A u 2 φ(S j 2 )) α > max k∈A u 1 φ(S k 2 ) φ(S k 1 ) ( A u 1 φ(S j 1 )) α (15)
Equivalent scheduling rule can be written for user 2. The proof of the theorem is given in [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF].

IV. HEURISTIC SOLUTION FOR FORECAST SCHEDULING

In this section, we propose a heuristic FS solution for an arbitrary number of users n which relies on FS with n = 2 users.

Consider an iterative algorithm for scheduling user(s) at time t with two main steps: (i) choose two users, i 1 and i 2 with the highest α-fair utility given by equations ( 16) and ( 17) where R j (t) is the mean past data rate (calculated over some time window) of the user j until time t; (ii) Select the best user (or both users) according to the FS rule for n = 2 for the time interval

[|1, T |].
It is recalled that FS for n = 2 can be computed using closed form formula or using a fast convex optimization solver. The algorithm is depicted in Table I. The user selection algorithm in Table I is not as optimal as the direct FS technique since it considers only the future SINR trajectories of the two best α-fair users. Interestingly, this approach provides significant gains with respect to Round Robin (RR) scheduler. It is recalled that at high mobility, the Proportional Fair (PF) scheduling converges to RR due to the short coherence time of the channel.

i 1 (t) = argmax j∈users φ(S t j ) R j (t) α (16) i 2 (t) = argmax j∈{users}-{i1(t)} φ(S t j ) R j (t) α , (17)
It is noted that alternatively and with a higher complexity, Step 2 in Table I can be performed using the FS for the two users selected in step 1 using a convex optimization solver. -Apply KKT resolution for the case of two users for i 1 and i 2 . Verify if time K exists (using (19) in [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF]). If it exists, use allocation (16-18) in [START_REF] Zaraoui | Analytical results for two users' forecast scheduling[END_REF]. Otherwise use [START_REF] Blondel | Recent advances in learning and control[END_REF], and select the user to schedule at time t accordingly (using the condition in Theorem 1) .

Step 4

-Update the mean past data rate for all users utilizing the received data rate at time t Step 5

-if t+1 ≤ T proceed to time t+1, T being the scheduling period Fig. 4. A road crossing a cell and the coverage area of a VSC.

V. NUMERICAL RESULTS

This Section presents numerical results for the FS heuristic obtained using a Long Term Evolution (LTE) network simulator written in Matlab. We consider a scenario with 10 vehicular users moving with a constant speed of 50km/h along a road crossing the cell. Close to the cell edge, the road traverses a coverage area of a Virtual Small Cell (VSC). It is recalled that the VSC is a remotely created small cell using a Large Scale Antenna System (LSAS) [START_REF] Tall | Virtual sectorization: design and self-optimization[END_REF] which is used to enhance the spatial SINR diversity along the trajectory (see Figure 4).

The spatial resolution of the REM is of 1m (which is feasible thanks to the interpolation of the GLM. The corresponding time interval at 50km/h is 70ms over which the SINR is considered constant. This time interval defines the time resolution of the FS during which a fixed allocation is applied, i.e. the same users are scheduled at time intervals defined by the technology (e.g. 1ms for LTE). The SINR experienced by the users along the road is presented in Figure 5. The first (stronger) and second maxima correspond to the transmissions at the vicinity of the VSC and the macro BS respectively. Simulation parameters are depicted in Table II.

In the following, we compare the TBUA solution with the (standard) FS solution of problem 3. The FS solution is solved using CVX (details below). The TBUA is solved using the heuristic solution described in Section IV with the analytical FS of Section III-C. We assume that the interfering neighboring cells have 10 percent of load (details on the interference model can be found in [START_REF] Zaaraoui | Forecast scheduling for mobile users[END_REF]). The interference conditions when using the REM will differ from those used to create the REM. We denote this as interference errors in the sequel. To show the robustness of the FS with respect to interference errors, we consider two cases: (i) a REM constructed with fully (100 percent) loaded interfering cells, which is denoted hereafter as "with interference errors" and (ii) a hypothetical REM that provides the exact SINR, namely 10 percent of loaded interfering cells, and is denoted as "without interfering errors".

Reference [START_REF] Zaaraoui | Forecast scheduling for mobile users[END_REF] explains that the scheduling decisions of the FS can be impacted to a certain extent by the interference errors, and the impact of these errors decreases with the increase in the scheduling periodicity. Moreover, the impact of the interference errors on the Mean User Throughput (MUT) is negligible. This behavior is explained by the high spatial dynamicity of SINR variations, the trends of which are well captured by the REM in spite of the interference errors. We consider the RR as a baseline for the different approaches. The rationale is that at high mobility considered here, the coherence time is too short to exploit fast fading in order to achieve opportunistic scheduling gain.

The objective functions of the optimization problems (2) and ( 3) are convex and can be solved using a convex optimization solver. The CVX Matlab library has been used (see [START_REF] Research | CVX: Matlab software for disciplined convex programming, version 2.0[END_REF] and [START_REF] Blondel | Recent advances in learning and control[END_REF]). The CVX solver verifies the convexity of the problem and solves it using SDPT3 or SeDuMi. SDPT3 implements the infeasible pathfollowing algorithms for solving conic programming problems whose constraint cone is a product of semidefinite cones. It uses a predictor-corrector primal-dual path-following method, with different types of search direction. SeDuMi is a linear/quadratic/semi-definite solver.

The MUT time evolution for the 10 users and for the three different FS solutions with 2s periodicity: the (standard) FS with and without interference errors, the TBUA heuristic with interference errors, and for the baseline RR are presented in Figure 6. One can observe oscillations that are correlated with the SINR variations (see Figure 5). Most of the time, the TBUA varies between the FSs and the RR scheduling. The curves for the FSs and TBUA solutions are most of the time above those of the RR scheduling. As stated above, the impact of interference errors is negligible, which is an important characteristic of the FS based solutions. Figure 7 compares the MUT gain with respect to the baseline RR for the ten users, averaged along the full trajectory crossing the cell, for different FS scheduling periodicities. For 2s periodicity, the MUT gain for the TBUA and FS is of 15 and 30 percent respectively, and it grows monotonically to 80 and 117 percent respectively for 20s periodicity. One can observe a close to linear monotonic growth of MUT gain till 15s scheduling priodicity, with non-significant improvement brought about the long term time-space diversity beyond this periodicity.

VI. CONCLUSION

The allocation of FS involves the solution of a convex optimization problem. The solution complexity scales up with the number of users and the scheduling duration. In view of reducing the computational complexity of the FS we have proposed a two step heuristic in which, for each scheduling period, we select the two users with the highest α-fair utility and then we apply to them the FS. It has been shown that for n = 2, the FS allocation can be written in a closed form, making the heuristic even more computational efficient. Numerical simulations show that the throughput gain achieved by the FS heuristic remains significant, even for short scheduling periodicity.

Fig. 1 .

 1 Fig. 1. Coverage map prediction using FRK.

Fig. 2 .

 2 Fig. 2. From SINR map to reachable bitrates.

Fig. 3 .

 3 Fig. 3. Scheduling using the best two α-fair heuristic.

Fig. 5 .

 5 Fig. 5. SINR for different users along the road.

  Forecast scheduling without interference error MUT in time Forecast scheduling with interference error MUT in time 2 best users algorithm MUT in time RR MUT in time

Fig. 6 .

 6 Fig.6. Time variation of users' throughput for a scheduling period of T=2s, averaged over a sliding window of 15 iterations.

Fig. 7 .

 7 Fig. 7. FS MUT gain with and without interference errors, and TBUA wrt the RR scheduler as a function of the scheduling period (in sec).

TABLE I TWO

 I BEST USERS ALGORITHM (TBUA) Algorithm of approximation with two best α-fair users Step 1 -Select the two best users i 1 and i 2 according to (16) and (17) at time t Step 2


Jérémy Vezinet

Olivier Julien Signal

Olivier Julien

Multi-Antenna GNSS Receiver for Space Launcher

Since 2011, he is a researcher and worked on multi-sensor hybridization research activities

INTRODUCTION

GNSS is today widely used in many navigation application: urban and terrestrial (for pedestrians, cyclists, automobiles or robots), maritime (for boats and ships), air navigation (for military and civil aircraft or UAVs) and spatial (for space or launch vehicles). Each of these domains can require a specific signal processing and navigation algorithm design based on the environment characteristics, the dynamic of the carrier or the requirement of the intended application.

For the particular case of space launchers, a significant variation in the attitude of the platform involves the use of multiple antennas located around the body in order to track a sufficient number of satellites during the flight. Several multi-antenna techniques have been addressed in the literature. Antenna diversity techniques are presented in [START_REF] Kim | Maximum Ratio Diversity Combining Receiver Using Single Radio Frequency Chain and Single Matched Filter[END_REF] and [START_REF] Akbari | Receiver Diversity Combining Using Evolutionary Algorithms in Rayleigh Fading Channel[END_REF] where a GNSS receiver uses combinations of signals from several antennas. Antenna diversity techniques offer a great improvement for radio link performance and uses redundancy of the received signal to increase the signal to noise ratio. Those techniques are usually employed to fight against signal fading. [START_REF] Grillenberger | An Integrated Diversity Switching GPS Receiver for Sounding Rocket Missions[END_REF] and [START_REF] Ebinuma | GPS Signal Tracking on Spinning Vehicles with Antenna Diversity Techniques[END_REF] also present antenna diversity techniques sounding rockets applications, vehicles that have high spinning movements, to keep tracking a sufficient number of satellite signals during the spinning motion. [START_REF] Puga | Multi-antenna GNSS Receiver Tracking Algorithm for Vehicles With Unconstrained Three-dimensional Motion[END_REF] presents an interesting approach by detailing a receiver that combines inputs from multiple antennas during signal tracking in order to keep a full sky visibility. The algorithm describes a fast-switching antenna selection process in order to ensure synchronization between the receiver and the incoming signal at every moment.

In addition to the desired multi-antenna architecture, the high dynamic motions of a launcher is a serious constraint for the design of a GNSS receiver. Indeed, acquisition and tracking of GNSS signals in presence of a non-null jerk is problematic since classic tracking loops have a limited capacity to follow frequency and phase variations induced by high dynamic motions. A launch vehicle is constantly submitted to high dynamic phases such as rotations or sudden accelerations (engines turn-on and turn-off) causing systematic losses of signal tracking.

One of the possible improvements for tracking loops used on GNSS receivers used in high dynamic vehicles can be the increase of the tracking loops' order. For example, [START_REF] Stephens | Controlled-Root Formulation for Digital Phase-Locked Loops[END_REF] presents a controlled-root formulation for digital PLL with coefficients allowing to model a fourth order PLL. Another improvement for tracking loop designs under high dynamic constraints consists in using aided PLL as depicted in [START_REF] Roncagliolo | GPS Digital Tracking Loops Design for High Dynamic Launching Vehicles[END_REF] and [START_REF] Roncagliolo | Optimized Carrier Tracking Loop Design for Real-Time High-Dynamics GNSS Receivers[END_REF]. In that structure, the FLL-aided PLL design allows tracking accelerations steps happening during some phases of flight of a launch vehicle.

The solution proposed in the paper depicts a multi-antenna GNSS solution for a launch vehicle. The described solution offers a multi-antenna, dual-frequency and multi-constellation capability and assess positioning performances for a high dynamic range vehicle. The proposed algorithm is an extended Kalman filter using code pseudorange and Doppler measurements for estimation of the user's position and velocity. The paper is organized as follows:

 The first section describes the antenna and receiver configuration considered during our simulations. Number and positions of antennas are presented just as the considered tracking parameters of the receiver.  The second section provides the description of the considered navigation filter.  The third section presents the antenna switching strategy employed when some signals goes from one antenna to another. In such a case, the tracking channel parameters of one antenna can be used for aiding the acquisition of the same signal through another antenna.  The fourth section details the algorithm performance assessment. Simulations in several configuration will be conducted: o A first basic configuration aims at using all the antenna measurements in the EKF algorithms. Therefore, in this scenario no lever arm will be considered. o The second configuration focuses on providing a navigation solution for each antenna. In this case the low number of visible satellites per antenna is an issue. o The last configuration is the more complex one, it involves the antenna switching process integrates all measurement (one per satellite) in one extended Kalman filter and includes lever-arm compensation.  The main conclusions of this paper, trade-off analysis and future work will be drawn in the last section.

ANTENNA AND RECEIVER CONFIGURATION

The system considered in this paper is composed of a GNSS receiver with three input antennas, mounted on a launch vehicle. Antennas are considered to be uniformly spaced around the cylindrical body of the launcher as illustrated in Figure 1. Three antennas are chosen in this configuration since their combined coverage should ensure that during every part of flight, all satellites in view are received whatever the attitude of the vehicle.

Figure 1. Spatial antenna distribution

The simulation presented in that paper use the GPS and Galileo constellations in the L1 and L5 frequency bands. Processing of the signals received by all antennas were based on the ENAC GNSS receiver simulator "GeneIQ". This receiver simulator works at the correlator output level. It thus employs very fine correlator output models including a large source of errors. For the purpose of this project, it has been adapted to be able to process signals from multiple antennas simultaneously and to take into account the attitude evolution of the platform. The position computation module has also been improved to take the multiple resulting measurements in consideration.

The navigation filter algorithm was executed at 100 Hz. The GNSS receiver simulator allows testing several tracking architectures but for the presented simulations the configuration is as illustrated in the Table 1. The RF filter bandwidth is assumed to be of 24 MHz (double-sided). The choice of the tracking loop parameters (PLL order, bandwidth, and discriminator) has been set to these values for the current paper based on typical parameter values in the state-of-the-art. However, a study for optimizing the choice of these parameters is planned.

Note that the system configuration is based on the use of 3 receivers, each associated to one antenna. Moreover, it is assumed that these 3 receivers share a common oscillator in order to have a unique system time.

NAVIGATION FILTER DESCRIPTION

The proposed navigation filter is an extended Kalman filter (EKF). The filter uses code pseudoranges and Doppler measurements so as to estimate the vehicle's position and velocity. By definition, a Kalman filter is a recursive estimator of the internal state of a linear stochastic system. The EKF formulation allows dealing with nonlinear systems such as navigation systems.

EKF State

In the proposed EKF, the state is described by the following 8 state vector:

𝑋 𝑘 = [𝑥 𝑥̇𝑦 𝑦̇𝑧 𝑧̇𝑏 𝐶 𝑑 𝐶 ] (1)
Where:  (𝑥, 𝑦, 𝑧) is the receiver's position vector in the ECEF system coordinates  (𝑥, 𝑦, 𝑧) is the receiver's velocity vector in the ECEF system coordinates  𝑏 𝐶 and 𝑑 𝐶 are respectively the clock bias and clock drift of the receiver Using the EKF, the estimation of the state vector 𝑋 𝑘 is divided in two steps:

 The time prediction: This step consists in providing an estimate of the state vector 𝑋 𝑘 (and of its covariance) based on its previous estimation 𝑋 𝑘-1 and on the state transition model. The estimate of the state vector through the time prediction step is called the a-priori estimate and denoted 𝑋 𝑘 𝑘 ⁄ -1 .

 The measurement update: This step consists in updating (or correcting) the a priori estimate 𝑋 𝑘 𝑘 ⁄ -1 (and its covariance) by using the measurement vector 𝑧 𝑘 . The estimate of the state vector through the measurement update step is called the a-posteriori estimate and denoted 𝑋 𝑘 𝑘 ⁄ .

State Transition Model

The state transition model of the EKF filter in the continuous time domain can be expressed as follows:

𝑋 ̇(𝑡) = 𝜑 ⋅ 𝑋(𝑡) + 𝑤(𝑡) (2)
Where:  𝑋 ̇ denotes the derivative of the state vector 𝑋  𝑤 is the additive process noise vector assumed centered, white and Gaussian  𝜑 is the continuous state transition matrix (Jacobian of the 𝑓 state transition function)

Then, in the discrete time domain, the transition state model becomes:

𝑋 𝑘 = 𝛷 ⋅ 𝑋 𝑘-1 + 𝑤 𝑘 (3)
Where:  𝑋 𝑘 denotes the state vector forward projection from the (𝑘 -1) 𝑡ℎ to the 𝑘 𝑡ℎ time epoch  𝛷 is the discrete state transition matrix and can be expressed as follows [START_REF] Shytermeja | Performance Evaluation of VDFLL Architecture for a Dual Constellation L1/E1 GNSS Receiver in Challenging Environments[END_REF]:

𝛷 = [ 𝐹 0 2×2 0 2×2 0 2×2 0 2×2 𝐹 0 2×2 0 2×2 0 2×2 0 2×2 𝐹 0 2×2 0 2×2 0 2×2 0 2×2 𝐹 ], 𝐹 = [ 1 Δ𝑇 0 1 ] (4)
 Δ𝑇 denotes the time interval between two consecutive estimations, representing the measurement update time of the central filter  𝑤 𝑘 is the discrete process noise vector. It is modeled as a Gaussian noise vector with zero mean and discrete covariance matrix 𝑄 𝑘

The process noise comes from two sources namely [START_REF] Shytermeja | Performance Evaluation of VDFLL Architecture for a Dual Constellation L1/E1 GNSS Receiver in Challenging Environments[END_REF], the user dynamic noise (constituted by the user's position and velocity terms) and the receiver's clock noise (local oscillator NCO noise).

The process noise is the first critical parameter affecting the performance of the filter. Therefore, an accurate tuning is required to fasten the EKF estimation convergence toward the true user state. Regarding the process noise covariance matrix, its discretetime form 𝑄 𝑘 is obtained from the continuous-time process noise matrix 𝑄. The 𝑄 matrix describes the uncertainty of the user's dynamic and is defined by:

𝑄 = 𝐸[𝑤 ⋅ 𝑤 𝑇 ] (5)
As depicted in [START_REF] Shytermeja | Performance Evaluation of VDFLL Architecture for a Dual Constellation L1/E1 GNSS Receiver in Challenging Environments[END_REF], 𝑄 covariance matrix contains:  The user's dynamic sensitive terms (𝜎 𝑥̇, 𝜎 𝑦̇, 𝜎 𝑧̇) , mainly driven by the modeled user dynamic  The receiver's oscillator noise terms including the oscillator's phase noise PSD, 𝜎 𝑏 , and the oscillator's frequency noise PSD, 𝜎 𝑑 , which by themselves depend on the Allan variance parameters ℎ0 and ℎ-2.

And can be expressed as follows:

𝑄 = [ 𝑄 𝑥 0 0 0 0 𝑄 𝑦 0 0 0 0 𝑄 𝑧 0 0 0 0 𝑄 𝑐𝑙𝑘 ] (6)
With:

 𝑄 𝑥 = [ 0 0 0 𝜎 𝑥2 ] (the same for 𝑄 𝑦 and 𝑄 𝑧 )

 𝑄 𝑐𝑙𝑘 = [ 𝜎 𝑏 2 0 0 𝜎 𝑑 2 ]
The discrete-time form of the process noise covariance sub-matrices 𝑄 𝑥 , 𝑄 𝑦 and 𝑄 𝑧 can then be expressed:

𝑄 𝑥,𝑘 = ∫ 𝛷 𝑥 ⋅ 𝑄 𝑥 ⋅ 𝛷 𝑥 𝑇 𝑡 𝑘-1 +Δ𝑇 𝑡 𝑘-1 = ∫ [ 1 Δ𝑇 0 1 ] ⋅ [ 0 0 0 𝜎 𝑥2 ] ⋅ [ 1 0 Δ𝑇 1 ] 𝑡 𝑘-1 +Δ𝑇 𝑡 𝑘-1 = [ Δ𝑇 3 3 ⁄ Δ𝑇 2 2 ⁄ Δ𝑇 2 2 ⁄ Δ𝑇 ] (7)
In an equivalent way, the receiver's clock process noise covariance discrete-time form can be expressed as follows:

𝑄 𝑐𝑙𝑘,𝑘 = [ 𝜎 𝑏 2 ⋅ Δ𝑇 + 𝜎 𝑑 2 ⋅ Δ𝑇 3 3 ⁄ 𝜎 𝑑 2 ⋅ Δ𝑇 2 2 ⁄ 𝜎 𝑑 2 ⋅ Δ𝑇 2 2 ⁄ 𝜎 𝑑 2 ⋅ Δ𝑇 ] (8)

Measurement Model

The measurement model describes the relation between the state of the system and the measurement vector. It is expressed as follows [START_REF] Shytermeja | Performance Evaluation of VDFLL Architecture for a Dual Constellation L1/E1 GNSS Receiver in Challenging Environments[END_REF]:

𝑧 𝑘 = ℎ(𝑋 𝑘 ) + 𝑣 𝑘 (9)
Where:  ℎ is the non-linear observation function that relates the measurement 𝑧 𝑘 to the state vector 𝑋 𝑘  𝑣 𝑘 is the measurement noise vector and is modelled as a zero-mean uncorrelated Gaussian noise process and independent to the process noise 𝑤𝑘  𝑧 𝑘 is the measurement vector and is composed of the code pseudoranges and the Doppler measurements Let 𝜌 𝑖 and 𝜌̇𝑖 respectively be the code pseudorange and Doppler measurements for a given satellite 𝑖. At each epoch 𝑘, these measurements are modeled as:

𝜌 𝑖 = ℎ 𝑝 (𝑋 𝑘 ) + 𝑛 𝑝,𝑘 (10)
Where:

ℎ 𝑝 (𝑋 𝑘 ) = 𝑟 𝑘 𝑖 + 𝑐 ⋅ 𝑏 𝐶 (11)
 𝑟 𝑘 𝑖 = √ (𝑥 𝑘 𝑖 -𝑥 𝑘 ) 2 + (𝑦 𝑘 𝑖 -𝑦 𝑘 ) 2 + (𝑧 𝑘 𝑖 -𝑧 𝑘 )
2 being the geometrical distance from the 𝑖 𝑡ℎ satellite to the user at epoch 𝑘  (𝑥 𝑖 , 𝑦 𝑖 , 𝑧 𝑖 ) being the coordinates of the 𝑖 𝑡ℎ satellite position in the ECEF frame  𝑛 𝑝 is the zero-mean Gaussian-distributed noise affecting the code pseudorange measurement And

𝜌̇𝑖 = ℎ 𝑑 (𝑋 𝑘 ) + 𝑛 𝑑,𝑘 (12)
Where: The measurement noise vector 𝑣 𝑘 is modelled as a zero-mean uncorrelated Gaussian noise process and is independent to the process noise 𝑤 𝑘 . Its covariance matrix 𝑅 𝑘 is defined by:

ℎ 𝑑 (𝑋 𝑘 ) = (
𝑅 𝑘 = 𝐸[𝑣 𝑘 ⋅ 𝑣 𝑘 𝑇 ] (14)

ANTENNA SWITCHING PROCESS

In order to deal with the high dynamic motion and the changes in attitude of a launch vehicle and ensure a sufficient number of received satellite signals, the multi-antenna structure has been considered in that study. Several solution adopting antenna diversity techniques have been presented in the literature ( [START_REF] Grillenberger | An Integrated Diversity Switching GPS Receiver for Sounding Rocket Missions[END_REF] and [START_REF] Ebinuma | GPS Signal Tracking on Spinning Vehicles with Antenna Diversity Techniques[END_REF]). In our proposed solution, we choose not to combine signals from different antennas but to help tracking the signal of a new satellite seen from a given antenna with the information obtained when tracking the same signal from another antenna. That strategy ensures a reduction of the acquisition time in case of fast spinning motion when losing and re-acquiring a signal periodically. Figure 1 shows the number of satellites seen by each antenna during the flight.

In a similar way, Figure 3, illustrates the number of satellite signals that are considered as locked and that can be used in the navigation filter. The event at 𝑡 = 150 𝑠 corresponds to a dropping phase during which the simulated trajectory is showing some singularities on the acceleration profile. The main consequence is the loss of lock of all signals. As the acquisition process was not studied in that paper, a perfect reacquisition is assumed a few seconds after that loss of lock.

Due to the proximity of the different antennas, it is very easy to initialize the DLL associated to a satellite newly seen by an antenna with the tracking parameters of the antenna already seeing that same antenna. However, this is not the case for the PLL. Indeed, the signals transmitted in the L1 and L5 bands by both GPS and Galileo constellations have wavelengths of around 19cm and 25cm. These values are smaller than the distance between antennas so the carrier phase of arriving signal to different antennas may be different. Moreover, the rotation motion of the vehicle also impacts in different ways the Doppler associated to each antenna for a given satellite. The Doppler difference between antennas will depend upon the rotation motion of the vehicle, the relative position of each antenna and the location of the satellite. As a consequence, the only way to take properly into account this phase/Doppler difference between antennas would be to have access to an accurate estimation of the vehicle attitude and motion. Otherwise, the rotation of the vehicle would have to be slow. In this paper, it is reminded that each antenna is linked to a specific receiver (each receiver sharing the same oscillator). The proposed solution is thus as follows: when antenna 1 is receiving a new signal, which is currently being tracked by antenna 2, the tracking parameters associated to antenna 2 are used for initializing the antenna 1 tracking parameter. More precisely, phase and frequency estimated on antenna 2 are transmitted for the initialization of the carrier tracking loop of the antenna 1. That process can be optimized in the future by simulating fast rotation motion and assessing the switching performances.

Figure 4 and Figure 5 show estimated C/N0 values for 2 satellites (PRN 8 and PRN37). For PRN 8, we see that antenna 3 is able to track the satellite until 𝑡 = 90𝑠. After that instant, antenna 1 is receiving the signal and antenna 2 is no longer tracking it. For PRN 37, there is a time interval during which antenna 1 and 3 are both tracking the satellite. During that phase, both antennas are able to track the signal and provide a measurement to the navigation filter. PLL associated to antenna 3 can provide to the PLL associated to antenna 1 the estimation of the frequency in order to start tracking the signal as soon as the satellite is in view.

PERFORMANCE ASSESSMENT

In order to assess navigation performance of this solution and to propose improvements regarding the antenna, receiver or tracking loops configuration, several sets of simulations were performed. These simulations are assuming the presence of thermal noise only during the trajectory. As mentioned earlier, it is important to understand that the simulator used works at the correlator output level and that for the purpose of this project, it has been adapted to be able to process signals from multiple antennas simultaneously and to take into account the attitude evolution of the platform. The pseudorange measurements used by the positioning filter are Iono-Free measurements based on L1/L5 measurements. It is also assumed that if the receiver loses the L1 or the L5 signal, the corresponding IF measurement is unavailable (there is no reversion mode towards a single frequency measurement).

A first batch of independent Kalman filters, in which each of them is allocated to a specific antenna is used to assess performances of individual antennas solutions. These solutions deal with estimating position and velocity of each antennas without using the multi-antenna structure. However, even in a dual constellation configuration, some antennas can be masked by the vehicle body during some parts of the flight (see Figure 1) and the navigation solution is then degraded. Then the main drawback of these independent antenna solutions is the low number of measurement that can turn in a diverging navigation solution as soon as the antenna face a reduced part of the sky. Figure 6 and Figure 7 illustrate the KF solution estimation error per antenna. The dot lines correspond to the 3𝜎 estimated covariance by the Kalman filters. As soon as the number of locked satellites drops (for antennas 2 in green and antenna 3 in red), the estimated covariance as well as the estimation error increase and does not provide a reliable solution. In the simulated configuration, it reflects the fact that antenna 1 is facing up during most of the flight while antennas 2 and 3 face down and thus have poor satellite visibility. As a comparison, Figure 8 and Figure 9 show a least square estimation of the position and velocity per antenna. As mentioned before, the main drawback of these individual solutions is the decrease of the number of visible satellites per antenna during some parts of the flight. A first simple multi-antenna solution that do not solve this visibility issue but that may improve the navigation performance is to average the three individual antenna solutions. Even if the averaging will be impacted by the degraded solutions per antenna in some parts of the flight (especially when the number of locked satellite is low), it can provide a good first combined navigation solution. The position and velocity estimation errors appears in Figure 10 and Figure 11 in pink.

The next multi-antenna solution tested for this study is a solution that consists in using all measurement gathered from the three antennas for computing a global navigation solution. That global navigation focuses in estimating the position and velocity of the gravity center of the launch vehicle. All measurement associated to the three antennas are combined in a main EKF and in a first configuration, there is no compensation of the lever arm between the antenna and the vehicle's center of gravity. The lever arm is not compensated as it is assumed that there is no access with sufficient accuracy to the attitude of the vehicle. The position and velocity estimation errors are showed in Figure 10 and Figure 11 in black. The least square position and velocity estimation error are illustrated in Figure 12 and Figure 13.

The solution in which the lever arm are compensated is also simulated. In this case, the lever arm expressed in the vehicle frame is assumed to be well known. As presented in equation ( 1), position and velocity in the Kalman filter are expressed in the ECEF frame. The lever arm is then expressed in the ECEF frame so as to compute the innovation vector with a lever arm compensated predicted position. That frame rotation is done using an estimation of the vehicle's attitude. The estimation is assumed to be provided by an inertial aid during that simulation. The improvement in terms of performance of the estimated position and velocity is illustrated in Figure 10 and Figure 11 in blue. It can be used so as to justify the need of an inertial aid in the final solution. The barycenter solution (in pink), as expected, drifts when one of the individual antenna solution drifts. This is due to the low number of satellite seen by antennas 2 and 3 after 100 seconds. That solution could be improved by detecting the reliable "per antenna" solution and discard the ones that are computed with a too low number of measurements.

The black plot shows the estimation errors for the solution that gather all measurements from all antennas but without lever arm compensation. A bias in the position estimation remains and could be an issue for large radius vehicles. In the current application a 2m radius has been assumed.

The blue plot, which is identical to the black one but that uses lever arm compensation shows that there is no longer bias on the position estimation. In the two last configurations, the Kalman filter can integrate a sufficient number of measurement in order to provide a reliable position and velocity estimation. For the solution with lever arm compensation, the position on each axis reaches a sub-meter accuracy.

CONCLUSION AND FUTURE WORK

The objective of this study is to propose a multi-antenna, multi-frequency and multi-constellation GNSS receivers. The solution is intended to be used for space launchers application and to deal with multiple constraints like the high dynamic of the vehicle or the masking of antennas by the vehicle. A position and velocity estimation algorithm using an extended Kalman Filter has been proposed integrating all measurement from the three antennas. A first set of results are provided based on simulations assessing the position and velocity estimation performance. An antenna switching strategy allowing to transmit estimated phase and frequency from an antenna to another has been used and will be improved in the future.

Those results offer a first analysis for refining the tracking loops architecture and parameters, the multi-antenna processing strategy and the navigation filter. Indeed, tracking parameters can be optimized so as to cope with the vehicle dynamics. Even if the increase of the carrier tracking loops bandwidth could have a significant effect on the noise at the measurement level and then on the navigation performances, it can help keeping lock of signals during high acceleration phases. Regarding the multiantenna strategy, even if a more complex process can help facing high rotation rates, the proposed way to initialize the PLL of an antenna with the estimated frequency from another antenna seems to be sufficient for the tested scenario. Finally, the navigation filter used in that solution (an 8-state EKF), could be improved by coupling the GNSS measurements to an IMU so as to estimate the fast changing attitude of the vehicle.

Figure 2 .Figure 3 .

 23 Figure 2. of visible satellites per antenna Figure 3. Number of locked satellites per antenna

Figure 4 .Figure 5 .

 45 Figure 4. Estimated C/N0 per antenna for PRN 8 Figure 5. Estimated C/N0 per antenna for PRN 37

Figure 6 .Figure 8 . 9 .

 689 Figure 6. KF Position estimation error per antenna Figure 7. KF Velocity estimation error per antenna

Figure 10 .Figure 12 .

 1012 Figure 10. Global KF Position estimation error Figure 11. Global KF Velocity estimation error

Table 1 . Code and Carrier tracking parameters L1/L5/E1/E5 Code Tracking Parameters

 1

	DLL Order

  𝑦̇𝑖, 𝑧̇𝑖) being the velocity of the 𝑖 𝑡ℎ satellite in the ECEF frame  𝑛 𝑑 is the zero-mean Gaussian-distributed noise affecting the Doppler measurement

	 (𝑎 𝑥 𝑖 , 𝑎 𝑦 𝑖 , 𝑎 𝑧 𝑖 ) denotes the line-of-sight unit vectors from the receiver to the 𝑖 𝑡ℎ satellite; 𝑎 𝑥 𝑖 = 𝑧 𝑘 𝑖 -𝑧 𝑘	𝑥 𝑘 𝑖 -𝑥 𝑘 𝑟 𝑘 𝑖 , 𝑎 𝑦 𝑖 =	𝑦 𝑘 𝑖 -𝑦 𝑘 𝑖 , 𝑎 𝑧 𝑖 = 𝑟 𝑘
	𝑟 𝑘 𝑖
	 (𝑥̇𝑖,
	𝑥̇𝑘 𝑖 -𝑥̇𝑘) ⋅ 𝑎 𝑥 𝑖 + (𝑦̇𝑘 𝑖 -𝑦̇𝑘) ⋅ 𝑎 𝑦 𝑖 + (𝑧̇𝑘 𝑖 -𝑧̇𝑘) ⋅ 𝑎 𝑧 𝑖 + 𝑐 ⋅ 𝑑 𝐶		(13)

ACKNOWLEDGMENTS

This works has been done in the framework of a French national funding called "FUI REGLO". We would like to thank the leader of this project SYNTONY GNSS for its support.


Simon Jones
email: jones5@rose-hulman.edu

Tonghua Zhang

Mathias Legrand

The computation of wavelet-Galerkin three-term connection coefficients on a bounded domain

Keywords:

Computation of triple product integrals involving Daubechies scaling functions may be necessary when using the wavelet-Galerkin method to solve differential equations involving nonlinearities or parameters with field variable dependence. Numerical algorithms for determining these triple product integrals, known as three-term connection coefficients, exist but tend to suffer from ill-conditioning. A more stable numerical solution algorithm is presented herein and shown to be both accurate and robust.

Introduction

The use of Daubechies wavelet families as Galerkin basis functions for solving differential equations is of growing interest [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Strang | Wavelets and Filter Banks[END_REF]. These oscillatory functions have compact support which allow sparse representation of complex responses on unbounded, bounded or periodic domains [START_REF] Amaratunga | Wavelet-Galerkin solution of boundary value problems[END_REF][START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Daubechies | Ten lecturers on wavelets[END_REF][START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF][START_REF] Pernot | A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis[END_REF][START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF][START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF]. For the discrete orthogonal wavelet-Galerkin method, Daubechies scaling functions are commonly used as the functional basis [START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF]. The Galerkin formulation for equations containing nonlinearities requiring the product of the field variable with itself or its derivative require the integration of a scaling function triple product. The Daubechies scaling functions cannot be defined explicitly making analytic integration intractable, and their fractal nature (i.e. discontinuities which are independent of scale) make numerical integration error prone [START_REF] Welstead | Fractal and wavelet image compression techniques[END_REF].

Innovative work by Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Zhang | Comments on 'The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], Latto et al. [START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF], and Romine et al. [START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF] provide algorithms to compute the exact solution to the three-term connection coefficients on a bounded domain. In each of these references the authors solve for the connection coefficients using a set of rank deficient scaling equations defined recursively using the two-scale definition of the scaling functions. The rank deficiency is filled by replacing a corresponding number of equations with theoretically independent moment equations, allowing determination of a unique set of connection coefficients. Due to the numerical error introduced during implementation of these algorithm, the scaling equations and moment equations are generally no longer independent. This has been found to lead to illconditioning of the system and calculation of erroneous connection coefficients. A novel algorithm is presented herein which can account for this numerical error by solving for a set of connection coefficients which satisfy all the constraining equations in a least-squares sense.

In Section 2 a brief review of Daubechies wavelet notation is included and key references which contain derivations of some necessary parameters are cited. Section 3 details the proposed method of computing the three-term connection coefficients, included an example calculation. The results are compared with existing coefficients found in the literature to validate the method. Conclusions are presented in Section 4.

Daubechies Wavelet Notation

The Daubechies scaling function is defined by a set of L filter coefficients p `W `2 OE0; L 1, where L is an even integer. The fundamental two-scale equation is defined as

.x/ D L 1 X `D0 p ` .2x `/ (1)
where .x/ is the scaling function with fundamental support over the finite intervals OE0; L 1. The filter coefficients p `are derived by imposing a number of constraints given by Daubechies [START_REF] Daubechies | Ten lecturers on wavelets[END_REF].

It is useful to define the n th derivative of the scaling function as

By amalgamating Eqs. ( 1) and ( 2) and accounting for the conditions in Ref [START_REF] Daubechies | Ten lecturers on wavelets[END_REF], it is possible to write [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] .

n/ .x/ D 2 n L 1 X `D0
p `.n/ .2x `/; n D 0; 1; : : : ; L=2 1:

It is also useful to define the inner product of the scaling function and its derivative over a bounded interval

n k .x/ D Z x 0 .y/ .n/ .y k/dy: (4)
The solutions n k .x/ are known as the two-term connection coefficients [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]; these coefficients are required in the next section when determining three-term connection coefficients. One algorithm for computing these two-term connection coefficients was derived by Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], with corrections presented by Zhang et al. [START_REF] Zhang | Comments on 'The computation of wavelet-Galerkin approximation on a bounded interval[END_REF].

Evaluation of the three-term connection coefficients

The three-term connection coefficients over a bounded domain are defined as follows [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] m;n j;k .x/ D Z x 0 .y/ .m/ .y j / .n/ .y k/dy [START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF] for 0 m; n .L=2 1/ and j; k; m; n; x 2 Z, with the following properties m;n j;k .x/ D 0 for jj j; jkj; or jj kj L 1

m;n j;k .x/ D 0 for x j; x k; or x 0 (7)

m;n j;k .x/ D m;n j;k .L 1/ for x j; x k; or x L 1: (8)
Substituting the two-scale relations ( 1) and (3) into Eq. ( 5) and performing a change of variable gives

m;n j;k .x/ D 2 mCn 1 L 1 X i a D0 L 1 X i b D0 L 1 X i c D0 p i a p i b p i c m;n 2j Ci b i a ;2kCi c i a .2x i a /: (9)
Accounting for the constraints provided by Eqs. ( 6) to ( 8), the scaling equations ( 9) can be written in matrix form as

2 .1 m n/ Q m;n .x/ D S Q m;n .x/ (10)
for x D 1; 2; : : : ; L 1, where S has entries compiled from summing the relevant triple products p i a p i b p i c as defined in Eq. ( 9). This implies the connection coefficient vector Q m;n .x/ belongs to the eigenspace corresponding to the eigenvalue 2 .1 m n/ from Eq. [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF]

q D .m C n C 1/ C 8 < : mCn X i D1 i if m C n L=2 mCn X i DL=2C1 3L 2 2i C .L C 2/L 8 if L=2 < m C n L 2: (14)
The eigenvectors corresponding to these q eigenvalues describe the solution space of the scaling equations [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF] for a given m and n. In fact, since the scaling equations depend only on the summation .m C n/ and not the specific derivatives, this set of eigenvectors gives the scaling equation solution space for all three-term connection coefficients whose derivatives sum to .m C n/ [START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF].

The unique solution for derivatives m and n is found by considering the set of moment equations which are derived in Ref [START_REF] Daubechies | Ten lecturers on wavelets[END_REF]; the derivation can be found in Refs. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Zhang | Comments on 'The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]:

X k k n m;n j;k .x/ D nŠ m j .x/ (15) X j j m m;n j;k .x/ D mŠ n k .x/: (16)
Thus Q m;n .x/ is uniquely described by the intersection of the scaling equation solution space with that of the moment equations. This implies the solution must be a linear combination of the q eigenvectors; the participation factors can be computed from the moment equations as detailed below.

Example calculation

Consider the specific case of m D 0; n D 1. Eq. ( 14) states q D 3 eigenvectors describe the solution space to the scaling equations, thus

Q m;n .x/ D c 1 u C c 2 v C c 3 w (17)
where u, v and w are the eigenvectors corresponding to the eigenvalue 2 .1 m n/ D 1 from Eq. ( 10); the constants c 1 , c 2 and c 3 are the respective participation factors to be determined. Substituting Eq. ( 17) into the moment equations results in 19) due to the constraints of Eqs. ( 6) to [START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF].

c 1 X k k n u C c 2 X k k n v C c 3 X k k n w D nŠ m j .x/ 8j 2 OE2 L; L 2 (18) c 1 X j j m u C c 2 X j j m v C c 3 X j j m w D mŠ n k .x/ 8k 2 OE2 L; L 2 (19)
These moment equations can be written in matrix form as

M 8 < : c 1 c 2 c 3 9 = ; D b ( 20
)
where M is a rectangular matrix of size 2g 3, comprised of the summation terms on the lefthand side of Eqs. ( 18) and ( 19), and b is a vector of length 2g composed of the respective righthand side terms. The participation factors are determined using a pseudoinverse

c D .M T M/ 1 M T b: (21)
This gives a robust "best fit" in a least-squares sense, which allows for any dependencies between the scaling and moment equations resulting from accumulated numerical error. The unique solution of the three-term connection coefficients can thus be found by substituting the participation factors found in Eq. ( 21) into Eq. ( 17). The algorithm is analogous for different values of m and n.

The accuracy of the computed Q m;n .x/ vector is quantified by first calculating the residual of each moment equation given in Eqs. ( 15) and ( 16); as Q m;n .x/ is a linear combination of the eigenvectors from Eq. ( 10), the scaling equations are automatically satisfied. To allow meaningful comparison of different derivative combinations, the residuals are normalized by the L2-norm of the righthand-side of the moment equation. The RMS value of this normalized residual vector is then calculated to give a scalar measure of the absolute error for a given n and m combination. This error measure for Daubechies scaling functions (L D 8) at all allowable derivative combinations is provided in Fig. 1. As shown, the error in Q m;n .x/ grows with higher derivatives, but even at n; m D 3 the error norm remains relatively low. As a comparison, the three-term connection coefficients for L D 6; m D 0; n D 1 published in Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] result in an RMS error norm of 1.55e-10, whereas for the current algorithm it is 1.82e-16. The three-term connection coefficients are tabulated in the Appendix for comparison.

Conclusions

Computation of the three-term connection coefficients is necessary when using the wavelet-Galerkin method to solve differential equations involving nonlinearities or parameters with variable dependence. Algorithms currently exist to solve for these coefficients but they have been found to suffer from ill-conditioning which can result in erroneous results. The current investigation introduces a novel solution algorithm which appears to be more numerically robust and comparatively more accurate than previously published algorithms.

A. Appendix

Table 1 lists the 125 three-term connection coefficients for L D 6, m D 0, n D 1, Q 0;1 .x/, computed using the algorithm described above.

Figure 1 .

 1 Figure 1. RMS of the residual error from all moment equations

  . The connection coefficient vector is of the form[START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Zhang | Comments on 'The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] unknowns for x 2 OE1; L 2 and 3L 2 9L C 7 unknowns for x D L 1. It can be shown the matrix S has q eigenvalues equal to 2 .1 m n/ , where[START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]

	Q	m;n .x/ D	h	m;n .1/; m;n .2/; ; m;n .L 1/	i T		(11)
		m;n .x/ D	8 < :h h	m;n x LC2 .x/; m;n x LC3 .x/; : : : ; m;n x 1 .x/ m;n 2 L .x/; m;n L 2 .x/ 3 L .x/; : : : ; m;n i T	i T	for x D 1; 2; ; L 2 for x D L 1	(12)
		m;n j .x/ D	8 < :h h	m;n j;x LC2 .x/; m;n j;x LC3 .x/; ; m;n j;x 1 .x/ i T m;n j; .x/; m;n j; .x/ j;C1 .x/; ; m;n i T	for x D 1; 2; : : : ; L 2 for x D L 1	(13)
	where D max.j C 2 L; 2 L/, D min.j C L 2; L 2/. The vector Q	m;n .x/ contains
	.L 2/ 3

Table 1 .

 1 Three-term connection coefficients for L D 6, m D 0, n D 1

	x	j	k	0;1 j;k .x/	x	j	k	0;1 j;k .x/	x	j	k	0;1 j;k .x/
	1	-3	-3	-7.15591411E-04	3	1	1	3.17107026E-01	5	-1	-2	8.74622766E-02
	1	-3	-2	1.35502680E-03	3	1	2	-1.33459578E-01	5	-1	-1	-2.48445901E-01
	1	-3	-1	7.87118799E-06	3	2	-1	5.29897817E-03	5	-1	0	3.13587753E-01
	1	-3	0	-6.14851011E-04	3	2	0	-1.15661098E-02	5	-1	1	-1.25272546E-01
	1	-2	-3	5.29897817E-03	3	2	1	2.54926579E-02	5	-1	2	-6.17732216E-03
	1	-2	-2	-1.15661098E-02	3	2	2	-1.93594670E-02	5	-1	3	3.29019694E-05
	1	-2	-1	2.54926579E-02	4	0	0	2.53068293E-08	5	0	-4	2.45926496E-06
	1	-2	0	-1.93594670E-02	4	0	1	4.96892129E-01	5	0	-3	2.37636417E-03
	1	-1	-3	-2.24262187E-02	4	0	2	3.29019100E-02	5	0	-2	9.35483714E-02
	1	-1	-2	6.06841495E-02	4	0	3	1.45424640E-03	5	0	-1	-6.27175506E-01
	1	-1	-1	-1.75635286E-01	4	1	0	-2.48445211E-01	5	0	0	3.61020792E-15
	1	-1	0	1.37823004E-01	4	1	1	3.13595282E-01	5	0	1	4.96891801E-01
	1	0	-3	8.23609549E-03	4	1	2	-1.25290082E-01	5	0	2	3.29026887E-02
	1	0	-2	1.96354133E-01	4	1	3	-6.13510292E-03	5	0	3	1.45229718E-03
	1	0	-1	-9.14074187E-01	4	2	0	-1.64544140E-02	5	0	4	1.52412730E-06
	1	0	0	7.09481500E-01	4	2	1	3.77788872E-02	5	1	-3	-4.45648096E-04
	2	-2	-2	-1.64544140E-02	4	2	2	-4.67007112E-02	5	1	-2	-2.07415149E-02
	2	-2	-1	3.77788872E-02	4	2	3	2.03167757E-02	5	1	-1	8.74622766E-02
	2	-2	0	-4.67007112E-02	4	3	0	-7.15591411E-04	5	1	0	-2.48445901E-01
	2	-2	1	2.03167757E-02	4	3	1	1.35502680E-03	5	1	1	3.13587753E-01
	2	-1	-2	8.73555499E-02	4	3	2	7.87118799E-06	5	1	2	-1.25272546E-01
	2	-1	-1	-2.49815835E-01	4	3	3	-6.14851011E-04	5	1	3	-6.17732216E-03
	2	-1	0	3.17107026E-01	5	-4	-4	-7.62063649E-07	5	1	4	3.29019694E-05
	2	-1	1	-1.33459578E-01	5	-4	-3	-4.46399112E-07	5	2	-2	1.33940712E-04
	2	0	-2	9.40317902E-02	5	-4	-2	1.94935669E-06	5	2	-1	4.92552163E-03
	2	0	-1	-6.20031558E-01	5	-4	-1	4.88738554E-07	5	2	0	-1.64513443E-02
	2	0	0	-1.91465368E-02	5	-4	0	-1.22963248E-06	5	2	1	3.78102693E-02
	2	0	1	5.42767481E-01	5	-3	-4	-3.24555703E-05	5	2	2	-4.67741857E-02
	2	1	-2	-2.24262187E-02	5	-3	-3	-7.26148591E-04	5	2	3	2.04916885E-02
	2	1	-1	6.06841495E-02	5	-3	-2	1.25180052E-03	5	2	4	-1.35890068E-04
	2	1	0	-1.75635286E-01	5	-3	-1	2.49826364E-04	5	3	-1	-3.24555703E-05
	2	1	1	1.37823004E-01	5	-3	0	-1.18818208E-03	5	3	0	-7.26148591E-04
	3	-1	-1	-2.48445211E-01	5	-3	1	4.45159357E-04	5	3	1	1.25180052E-03
	3	-1	0	3.13595282E-01	5	-2	-4	1.33940712E-04	5	3	2	2.49826364E-04
	3	-1	1	-1.25290082E-01	5	-2	-3	4.92552163E-03	5	3	3	-1.18818208E-03
	3	-1	2	-6.13510292E-03	5	-2	-2	-1.64513443E-02	5	3	4	4.45159357E-04
	3	0	-1	-6.27152133E-01	5	-2	-1	3.78102693E-02	5	4	0	-7.62063649E-07
	3	0	0	2.88213076E-04	5	-2	0	-4.67741857E-02	5	4	1	-4.46399112E-07
	3	0	1	4.96129851E-01	5	-2	1	2.04916885E-02	5	4	2	1.94935669E-06
	3	0	2	3.48068749E-02	5	-2	2	-1.35890068E-04	5	4	3	4.88738554E-07
	3	1	-1	8.73555499E-02	5	-1	-4	-4.45648096E-04	5	4	4	-1.22963248E-06
	3	1	0	-2.49815835E-01	5	-1	-3	-2.07415149E-02