N Kacem

COMPLEX DYNAMICS AND MULTISTABILITY IN NONLINEAR RESONANT NANOSYSTEMS BEYOND THE DUFFING CRITICAL AMPLITUDE

Keywords: Nonlinear Dynamics, Nanotechnology, Multistability, Complex Dynamics

Analytical multi-physics models which include main sources of nonlinearities for nanoresonators electrostatically actuated are developed in order to assess complex dynamics in nanosystems beyond the Duffing critical amplitude. In particular, multistability is investigated for doubly clamped beams and cantilevers. The bifurcation topology of a particular multistable behavior (up to five amplitudes for a given frequency) is parametrically identified and experimentally validated.

INTRODUCTION

Micro and nanoelectromechanical (MEMS/NEMS) devices have been the subject of extensive research for a number of years and have generated much excitement as their use in commercial applications has increased. Indeed, MEMS technology has opened up a wide variety of potential applications not only in the inertial measurement sector, but also spanning areas such as communications (filters, relays, oscillators, LC passives, optical switches), biomedicine (point-of-care medical instrumentation, microarrays for DNA detection and high throughput screening of drug targets, immunoassays, in-vitro characterization of molecular interactions), computer peripherals (memory, new I/O interfaces, read-write heads for magnetic disks) and other miscellaneous areas such as in projection displays, gas detection and mass flow detection.

NEMS are the natural successor to MEMS as the size of the devices is scaled down to the submicron domain. This transition is well adapted with the resonant sensing technique for a large panel of applications. One reason of down scaling resonant sensors to the NEMS size is the ability to detect very small physical quantities by increasing their sensitivity [START_REF] Ekinci | Nanoelectromechanical systems[END_REF]. In particular, NEMS have been proposed for use in ultrasensitive mass detection [START_REF] Ekinci | Ultrasensitive nanoelectromechanical mass detection[END_REF][START_REF] Jensen | An atomic-resolution nanomechanical mass sensor[END_REF], radio frequency (RF) signal processing [START_REF] Nguyen | Micromechanical components for miniaturized low-power communications[END_REF][START_REF] -C. Nguyen | Tunable, switchable, high-q vhf microelectromechanical bandpass filters[END_REF], and as a model system for exploring quantum phenomena in macroscopic systems [START_REF] Cho | Physics -researchers race to put the quantum into mechanics[END_REF][START_REF] Lahaye | Approaching the quantum limit of a nanomechanical resonator[END_REF].

Unfortunately, the nonlinear regime for nanomechanical resonators is easily reachable, so that the useful linear dynamic range of the smallest NEMS devices is severely limited. In fact, many applications we are hoping for in the near future will involve operation in the nonlinear regime, where the response to the stimulus is suppressed and frequency is pulled away from the original resonant frequency. Actually, it is a challenge to achieve large-amplitude motion of NEMS resonators without deteriorating their frequency stability [START_REF] Feng | Phase noise and frequency stability of very-high frequency silicon nanowire nanomechanical resonators[END_REF]. The relative frequency noise spectral density [START_REF] Robins | Phase Noise in Signal Sources[END_REF] of a NEMS resonator is given by:

S f = 1 2Q 2 S x P 0 (1)
where S x is the displacement spectral density and P 0 is the displacement carrier power, ie the RMS drive amplitude of the resonator 1 2 A 2 . Remarkably, driving the resonator at large oscillation amplitude leads to better SNR and, thus, simplifies the design of the electronic feedback loop. However, doing so in the nonlinear regime reduces the sensor performances since the frequency instability of a nonlinear resonator is proportional to its oscillation amplitude. Moreover, even when NEMS resonators are used as oscillators in closed-loop, a large part of noise mixing [START_REF] Roessig | Nonlinear mixing in surface-micromachined tuning fork oscillators[END_REF][START_REF] Kaajakari | Phase noise in capacitively coupled micromechanical oscillators[END_REF] due to nonlinearities drastically reduces their dynamic range and alters their detection limit.

This paper is about investigating such limitations beyond the Duffing critical amplitude [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF][START_REF] Kacem | Dynamic range enhancement of nonlinear nanomechanical resonant cantilevers for highly sensitive nems gas/mass sensor applications[END_REF], based on the nonlinear dynamics of nanomechanical resonators. In particular, multistability is tracked analytically and experimentally for clamped-clamped nanobeams and nanocantilevers showing a complex multistable dynamics with up to five amplitudes for a given frequency.

MODELS

Clamped-clamped beam resonators

A clamped-clamped nanobeam is considered (Figure 1) subject to a viscous damping with coefficient c per unit length and actuated by an electric load v(t) = V dc + V ac cos( Ωt ), where V dc is the DC polarization voltage, V ac is the amplitude of the applied AC voltage, and Ω is the excitation frequency.

Equation of motion

The equation of motion that governs the transverse deflection w(x, t) is written as:

EI ∂ 4 w(x, t) ∂ x4 + ρbh ∂ 2 w(x, t) ∂ t2 + c ∂ w(x, t) ∂ t -   Ñ + Ebh 2l l 0 ∂ w(x, t) ∂ x 2 dx   ∂ 2 w(x, t) ∂ x2 = 1 2 ε 0 b V dc + V ac cos( Ωt ) 2 (g -w(x, t)) 2 (2)
where x is the position along the nanobeam length, E and I are the Young's modulus and moment of inertia of the cross section. Ñ is the applied tensile axial force due to the residual stress on the silicon or the effect of the measurand, t is time, ρ is the material density, h is the nanobeam thickness, g is the capacitor gap width, and ε 0 is the dielectric constant of the gap medium. The last term in Equation ( 2) represents an approximation of the electric force assuming a complete overlap of the area of the nanobeam and the stationary electrode. The boundary conditions of Equation ( 2) are:

w(0, t) = w(l, t) = ∂ w ∂ x (0, t) = ∂ w ∂ x (l, t) = 0 (3)

Normalization

For convenience and equations simplicity, we introduce the nondimensional variables:

w = w g , x = x l , t = t τ (4)
where τ = 2l 2 h 3ρ E . The substitution of Equation (4) into Equations ( 2) and (3) leads to:

∂ 4 w ∂x 4 + ∂ 2 w ∂t 2 + c ∂w ∂t -   N + α 1 1 0 ∂w ∂x 2 dx   ∂ 2 w ∂x 2 = α 2 [V dc + V ac cos(Ωt)] 2 (1 -w) 2 (5)
w(0, t) = w(1, t) = ∂w ∂x (0, t) = ∂w ∂x (1, t) = 0 (6)
The parameters appearing in Equation ( 5) are:

c = cl 4 EIτ N = Ñ l 2 EI α 1 = 6 g h 2 α 2 = 6 ε 0 l 4 Eh 3 g 3 Ω = Ωτ (7)

Solving

The total beam displacement w(x, t) is a sum of a static dc displacement w s (x) and a timevarying ac displacement w d (x, t). However, for our devices, it is easy to check that the static deflexion is negligible. Typically, the measured quality factors Q are in the range of 10 4 -5.10 4 and the V dc ≤ 200V ac. Thus, the ratio between the static and the dynamic deflection is:

w s (x) w d (x, t) ≈ V dc 2Q.V ac ≤ 1% (8)
A reduced-order model is generated by modal decomposition transforming Equation ( 5) into a finite-degree-of-freedom system consisting of ordinary differential equations in time. We use the undamped linear mode shapes of the straight nanobeam as basis functions in the Galerkin procedure. To this end, we express the deflection as :

w(x, t) = Nm k=1 a k (t)φ k (x) (9)
where N m is the number of modes retained in the solution, a k (t) is the k th generalized coordinate and φ k (x) is the k th linear undamped mode shape of the straight nanobeam, normalized such that

1 0 φ k φ j = δ kj where δ kj = 0 if k = j and δ kj = 1 if k = j.
The modal basis is formed of the eigenmodes of a linear undamped straight nanobeam. The latter are the solutions of the following equation:

d 4 φ k (x) dx 4 = λ 4 k φ k (x) (10)
We multiply Equation ( 5) by φ k (x) (1 -w) 2 in order to include the complete contribution of the nonlinear electrostatic forces in the resonator dynamics without approximation. Next, we substitute Equation ( 9) into the resulting equation, use Equation [START_REF] Roessig | Nonlinear mixing in surface-micromachined tuning fork oscillators[END_REF] to eliminate d 4 φ k (x) dx 4 , integrate the outcome from x = 0 to 1, and obtain a system of coupled ordinary differential equations in time.

Since the first mode is the dominant mode of the system and the higher modes are negligible, it is enough to consider n = 1 [START_REF] Kacem | Computational and quasi-analytical models for non-linear vibrations of resonant mems and nems sensors[END_REF] and obtain: ä1 + (500.564 + 12.3N )a 1 + (927

+ 28N + 151α 1 ) a 1 3 + 347α 1 a 1 5 +(1330.9 + 38.3N )a 1 2 + 471α 1 a 1 4 + 2.66c 1 a 1 ȧ1 + 1.85c 1 a 1 2 ȧ1 +c 1 ȧ1 + 2.66a 1 ä1 + 1.85a 1 2 ä1 = - 8 3π α 2 [V dc + V ac cos(Ωt)] 2 (11)
In addition, since near-resonant behavior is the principal operating regime of the proposed system, a detuning parameter, σ is introduced, as given by:

Ω = ω 1 + εσ (12)
Separating the resulting equations and averaging them over the period 2π Ω in the t-domain results in the system's averaged equations, in terms of amplitude and phase, which are given by:

