Marco Fedi
email: marco.fedi@istruzione.it

Hydrodynamics of the dark superfluid: IV. Cosmological model

Keywords: Standard model of cosmology, cosmological redshift, quantum cosmology, superfluid vacuum, light propagation, dark energy PACS 98.80.Es, 98.62.Py, 42.25.Bs, 95.36.+x, 95.35.+d

Introduction

The standard FLRW ΛCDM big-bang model presents, as known, many unsolved criticalities, as also recently remarked by López-Corredoira [START_REF] López-Corredoira | Tests and problems of the standard model in cosmology[END_REF]. Cosmic inflation is for instance a necessary complementary theory to the model but itself not justified resorting to current physics. These problems range from the issues of homogeneity, flatness and magnetic monopoles (as we know, artificially solved by introducing inflation) up to baryon asymmetry, galaxy structure formation, dwarf galaxies scarcity, necessity of accelerated expansion, stable density of vacuum energy despite expansion etc. Nevertheless, the general acceptance of the standard model and the numerous efforts to make it work properly, are due to the fact that we have currently found a single valid explanation to the redshift observed by E. Hubble, that of an expanding universe. However, after decades, the feeling is now to insist on a path which is maybe not correct. As Hossenfelder has recently stated [START_REF] Hossenfelder | Science needs reason to be trusted[END_REF], "science needs reason to be trusted". We therefore believe that we need to reinterpretate the cosmological redshift, by starting from a better understanding of the physical vacuum and of the nature of light, as previously discussed in [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF][START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF]. The results are encouraging. Since the universe is filled at 95% with dark energy and dark matter (the dark superfluid (DS)), light has to propagate through them rather than through a real, Newtonian vacuum. A foaming quantum vacuum where particle-antiparticle pairs continuously appear and annihilate is perhaps nothing but the stochastic quantum hydrodynamic oscillation in the order parameter of the DS, in which vortex-antivortex pairs form and annihilate [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF], as also observed in liquid helium when cooled under a critical temperature and in polaritons superfluids [START_REF] Boulier | Lattices of quantized vortices in polariton superfluids[END_REF], even at room temperatures [START_REF] Lerario | Room-temperature superfluidity in a polariton condensate[END_REF]. We have speculated that a photon too might be a quantum hydrodynamic perturbation of the DS, a phonon propagating through it, since we know that phonons also propagate through fluids and superfluids, even as transversal waves [START_REF] Bolmatov | The phonon theory of liquid thermodynamics[END_REF] and possess wave-particle duality such as photons. Furthermore, the speed of sound in dilatant fluids cannot be surpassed, thus, by hypothesizing a dilatant behavior of the DS (due to its granular, quantum nature) under shear stress occurring in a relativistic regime, we would also explain the reason why the speed of light cannot be surpassed [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF]. Once it has transiently passed to a solid state, only sound, i.e. light in this case, can still propagate through it and this would pave the way to merging special relativity with quantum physics, specifically with quantum hydrodynamics (see chart and details in [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF], Sect. 4). Once taken into account the existence of the DS, the analogies photon-phonon were too many to make us desist from re-investigating the nature of light. In this case the speed of light would be dependent on the density of the DS. By taking into account this dependence, we obtain an exact solution to the complete data set concerning the cosmological redshift with a single equation, which indeed describes the speed of light through the DS. This without resorting to expanding and accelerating universe and obtaining a simpler and effective cosmological model (see table in the annex for a model comparison). Here the energy density of free space coincides with that of dark energy (whose value is still debated, see Sect. 2) that permeates the whole cosmos and we believe it possess, along with dark matter, superfluid features, analogous to those of a Bose-Einstein condensate [START_REF] Das | Dark matter and dark energy from Bose-Einstein condensate[END_REF], at the average temperature of the cosmic microwave background (CMB). As known, in cosmology it is indeed useful to treat the universe as a fluid, which some authors let coincide with dark energy [START_REF] Amendola | Dark Energy. Theory and observations[END_REF][START_REF] Huang | Dark energy and dark matter in a superfluid universe[END_REF][START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF], if a single-fluid model is taken into account. Exactly as Maxwell derived his equation for the speed of light from the density and the elasticity of the hypothesized ether, translated into the magnetic permeability and permittivity of vacuum, we use in our computations a formula for the speed of light based on its propagation through the D-S. In superfluid quantum gravity (SQG) [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF], the ether wind --or, to be precise, the flow of DS due to Bernoulli pressure exerted by massive particles being quantum vortices in the DS [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF][START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF] --corresponds to the gravitational field, as previously hypothesized also by other authors [START_REF] Cahill | Gravity as quantum foam in-flow[END_REF][START_REF] Kirkwood | The Physical Basis of Gravitation[END_REF]. In this case, the premise of the Michelson-Morley test (relative motion Earth-ether) [START_REF] Michelson | On the Relative Motion of the Earth and the Luminiferous Ether[END_REF] would have been inappropriate and no ether could be detected. The negative result of that test has forced us to reject the existence of a dark medium permeating the physical vacuum and this despite growing evidences from quantum physics (quantum vacuum) and cosmology (dark energy) about the non-existence of a real, Newtonian vacuum [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF][START_REF] Huang | A Superfluid Universe[END_REF][START_REF] Sbitnev | Physical Vacuum is a Special Superfluid Medium[END_REF][START_REF] Sbitnev | Dark matter is a manifestation of the vacuum Bose-Einstein condensate[END_REF][START_REF] Volovik | The Universe in a helium droplet[END_REF][START_REF] Das | Dark matter and dark energy from Bose-Einstein condensate[END_REF] as also investigated in recent works [START_REF] Amelino-Camelia | In vacuo dispersion features for gamma-ray-burst neutrinos and photons[END_REF][START_REF] Mignani | Evidence for vacuum birefringence from the first optical-polarimetry measurement of the isolated neutron star RX J1856.5-3754[END_REF], and on the contrary about the existence of a dark medium owning superfluid features which photons might propagate through as quasi-particles [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF], i.e. as quantized waves which possess a momentum. With this different premise, let us come back to the interpretation of the cosmological redshift and we will see that an exact law for the observed redshift (including the distant supernovae) may be obtained from the role of density on light propagation, without critical side-effects and apparently solving various other issues of modern cosmology, as summarized in the annex.

Structure of the universe

Here we assume that the dark superfluid be distributed in a spherical volume and that its density fades out in distance, up to be zero at its boundary. In the present model, the universe is therefore a bubble of dark superfluid showing a 1/r density decrease (see Fig. 1). By assuming that this bubble is not expanding, we obtain in Sect. 3 the correct chart for the cosmological redshift. Outside of this bubble nothing, nor light can propagate. By approaching the edge of the universe, light would progressively lose energy, redshifting up to vanish. In the model the spherical volume filled with dark superfluid is also spinning and we therefore agree with Longo [START_REF] Longo | Detection of a dipole in the handedness of spiral galaxies with redshifts z ∼ 0.04[END_REF], who hypothesizes a rotating primordial universe. A spinning bubble (at least previously spinning) is not necessary for our redshift computations but helps to solve other well-known issues. As for laboratory superfluids, where single-handed vortices emerge from a rotating superfluid (Fig. 1), such a rotation would explain baryogenesis without antimatter, solving the baryon asymmetry and also explaining the statistical one-handed bias of spiral galaxies [START_REF] Longo | Detection of a dipole in the handedness of spiral galaxies with redshifts z ∼ 0.04[END_REF] and the left-handed bias of life molecules [START_REF] Tasker | Why is life left-handed? The answer is in the stars[END_REF]. About the reasons for considering massive fundamental particles as quantized vortices in the DS see [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF]. As regards baryon asymmetry, the decay of B mesons is indeed predicted to favor the production of matter over antimatter but not enough to explain the huge preponderance of matter in the present universe. Recent tests at the LHCb [START_REF]The LHCb collaboration: Measurement of matterantimatter differences in beauty baryon decays[END_REF] showed that baryons too seem to violate the CP symmetry by following different decay paths but stronger data are necessary and this phenomenon is likely too rare in nature. Thus, with the present knowledge, baryon asymmetry remains unsolved in the standard ΛCDM big-bang model. On the contrary, Fig. 1 shows how a rotating superfluid could generate only same-handed particles, explaining the observable absence of antimatter in the universe.

The problem of vacuum energy density