Ȧ = - 1 2 ǫξ 0 A - 1 8 ǫξ 2 A 3 + 1 2 ǫ κ ω 1 sin β + O(ε 2 ) (13)
A β = -Aσǫ + 3 8 ǫ χ 3 ω 1 A 3 + 5 16 ǫ χ 5 ω 1 A 5 - 7 10 ǫω 1 A 3 + 1 2 ǫ κ ω 1 cos β + O(ε 2 ) ( 14
)
where ω 1 = √ 500.564 + 12.3N , ξ 0 = c 1 , ξ 2 = 1.85c 1 , χ 3 = 927 + 28N + 151α 1 , χ 5 = 347α 1 and κ = 16 3π α 2 V acV dc.

Cantilevers

In order to develop a model for nanocantilever beams, a slender uniform flexible beam is considered as shown in Figure 2. The beam is initially straight and it is clamped at one end and free at the other end, subject to viscous damping with a coefficient c per unit length and actuated by an electric load v(t) = V dc + V ac cos( Ωt ), where V dc is the DC polarization voltage, V ac is the amplitude of the applied AC voltage, t is time and Ω is the excitation frequency. In addition, the beam follows the Euler-Bernoulli beam theory, where shear deformation and rotary inertia terms are negligible. In order to derive the nonlinear equation of motion describing the flexural vibration of a cantilever beam electrostatically actuated, we follow a variational approach, based on the extended Hamilton principle and used by Crespo da Silva and Glynn [START_REF] Crespo Da Silva | Nonlinear flexural-flexural-torsional dynamics of inextensional beams. i: Equations of motion[END_REF][START_REF] Crespo Da Silva | Nonlinear flexural-flexural-torsional dynamics of inextensional beams. ii: Forced motions[END_REF] and Crespo da Silva [START_REF] Crespo Da | Non-linear flexural-flexural-torsional-extensional dynamics of beams. i: Formulation[END_REF][START_REF] Crespo Da | Non-linear flexural-flexural-torsional-extensional dynamics of beams. ii: Response analysis[END_REF].

EI

w′′′′ + w′ w′ w′′

′ ′ + ρbh ẅ + c ˙w = - 1 2 ρbh w′ x l ∂ 2 ∂ t2 s 1 0 ( w′ ) 2 ds 2 ds 1 ′ + 1 2 ε 0 b V dc + V ac cos( Ωt ) 2 (g -w) 2 (15)
where primes and dots denote respectively the partial differentiation with respect to the arclength x and to the time t. w is the beam bending deflection, E and I are the Young's modulus and geometrical moment of inertia of the cross section. l and b are the length and width of the nanobeam, ρ is the material density, h is the nanobeam thickness in the direction of vibration, g is the capacitor gap width, and ε 0 is the dielectric constant of the gap medium. The first term in the left-hand side of Equation ( 15) is due to the nonlinear expression for the curvature of the beam, while the first term in the right-hand side, which involves a double time derivative, is the nonlinear inertial term. The boundary conditions are:

w(0, t) = w′ (0, t) = w′′ (l, t) = w′′′ (l, t) = 0 (16)
Substituting Equation (4) into Equations ( 15) and ( 16), we obtain:

w iv + ẅ + c ẇ + δ 1 w ′ w ′ w ′′ ′ ′ = -δ 2 w ′ x 1 ∂ 2 ∂t 2 x 1 0 (w ′ ) 2 dx 2 dx 1 ′ + δ 3 V dc V ac 1 + Vac V dc cos(Ωt) 2 (1 -w) 2 (17)
w(0, t) = w ′ (0, t) = w ′′ (1, t) = w ′′′ (1, t) = 0 ( 18
)
where primes and dots denote respectively the partial differentiation with respect to the dimensionless arclength x and to the dimensionless time t. The parameters appearing in Equations ( 17) are:

c = cl 4 EIτ , δ 1 = g l 2 , δ 2 = 1 2 g l 2 δ 3 = 6V ac V dc ε 0 l 4 Eh 3 g 3 , Ω = Ωτ (19)
The electrostatic force in Equation ( 17) is expanded in a fifth order Taylor series, in order to capture 5 possible amplitudes for a given frequency in the mixed behavior [START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF]. Then, Equation ( 9) is substituted into the resulting equation, Equation ( 10) is used to eliminate d 4 φ k (x) dx 4 , and the outcome is multiplied by φ k and integrated from x = 0 to 1 for k ∈ [1, n] ∩ N . Thus, a system of coupled ordinary differential equations in time is obtained.

In addition, since near-resonant behavior is the principal operating regime of the proposed system, a detuning parameter σ is introduced, as given by:

Ω = ω 1 + εσ ( 20
)
where

ω 1 = λ 2 1 -Vac V dc δ 3 -2 V dc
Vac δ 3 and ε is the small nondimensional bookkeeping parameter. Separating the resulting equations and averaging them over the period 2π Ω in the t-domain results in the system's averaged equations in terms of amplitude A and phase β given by:

Ȧ = εδ 3 sin β ω 1 0.78 + 1.11A 2 + 2.44A 4 - εc 2 A + O(ε 2 ) (21) β = εσ -ε δ 3 cos β ω 1 0.78 A + 3.32A + 12.19A 3 + ε A 4 δ 3 ω 1 12.51 V dc V ac -6.25 V ac V dc -ε 15.16A 2 δ 1 ω 1 + εδ 3 A 2 ω 1 1.76 V ac V dc + 3.52 V dc V ac + 2.3εA 2 δ 2 ω 1 + O(ε 2 ) (22)
The steady-state motions occur when Ȧ = β = 0, which corresponds to the singular points of Equations ( 13) and ( 14) or ( 21) and ( 22). Thus, the frequency response equation can be written in its parametric form A = K 1 (β), Ω = K 2 (β) as a function of the phase β. This set of two equations is easily implementable in Matlab or Mathematica. This ability makes the model suitable for NEMS designers as a quick tool to optimize the resonant sensors performance. The plots of Figure 3 were carried out for a resonant nano-cantilever with the following set of parameters: l = 12.5µm, h = 300nm, b = 500nm and V dc = 50V ac . g and V ac were used for parametric studies. This analytical model enables the capture of all the nonlinear regimes in the resonator dynamics and describes the competition between the mechanical hardening and the electrostatic softening behaviors.

Practically, when h/g << 1 and for a high quality factor, the dynamics is dominated by the hardening nonlinearities and in the opposite case, the frequency response is nonlinearly softening. Between these two configurations (h/g = 0.15 for instance), a mixed hardeningsoftening behavior is inescapable. This leads to a multistability with five possible amplitudes for a given frequency. It is highly unstable and difficult to control for NEMS designers.

EXPERIMENTAL VALIDATION

The device consists in a silicon doubly clamped beam electrostatically in plane actuated and detected using two electrodes which allows for 2 port electric measurements. The resonator is 200 µm long, 4 µm thick, 2 µm wide, the actuation gap is 1 µm, and the detection gap 750 nm. The device was placed in a vacuum chamber, and the electrical measurements were performed at room temperature using a low noise lock-in amplifier. The drive voltage is V ac = 0.6 V and the beam is polarized with V dc = 10 V . Figure 4 shows the frequency response of the device, with up-and down-sweeps. The quality factor obtained with this polarization voltage and in a linear regime is 4000. The Duffing critical amplitude is then A c = 53 nm, i.e. V c = 25 µV . The peak obtained is then far beyond A c , up to 70% of the gap. The frequency response shows 4 bifurcation points noted B i , at which jumps J i occur to destination points D i on stable branches, according to the direction of the sweep : as the frequency is swept up from f 0 , the output voltage follows the path labelled f 0 -B 4 -D 4 -B 2 -D 2 -f 1 , and as it is swept down from f 1 , the path f 1 -B 1 -D 1 -B 3 -D 3 -f 0 is followed. In the presence of the 3 other bifurcation points, the point B 4 may be called the mixed behavior initiation point. It is highly unstable: it is located at relatively high amplitude (i.e. in a state of high potential energy) as opposed to point B 1 or as opposed to a typical softening behavior, and the state variables have to jump to a destination stable branch at even higher amplitude.

The parametric analytical frequency response is superposed to the experimental points in Figure 4. The model shows an excellent agreement: the 4 bifurcation points exist and are well located, and the stable branches coincide very well with the measurement. This confirms the performance and the accuracy of the model at high amplitude in the nonlinear regime.

CONCLUSIONS

Analytical models have been developed to quantitatively assess the nonlinear dynamics of NEMS resonators. These models include the main sources of nonlinearities (electrostatic, geometric and inertial) and they are based on the modal decomposition (the Galerkin discretization) combined with a perturbation technique (the averaging method). The shape of the model output (two parametric equations) has the advantage to be simple and easy to implement for resonant NEMS designers. Particularly, the importance of the fifth order nonlinearities has been demonstrated through the analytical as well as experimental identification of the mixed behavior.

Figure 1 :

 1 Figure 1: Schema of an electrostatically actuated clamped-clamped nanobeam

Figure 2 :

 2 Figure 2: Schema of an electrostatically actuated nanocantilever

Figure 3 :

 3 Figure 3: Analytical forced frequency responses for Q = 10 4 and several values of g and V ac . W max is the beam displacement at its free end normalized by the gap g. The different bifurcation points are {B 1 , B 2 , B 3 , B 4 }, the point B 4 characterizes the initiation of the mixed behavior.

Figure 4 :

 4 Figure 4: Analytical and experimental frequency curves showing a mixed behavior and the followed paths respectively in a sweep up frequency f 0-B 4 -D 4 -B 2 -D 2 -f 1 and a sweep down frequency f 1 -B 1 -D 1 -B 3 -D 3 -f 0 . {J 1 , J 2 , J 3 , J 4 }are the four jumps cauterizing a typical mixed behavior of NEMS resonators, {B 1 , B 2 , B 3 , B 4 } are the different bifurcation points and {D 1 , D 2 , D 3 , D 4 } are the destination points after jumps. The two branches [B 1 , B 2 ] and [B 3 , B 4 ] in dashed lines are unstable.

N. Kacem


Mathieu Huard
email: mathieu.huard@cea.fr

Rida T Farouki
email: farouki@ucdavis.edu

Nathalie Sprinsky

Luc Biard
email: luc.biard@imag.fr

Nathalie Sprynski
email: nathalie.sprynski@cea.fr

C 2 interpolation of spatial data subject to arc-length constraints using Pythagorean-hodograph quintic splines

published or not. The documents may come

C2 interpolation of spatial data subject to arc-length constraints using Pythagorean-hodograph quintic splines

Introduction

The Leti1 department of the CEA2 in Grenoble, France has recently developed the Morphosense, a flexible ribbon-like device instrumented with micro-accelerometers and micro-magnetometers according to a known distribution along its length [START_REF] David | Method and device for acquisition of a geometric shape[END_REF]. Each sensor knot provides measurements that can be used to completely determine the spatial orientation of the knot. From this orientation data, and the known distances between consecutive sensors, it is desired to reconstruct the geometrical shape assumed by the Morphosense when placed in contact with a physical object (see Figure 1).

It should be emphasized that the Morphosense shape recovery is not simply a classical Hermite or spline interpolation problem, since the absolute (or even relative) spatial locations of the sensor knots are not known. Instead, only the curvilinear distances between successive sensors are known. Thus, the shape recovery problem amounts to constructing a spline curve that interpolates a given sequence of Serret-Frenet frames, with precribed arc lengths for each segment between such frames. This problem has been previously studied [START_REF] Huard | Serret-Frenet frame interpolation under length constraints for surface reconstruction[END_REF][START_REF] Huard | Reconstruction of quasi developable surfaces from ribbon curves[END_REF][START_REF] Sprynski | Reconstruction de courbes et surfaces à partir de données tangentielles[END_REF] assuming an arc-length parameterization, which allows control over the curve segment lengths between sensors. The proposed methods aim to interpolate the tangent indicatrix of the curve on the unit sphere. Since they rely on numerical integration, however, such methods do not yield analytic definitions of the reconstructed curve.

A new global approach to the shape reconstruction is presented herein, based on the polynomial Pythagorean-hodograph (PH) curves [START_REF] Farouki | Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable[END_REF]. The proposed solutions are spatial C 2 PH quintic splines, with segments that match the arc-length constraints. Specifically, two spatial C 2 PH quintic spline interpolation schemes, subject to segment arc-length constraints, are developed. The first scheme interpolates a given sequence of points in R3 , while the second interpolates given derivative vectors (providing an approximate solution to the Morphosense shape reconstruction problem).

The distinctive feature of polynomial PH curves is that their parametric speed, i.e., the derivative of arc length with respect to the curve parameter) is a polynomial -rather than the square root of a polynomial -in the curve parameter. This special algebraic structure of PH curves offers many useful computational advantages over "ordinary" polynomial parametric curves. For example, the PH curves admit rational offsets and exact arc-length computation, and are well-suited to real-time CNC interpolator algorithms, used to drive computer-controlled machines along curved paths with a prescribed speed variation. A comprehensive treatment may be found in [START_REF] Farouki | Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable[END_REF].

The construction of planar C 2 PH spline interpolants involves the solution of a system of N quadratic equations in N complex unknowns [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF][START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF]. In the construction of spatial C 2 PH splines, however, there are more freedoms 3 than are required to satisfy the interpolation conditions [START_REF] Farouki | Spatial C 2 PH quintic splines[END_REF][START_REF] Farouki | Design of C 2 spatial Pythagorean-hodograph quintic splines by control polygons[END_REF]. A method was presented in [START_REF] Farouki | Design of C 2 spatial Pythagorean-hodograph quintic splines by control polygons[END_REF] to fix the residual freedoms by optimizing a shape measure for the interpolant. In the present context, the residual freedoms will be exploited to achieve satisfaction of the prescribed arc-length constraints for successive spline segments, a physically and geometrically important aspect of the Morphosense shape reconstruction problem.

The remainder of this paper is organized as follows. Section 2 first presents some preliminaries concerning the definition and properties of spatial PH curves. The system of equations that defines a spatial C 2 PH quintic spline is then introduced in Section 3. The interpolation of a set of points under segment arc-length constraints is discussed in Section 4, and the interpolation of derivatives under segment arc-length constraints is addressed in Section 5. The methods are illustrated using computed examples in both cases. Finally, Section 6 summarizes the capabilities and limitations of these methods, and briefly discusses possible enhancements and extensions.

Spatial Pythagorean-hodograph curves

A polynomial PH curve r(t) = (x(t), y(t), z(t)) in R 3 is characterized by the fact that its hodograph or derivative r ′ (t) = (x ′ (t), y ′ (t), z ′ (t)) satisfies the Pythagorean relation

x ′2 (t) 2 + y ′2 (t) + z ′2 (t) = σ 2 (t) (1)
for some polynomial σ(t). The polynomial σ(t) = |r ′ (t)| is the parametric speed of the curve -i.e., the rate of change ds/dt of its arc length s with respect to the curve parameter t. The Pythagorean condition ( 1) is satisfied if and only if the hodograph components can be expressed in terms of four polynomials u(t), v(t), p(t), q(t) in the form

x ′ (t) := u 2 (t) + v 2 (t) -p 2 (t) -q 2 (t) , y ′ (t) := 2 [ u(t)q(t) + v(t)p(t) ] , ( 2
)
z ′ (t) := 2 [ v(t)q(t) -u(t)p(t) ] ,
and the parametric speed is then

σ(t) := u 2 (t) + v 2 (t) + p 2 (t) + q 2 (t) . (3)
The quaternion and Hopf map forms [START_REF] Choi | Clifford algebra, spin representation, and rational parameterization of curves and surfaces[END_REF] are two alternative (equivalent) formulations of spatial PH curves that embody the structure (1)- [START_REF] David | Method and device for acquisition of a geometric shape[END_REF]. The Hopf map form will be adopted here. This generates the form (2) from a pair of complex 4 polynomials

α(t) := u(t) + i v(t) , β(t) := q(t) + i p(t) (4)
through the expression