From measurements of the WMAP probe, the energy density of vacuum (we assume dark energy) results to be ∼ 6 × 10 -9 J/m 3 . However, in our model, as long as quantum vortices do not form in the DS [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF][START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF], dark energy does not exert gravitational pull. This is clear also considering the fact that in a detected flat universe (Ω = 1) accelerated expansion cannot occur (as on the contrary necessary in the standard model). Thus, assuming that dark energy at rest (not hydrodynamically perturbed [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF][START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF]) does not exert gravitational pull, the measurement of the WMAP probe would not be correct. It appears to be first necessary to better theorize what dark energy is and how it behaves, to decide its role in gravitation. And it is also necessary to understand the quantum foundation of gravity, which could exactly bring Fig. 1 Above: in this superfluid model the universe, here imaginarily observed from an outer point of view, appears as a non-expanding, spinning bubble of dark superfluid, whose density fades out from the center as 1/r. It is not excluded that other universes may exist next to it but light (and any other kind of information) could not travel between them. The decreasing density would explain the redshift of distant supernovae and its rotation may justify the baryogenesis without antimatter and the formation of the indirectly observed cosmic web of dark matter [START_REF] Eckert | Warmhot baryons comprise 5-10 per cent of filaments in the cosmic web[END_REF] (the scaffolding of the universe) as a hydrodynamic perturbation of the DS (Fig. 2). Below: the formation of single-handed quantum vortices [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF] in a rotating superfluid as the probable reason for baryon asymmetry and for the single-handed bias observed in spiral galaxies [START_REF] Longo | Detection of a dipole in the handedness of spiral galaxies with redshifts z ∼ 0.04[END_REF] and organic molecules [START_REF] Tasker | Why is life left-handed? The answer is in the stars[END_REF]. Fig. 2 Left: (a) Metal atoms trapped in superfluid helium vortices highlight a structure of vortex filaments [START_REF] Moroshkin | Nanowire formation by gold nano-fragment coalescence on quantized vortices in He II[END_REF]; (b) galactic filaments of dark matter [START_REF] Eckert | Warmhot baryons comprise 5-10 per cent of filaments in the cosmic web[END_REF] which baryon matter aggregates on forming stars and galaxies. Here the relationship between dark energy and dark matter is the same existing between superfluid helium and the vortex filaments which manifest in it, i.e. we suppose dark matter is a hydrodynamic manifestation of superfluid dark energy [START_REF] Huang | Dark energy and dark matter in a superfluid universe[END_REF][START_REF] Sbitnev | Dark matter is a manifestation of the vacuum Bose-Einstein condensate[END_REF]. into play dark energy [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF]. With the DS at rest only pressure is exerted (from its energy density, J/m 3 ⇒ Pa) which balances the gravitation of baryon and dark matter. Since in a superfluid a quantum vortex is a topological defect where the field density drops to zero in a region having radius ξ (which refers to the vortex healing length, see [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF][START_REF] Huang | A Superfluid Universe[END_REF]), it is very likely in our model that the value for vacuum's energy density could be extremely high and finite, surely closer to the calculations of quantum field theory than to those of the WMAP spacecraft. If we imagine fundamental massive par-ticles as bubbles in a fluid (also in a bubble, as in a vortex, the fluid density drops to zero) it is easy to understand how vacuum can possess a much higher density than ordinary matter, even of that of a neutron star, since the neutron is represented in our superfluid approach [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF] as a three-vortex system (three quarks) and in each vortex we have field density |ψ| 2 = 0. Baryon matter is not more dense than the DS but is rather a small part of it possessing great kinetic energy (i.e. spin, as a particle's total internal kinetic energy [START_REF] Esposito | On the Role of Spin in Quantum Mechanics[END_REF][START_REF] Recami | [END_REF][START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF]). We perceive the vortex and its action through the medium (gravity, electromagnetism), not the medium directly, since it interacts with us through its topological defects (particles) and the forces they generate. This is in our opinion the reason why we cannot directly detect dark energy. We have to search for its manifestations and the most evident may be light itself [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF]. Thus, the reason why we cannot detect the gravitational pull of this enormous amount of dark energy at rest is because no gravity is exerted unless stable (not the ephemeral vacuum's fluctuations) vortices arise. By way of comparison, we could say that we cannot detect the air around us unless it is perturbed. And if it is not, this does not mean we live in a vacuum. So, the question with dark energy is "what we expect to observe". We think that the universe, hence the observer himself [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF], are quantum hydrodynamic manifestations of the DS. We could be what we are trying to observe.

Exact law for cosmological redshift without expansion

We use here the model of photon [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF] as a quasi-particle (transverse phonon, respecting its wave-particle duality) which acoustically propagates through the DS. Superfluid light in bulk nonlinear media has been described by Carusotto [START_REF] Carusotto | Superfluid light in bulk nonlinear media[END_REF] and a superfluid behavior of photons has been observed in exciton-polariton condensates [START_REF] Leboeuf | Superfluid motion of light[END_REF][START_REF] Dominici | Real-space collapse of a polariton condensate[END_REF]. Also the DS would be a sort of dark Bose-Einstein condensate [START_REF] Das | Dark matter and dark energy from Bose-Einstein condensate[END_REF], at the temperature of the CMB radiation 2 . If light is the acoustic excitation of the DS, a sound that we "hear" with our eyes, we can therefore use the formula for the speed of sound in a fluid medium, v s = K/ρ, expressing the ratio bulk modulus to density, and resort to the isentropic compressibility, β = 1/K, to define the speed of light as c = 1

β d ρ d (1)
2 In our model the CMB radiation is not the residual radiation from recombination (no Big Bang) but the intrinsic average temperature (∼ 2.725K) of the dark superfluid at rest, similar to that of other superfluids (∼ 2.172K for phase transition in 4 He). Indeed, in a rotating superfluid universe (Fig. 1) baryon matter can form without Big Bang and without producing antimatter. This helps to solve baryon asymmetry.

as also in [START_REF] Amendola | Dark Energy. Theory and observations[END_REF]. As expected, we see that the speed of light is inversely proportional to the density of the DS. Let us plot (Fig. 3) c = P d /ρ d with P d = P d 0 + kρ d , where P d 0 and k are arbitrary constants, to observe how the speed of light varies with density. Since we assumed that density decreases moving away from the center of the universe (Fig. 1), the speed of light will increase in distance equally in all directions from a central point of observation. The most distant a galaxy or a supernova, the more redshift its light will undergo at detection. Here it is important to specify that once light reaches our instruments on Earth or in the space, it settles on the local value of speed, i.e. the known c = 299972458 m/s, the value allowed by the local density (where local might be approximatively extended to beyond the solar system, considering the cosmic distances involved in the redshift measurement). A frequency deviation from that at the source has therefore occurred (due to the lower velocity in our frame of reference), which is detected as a redshift (Fig. 4) and this also applies to distant supernovae, as ascertainable from Fig. 5. In short, the curve of the observed redshift is nothing but the reciprocal of the speed of light and this is the core of the present model. Also in Hubble's law, z = H 0 D/c, the 1/c dependence appears, albeit there c is a constant. And this is the problem, since the law is in that case linear and does not fit the observations, being not a suitable mathematical law for the cosmological redshift. In our formula the variation of distance is implicit in the varying (1/r) value of density and the Hubble parameter disappears, substituted if anything by a sim- ple coefficient that also eliminate the units or we can otherwise measure the redshift in [s/m]. Abandoning the expanding Big Bang universe, the use of FLRW metric becomes invalid and we explain the observations gaining a dramatic simplification in the mathematics. Eventually, to conclude our speculative reasoning, it is interesting to see what happens if the point of observation (Earth) is not exactly in the center of the bubble-universe. Fig. 6 shows that this fact would add deviations from the predicted redshift curve, since the measurement would not be anisotropic anymore, as regards the chosen direction of observation. For instance, observing along r, would yield different redshift data for some equidistant galaxies, e.g. those located at the equidistant points D and K. This can in part explain the deviation from the theoretical curve which, for some objects, is much higher than expected.

Galaxy formation

Structure formation at the level of single galaxies is another weak spot of the standard model. Starting from a conjecture by Lathrop [START_REF] Huang | A Superfluid Universe[END_REF]31] we think that galaxies have formed from Fig. 6 deviations of the observed redshift from the predicted curve (Fig. 4,5) due to the position of the point of observation (Earth) in the universe. Concentric grey circles express a discretized variation of density in distance and white circles represent equidistant points from the Earth (⊕). We notice that equidistant galaxies (e.g. B and M) may show a different redshift. This could explain part of the observed deviations (Fig. 5). However, taking into account all the possible different directions (observing for instance along s, t, u) the predicted curve is on average respected, with more or less deviation, assuming the observer is not excessively decentralized toward the edge of the universe. This means that our galaxy is relatively central and that the observation of the sky from peripheral regions of the universe would be strongly anisotropic. An analysis of redshift deviations to verify the hypothesis illustrated in Fig. 6 is suggested and could yield an approximate position of the Earth in the superfluid, finite, cosmos). Fig. 7 Structure formation of a spiral galaxy (here NGC 1365) in a superfluid universe. Elaboration from an idea by Lathrop [31,[START_REF] Huang | A Superfluid Universe[END_REF], compatible with the superfluid model. A giant vortex of dark matter in the DS attracts dust (also gravity is here considered as a quantum hydrodynamic phenomenon in the DS, see SQG [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF]) which clumps and accretes on the vortex (we say due to Bernoulli pressure [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF]25,[START_REF] Huang | A Superfluid Universe[END_REF]). Because of the gravitational interaction, baryon clumps attract each other toward the center of the vortex, forming a black hole.

gas and dust accretion onto large vortices of dark matter, due to gravitational interaction (here interpreted as the action of Bernoulli pressure [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF] therefore again a hydrodynamic phenomenon in the DS). This framework is illustrated in Fig. 7. The rotation of a galaxy would directly derive from that of the original vortex. The accreted clumps of baryon matter attract each other, moving toward the center of the ring. This originates a black hole, believed to be present at the center of galaxies. The way clumps are distributed on the ring would determine the shape of the galaxy. It is interesting to note that since a uniform distribution is unlikely, galaxies have generally assumed a spiral shape, while elliptical galaxies derive from merged galaxies. Irregular galaxies would have formed due to baryon and dark matter clustering without a "pilot-vortex". About the anomalous redshift deviations from galaxies in the Virgo cluster, Arp [START_REF] Arp | What is the mean redshift of the Virgo cluster?[END_REF] thinks that there seem to be no alternative explanation other than a dependence on morphological galaxy types, where Sc galaxies show the largest redshift and the Sb the smallest value. Within our model, we argue it may depend on the different distribution and concentration of dark matter in galaxies depending on their shape. Also in this case, different density in the DS would therefore affect the equation of light (2).

Conclusion

As summarized in the table (see annex), our investigation shows that a cosmological model assuming a (once?) rotating, superfluid, approximatively spherical, universe whose density fades out as 1/r and where light corresponds to its acoustic excitation [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF] is simpler and at first glance more effective than the FLRW ΛCDM big bang model, to justify the ontogeny of our universe and the observed redshift according to the available data set [START_REF] Amanullah | Spectra and HST light curves of six type Ia supernovae at 0.511¡z¡1.12 and the Union2 compilation[END_REF], without the necessity of complementary theories. Just taking into account the dependence of the speed of light on the density of the DS, we obtain a law for the redshift, without resorting to an expanding universe and to the mathematics of general relativity (which also has to get along with quantum theory and the framework becomes increasingly critical). On the contrary, this superfluid model predicts an approximately static universe, where galaxies have not formed from primordial fluctuations which have then expanded but in loco, due to quantum hydrodynamic perturbations of the DS. We deduce that the DS has existed ever since [START_REF] Ali | Cosmology from quantum potential[END_REF], letting it obey the principle of energy conservation, and its history is maybe older than the universe itself, without an adimensional singularity containing all the energy-mass of the universe which, at some point, decided to instantly swell and then to continue its (accelerated) expansion.

A superfluid model seems able to also justify the indirectly observed cosmic web of dark matter (Fig. 2), baryon asymmetry (Fig. 1), the statistically relevant single-handed bias of spiral galaxies [START_REF] Longo | Detection of a dipole in the handedness of spiral galaxies with redshifts z ∼ 0.04[END_REF] and that of organic molecules [START_REF] Tasker | Why is life left-handed? The answer is in the stars[END_REF], the small amount of dwarf galaxies, the structure formation of single galaxies (Fig. 7) and the reason why the speed of light cannot be surpassed [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF]. Dark energy theory should be improved in that dark energy would not exert gravitational pull (nor a supposed "negative gravity") as long as vortices (massive particles) do not form in it [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF][START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF]: at rest, dark energy only exerts repulsive pressure (as energy density) able to balance the gravitation of baryon matter and to avoid a gravitational collapse. The present model suggests that quantum cosmology does not need to reconcile general relativity with quantum physics. Even because such a reconciliation seems to lead us anyway to the superfluid approach, where curved spacetime gives way to the hydrodynamics of the DS [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF]. based on dark matter evolution model the small number of dwarf galaxies cannot be explained not present; the dark matter web on which galaxies form arises from the hydrodynamic perturbation of the cosmic superfluid, with no prevalence of dwarf galaxies baryogenesis due to the Big bang and linked to its criticalities due to the rotation of the cosmic superfluid massive fundamental particles arise and possess spin; this also solves several issues of particle physics [START_REF] Fedi | Hydrodynamics of the dark superfluid: I. Genesis of fundamental particles[END_REF] baryon asymmetry and one-handed bias not explained cosmic superfluid spinning in a single direction justifies the absence of antimatter in the universe, the statistical one-handed bias of spiral galaxies [START_REF] Longo | Detection of a dipole in the handedness of spiral galaxies with redshifts z ∼ 0.04[END_REF], and the left-handed bias of life-molecules [START_REF] Tasker | Why is life left-handed? The answer is in the stars[END_REF] redshift of distant supernovae not explained by linear Hubble's law. Accelerated expansion necessary. explained by a single equation, no accessory hypotheses needed galaxy structure formation not explained from gravitatonal accretion of dust onto great vortices in the cosmic superfluid (according to Lathrop). This explains the shape of galaxies, their rotation and the presence of a black hole at their center future evolution of the universe big crunch or thermal death depending on ρc approximatively static universe where the gravitation of baryon and dark matter [START_REF] Fedi | Hydrodynamics of the Dark Superfluid: III. Superfluid Quantum Gravity[END_REF] is balanced by the internal pressure of dark energy speed of light as upper limit to acceleration not explained due to the dilatancy of the cosmic superfluid under shear stress occurring in relativistic regime (v → c) [START_REF] Fedi | Hydrodynamics of the dark superfluid: II. photon-phonon analogy[END_REF] nature of light problem light as an odd wave without propagation medium which is also a particle and has a constant speed everywhere in the universe light corresponds to transverse phonons through the quasi-lattice of the dark superfluid, i.e. to the acoustic excitation of dark energy. Actually, a photon is therefore a quasi-particle (phonon) still owning wave-particle duality (due to possessing quantized momentum). The speed of light depends on the density of the DS and is therefore not constant in different parts of the universe. The speed increases with the radius of the spherical universe, following the decrease in density.

A -Model comparison

Fig. 3

 3 Fig. 3 Speed of light-density plot. The scale has been set by way of example acting on the arbitrary constants P d 0 and k, in the interval between ρ d -→ 0 and ρ d = 1.

Fig. 4

 4 Fig. 4 Redshift caused by decreasing speed of light due to the density of the DS fading out with distance. As the curve suggests, the known value for c may be valid with good approximation only within cosmologically short distances.

Fig. 5

 5 Fig. 5 Distance modulus-redshift plot. Data including the supernova Union2 compilation [17].


THE SWITCHING SYSTEMS THEORY

The optimal station keeping control problem of a geostationary satellite equipped with electric thrusters is recast as a switched system commutation times optimisation problem, considering that a system with a bang-bang control profile is composed by several subsystems, one for which all the controls are off and the other ones one for which one of the controls is on. In order to optimise the commutation times, the optimal firing sequence has to be known in advance. This sequence is provided by a two-step decomposition technique and the proposed method can be interpreted as a third step. Simulation results on a realistic example validate the benefit of this third optimisation step on the control sequence fuel consumption.

Introduction

Spacecraft orbiting the Earth on Geostationary Earth Orbits (GEO) undergo orbital disturbing forces, resulting in a natural drift outside their operational station keeping (SK) windows (a rectangular box of a given geographical longitude and latitude range). It is therefore mandatory to design an accurate SK control strategy, in order not to let the spacecraft operating conditions deteriorate.

The usual spacecraft propulsion system is composed of chemical thrusters (see Soop 12 or Sidi [START_REF] Sidi | Spacecraft Dynamics and Control[END_REF] ). However, the idea of using electric propulsion has been proposed in the sixties (see for instance the references Barret 2 and Hunziker [START_REF] Hunziker | Low-Thrust Station Keeping Guidance for a 24-Hour Satellite[END_REF] ), followed by some theoretical developments in the eighties (as for example the development conducted by Anzel 1 and Eckstein [START_REF] Eckstein | Optimal Station Keeping by Electric Propulsion With Thrust Operation Constraints[END_REF] ). Despite the difficult operational constraints arising from the use of the electric propulsion (disjunctive thrusts, minimum elapsed time between two thrusts, on-off profile of the thrusters), this kind of propulsion is a viable alternative to the commonly used chemical one thanks to its high specific impulse. Indeed, this particular feature naturally imply fuel consumption savings enabling to increase the spacecraft longevity.

The minimum-fuel station keeping problem may be mathematically expressed as an optimal control problem (OCP). Two types of methods may be used to solve optimal SK control problems. The collocation based direct methods rely on a discretisation of the state and the control vectors in order to create a non linear programming problem (see Betts 4 or Hull [START_REF] Hull | Conversion of optimal control problems into parameter optimization problems[END_REF] ) whose solution is a discretized approximation of the optimal solution. Indirect methods are based on the Pontryagin Maximum Principle (PMP). Despite the difficulty of handling the operational constraints of the electric propulsion within the PMP framework, Gazzino et al. [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] set up a decomposition method to solve the optimal station keeping control problem taking the operational constraints into account. However, due to the great number of tuning parameters and to the complexity of the problem, the obtained local solutions appear to be highly sensitive to the initialisation process, thus resulting in some convergence issues in some cases or in a possible optimality deficiency.

If the operational constraints are respected, the spacecraft trajectory is composed of thrusting arcs separated by coasting arcs. Thus, the control profile switches from a time interval where all thrusters are off to a time interval for which one thruster is on, and vice-versa. Therefore, the system can be viewed as a switched system composed by one subsystem per thruster and one subsystem to describe the case where all thrusters are off. Xu and Antsaklis [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] proposes a method in order to compute the optimal switching times of switched systems thanks to a time change of coordinates. The idea of this paper is to use the method of Xu and Antsaklis [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] to optimise the switching times for the optimal station keeping control problem. As this method requires to know in advance the optimal firing sequence of the thrusters, the two-step decomposition technique of Gazzino et al. [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] is first used as an initial guess. The proposed optimisation step can hence be considered as a third step of the decomposition method. The efficiency of the proposed technique for improving the results of Gazzino et al. [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] is illustrated on a realistic numerical example.

Problem statement 2.1 Geostationary Spacecraft Dynamics

The motion of a spacecraft orbiting the Earth on a Geostationary Earth Orbit (GEO) can be described with the equinoctial orbital elements as defined in Battin: 3

x eoe = [a e x e y i x i y MΘ ] t ∈ R 6 , ( 1
)
where a is the semi-major axis, (e x , e y ) the eccentricity vector, (i x , i y ) the inclination vector, MΘ = ω + Ω + M -Θ is the mean longitude where Ω is the right ascension of the ascending node, ω is the perigee's argument, M is the mean anomaly and Θ(t) is the right ascension of the Greenwich meridian. The dynamics of the spacecraft may be represented by the following non linear state-space model:

dx eoe dt = f L (x eoe , t) + f G (x eoe , t)u. (2)
f L ∈ R 6 is the Lagrange contribution part of the external disturbing forces. For a GEO spacecraft, these disturbing forces are the non-spherical part of the Earth gravitational potential that mainly affects the mean longitude MΘ , the solar radiation pressure (SRP) that mainly affects the eccentricity vector (e x , e y ) and the gravitational attraction of the sun and the moon that mainly affects the inclination vector (i x , i y ). These forces are described by the CNES ORANGE model (cf. Campan et Brousse [START_REF] Campan | ORANGE: Orbital analytical model for geosynchronous satellite[END_REF] ). f G ∈ R 6×3 is the Gauss contribution part for the disturbing forces that do not derive from a potential. u = u R u T u N t ∈ R 3 is the control vector expressed in the local orbital RT N frame (also written RS W) defined in Soop [START_REF] Soop | Handbook of Geostationary Orbits[END_REF] by:

• N is the unit vector along the kinetic momentum;

• R is the unit vector along the direction Earth's center -satellite;

• T completes the right-handed orthogonal direct basis.

In order to deal with the station keeping problem, the relative state of the satellite with respect to the station keeping state is defined as:

x sk = a sk 0 0 0 0 MΘ sk t

where a sk is the synchronous semi-major axis and MΘ sk is the station mean longitude.

The relative dynamics equations are developed by linearization of Equation ( 2) about the station keeping point (3). By denoting x = x eoex sk the relative state model for the SK problem reads:

dx dt = A(t)x + D(t) + B(t)u, (4)
where the matrices A ∈ R 6×6 , B ∈ R 6×3 , C ∈ R 3×6 and D ∈ R 6 are defined as follows:

A(t) = ∂ ∂x eoe ( f L (x eoe (t), t)) x eoe =x sk , (5)
B(t) = f G (x sk , t), (6)
D(t) = f L (x sk , t), (7)
The geographical coordinates of the satellite:

y eoe = T (x eoe , t)x eoe , (8)
are of interest because the station keeping problem consists in constraining them in the vicinity of the station position y sk = r sk 0 λ sk t where r sk is the synchronous radius and λ sk is the station keeping geographical longitude. The relative geographical position with respect to the station-keeping position is denoted by:

y = y eoe -y sk = T (x sk , t)x = C(t)x, (9)
and is obtained by linearising Equation (8).

After the linearisation, the station keeping requirement on the latitude and the longitude of the spacecraft is expressed by constraining the relative geographical position with respect to the centre of the SK window:

|[0 1 0]C(t)x(t)| δ and |[0 0 1]C(t)x(t)| ∀t ∈ [0, T ], (10)
with δ being the half width of the SK window in the latitude and longitude directions. A trajectory satisfying the SK constraints ( 10) is called SK-feasible.

Electric Thrusting System

The considered satellite is equipped with four electric thrusters mounted on the anti-nadir face, each of them having an orientation defined by a cant angle θ and a slew angle α. These angles define the North-East, the South-East, the North-West and the South-West directions of thrust. The satellite dynamics can be expressed considering the four thrusts provided by the four engines as control variables. The control u(t) expressed in the local orbital frame is a linear combination of the four thrusts such that u = ΓF, where [START_REF] Betts | Survey of Numerical Methods for Trajectory Optimization[END_REF] . The thrust direction matrices Γ j ∈ R 3 are defined such that:

Γ = [Γ 1 | Γ 2 | Γ 3 | Γ 4 ] ∈ R 3×4 and F = [F 1 F 2 F 3 F 4 ] t ∈ [0, F max ]
Γ j = 1 m -sin θ j cos α j -sin θ j sin α j -cos θ j t , (11)
where the cant and slew angles θ j and α j respectively are defined exactly as in Anzel. [START_REF] Anzel | Controlling a stationary orbit using electric propulsion[END_REF] For the sake of simplicity, the thrust vector is normalized by the maximum level of thrust F max . It is thus possible to write F = F max F with F ∈ [0, 1] 4 . As the thrust are on-off, the thrust profile is modelled as a rectangular signal that is parametrized by the date t i, j corresponding to the middle instant of the thrust and by its half width duration denoted ∆t i, j :

F i (t) = F max P i j=1 RectangleFunction(t, t i, j , ∆t i, j ). ( 12
)
where P i is the number of thrusts of thruster i. Some technological operational constraints that are inherently related to the actuation system have to be taken into account:

(a) thrusters cannot have simultaneous thrusts;

(b) a thrust must last at least T l : 2∆t i, j T l ;

(c) two successive thrusts of a given thruster must be separated of an interval of latency equal to T s ;

(d) two thrusts of two different thrusters must be separated by an interval of latency equal to T d .

The constraint for the time latency between the thrust k of thruster i and the thrust l of thruster j is thus mathematically expressed as:

|t i,k -t j,l | -(∆t i,k + ∆t j,l ) K i, j , (13)
for k = 1 . . . P i and l = 1 . . . P j , where K i, j = T s if i = j (constraint (c)) and K i, j = T d otherwise (constraint (d)). In addition, some other convenient constraints are enforced in order to make sure that the thrusts will not begin before 0 and end after T : t i, j -∆t i, j 0 and t i, j + ∆t i, j T.

A trajectory satisfying the operational constraints for the electric propulsion system is called operationally feasible.

Statement of the Station Keeping Optimal Control Problem

The main goal of the station keeping system is to maintain the longitude and the latitude of the satellite in a box defined by its size δ by acting on the orbital parameter via the 4 thrusters while optimising the fuel consumption. The associated optimal control problem is in general defined over a fixed horizon for the computation of optimal open loop control laws. In this context, optimality means that a minimum fuel-solution is looked for to extend the operational life time of the satellite. Therefore, performing minimum-fuel station keeping is minimizing the performance index:

J = 2F max 4 thruster i=1 ||Γ i || 1 P i j=1 ∆t i, j , (15)
Considering all the constraints described above, the minimum-fuel SK problem to solve may be summarized as the following optimal control problem:

Problem 1 Find the sequence of dates {t i, j } and durations {∆t i, j }, for i = 1 . . . 4, j = 1 . . . P i (with P i fixed) solutions of the minimization problem:

min t i, j ,∆t i, j J = 4 i=1 P i j=1 ∆t i, j , s.t.                          ẋ(t) = A(t) x(t) + D(t) + B(t)ΓF max F(t), x(0) = 0, |[0 1 0]C(t)x(t)| δ, |[0 0 1]C(t)x(t)| δ, 2∆t i, j T l , t i, j -∆t i, j 0, t i, j + ∆t i, j T, |t i,k -t j,l | -(∆t i,k + ∆t j,l ) K i, j , (16)
• 3. Problem resolution in the Switching Systems Framework

Transformation of the Optimal Control Problem

The analysis of the geostationary orbital perturbations done in the existing literature (see e.g. Campan and Bousse, [START_REF] Campan | ORANGE: Orbital analytical model for geosynchronous satellite[END_REF] Sidi 11 or Soop [START_REF] Soop | Handbook of Geostationary Orbits[END_REF] ) allow to determine some rule of thumb for geostationary SK strategies. The North-South out-ofplane effect of the Sun and the Moon attraction are the most predominant forces that have to be corrected each half orbit, once in the South direction and half an orbit later in the North direction. According to this rule of thumb, the control scheme that should be applied is such that the two North thrusters are firing simultaneously, as well as the two South thrusters of the considered propulsion system. Nevertheless, this resulting control sequence does not fulfil the operational constraint (a).

In the reference Gazzino et al., [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] a two-step decomposition method have been used in order to solve the station keeping problem while satisfying the operational constraints (a)-(d). However this decomposition method suffers from the drawback of necessitating a fine tuning of a lot of parameters. The aim of the proposed contribution is to add a third step that can help to simplify this parameter tuning step.

From the method of Gazzino et al., [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] the control sequence that has been created is a combination of coast arcs, of firing arcs for the North-East thruster, of firing arcs for the South-East thruster, of firing arcs of the South-East thruster and of firing arcs for the South-West thruster, that is to say, a combination of the five subsystems:

ẋ(t) =                              A(t)x(t) + D(t), A(t)x(t) + D(t) + B(t)ΓF max 1 0 0 0 T , A(t)x(t) + D(t) + B(t)ΓF max 0 1 0 0 T , A(t)x(t) + D(t) + B(t)ΓF max 0 0 1 0 T , A(t)x(t) + D(t) + B(t)ΓF max 0 0 0 1 T . ( 17
)
By denoting U i a R 4 vector with a 1 at the i th position, and by extension U 0 the R 4 zero vector, all the dynamics functions:

f i (x(t)) = A(t)x(t) + D(t) + B(t)ΓF max U i , i ∈ {0, . . . , 4}, (18)
can represent a subsystem of the overall system. Hence, the time horizon optimisation [0, T ] can be splitted into K intervals T k for which the control vector V k is one of the five admissible control vectors U i , i ∈ {1, . . . , 4}. With these considerations, all the dynamics functions f k can represent a subsystem of the overall system. Hence, the commutation from V k to V k+1 between the intervals T k and T k+1 at commutation time t k can be viewed as a commutation between two of the five subsystems. In a switched system framework, the system dynamics can be written as:

for t ∈ T k , ∃i k ∈ {0, . . . , 4}, ẋ(t) = f i k (x(t)), (19)
where the following equality holds :

V k = U i k .
In accordance with the operational constraints (c)-(d) two thrusting arcs must be separated by a coasting arc. Thus the control sequence has the following structure:

U 0 , V j 1 , U 0 , V j 2 , U 0 , . . . , U 0 , V j P , U 0 ( 20
)
where V j k = U j with j ∈ {1, . . . , 4} and P is the optimal number of thrusts. There are therefore 2P switches and 2P + 1 intervals on which the control is constant. It is possible to write T = t 2P+1 as if the final time were the last commutation time. In order to satisfy the operational constraints (c) and (d), the length of the firing arcs must respect the constraints:

• t 2k+1t 2k T s if the coasting arc lies between two firing arcs of the same thruster,

• t 2k+1t 2k T d if the coasting arc lies between two firing arcs of two different thrusters.

The operational constraint (b) is then written as:

t 2k -t 2k-1 T l . (21)
With this reformulation using the five subsystems of the overall system, the rectangular signal of the thrust can be parametrised by its start and end times as:

F i (t) = F max P i j=1
RectangleFunction(t, t j,on , t j,off ).

(

) 22
The performance index is thus rewritten as:

J = F max P i=1 ||Γ i || 1 (t 2i -t 2i-1 ), (23)
since odd intervals are coasting arcs and the even intervals are firing arcs. The Optimal Control Problem 1 is thus rewritten in the following form:

Problem 2 Find the optimal switching sequence {t k }, k ∈ {1, . . . , 2P}, with P fixed, solution of the minimisation problem:

min {t k } J = F max P i=1 (t 2i -t 2i-1 ), s.t.                                ẋ(t) = f i k (x)(t), |[0 1 0]C(t)x(t)| δ, |[0 0 1]C(t)x(t)| δ, x(0) = 0, t 2k -t 2k-1 T l , t 2k+1 -t 2k T g , (24)
with T g = T s for the firing arcs of the same thruster and T g = T d for the firing arcs of two different thrusters.

•

The Problem 2 is an OCP for which the decision variables are the commutation times. It has been partly solved in Gazzino et al. [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] by a numerical two-step decomposition approach but with no guarantee of optimality of the firing times thus obtained. The first step of this method consists in removing the operational constraints and in solving the simplified SK OCP with a direct method initialised by an indirect method. The trajectory obtained after this first step is in general SK-feasible but not operationally feasible. In order to recover rectangle functions from the continuous control profile, a threshold parameter ς has to be chosen carefully. The second step aims at creating an operationally feasible trajectory from the one computed at step 1. To this end, two equivalence schemes have been created, the first one being called the Consumption Based Equivalence (CBE) and the second one called the Effect Based Equivalence (EBE) (see Gazino et al 8 for details and the significance of these appellations). The SK-feasible and operationally feasible trajectory resulting from step 2 is very sensitive to the value of the threshold parameters ς. Therefore, our objective is to improve this preliminary solution by using the method of Xu and Antsaklis [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] based on the switched systems theory.

Optimisation of the Commutation Times with the Switched Systems Theory

The reference Xu and Antsaklis [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] presents a technique to optimise the commutation times between the different subsystems of a compound system that is used in order to optimise the switching times between the coasting and the firing arcs. This method resorts to a time change of coordinates and a parametrisation of the commutation times. However it is mandatory to assume that the order of the sequences of the active subsystems is a feasible one, obtained by application of the result of Gazzino et al. [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] Seeking this on and off times is equivalent to find an optimal subsystems commutation sequence.

A time change of coordinates is used to parametrise the switching times:

t = t k + ∆ k (τ -k) if t ∈ [t k , t k+1 ], (25)
with ∆ k = t k+1t k . The system dynamics is rewritten with the new time variable:

∂x(τ) ∂τ =                  ∆ 2k-2 [A(t(τ))x(τ) + D(t(τ))] if τ ∈ [2k -2, 2k -1], ∆ 2k-1 [A(t(τ))x(τ) + D(t(τ)) + B(t(τ))ΓF max U k ] if τ ∈ [2k -1, 2k],
. . .

∆ 2P [A(t(τ))x(τ) + D(t(τ))] if τ ∈ [2P, 2P + 1], (26)
for k = 1, . . . , P. The state vector can be now considered as a function of the new time variable τ and of the switching times t k : x = x(τ, {t k }).

In order to handle state constraints, in the reference Gazzino et al., [START_REF] Gazzino | A Minimum-Fuel Fixed-Time Low-Thrust Rendezvous Solved with the Switching Systems Theory[END_REF] where a fuel-optimal rendezvous problem has been solved with the Xu and Antsaklis [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] technique, a target term has been added to the objective function. However, as the station keeping constraint are expressed as trajectory state constraints, this constraint is handled in the objective function by adding a penalisation term:

ψ(t,x(t)) = 1 2 [C 2 (t)x(t) -δ] 2 sign C 2 (t)x(t) -δ + 1 + 1 2 [C 2 (t)x(t) + δ] 2 sign -C 2 (t)x(t) -δ + 1 + 1 2 [C 3 (t)x(t) -δ] 2 sign C 3 (t)x(t) -δ + 1 + 1 2 [C 3 (t)x(t) + δ] 2 sign -C 3 (t)x(t) -δ + 1 , (27)
where

C 2 (t) = [0 1 0]C(t) and C 3 (t) = [0 0 1]C(t).
Applying the temporal change of coordinates, the new objective function taking the station keeping window penalisation term into account reads:

J(x, {t k }) = F max P i=1 (t 2i -t 2i-1 ) + µ 2P+1 k=1 (t k -t k-1 ) k k-1 ψ(t(τ), x)dτ, (28)
with µ a parameter that has to be chosen. The station OCP reads thus in the switching systems framework:

Problem 3 Find the optimal, switching sequence {t k }, k ∈ {1, . . . , 2P}, with P fixed, solution of the minimisation problem:

min x,{t k } J(x, {t k }) = F max P i=1 (t 2i -t 2i-1 ) + µ 2P+1 k=1 (t k -t k-1 ) k k-1 ψ(t(τ), x)dτ, such that                          ẋ(τ) =        ∆ 2k-2 [A(t(τ))x(τ) + D(t(τ))] if τ ∈ [2k -2, 2k -1], ∆ 2k-1 [A(t(τ))x(τ) + D(t(τ)) + B(t(τ))ΓF max U k ] if τ ∈ [2k -1, 2k] x(0) = 0, t 2k -t 2k-1 T l , t 2k+1 -t 2k T g , (29)
with T g = T s for the firing arcs of the same thruster and T g = T d for the firing arcs of two different thrusters.

•

As the switching times determine the structure of the control profile, the state trajectory is recovered by propagation of the system dynamics from the initial condition.

Problem 3 can be solved with a descent algorithm. To this end, it is necessary to compute the derivative of the performance index in order to obtain the descent direction. The computation of derivative requires to distinguish the odd and the even commutation times:

dJ dt 2l-1 = -F max + µ 2l-1 2l-2 ψ(t(τ), x(τ)dτ + µ(t 2l-1 -t 2l-2 ) 2l-1 2l-2 (τ -2l + 2) ∂ψ ∂t (t(τ), x(τ)) + ∂ψ ∂x t (t(τ), x(τ)) ∂x ∂t 2l-1 (τ) dτ -µ 2l 2l-1 ψ(t(τ), x(τ)dτ + µ(t 2l -t 2l-1 ) 2l 2l-1 (-τ + 2l + 2) ∂ψ ∂t (t(τ), x(τ)) + ∂ψ ∂x t (t(τ), x(τ)) ∂x ∂t 2l-1 (τ) dτ, (30)
dJ dt 2l = F max -µ 2l+1 2l ψ(t(τ), x(τ)dτ + µ(t 2l+1 -t 2l ) 2l+1 2l (-τ + 2l + 1) ∂ψ ∂t (t(τ), x(τ)) + ∂ψ ∂x t (t(τ), x(τ)) ∂x ∂t 2l-1 (τ) dτ + µ 2l 2l-1 ψ(t(τ), x(τ)dτ
+ µ(t 2l -t 2l-1 ) 2l 2l-1 (τ -2l + 1) ∂ψ ∂t (t(τ), x(τ)) + ∂ψ ∂x t (t(τ), x(τ)) ∂x ∂t 2l-1 (τ) dτ. (31)
The derivatives of the state vector with respect to the commutation times verify the following relation:

d dτ ∂x(τ) ∂t k = ∂ ∂t k dx(τ) dτ . (32)
It is thus possible to derive a dynamic equation for ∂x(τ) ∂t k by differentiation of the dynamics equation of the state vector x(τ) with respect to the commutation time t k . Hence the function τ → ∂x(τ) ∂t k is obtained by integrating the following dynamics equations:

d dτ ∂x(τ) ∂t 2l-1 =                                                                            (t 2k-1 -t 2k-2 )A(t(τ)) ∂x(τ) ∂t 2l-1 if τ ∈ [2k -2, 2k -1], k l (t 2k -t 2k-1 )A(t(τ)) ∂x(τ) ∂t 2l-1 if τ ∈ [2k -1, 2k], k l A(t(τ))x(τ) + D(t(τ)) + (t 2l-1 -t 2l-2 ) (τ -2l + 2) ∂A(t) ∂t x(τ) + ∂D(t) ∂t +A(t(τ)) ∂x(τ) ∂t 2l-1 if τ ∈ [2l -2, 2l -1], -[A(t(τ))x(τ) + D(t(τ)) + B(t(τ))U l ] + (t 2l -t 2l-1 ) (-τ + 2l + 2) ∂A(t) ∂t x(τ) + ∂D(t) ∂t + ∂B(t) ∂t U l +A(t(τ)) ∂x(τ) ∂t 2l-1 if τ ∈ [2l -1, 2l], (33)
and

d dτ ∂x(τ) ∂t 2l =                                                                            (t 2k-1 -t 2k-2 )A(t(τ)) ∂x(τ) ∂t 2l if τ ∈ [2k -2, 2k -1], k l (t 2k -t 2k-1 )A(t(τ)) ∂x(τ) ∂t 2l if τ ∈ [2k -1, 2k], k l -[A(t(τ))x(τ) + D(t(τ))] + (t 2l+1 -t 2l ) (-τ + 2l + 1) ∂A(t) ∂t x(τ) + ∂D(t) ∂t +A(t(τ)) ∂x(τ) ∂t 2l if τ ∈ [2l, 2l + 1], A(t(τ))x(τ) + D(t(τ)) + B(t(τ))U l + (t 2l -t 2l-1 ) (τ -2l + 1) ∂A(t) ∂t x(τ) + ∂D(t) ∂t + ∂B(t) ∂t U l +A(t(τ)) ∂x(τ) ∂t 2l if τ ∈ [2l -1, 2l]. (34)
The derivative of the penalisation function with respect to the state vector is:

∂ψ ∂x (t, x(τ)) = C 2 (t) [C 2 (t)x(t) -δ] 2 sign C 2 (t)x(t) -δ + 1 + C 2 (t) [C 2 (t)x(t) + δ] 2 sign -C 2 (t)x(t) -δ + 1 + C 3 (t) [C 3 (t)x(t) -δ] 2 sign C 3 (t)x(t) -δ + 1 + C 3 (t) [C 3 (t)x(t) + δ] 2 sign -C 3 (t)x(t) -δ + 1 , (35)
and its derivative with respect to time is:

∂ψ ∂x (t, x(τ)) = dC 2 (t) dt x [C 2 (t)x(t) -δ] 2 sign C 2 (t)x(t) -δ + 1 + dC 2 (t) dt x [C 2 (t)x(t) + δ] 2 sign -C 2 (t)x(t) -δ + 1 + dC 3 (t) dt x [C 3 (t)x(t) -δ] 2 sign C 3 (t)x(t) -δ + 1 + dC 3 (t) dt x [C 3 (t)x(t) + δ] 2 sign -C 3 (t)x(t) -δ + 1 , (36)
The SK problem can therefore be solved according to the following algorithm:

Numerical Results

The proposed commutation times optimisation technique is applied for the station keeping of a 4850 kg satellite equipped with the electric propulsion system presented in Section 2.2. Two cases are studied. For the first one, we were able to correctly tune the threshold parameter ς so that the trajectory obtained after the step 2 is SK-feasible and operationally feasible. The proposed technique optimises the commutation times in order to decrease the consumption of the Sk strategy.

For the second case, the parameter ς was harder to determine, resulting in a trajectory operationally feasible but not SK-feasible. The proposed technique manages to optimise the commutation times so that the ensuing trajectory is SK and operationally feasible while reducing also drastically the obtained consumption..

Commutation Time Optimisation for a SK-Feasible Trajectory

The two-step decomposition technique of Gazzino et al. [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] is used first in order to obtain a SK and operationally feasible trajectory. Solving the problem without the operational constraints (step 1) shows an optimal control profile for which the North and South thrusters fire simultaneously, which is also the optimal configuration according to the station keeping literature. Therefore, the step 2 is mandatory to get an operationally feasible control profile.. It was possible to choose ς = 0.8 for the CBE and ς = 0.9 for the EBE. The µ coefficient of Equation ( 28) is chosen to be 1.10 7 .

Figures 1 and2 depict the control profiles obtained by the proposed method starting respectively from the CBE scheme and from the EBE scheme. Figures 3 and4 exhibit the trajectories resulting from these control profiles in the latitude-longitude plane. The two equivalence schemes lead to different consumption values (see Table 1), and the proposed optimisation of the switching times leads to a substantial fuel saving. In general, a high value of the threshold parameter leads to a smaller fuel consumption. 1: Consumptions for the two-step and the three-step solutions.

Commutation Time Optimisation for a non SK-Feasible Trajectory

The CBE and EBE schemes enforce the operational constraints but do not tackle the SK constraints. Therefore the threshold parameter ς has to be precisely adjusted so that a SK-feasible trajectory can be computed at the end of step 2. Nevertheless, it may happen that the threshold parameter ς is so hard to determine that the computed middle times and thrust half durations lead to a non SK-feasible trajectory. The proposed technique based on the switched systems manages to optimise the commutation time so that the trajectory at the end of step 3 is both operationally feasible and SK-feasible.

Figure 6 depict the state trajectory at the end of step 2 with ς = 0.9 computed with the CBE scheme. In this case, the SK constraints are not respected. Despite this, the step 3 enforces the SK constraints while making sure that the operational constraints are still respected. Figure 5 displays the corresponding control profiles. Table 2 shows the decrease in consumption realised by step 3. Note that for a non SK-feasible trajectory computed with the EBE scheme, the proposed computation times optimisation technique manages as well to enforce the station keeping constraints. -: solution of the consumption based equivalence with threshold 0.9, --: solution of the switched systems optimisation.

Equivalence CBE scheme Two-step 12.37 m/s solution Threshold value 0.9 Three-step 0.31 m/s solution Table 2: Consumptions for the two-step and the three-step solutions in the case were the SK constraints are not respected at the end of step 2.

Conclusion

In this paper, an optimisation of the commutation times originating in the theory of the switched systems has been used in order to optimise the on-off firing sequence for the minimum-fuel station keeping of a satellite equipped with electric thrusters. In fact, the system can be naturally decomposed into several subsystems, one per thruster whose control is on and one corresponding to a coasting arc. As the proposed method requires to know beforehand a feasible firing sequence, a two-step decomposition technique has been used as an initial guess for the optimal sequence of firing arcs.. The optimisation of the switching times allow to overcome the drawbacks of the previous steps.

  Solving the Minimum-Fuel Low-Thrust Geostationary Station Keeping Problem via the Switching Systems Theory C. Gazzino † D. Arzelier , L. Cerri * , D. Losa ‡ , C. Louembet and C. Pittet * CNRS, LAAS-CNRS, Université de Toulouse 7 avenue du Colonel Roche, 31031 Toulouse cedex 4, France

Algorithm 1

 1 Optimise the commutation times {t k } Require: x(0) the initial position Ensure: the trajectory is operationally and SK-feasible {Step 1: simple SK problem} 1: Remove the operational constraints (a)-(d) from Problem 1 2: u i (t) i=1,...,4,t∈[0,T ] ← the optimal control for the SK-feasible trajectory computed by an indirect method initialised by a direct method (see Gazzino et al. 8 ) 3: ς ← threshold parameter to transform the continuous control profile in a rectangle profile 4: {t i,k , ∆t i,k } ← the middle times t i,k and the half durations ∆t i,k of the thrusts, {Step 2: equivalence schemes} 5: {t i,k , ∆t i,k } ← the middle times and half durations computed by the CBE or the EBE. {Step 3: commutation times optimisation with an interior point algorithm} 6: {t k } ← the commutation times computed from {t i,k , ∆t i,k } 7: {t (0) k } ← {t k } 8: {t k } ← the solution of Problem 3 with an interior point algorithm

Figure 1 :

 1 Figure 1: Control profiles: -: solution of the consumption based equivalence with threshold 0.8, --: solution of the switched systems optimisation.

Figure 2 :

 2 Figure 2: Control profiles: -: solution of the consumption based equivalence with threshold 0.8, --: solution of the switched systems optimisation.

Figure 5 :Figure 6 :

 56 Figure 5: Control profiles: -: solution of the consumption based equivalence with threshold 0.9, --: solution of the switched systems optimisation.

Table

	0.06
	0.04
	0.02
	Latitude ( •	0
	-0.02
	-0.04
	-0.06
		117.94	117.96	117.98	118	118.02	118.04
						Longitude ( • )
					Equivalence scheme	CBE	EBE
	0.06		Two-step solution	8.88 m/s 1.72 m/s
					Threshold value	0.8	0.9
	0.04		Three-step solution	0.45 m/s 0.31 m/s
	0.02
	Latitude ( •	0
	-0.02
	-0.04
	-0.06
	117.94	117.96	117.98	118	118.02	118.04	118.06
						Longitude ( • )

) Figure 3: Latitude vs. longitude: -: solution of the consumption based equivalence with threshold 0.8, --: solution of the switched systems optimisation.

)

Figure 4: Latitude vs. longitude: -: solution of the effect based equivalence with threshold 0.9, --: solution of the switched systems optimisation.


Baisong Chen

Zhuo Fu

Yuchun Pan

Jihua Wang
email: wangjh@nercita.org.cn

Zhixuan Zeng

Single Leaf Area Measurement Using Digital Camera Image

Keywords: Digital camera, Leaf area, Area measurement, Image processing

Leaf area index is one of the most important parameters in ecological and environmental studies. This paper presents a method for single leaf area measurement based on the counting of the leaf pixels in digital leaf image. Initially, the target leaf is put on a piece of white paper on which four endpoints of an equal-length and orthogonal cross are printed. A background color threshold is then set to separate the white background pixels from the leaf pixels. After the removal of the background pixels, the leaf pixels are left; and the ratio of the number of leaf pixels to the number of pixels of the reference facet formed by the four control endpoints is equal to that of their areas. The single leaf area can then be computed easily based on this proportional relation. Analysis and experimental results indicate that the proposed method is an efficient and precise method for single leaf or leaf-like object area measurement.

Introduction

Leaf is a key functional organ of photosynthesis and transpiration, the size of left area impacts greatly on the physiological function of the plants. LAI (Leaf Area Index) is also an important input parameter of some plant, ecology and environment models [START_REF] Chen | etal: Derivation and Validation of Canada-Wide Coarse-Resolution Leaf Area Index Maps Using High-Resolution Satellite Imagery and Ground Measurements[END_REF][START_REF] Jongschaap | Run-Time Calibration of Simulation Models by Integrating Remote Sensing Estimates of Leaf Area Index and Canopy Nitrogen[END_REF], thus fast and accurate estimation of leaf area is significant. The currently available approaches for single leaf area measurement are mainly including: using leaf area measuring tools, grids approximation, weighing [START_REF] Diao | Single Leaf Area Estimation Models Based on Leaf Weight of Eucalyptus in Southern China[END_REF], using planimeter [START_REF] Gist | An Optical Planimeter for Leaf Area Determination[END_REF], displacement volume computation [START_REF] Chen | Leaf Area Index of Boreal Forests: Theory, Techniques, and Measurements[END_REF], regression models [START_REF] Williams | Nondestructive Leaf Area Estimation of 'Niagara' and 'Dechaunac' Grapevines[END_REF][START_REF] Vertessy | Relationships between Stem Diameter, Sapwood Area, Leaf Area and Transpiration in a Young Mountain Ash Forest[END_REF][START_REF] Kumar | Calibration and Validation of Regression Model for Non-Destructive Leaf Area Estimation of Saffron (Crocus Sativus L.)[END_REF], and image processing [START_REF] Xiao | A Simple Non-Destructive Method to Measure Leaf Area Using Digital Camera and Photoshop Software[END_REF][START_REF] Caldas | Measurement of Leaf Area With a Hand-Scanner Linked to a Microcomputer[END_REF][START_REF] Li | Extraction of Leaf Area Index of Wheat Based on Image Processing Technique[END_REF], and so on. Among them, using leaf area measuring tools or planimeter is highly accurate, however, with high cost of equipments as well. And other methods such as the grids estimation, weighing, displacement computation and regression models are available for different kinds of leaves and low in cost, but their accuracy is low too. Considering the cost of tools, applicability, accuracy, efficiency, and the convenience for outdoor measurement, the image processing method can be a considerable method. Based on this, this paper proposes a method for computing single leaf area by separating the background and counting the leaf pixels.

2

Materials and method

Leaves for experiment and their areas

Thirty pieces of Japan Euonymus leaves with different shapes were selected as experimental materials. They were numbered and put on a 10×10cm white referent paper, respectively. And then they were scanned with a high resolution, respectively. The number of pixels N B in the entire 10×10cm region was then counted in Photoshop, and the number of non-background pixels N L were then computed by removing the white background pixels. Because the scanner images objects in a progressive and orthographic projection manner, the proportion of the area of leaf S L and the area of background S B will equate to that of the number of leaf pixels and the number of background pixels. (

Using equation ( 1), we can compute the actual area of every leaf respectively, and their areas are then considered as true values in the validation section.

Method and procedures

The principle of this proposed method is the same as that of the above method for measuring leaf area by a scanner, that is, using the proportions of the area and the number of pixels. For the convenience of taking pictures on field, we use a common digital camera rather then a scanner. And for the convenience of program processing and the insurance of accuracy, the referent paper is replaced by four endpoints of an equal-length and orthogonal cross printed on a piece of paper. So the main procedures are as follows:

1) Four endpoints of an equal-length (10×10cm) and orthogonal cross are initially printed on the middle of a piece of common A4 paper. Leaves are then put onto this paper respectively, unwrapped with its centre lapping on that of the cross, leaving four endpoints uncovered. If the leaf cannot be unwrapped naturally, a cut or clip is then needed.

2) Choose suitable settings and parameters for the camera, and then take pictures of leaves respectively, bearing in mind that the primary optical axis should be perpendicular to the background paper with its focus aiming at the centre of the leaf. This ensures the correctness of above equation [START_REF] Chen | etal: Derivation and Validation of Canada-Wide Coarse-Resolution Leaf Area Index Maps Using High-Resolution Satellite Imagery and Ground Measurements[END_REF].

3) After the above image acquisition, we obtain images with white background, green leaf, and four black dots. Initially, we find the coordinates of the four endpoints on an image. And from these coordinates, we can compute the distances and included angle of the four endpoints. Because we had 4 endpoints of an equal-length and orthogonal cross printed on the paper, then their corresponding points in an image should also be equal-length and orthogonal, if the positions of camera and leaf are right when the picture is taken. Thus, the difference of the in-image distance and the actual distance should not be larger than a threshold, say, 1% of the longer distance. Also, the difference of the included angle should be less than one degree. If they do not meet these conditions, then the positions of the camera and the leaf may be wrong during taking picture, for example, the camera would be in a gradient position. If this is the case, then, equation ( 1) would no longer hold in such images, and thus they need to have picture retaken. If the image is available, a further operation can then be performed to remove the background. The background is white (RGB: 255, 255, 255), but it will be slightly changed due to the color error of the CCD, thus we can set a color threshold, say RGB greater than 230, to filter those background pixels. After the background pixels are counted, there are pixels of leaf and those four endpoints left. And the four points can be very tiny, so we can ignore their area. This means these non-background pixels can be considered as leaf pixels. Thus, we can compute the numbers of background and leaf pixels, and the area of background (10×10cm). Finally, the area of the interested leaf can be computed via equation (1).

Flowchart for the image processing

Based on the above descriptions, we can conclude the main image processing procedures as: 1) coordinates acquisition, 2) distances and included angle computation, 3) error detection, 4) separation of background, 5) left pixels count, 6) leaf area computation. And then the flowchart for programming could be as follows.