r ′ (t) := ( |α(t)| 2 -|β(t)| 2 , 2 Re(α(t)β(t)) , 2 Im(α(t)β(t) ) (5)
(where z denotes the complex conjugate of z), and the parametric speed is then

σ(t) := |α(t)| 2 + |β(t)| 2 . ( 6
)
The simplest non-trivial spatial PH curves are the cubics, generated by choosing linear complex polynomials for (4), substituting in (5), and integrating. Since the spatial PH cubics have insufficient shape flexibility for C 2 spline interpolation, however, we employ PH quintics. These are defined in terms of two quadratic complex polynomials, specified in Bernstein form as

α(t) := α 0 (1 -t) 2 + α 1 2(1 -t)t + α 2 t 2 , β(t) := β 0 (1 -t) 2 + β 1 2(1 -t)t + β 2 t 2 .
Writing a point (x, y, z) ∈ R 3 in the form (x, y + i z), the control points q k = (x k , y k + i z k ) for the Bézier form

r(t) := 5 k=0 q k 5 k (1 -t) 5-k t k
of the PH quintic generated by substituting α(t), β(t) into (5) and integrating can be written as

q 1 := q 0 + 1 5 (|α 0 | 2 -|β 0 | 2 , 2α 0 β 0 ) , q 2 := q 1 + 1 5 (Re(α 0 α 1 -β 0 β 1 ), α 0 β 1 + α 1 β 0 ) , q 3 := q 2 + 1 15 (2(|α 1 | 2 -|β 1 | 2 ) + Re(α 2 α 0 -β 2 β 0 ), 4α 1 β 1 + α 0 β 2 + α 2 β 0 ) , ( 7
)
q 4 := q 3 + 1 5 (Re(α 1 α 2 -β 1 β 2 ), α 1 β 2 + α 2 β 1 )
,

q 5 := q 4 + 1 5 (|α 2 | 2 -|β 2 | 2 , 2α 2 β 2 ) ,
where q 0 = (x 0 , y 0 + i z 0 ) is a freely-chosen integration constant.

Remark 1. If α 0 -2α 1 + α 2 = β 0 -2β 1 + β 2 = 0, then α(t)
and β(t) are actually degree-elevated linear polynomials, and the resulting PH curve is a cubic rather than a quintic.

C 2 PH quintic spline equations

For brevity, we restrict our attention here to uniformly-parameterized splines -i.e., the parameter domain for each spline segment r i (t) is t ∈ [ 0, 1 ]. This is appropriate to the Morphosense shape reconstruction problem, since the sensors typically have a uniform distribution along its length. Details of the generalization to non-uniform parameterizations may be found in [START_REF] Farouki | Spatial C 2 PH quintic splines[END_REF]. Given a set of complex coefficients (a i , b i ) for i = 0, . . . , N +1, consider the N quintic PH curves r i (t), t ∈ [ 0, 1 ] generated for i = 1, . . . , N by the complex quadratic polynomials

α i (t) := 1 2 (a i-1 + a i ) (1 -t) 2 + a i 2(1 -t)t + 1 2 (a i + a i+1 ) t 2 , β i (t) := 1 2 (b i-1 + b i ) (1 -t) 2 + b i 2(1 -t)t + 1 2 (b i + b i+1 ) t 2 . ( 8
)
Adoption of this form automatically guarantees continuity of first and second derivatives between consecutive PH spline segments r i (t) and r i+1 (t). Specifically, one can verify [START_REF] Farouki | Design of C 2 spatial Pythagorean-hodograph quintic splines by control polygons[END_REF] that

r ′ i (1) = r ′ i+1 (0) = ( 1 4 (|a i + a i+1 | 2 -|b i + b i+1 | 2 ) , 1 2 (a i + a i+1 )(b i + b i+1 ) ) , r ′′ i (1) = r ′′ i+1 (0) = ( |a i+1 | 2 -|a i | 2 + |b i+1 | 2 -|b i | 2 , 2 (a i+1 b i+1 -a i b i ) ) . ( 9
)
Note that, although C 1 and C 2 continuity is thereby ensured, one must still enforce C 0 continuity. Integration of the hodograph r ′ i (t) allows the initial point r i (0) of each segment to be freely chosen, but satisfaction of the C 0 condition r i (1) = r i+1 (0) with the succeeding segment must be achieved by means of an appropriate scheme (see Sections 4 and 5 below).

From [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF], one can verify that the total arc length

l i := 1 0 |α i (t)| 2 + |β i (t)| 2 dt ( 10
)
of the spline segment r i (t) is defined in terms of the complex coefficients a i , b i by the expression

120 l i = 6 (|a i-1 | 2 + |b i-1 | 2 ) + 54 (|a i | 2 + |b i | 2 ) + 6 (|a i+1 | 2 + |b i+1 | 2 ) + 13 (a i-1 a i + b i-1 b i + a i-1 a i + b i-1 b i + a i a i+1 + b i b i+1 + a i a i+1 + b i b i+1 ) + a i-1 a i+1 + b i-1 b i+1 + a i-1 a i+1 + b i-1 b i+1 . ( 11
)
Note that, although l i is expressed in terms of complex numbers, it has a real value.

Interpolation of points under length constraints

Consider first the problem of interpolating a sequence of points by a spatial C 2 PH quintic spline, with prescribed arc lengths for each spline segment. Specifically, given points p 0 , . . . , p N ∈ R 3 and segment arc lengths l 1 , . . . , l N the problem is to identify a set of PH quintic segments r 1 (t), . . . , r N (t) defined on t ∈ [ 0, 1 ] such that

r i (1) = r i+1 (0) = p i , i = 1, . . . , N -1 and 1 0 |r ′ i (t)| dt = l i , i = 1, . . . , N
with r 1 (0) = p 0 and r N (1) = p N . By using the polynomials ( 8), the first and second derivative continuity conditions ( 9) are automatically satisfied.

Interpolation of points

Each segment r i (t) for i = 1, . . . , N is required to satisfy r i (0) = p i-1 and r i (1) = p i . On integrating the hodograph r ′ i (t), one can choose r i (0) = p i-1 as the integration constant. Achieving r i (1) = p i then corresponds to satisfying the equation

1 0 r ′ i (t) dt = p i -p i-1 := ∆p i ( 12
)
for i = 1, .., N . Writing ∆p i = (∆x i , ∆y i + i ∆z i ) and using [START_REF] Farouki | Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable[END_REF], this equation is equivalent to

1 0 |α i (t)| 2 -|β i (t)| 2 dt = ∆x i , 1 0 2 α i (t)β i (t) dt = ∆y i + i ∆z i , ( 13
)
which can be expressed in terms of the complex numbers a i , b i as

6 (|a i-1 | 2 -|b i-1 | 2 ) + 54 (|a i | 2 -|b i | 2 ) + 6 (|a i+1 | 2 -|b i+1 | 2 ) + 13 (a i-1 a i -b i-1 b i + a i-1 a i -b i-1 b i + a i a i+1 -b i b i+1 + a i a i+1 -b i b i+1 ) + a i-1 a i+1 -b i-1 b i+1 + a i-1 a i+1 -b i-1 b i+1 = 120 ∆x i , ( 14
)
6 a i-1 b i-1 + 54 a i b i + 6 a i+1 b i+1 + 13 (a i-1 b i + a i+1 b i + a i b i-1 + a i b i+1 ) + a i-1 b i+1 + a i+1 b i-1 = 60 (∆y i + i ∆z i ) , ( 15
)
for i = 1, . . . , N . Thus, since equations ( 14) are real and equations (15) are complex, interpolation of the points p 0 , . . . , p N is achieved through the satisfaction 3N scalar conditions on the complex unknowns a 0 , . . . , a N +1 and b 0 , . . . , b N +1 .

Satisfaction of length constraints

From [START_REF] Huard | Serret-Frenet frame interpolation under length constraints for surface reconstruction[END_REF], satisfaction of the arc length constraints yields the equations

6 (|a i-1 | 2 + |b i-1 | 2 ) + 54 (|a i | 2 + |b i | 2 ) + 6 (|a i+1 | 2 + |b i+1 | 2 ) + 13 (a i-1 a i + b i-1 b i + a i-1 a i + b i-1 b i + a i a i+1 + b i b i+1 + a i a i+1 + b i b i+1 ) + a i-1 a i+1 + b i-1 b i+1 + a i-1 a i+1 + b i-1 b i+1 = 120 l i ( 16
)
for i = 1, . . . , N . Comparing equation [START_REF] Farouki | Design of C 2 spatial Pythagorean-hodograph quintic splines by control polygons[END_REF] and the first equation in [START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF], it is apparent that (16) differs from ( 14) only in the replacement of all minus signs with plus signs on the left-hand side.

The real equations ( 14) and ( 16) can thus be simplified by replacing them with one-half of their sum and their difference, namely

6 |a i-1 | 2 + 54 |a i | 2 + 6 |a i+1 | 2 + 13 (a i-1 a i + a i-1 a i + a i a i+1 + a i a i+1 ) + a i-1 a i+1 + a i-1 a i+1 = 60 (l i + ∆x i ) , ( 17
)
6 |b i-1 | 2 + 54 |b i | 2 + 6 |b i+1 | 2 + 13 (b i-1 b i + b i-1 b i + b i b i+1 + b i b i+1 ) + b i-1 b i+1 + b i-1 b i+1 = 60 (l i -∆x i ) . ( 18
)
Interpolating the points p 0 , . . . , p N with prescribed segment arc lengths l 1 , . . . , l N thus amounts to solving the system defined by the N complex equations (15) and 2N real equations ( 17)-( 18) for i = 1, . . . , N . Altogether, these amount to imposing 4N scalar conditions on the 4N + 8 scalar degrees of freedom embodied in the complex values a 0 , . . . , a N +1 and b 0 , . . . , b N +1 .

Shape optimization procedure

In classical C 2 cubic spline interpolation problems, end conditions must be imposed on the first and last segment, in order to match the number of equations and free unknowns. The commonly-used end conditions are (i) quadratic end spans, involving degree reduction of r 1 (t) and r N (t); (ii) the "not-a-knot" condition, in which r 1 (t), r 2 (t) and r N -1 (t), r N (t) meet with C 3 continuity; and (iii) periodic end conditions, appropriate to smooth closed curves [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF].

Since the present focus is on open curves, analogs to end conditions (i) or (ii) must be considered for C 2 PH quintic splines. The "not-a-knot" condition is difficult to formulate with PH curves, but a degree reduction of r 1 (t) and r N (t) from quintic to cubic can be easily implemented. It occurs when α 1 (t), β 1 (t) and α N (t), β N (t) are actually degree-elevated linear (rather than true quadratic) polynomials. By Remark 1, this occurs when

a 0 -2 a 1 + a 2 = b 0 -2 b 1 + b 2 = 0 , a N +1 -2 a N + a N -1 = b N +1 -2 b N + b N -1 = 0 .
These additional conditions provide a system of 2N + 4 equations in the 2N + 4 complex variables a 0 , . . . , a N +1 and b 0 , . . . , b N +1 . It should be noted that these equations are not algebraic, because the variables and their conjugates appear in them. They can be recast as a system of 4N + 8 real equations in 4N + 8 real variables (the real and imaginary parts of a 0 , . . . , a N +1 and b 0 , . . . , b N +1 ), but there is no guarantee of real solutions, because the equations are non-linear.

In fact, preliminary experiments show that in most cases there are no solutions for which end segments with prescribed lengths are degree-elevated cubics (primarily due to the length constraint). Instead of end conditions, an approach based on shape optimization is used here to fix the residual degrees of freedom, similar to that employed in [START_REF] Farouki | Design of C 2 spatial Pythagorean-hodograph quintic splines by control polygons[END_REF]. The idea is to rely on the good shape properties and simpler nature of PH cubic curves. As in [START_REF] Farouki | Design of C 2 spatial Pythagorean-hodograph quintic splines by control polygons[END_REF], minimization of the global energy function

F (a 0 , b 0 , . . . , a N +1 , b N +1 ) := N i=1 |a i-1 -2a i + a i+1 | 2 + |b i-1 -2b i + b i+1 | 2 (19)
will be considered. This can be regarded as a measure of the "distance" of the constructed C 2 PH quintic spline from a PH cubic spline -in particular, minimization of (19) will exactly reproduce a PH cubic spline if the given points and segment arc lengths are compatible with such a curve.

The problem is then one of optimization of the quadratic function ( 19), subject to the quadratic constraints defined by ( 15) and ( 17)-( 18) for i = 1, . . . , N . In order to ensure rapid convergence, it is critical to have a good initial estimate for a 0 , . . . , a N +1 and b 0 , . . . , b N +1 .

Initialization of coefficients

To obtain initial values for the coefficients a i , b i the equations ( 15) and ( 17)-( 18) are first simplified by setting

∆a i = a i -a i-1 , ∆b i = b i -b i-1 for i = 1, . . . , N , substituting a i-1 = a i -∆a i , a i+1 = a i + ∆a i+1 , b i-1 = b i -∆b i , b i+1 = b i + ∆b i+1 ,
into those equations, and retaining only zero-order terms in ∆a i , ∆b i . The equations ( 15) and ( 17)-( 18) then yield

a i b i = 1 2 (∆y i + i ∆z i ) , |a i | 2 = 1 2 (l i + ∆x i ) , |b i | 2 = 1 2 (l i -∆x i ) . ( 20
)
In order to ensure satisfaction of the compatibility condition

|a i b i | 2 = |a i | 2 |b i | 2 ,
equations (20) are then modified to

a i b i = 1 2 γ i (∆y i + i ∆z i ) , |a i | 2 = 1 2 γ i (∆p i + ∆x i ) , |b i | 2 = 1 2 γ i (∆p i -∆x i ) , ( 21
)
where ∆p i := |∆p i | and γ i := l i /∆p i ≥ 1. The solutions of equations (21) have the form

a i = 1 2 γ i (∆p i + ∆x i ) exp(i (θ i + φ i )) , b i = 1 2 γ i (∆p i -∆x i ) exp(i φ i ) , ( 22
)
where φ i is a free parameter, and θ i is defined by

(cos θ i , sin θ i ) := (∆y i , ∆z i ) (∆y i ) 2 + (∆z i ) 2 .
To determine the free parameters φ 1 , . . . , φ N we first set φ 1 = 0, and then obtain φ i from φ i-1 for i = 2, . . . , N by minimizing |∆a i | 2 + |∆b i | 2 with respect to it. We have

|∆a i | 2 + |∆b i | 2 = γ i-1 ∆p i-1 + γ i ∆p i -γ i-1 γ i (∆p i-1 -∆x i-1 )(∆p i -∆x i ) cos(∆φ i ) -γ i-1 γ i (∆p i-1 + ∆x i-1 )(∆p i + ∆x i ) cos(∆θ i + ∆φ i ) ,
where ∆θ i = θ i -θ i-1 and ∆φ i = φ i -φ i-1 . Setting the derivative of this expression with respect to ∆φ i equal to zero yields

tan ∆φ i = -sin ∆θ i cos ∆θ i + (∆p i-1 -∆x i-1 )(∆p i -∆x i ) (∆p i-1 + ∆x i-1 )(∆p i + ∆x i )
.

This identifies two values ∆φ i ∈ (-π, +π ] corresponding to minimum and maximum values of

|∆a i | 2 + |∆b i | 2
, and we must distinguish the former from the latter by evaluating this expression. Obtaining φ i from ∆φ i , the values a i and b i are then completely determined from (22). The six data points correspond to parameter values t = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 on the circular helix r(t) = (cos 2t, sin 2t, 0.3 t), with exact helix segment arc lengths used for the initialization scheme.

Computational procedure

The examples presented below were all computed using the interior-point algorithm of the MATLAB "fmincon" constrained-optimization function. The objective function ( 19) is quadratic and convex in the variables a 0 , b 0 , . . . , a N +1 , b N +1 . The constraints (15) and ( 17)-( 18) are also quadratic, but non-convex. The efficiency of the method is enhanced by exploiting the fact that the Jacobian and Hessian matrices for the cost function and the constraints are sparse and have simple, closed-form expressions. All of the calculations were performed on a PC with a 3 GHz cpu.

It should be noted that solutions may not exist for arbitrarily-chosen data. In particular, when the points p 0 , . . . , p N are sampled from a smooth analytic curve, the experiments indicate that the specified segment arc lengths l 1 , . . . , l N should be "close" to those of the analytic curve, to ensure existence of a solution. Because of the non-linear nature of the problem, it is difficult to formulate precise a priori "closeness" criteria to guarantee the existence of a solution. The segment arc lengths should therefore be regarded as parameters that can be interactively adjusted by the user to ensure the existence of an interpolant with desirable shape properties.

For admissible points and arc-length constraints, the optimization process is found to be robust and rapidly convergent to machine precision, so the prescribed points and arc lengths between them are precisely matched. Although convergence is often observed without invoking the initialization step described in Section 4.4, this step is strongly preferred because (i) convergence is much faster; and (ii) the likelihood of converging to a non-optimal local minimum is greatly reduced.

Numerical examples

Two kinds of applications are considered here, with corresponding initial data sets.

(a) In the first application context, the data (points p i and segment lengths l i ) come from physical measurements, and correspond to a real curve. The data are therefore physically and geometrically coherent, and the interpolation scheme corresponds to a reconstruction problem that aims to find a solution close to the physical curve. Such initial data may, for instance, be obtained by a physical shape sensing system similar to the Morphosense ribbon. In this context, two key observations arise from an analysis of the computed examples.

• The C 2 PH quintic spline interpolation under segment length constraints using exact data allows more faithful curve reconstruction than classical interpolation methods. As is evident in Figure 3, the method is a good alternative to a costly refinement of the sampling: when prescribing lengths close to the actual lengths of the sampled curve, the method provides a C 2 spline curve closer to the original curve than the natural C 2 cubic spline. • If the data are intrinsically constrained by a strong underlying structure (the Morphosense, for instance, follows geodesic curves of the surface on which it is placed), the method gives good results when using the actual lengths between the knots. However, when the lengths are modified even slightly, the reconstructed curves may exhibit undesired oscillations that are often difficult to control. Figure 4 illustrates this in the case of the data points t = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 uniformly sampled from the circular helix r(t) = (cos 2t, sin 2t, 0.3 t), using both the exact helix arc lengths, and arc lengths deviating by ±1% from the exact values. Table 1 lists, for this example, the number of iterations and cpu time required for convergence of the optimization process to machine precision, with and without use of the initialization scheme described in Section 4.4. Although the optimization typically converges with even "random" choices for the initial values, it is apparent that the method of Section 4.4 can substantially accelerate the convergence (and is more likely to identify the global minimum). Table 2 illustrates how the number of iterations and cpu time required for convergence to machine precision vary with the number of data points uniformly sampled from the circular helix, using the initialization scheme of Section 4.4.

(b) In a different context, one may wish to only interpolate the points p i and exploit the segment length constraints l i as "tension" or "design" parameters. In this case, the lengths l i must be chosen carefully. Obviously, the lengths must satisfy the lower bounds l i ≥ |∆p i |, and excessive lengths can introduce undesired behavior (such as loops). Figure 5 shows a case in which different lengths are chosen for the second segment, resulting in two radically different solutions for the curve (even with initialization without initialization iterations cpu time iterations cpu time Table 2: Dependence of number of required iterations and total cpu times on the number of sampled data points for the circular helix example, using the initialization scheme described in Section 4.4.

l i = 0.
Figure 5: Two views of an example interpolating points p 0 = (1, 1, 0), p 1 = (0, 4, 3), p 2 = (2, 4, 2), p 3 = (3, 3, 3) with equal lengths l 1 = 4.9 and l 3 = 2.0 for the first and last segment but the different values l 2 = 2.9 (dashed curve) and l 2 = 2.5 (solid curve) for the middle segment (a 16% difference).

for adjacent segments of fixed length). Figure 6 shows another case, in which three different PH quintic splines interpolating the points

p 0 = (1, 1, 0) , p 1 = (0, 4, 3) , p 2 = (2, 1, 1) , p 3 = (2, 4, 2) , p 4 = (3, 3, 3)
are constructed, with fixed segment lengths l 2 = 5.0, l 3 = 3.9, l 4 = 2.75, and three different values for l 1 . This example illustrates the use of the segment length as a design parameter. The number of iterations and total cpu time required for convergence to machine precision are listed in Table 3, with and without the initialization procedure described in Section 4.4. Table 3: Number of iterations for convergence to machine precision and compute times required for the example in Figure 6, with and without use of the initialization scheme described in Section 4.4.

It should be noted that solutions may not exist when the prescribed segment arc lengths l i are large compared to the spacings ∆p i between the points -and, if they do exist, they are generally of undesirable shape. Figure 7 illustrates this in the case of the four equidistant points

p 0 = (0, 0, 0) , p 1 = (0, 1, 0) , p 2 = (1, 1, 0) , p 3 = (1, 1, 1) , ( 23
)
with one or more of the ratios γ i = l i /∆p i being large compared to unity. For comparison, the case (γ 1 , γ 2 , γ 3 ) = (1.15, 1.08, 1.15) is also shown in Figure 7 -this corresponds to the solution of the optimization problem when the length constraints are relaxed. Finally, the existence of solutions that correspond to local minima of the objective function is illustrated in Figure 8. The four data points for this example are the same as in Figure 5, and the prescribed segment lengths are l 1 = 9.25, l 2 = 4.75, l 3 = 3.70. When the initialization scheme of Section 4.4 is employed, the optimization converges to the global minimum, but may otherwise converge to one of the solutions representing a local minimum.

C 2 interpolation of derivatives under length constraints

We consider now the interpolation of derivative vectors d 0 , . . . , d N (rather than points p 0 , . . . , p N ) by C 2 PH quintic splines, subject to prescribed lengths l 1 , . . . , l N for each segment.

Interpolation of derivatives

Given derivative vectors d 0 , . . . , d N and prescribed segment lengths l 1 , . . . , l N we wish to compute a C 2 PH quintic spline, as defined in Section 3, with segments r i (t) for i = 1, . . . , N such that (a) r i (1) = r i+1 (0

) for i = 1, . . . , N -1 , (b) r ′ 1 (0) = d 0 , r ′ N (1) = d N , r ′ i (1) = r ′ i+1 (0) = d i for i = 1, . . . , N -1 , (c) 1 0 |r ′ i (t)| dt = l i for i = 1, . . . , N .
As with the previous method, we note that condition (a), together with the constructions introduced in Section 3, automatically generates a C 2 PH quintic spline. Each of the conditions (a)-(c) is considered in greater detail below.

(a) This condition simply states that the spline must have positional continuity, i.e., each segment r i+1 (t) must begin where the previous segment r i (t) ends. This is achieved by taking r i+1 (0) = r i (1) as the integration constant upon integrating r ′ i+1 (t). Thus, only the first integration constant r 1 (0) is free, corresponding to a translational freedom of the solution.

(b) Since PH curves are generated by integrating their hodograph, this condition can be formulated directly from equation ( 9) as

1 4 (|a i + a i+1 | 2 -|b i + b i+1 | 2 ) = dx i , 1 2 (a i + a i+1 )(b i + b i+1 ) = dy i + i dz i , (24)
where d i = (dx i , dy i , dz i ) for i = 0, . . . , N . As with the previous method, the first of equations ( 24) is real, while the second is complex. The relations (24) thus amount to 3(N + 1) scalar equations (three for each interpolated derivative vector).

(c) As in the previous method, the N equations (16) expressing the length constraints are imposed.

The system defined by ( 16) and (24) comprises 4N +3 scalar equations, to be satisfied by the 4N +8 scalar freedoms embodied in the complex coefficients a 0 , . . . , a N +1 and b 0 , . . . , b N +1 . This problem is therefore somewhat more constrained than that of Section 4 -three additional constraints arise from the interpolation of the initial derivative d 0 . In the previous method, interpolation of the first point was achieved by choice of the integration constant, but this is not possible for interpolation of the initial derivative vector in the present context. As in Section 4, the spline computation will be based on minimization of the function (19) under the constraints ( 16) and (24).

Initialization of coefficients

The initialization method is analogous to that employed in Section 4. Substituting a i-1 = a i -∆a i , a i+1 = a i +∆a i+1 , and b 16) and (24), and retaining only the zero order terms, we obtain the simplified equation system

i-1 = b i -∆b i , b i+1 = b i +∆b i+1 in equations (
|a i | 2 + |b i | 2 = l i , |a i | 2 -|b i | 2 = dx i , a i b i = 1 2 (dy i + i dz i ) . ( 25
)
By combining the first and second equations, these can be re-formulated as

|a i | 2 = 1 2 (l i + dx i ) , |b i | 2 = 1 2 (l i -dx i ) , a i b i = 1 2 (dy i + i dz i ) . ( 26
)
This system is similar to the simplified equations (20) for the previous method, and can be viewed as a "continuous" version of the latter, with the the derivatives (dx i , dy i , dz i ) replacing the differences (∆x i , ∆y i , ∆z i ). Figure 9 shows an example curve generated by the initial approximation.

Numerical examples

The computational procedure employed in the present context is, with appropriate modifications, as described in Section 4.5. We consider again both curve reconstruction and curve design contexts. The derivative interpolation scheme described above was motivated by the Morphosense shape reconstruction problem. However, the Morphosense provides unit orientation vectors, rather than derivatives of non-unit magnitude. Several methods are currently being investigated to address this issue. Figure 10 Figure 11 illustrates a reconstruction from arc lengths and unit tangent vectors sampled from the circular helix r(t) = (cos 2t, sin 2t, 0.3 t) at parameter values t = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 (compare with Figure 4). The unit tangents are scaled by the segment arc lengths as in the preceding example. As with the point interpolation method of Section 4, the shape is seen to be quite sensitive to the prescribed length constraints. If the specified segment lengths differ significantly from those of the analytical curve, the interpolant exhibits strong oscillations. Finally, Figure 14 illustrates the behavior of interpolants to prescribed derivatives and segment arc lengths in the case of "large" lengths. As a reference case, the arc lengths are specified as

Figure 1 :

 1 Figure 1: Discrete geometric data furnished by the Morphosense ribbon for shape reconstruction.

  Figure 2 compares the curves generated by the coefficients a 0 , . . . , a N +1 and b 0 , . . . , b N +1 obtained from the above initialization scheme, with that generated by a random choice of the coefficients. The former is clearly much closer to interpolating the prescribed data points and segment arc lengths, resulting in much faster convergence to the desired interpolant -see Section 4.6 below.

Figure 2 :

 2 Figure 2: PH quintic spline curves defined by coefficients a 0 , . . . , a N +1 and b 0 , . . . , b N +1 generated by the initialization scheme of Section 4.4 (left) and by a set of randomly-chosen coefficients (right).The six data points correspond to parameter values t = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 on the circular helix r(t) = (cos 2t, sin 2t, 0.3 t), with exact helix segment arc lengths used for the initialization scheme.

Figure 3 :

 3 Figure 3: Approximation of a sine curve (solid) by a C 2 cubic spline (dotted) interpolating discrete points, and by a C 2 PH quintic spline (dashed) interpolating both points and segment arc lengths. The data points correspond to the parameter values t = 0, 1 3 π, 2 3 π, π on the curve r(t) = (t, sin t, 0).

Figure 4 :

 4 Figure 4: (a) A helix with known segment lengths L i between sample points p 0 , . . . , p 5 . Also shown are reconstructions interpolating the given points with segment lengths (b) l i = 0.99 L i ; (c) l i = L i ; and (d) l i = 1.01 L i illustrated as (x, z) and (x, y) plots. (e) Superposition of the reconstructions.

Figure 6 :

 6 Figure 6: Examples of interpolants to fixed points and equal segment lengths l 2 , l 3 , l 4 but different l 1 values. Upper left curve: l 1 = 6.3. Upper right curve: l 1 is increased by 15% to 7.24. Lower left curve: l 1 is decreased by 15% to 5.35. Lower right: all three of the interpolants are superimposed.

Figure 7 :

 7 Figure 7: Examples of C 2 PH quintic splines interpolating the points (23) with prescribed segment lengths l i that substantially exceed the spacings ∆p i between successive points (dashed curves). The ratios γ i = l i /∆p i have the values (γ 1 , γ 2 , γ 3 ) = (1.1, 2, 10) on the left; (5, 2, 5) in the center; and (10, 3, 10) on the right. The solid "reference" curve corresponds to (γ 1 , γ 2 , γ 3 ) = (1.15, 1.08, 1.15).

Figure 8 :

 8 Figure 8: Interpolants corresponding to the global minimum (solid curve) and local minima (dashed and dotted curves) of the objective function. The data points are as in Figure 5 with segment lengths (l 1 , l 2 , l 3 ) = (9.25, 4.75, 3.70). The initialization scheme ensures convergence to the global minimum.

Figure 9 :

 9 Figure 9: Curve determined by the initial values for derivative and arc length data from a helix.

  shows preliminary results obtained from noisy data measured by a Morphosense. For this example, the normalized derivatives obtained from the sensors are d 0 = (0.3862, -0.0315, 0.9219) , d 1 = (0.3775, -0.1613, 0.9118) , d 2 = (0.4077, -0.1947, 0.8921) , d 3 = (0.4444, -0.3071, 0.8415) , d 4 = (0.3224, -0.4485, 0.8336) , d 5 = (0.3913, -0.4889, 0.7797) , d 6 = (0.2771, -0.5945, 0.7548) , d 7 = (0.1710, -0.7053, 0.6880) , and the measured distances between sensors are l 1 = 170.05, l 2 = 169.50, l 3 = 168.26, l 4 = 169.33, l 5 = 169.53, l 6 = 168.82, l 7 = 170.10. For the reconstruction, the vectors d i are scaled by 1 2 (l i +l i+1 ) for i = 1, . . . , 6 while d 0 and d 7 are scaled by l 1 and l 7 respectively. With initialization, the solution converged in 31 iterations and 0.28 seconds (260 iterations and 2.25 seconds without initialization). These convergence figures, with initialization, are representative of all subsequent examples.

Figure 10 :

 10 Figure 10: Acquisition and reconstruction of a curve on a cone with the Morphosense ribbon. The analytical curve (solid) is acquired by the ribbon and reconstructed as a PH quintic spline (dashed).

Figure 11 :

 11 Figure 11: (a) A helix with known segment lengths L i between sampled derivatives d 0 , . . . , d 5 . Reconstructions interpolating the given derivatives with segment lengths (b) l i = 0.99 L i ; (c) l i = L i ; and (d) l i = 1.01 L i are illustrated as (x, z) and (x, y) plots. (e) Superposition of the reconstructions.

Figure 12 :

 12 Figure 12: An example with different lengths l 1 for the first segment. For the dashed curve, the first segment has length l 1 = 2, and this is increased by 1% for each of the successive (solid) curves.

Figures 12

 12 Figures 12 and 13 show examples of reconstructions from arbitrarily-chosen derivative vectors and arc lengths, illustrating again the strong influence of the length constraints on the curve shape. In Figure12, the derivatives and segment arc lengths ared 0 = (2, 1, 1), d 1 = (1, 1, 1), d 2 = (0, 1, 0), d 3 = (1, 1 2 , 12) and l 1 = 2.0, l 2 = 1.4, l 3 = 1.0, with l 1 being increased by 1% for successive curves. In the Figure13example, we choose d 0 = (-2, -1, -1),d 1 = (1, 1 2 , 1 2 ), d 2 = (0, 1, 1 2 ), d 3 = (0, -1, -12 ) and l 1 = 1.2, l 2 = 1.4, l 3 = 0.7, with l 2 being increased by 1% for successive curves. The segment lengths evidently serve as useful tension parameters, facilitating manipulation of the curve shape in geometric design applications.

Figure 13 :

 13 Figure 13: An example with different lengths l 2 for the second segment. For the dashed curve, the second segment has length l 2 = 1.4. This increases by 1% for each of the successive (solid) curves.

Table 1 :

 1 Number of iterations for convergence to machine precision and compute times required for the helix example in Figure4, with and without the initialization procedure described in Section 4.4.

	99 L i	30	0.24 sec.	130	0.85 sec.
	l i = 1.00 L i	49	0.29 sec.	135	0.67 sec.
	l i = 1.01 L i	24	0.17 sec.	124	0.73 sec.
		data points iterations cpu time
		4	15	0.12 sec.
		6	37	0.33 sec.
		8	45	0.47 sec.
		10	62	0.80 sec.
		15	69	1.03 sec.
		20	141	3.14 sec.

Laboratoire d'électronique des technologies de l'information.

Commissariat à l'Energie Atomique et aux Energies Alternatives.

These extra freedoms are necessary to ensure a rotation-invariant representation[START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF][START_REF] Farouki | Hermite interpolation by rotation-invariant spatial Pythagorean-hodograph curves[END_REF].

Bold characters will be used to denote both complex numbers and points/vectors in R 3 -the meaning should be clear from the context.

choice of derivative magnitudes was briefly discussed in Section 5.3). C 2 PH quintic splines cannot interpolate Serret-Frenet frames, but one may come closer to this requirement by relaxing from C 2 to G 2 continuity. Further investigations along these lines will be reported in due course.