Fig. 1. The flowchart for the image processing of single leaf area measurement

Here are some details need to notice during programming: 1) after four endpoints are picked, two distances can be computed, then the included angle can be calculated by: cosα = (vector1 • vector2) / (|vector1| * |vector2|). Again, the differences of distances and included angle should not be exceeding, otherwise the program should give an alert and terminate; 2) the RGB threshold for background separation can be different because of different lighting situations and camera settings, based on many experiments, we recommend RGB (230, 230, 230) as a threshold, pixels with R, G, and B values all greater than this can be considered as background; 3) if the differences of distance and included angle are not exceeding, then the average number of pixels should be used to compute the number of pixels of the 10×10cm square; 4) 1) is just one side area of the leaf area, if area of double sides is wanted, a multiplication with two is then needed.

Results and discussions

For the convenience of program processing, based on the above method, we developed a program within VC 6.0 and its MFC and an open JPEG lib called JpegLib. To test the accuracy and efficiency of this method, all the selected leaves were measured by this program, then, we computed the RMSE for every single leaf via the following formula:

å = - = n i t m S S n RMSE 1 2 ) ( 1 . ( 2
)
Where n is total number of leaves, S m is the area of a numbered leaf measured by program, S t is its corresponding true value of area. To find a best resolution for this type of leaf, we took pictures of these leaves with different resolutions, and then we computed the RMSE and recorded the average processing time for every image, results are showed in the Table 1 below. As indicated by Table 1, the RMSE remains the same when the resolution is higher than 1600×1200; instead, higher resolution image demands larger storage space and longer processing time. A Jpeg decoding is needed in every time when the program reads the image file, also, a full traverse of all pixels is needed when the program tries to count the amounts of the background and leaf pixels. Thus, there is no benefit to using a resolution higher than 1600×1200 for leaves like Japan Euonymus's. Therefore, the best image-taking resolution is 1600×1200.

To test the stability of this method, we took 20 pictures for a single leaf using the best resolution, and then we computed the area for every picture, and based on the results, we got the following basic statistics in Table 2.

As shown in Table 2, the coefficient of variation and mean relative error are both small, indicating that the areas are mostly approximate, and this method is highly accurate and stable.

The principle of this proposed method is quite similar to that of the reference [START_REF] Xiao | A Simple Non-Destructive Method to Measure Leaf Area Using Digital Camera and Photoshop Software[END_REF]. To compare these two methods, we took one picture for every selected leaf at the best resolution, and then we computed areas using both methods, respectively. The results are shown in Fig. 2 below.

Fig. 2. The comparison of measurements of 30 pieces of leaves using two different methods

As shown in Fig. 2, the results of the two methods are quite similar, with high correlation too. But, there is still an improvement in accuracy in our method. Because of the virtual equal-length and orthogonal cross, the program can determine if the positions of camera and the leaf are correct when the picture is taken. This ensures the correctness of equation ( 1) and thus avoids the gross error and improves the accuracy.

Our method distinguishes the leaf pixels from the background ones by means of their color difference. This brings additional benefits such as: automatic removal of the holes and the sear or ill spots. This method can also be used to measure areas of some leaf-like slices such as paper money and so on.

Conclusion

This paper proposes a method for single leaf area measurement by means of counting of the background and leaf pixels. In this method, a virtual equal-length and orthogonal cross is adopted to avoid the potential gross error and improve the accuracy of leaf area measurement. Experiments indicate that the best image resolution for moderate or small leaves is 1600×1200. This method is also available for the area measurement for any leaf-like slice.

  by equation (

Table 1 .

 1 The RMSE and processing time for single leaf area calculation at different image

	resolutions
	Resolution	RMSE (cm 2 )	Average processing time (ms)
	640×480	0.07	23
	1024×768	0.04	78
	1600×1200	0.02	221
	2048×1536	0.02	617
	2272×1704	0.02	1387
	Configurations: CPU Pentium D3.2G, 2G DDR2 800RAM, WinXP;
	The best resolution would vary for different objects.

Table 2 .

 2 The statistical results of 20 area measurements for a single leaf at resolution of 1600×1200

	True area(cm 2 )	Min(cm 2 ) Max (cm 2 )	Mean(cm 2 )	Standard deviation(cm 2 )	Coefficient of variation	Mean relative error
	5.26	5.16	5.35	5.27	0.06	1.1%	0.2%
	Note: Mean relative error = |measured area -true area| / true area × 100%.

Acknowledgments. This work was supported by the National Key Technology R&D Program of China (No. 2008BAJ08 B03). The authors would also like to thank the anonymous reviewers for their constructive comments and suggestions.