l 1 = l 2 = l 3 = 1 and the derivatives are defined by the unit vectors d 0 = (1, 0, 0) , d 1 = (0, 1, 0) , d 2 = (0, 0, 1) , d 3 = (1, 0, 0) .

(27)

When the segment lengths are increased to l 1 , l 2 , l 3 = 3 (without changing d 0 , d 1 , d 2 , d 3 ) the overall size of the curve increases, but the individual segments become very "flat" because of the relatively small derivative magnitudes. However, if d 0 , d 1 , d 2 , d 3 are also scaled by a factor of 3, the resulting interpolant is simply a uniformly-scaled instance of the reference curve. This behavior, apparent in Figure 14, illustrates the importance of using commensurate derivative magnitudes and segment arc lengths to ensure interpolants of good shape quality.

Conclusion

Two novel C 2 curve interpolation schemes, motivated by applications to shape reconstruction from discrete micro-electronic sensor data, have been proposed. A distinctive feature of these methods is the satisfaction of prescribed segment length constraints, based on the special algebraic structure of spatial Pythagorean-hodograph (PH) curves. This feature is very important in contexts where the precise spacing between sensors is known, and also provides useful shape or tension parameters for manipulating curves in geometric design contexts. It is expected that many other physical or geometrical modelling applications can also benefit from this capability. The interpolation schemes involve the optimization of a quadratic cost function, subject to a system of quadratic constraints, and methods to identify a good starting approximation ensure an efficient and robust iterative solution of this problem. Computed examples show that the prescribed segment lengths strongly influence the resulting curve shapes. In particular, for the problem of curve reconstruction from a smooth physical shape, specifying the precise segment lengths significantly improves the accuracy of the reconstruction -especially in the case of sparse sampling.

The schemes presented herein constitute preliminary steps toward solving the Morphosense shape reconstruction problem, introduced in Section 1. They incorporate the important feature of reproducing the precise curvilinear distances between sensor locations, but the method of Section 5 interpolates parametric derivative vectors, rather than Serret-Frenet frames with unit vectors (the


Manabu Tsukada
email: manabu.tsukada@inria.fr

José Santa
email: jose.santa@cud.upct.es

Satoshi Matsuura
email: matsuura@is.naist.jp

Thierry Ernst
email: thierry.ernst@mines-paristech.fr

Kazutoshi Fujikawa
email: fujikawa@itc.naist.jp

AnaVANET: an experiment and visualization tool for vehicular networks

Keywords: Vehicular Ad-hoc Networks, Experimental Evaluation, Wireless Multihop Communication, Network Mobility, Visualization Tool

The experimental evaluation of wireless and mobile networks is a challenge that rarely substitutes simulation in research works. This statement is even more evident in vehicular communications, due to the equipment and effort needed to obtain significant and realistic results. In this paper, key issues in vehicular experimental evaluation are analyzed by an evaluation tool called AnaVANET, especially designed for assessing the performance of vehicular networks. This software processes the output of well-known testing tools such as ping or iperf, together with navigation information, to generate geo-aware performance figures of merit both in numeric and graphical forms. Its main analysis capabilities are used to validate the good performance in terms of delay, packet delivery ratio and throughput of NEMO, when using a road-side segment based on IPv6 GeoNetworking.

Introduction

Vehicular networks are essential for Intelligent Transportation Systems (ITS) to optimize the road traffic and achieve safe, efficient and comfortable human mobility. Essentially, there are two main communication paradigms in vehicular communications, vehicle to vehicle (V2V) and vehicle to infrastructure (V2I), depending on whether the communication is performed directly between vehicles or using nodes locally or remotely installed on the road infrastructure.

When the V2V paradigm is considered, the research field is commonly called Vehicular Ad-hoc Networks, or VANET. Although there are a lot of works related to VANET applications and basic research at physical, MAC and network layers, there is a significant lack of real evaluation analysis in this field, due to cost and effort implications. A number of experimentation works and supporting tools should be improved in the short term, in order to give real evidences to car manufacturers and road operators of the benefits of vehicular communications.

Conventional network measurement tools (e.g. iperf, ping or traceroute) assume fixed networks and assess network performances in a end-to-end basis. However, under dynamic network conditions such as in the vehicular networks case, it is difficult to measure detailed status of networks by using solely these tools, because vehicles are always changing their location and the performance of wireless channels fluctuates. In order to solve these issues, we have developed a packet analysis and visualization tool called AnaVANET1 , which considers the peculiarities of the vehicular environment for providing an exhaustive evaluation software for outdoor scenarios. Both V2V and V2I networks can be efficiently analyzed, thanks to the integrated features for collecting results, post-processing data, generate graphical figures of merit and, finally, publish the results in a dedicated web site (if desired).

The rest of the paper is organized as follows. Section 2 introduces the readers about the network layer protocols in vehicular networks experimentation. Then, the issues and requirements for evaluating vehicular networks are listed in Section 3. The evaluation methodology desired in this frame is described in Section 4 and, as a result of our analysis, the design and implementation of the AnaVANET evaluation tool is detailed in Section 5, together with a reference evaluation of a network testbed using the tool in Section 6. Finally, Section 7 concludes the paper summarizing the main results and addressing future works.

Network Protocols in Vehicular Networks

Network protocols in vehicular networks can be classified in infrastructureless scenarios, i.e. V2V, and infrastructure-based scenarios, i.e. V2I, as showed in Fig. 1. The infrastructure-less scenario is well-known by the research area of VANET or Mobile Ad-hoc Networks (MANET). These approaches are designed to enable wireless communications in dynamic topologies without any infrastructure. Routing protocols here are further classified as topology-based and position-based routing protocols. Topology-based protocols were divided into two main branches by the IETF MANET working group: reactive, where nodes periodically exchange messages to create routes (e.g. AODV [START_REF] Perkins | Ad hoc On-Demand Distance Vector (AODV) Routing[END_REF]), and proactive, in which control messages are exchanged on demand when it is necessary to reach a particular node (e.g. OLSR [START_REF] Clausen | Optimized Link State Routing Protocol (OLSR)[END_REF]).

Unlike topology based routing, position based routing does not need to maintain part of the network structure in order to forward packets towards the destination node. When routing packets based on position, nodes forward the packets with the aim of reaching the nodes within a geographical location. Thus, position based routing can eliminate the problem that appears in topology based protocols when routes become quickly unavailable in high mobility scenarios. In Greedy Perimeter Stateless Routing (GPSR) [START_REF] Karp | Gpsr: Greedy perimeter stateless routing for wireless networks[END_REF], for instance, the intermediate nodes make a decision based on the destination position and neighbor positions. The Car-to-Car Communication Consortium (C2CC) also specified the C2CNet protocol, which was later enhanced by the GeoNet project to support IPv6. Within the ITS standardization domain, GeoNetworking [START_REF]Intelligent Transport Systems (ITS); Vehicular Communications; Part 4: Geographical Addressing and Forwarding for Point-to-Point and Point-to-Multipoint Communications; Sub-part 1: Media-Independent Functionality[END_REF] is being completed by ETSI at the moment, integrating several geo-aware strategies to better route packets in vehicular networks.

On the other side, infrastructure-based protocols have been focused on the global connectivity of nodes to the Internet. Mobile IPv6 [START_REF] Perkins | Mobility Support in IPv6. RFC 6275[END_REF] solved the mobility problem for mobile hosts and, later, Network Mobility Basic support (NEMO) [START_REF] Devarapalli | Network Mobility (NEMO) Basic Support Protocol[END_REF] provided a solution for the mobility of a whole network (e.g. a vehicle or bus), which has been recommended by the ISO TC204 WG16 to achieve Internet mobility for vehicles.

Issues and Requirements for VANET Evaluation

Using multi-hop and dynamic routing strategies presents a challenge in the evaluation of vehicular networks. Common end-to-end evaluation tools such as ping6 and iperf are useless to track the effect of route change, because they are unaware of the path taken during a communication test. An additional lack of these tools is the possibility to measure the performance of hop-by-hop links, since the study is carried out end-to-end. Also, geographical and external factors such as nodes position, distance between nodes or obstacles are not linked with network performance figures of merit.

With the aim of summarizing these main requirements when evaluating multi-hop vehicular networks, the next needs are found essential by the software tools used in experimental campaigns for evaluating both V2V and V2I:

Path detection The topology of a vehicular network with dynamic routing changes frequently as vehicles move, and the communication path is changed accordingly. Thus, the tool should take note of the communication path used in every moment. Communication performance in links Once the communication path is tracked, the tool should measure the performance in a link-by-link as well as end-toend basis.

Geographical awareness

The network performance in a link depends on various geographical factors, such as the distance between the nodes, the movement speed and direction, and the existence of obstacles in the communication link. Thus, the evaluation tool should take the above geographical factors into account. Intuitive visualization Performance figures of merit and environmental information should be shown together in a synchronized way, and the spatiotemporal data series should be available in post process to play them at different speeds, stop when desired, or replayed freely as he or she wants. Independence from network protocols Given the various network layer protocols in vehicular communications, the evaluation tool should be independent from the one chosen. Independent from devices Since the configuration of vehicle and infrastructure devices may differ, the evaluation tool should not rely on any specific device functionality. Adaptation to various scenarios The software evaluation tool should accommodate to all possible communication scenarios (moving or static, urban or highway, etc.). Easiness for data collection Since a lot of experiments could be needed in a extensive campaign, the easiness of gathering data and deploying software modules in devices is essential.

As it is later described, the evaluation tool presented in this work (Ana-VANET) copes with the previous requirements.

Evaluation Methodology

The evaluation goals are to analyze which testing conditions affect which data flows or network protocols. For achieving this end it is necessary to design a proper evaluation methodology. Within it we should consider the tendency of results by repeating tests with the same settings or varying parameters under study, such as the network protocol, the mobility of nodes or the data volume. The overall analysis should be supported by a proper evaluation tool, such as the later presented AnaVANET. This section details both the testing conditions and the possible routing protocols to consider, as it is summarized in Fig. 2, by introducing the concept and presenting our real use case for testing the performance of NEMO over IPv6 GeoNetworking.

Testing Conditions

Testbed Platform The testbed used for the evaluation of a network architecture should be carefully chosen to implement most relevant nodes in real software and hardware. In vehicular communications, this is extremely important, since a good deployment could be needed in case of testing V2V multi-hop networks. Fig. 2. Evaluation Methodology

In our particular case, the testbed comprises a set of four vehicles and two roadside stations, as illustrated in Fig. 3. Each vehicle is equipped with a mobile router (MR), with at least two interfaces: an Ethernet link to connect mobile network nodes (MNNs) within the in-vehicle network, and a wireless adapter in ad-hoc mode used for both V2V and V2I communications. On the roadside, access routers (ARs) are prepared to be fixed on the top of a building or any other elevated point near the road. Each one provides two interfaces: an Ethernet link for a wired Internet access, and a wireless adapter in ad-hoc mode to connect with vehicles in the surroundings. At a backend point in the Internet, a home agent (HA) is installed to support Internet mobility of MRs by using NEMO.

Among the various testbed conditions, the hardware specification (CPU, memory, etc), antenna and wireless settings are important factors for the evaluation, since they will highly affect the results. In our case, MRs are Alix3d3 embedded boxes provided with a Linux 2.6.29.6 kernel. Each MR has a mini-pci wireless card Atheros AR5414 802.11 a/b/g Rev 0, and an antenna 2.4GHz 9dBi indoor OMNI RP-SMA6 is used. The frequency used has been 2.422Ghz and the data rate has been fixed to 6 Mbits/s.

Testing Scenarios Fixing the evaluation scenarios beforehand is essential in the planning of a testing campaign. In general, the main factors that determine the possible scenarios are:

Mobility Static scenarios can be chosen to test the network operation in a controlled way, but also dynamic ones can be used in a realistic evaluation. A dynamic scenario is considered in our case. Location The place in which the tests are carried out impacts on the network performance, due to signal propagation blockage issues above all. In our case a semi-urban scenario is used within the INRIA-Rocquencourt installations.

Number of vehicles

The number of hops between the source and the destination vehicles affect the communication delay and the higher probability of packet looses, due to route changes or MAC transmission issues. Up to four vehicles are considered in our case.

As summarized in Fig. 4, testing scenarios have been divided into urban and highway; mobility has been set to static, urban-like speed, and high speed.

Data Flows and Performance Indicators

A number or protocols and data flows can be set for evaluations, however, only the most representative and more used in the literature should be considered to study concrete performance indicators. For instance, in our case UDP, TCP and ICMPv6 are used to measure the network performance between two communication end-nodes (MNN to MNN) mounted within two vehicles:

UDP is a connection-less unidirectional transmission flow. The traffic is generated by iperf in our case. It is considered that with UDP the performance indicators under consideration can be the packet delivery ratio (PDR), throughput and jitter. TCP is a connection-oriented bidirectional transmission flow. This traffic is also generated by iperf in our case. The performance indicator under consideration here has been the maximum throughput. ICMPv6 is a bi-directional transmission flow. The traffic is generated by ping6 in our case. The performance indicator under consideration can be the road trip delay time (RTT) and PDR.

System Design and Implementation of AnaVANET

AnaVANET (initially standing for Analyzer of VANET) is an evaluation tool implemented in Java to assess the performance of vehicular networks. It takes as input the logs generated by the iperf, tcpdump and/or ping6, together with navigation information in NMEA format, to compute the next performance metrics: network throughput, delay, jitter, hop count and list of intermediate nodes in the communication path, PDR end-to-end and hop-by-hop, speed, and instantaneous position.

AnaVANET is put in the context of the evaluation scenario described in the previous section in Fig. 5 , showing also the main inputs and outputs of the tool. The sender MNN (left most vehicle) is in charge of generating data traffic, and both the sender and the receiver (right most vehicle) MNNs record a high level log, according to the application used to generate network traffic (iperf and ping6 for the moment). All MRs record information about forwarded data packets by means of the tcpdump tool, and log the vehicle position continuously. All this data is post-processed by the AnaVANET core software and then analyzed. The tool traces all the data packets transmitted from the sender node to detect packet losses and calculate statistics for each link and end-to-end, and then merge all these per-hop information with transport level statistics of the traffic generator. As a result, AnaVANET outputs an XML file with statistics on a one-second basis, and a packet trace file with the path followed by each data packet. Once generated, performance metrics can be graphically showed through plots generated by gnuplot and a website where all tests are available. The screenshot of the website is shown in left bottom of Fig. 5. Accessing the website one can replay the tests on a map to see momentary figures of merit.

On the map, the position and movement of the vehicle are depicted with the speed of each vehicle and the distance between them. The transferred data size, bandwidth, packet loss rate, RTT and jitter, for each link and end-to-end are displayed. The network performance is visualized by the width of links and the colors used to draw them.

Evaluation of NEMO over IPvGeoNetworking

Early versions of AnaVANET were designed for evaluating infrastructure less network protocols, as used in our previous works for analyzing OLSR in vehicular environments [START_REF] Santa | Assessment of vanet multi-hop routing over an experimental platform[END_REF] and later tests of IPv6 over C2CNet [START_REF] Tsukada | Experimental evaluation for IPv6 over VANET geographic routing[END_REF] in the FP7 GeoNet project.

The current version of AnaVANET can also analyze infrastructure-based network protocols such as NEMO. In this section, we report a summary of the results collected in the evaluation of NEMO over IPv6 GeoNetworking when a vehicle connects with a node located in the Internet using two roadside units as access routers. The umip.org2 implementation of NEMO is used and the cargeo6.org3 software is used for IPv6 GeoNetworking. ICMPv6 and UDP evaluations in handover scenarios were performed at INRIA Paris-Rocquencourt campus with the two ARs previously presented in the testbed description. The speed of the vehicle was limited to less than 15 km/h, like in a low mobility urban scenario. The reader can directly click in Fig. 6 -Fig. 9 to see the correspondent result in the AnaVANET web viewer, to further perceive the details of the gathered results.

ICMPv6 echo requests (64 bytes) are sent from the MNN to a common computer located in the wired network twice in a second, which replies with ICMPv6 echo replies. The results collected in the ICMPv6 tests are plotted in Fig. 6. The lower part shows the itinerary of the vehicle and the locations of AR1 and AR2 on the map, whereas the upper part shows the RTT, the packet loss and the result of the mobility signaling. The X-axis and the Y-axis of the upper part are the latitude and the longitude of the vehicle, corresponding to the road stretch indicated in the lower part of the figure. When either the request or the reply is lost, the RTT is marked with a zero value and, at the same time, a packet loss is indicated. A binding registration success is plotted when the NEMO binding update (BU) and the corresponding binding acknowledgment (BA) are successfully processed. On the contrary, if either of them is lost, a binding registration fail is plotted at the position. Fig. 7 shows the same results of the test, but referred to the test time. The upper graph shows the RTT and the distance to the two ARs; the middle one shows the PDR obtained with the two ARs; and, finally, the lower plot shows the status of the NEMO signaling. A NEMO success means that the binding registration has been successfully performed, and a fail indicates that either the BU or the BA has been lost.

The results collected in the UDP tests are plotted in Fig. 8. UDP packets are sent from the MNN to the wired node at a rate of 1 Mbps and a length of 1250 bytes. The lower part of the figure shows the itinerary of the vehicle, and the upper part corresponds to the PDR obtained with the ARs and the binding registration results, as in the previous case. The road stretch is the same one used above, but the vehicle moves on the contrary direction in this case.

In the time-mapped results showed in Fig. 9, the upper graph shows the UDP throughput from the MNN to the wired node, the middle part shows the PDR to the two ARs, and the lower plots the status of the NEMO signaling. Success of NEMO status means that the binding registration is successfully performed and Fail means that either the BU or the BA is lost. The paper has presented the peculiarities of evaluating vehicular networks experimentally, through presenting the most used protocols and detailing the needs of the software tools to be used for this task. After that, the importance of the testing methodology is described, and a reference design of a vehicular network evaluation is used to exemplify it. The testbed design and implementation, testing scenarios, routing protocols and data flows, are found essential to be fixed beforehand to avoid improvisation during the testing campaign. The AnaVANET platform is then presented as an efficient evaluation software to process the data gathered by common testing tools, and then generate lots of performance indicators of the trials. The capabilities of AnaVANET are exploited in a novel evaluation of NEMO over IPv6 GeoNetworking, using the tool to gather RTT, PDR and channel throughput information. The results reveal that mobile IPv6 connectivity can be maintained in a V2I case using GeoNetworking over WiFi to pass NEMO IPv6 traffic between vehicles and infrastructure.

Our future work includes, first, a link layer extension of the system to analyze the channel quality (RSSI), load ratio and coverage map. Second, it is considered the support for multicast data flows, since it is essential for the dissemination of events in vehicular networks. Third, we plan to evaluate a real application developed for cooperative ITS.

Fig. 1 .

 1 Fig. 1. Network Protocols in vehicular networks

Fig. 4 .

 4 Fig. 3. Network Configuration

Fig. 5 .

 5 Fig. 5. Overview of AnaVANET

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Map-based RTT, Packet Losses and Mobility Signaling of ICMP evaluation in a handover scenario

http://anavanet.net/

http://umip.org

http://www.cargeo6.org

Acknowledgment

This work has been sponsored by the European 7th FP, through the ITSSv6 (contract 270519), FOTsis (contract 270447) and GEN6 (contract 297239) projects, and the Spanish Ministry of Science and Innovation, through the Walkie-Talkie project (TIN2011-27543-C03).