Minjia Shi

Denis S Krotov
email: krotov@math.nsc.ru

Patrick Solé
email: sole@enst.fr

On the coset graph construction of distance-regular graphs *

Keywords:

We show that no more new distance-regular graphs in the tables of the book of (Brouwer, Cohen, Neumaier, 1989) can be produced by using the coset graph of additive completely regular codes over finite fields.

Introduction

Since the times of Delsarte [?], an important tool to construct distanceregular graphs has been the coset graph of a completely regular code. Many examples can be found in Chapter 11 of [?], involving notably the Golay codes and the Kasami codes (for a new description of these codes see e.g. [?]). More recently, a new distance-regular graph related to the dodecacode, an important additive F 4 -code of length 12 was constructed in [?]. This solved both a thirty year old problem in graph theory [?, Ch. 14], and a forty-five year old problem in coding theory [?].

• v = 343 = 7 3 , {48, 35, 9; 1, 7, 40} (Corollary ??);

• v = 625 = 5 [START_REF] Brouwer | Distance-Regular Graphs[END_REF] , {52, 42, 16; 1, 6, 28} (Section ??);

• v = 729 = 3 [START_REF] Krotov | On weight distributions of perfect colorings and completely regular codes[END_REF] , {56, 42, 20; 1, 6, 28} (Section ??);

• v = 512 = 2 [START_REF] Maruta | The nonexistence of some ternary linear codes of dimension 6[END_REF] , {63, 48, 10; 1, 8, 54} (Lemma ??);

• v = 729 = 3 [START_REF] Krotov | On weight distributions of perfect colorings and completely regular codes[END_REF] , {80, 63, 11; 1, 9, 70} (Corollary ??);

• v = 729 = 3 [START_REF] Krotov | On weight distributions of perfect colorings and completely regular codes[END_REF] , {104, 66, 8; 1, 12, 88}; a graph does not exist [?].

In this note we discuss the possibility to construct distance-regular graphs with the parameters above as coset graphs of linear completely regular codes in Hamming graph H(n, q), where q is prime. (Note that the existence of a completely regular code in H(n, q) where q = p s for prime p implies the existence of a completely regular code in H(n(q -1)/(p-1), p) with the same intersection matrix, because H(n(q -1)/(p -1), p) covers H(n, q), see also Lemma ?? below.)

Definitions and basic facts

Graphs

All graphs in this note are finite, undirected, connected, without loops or multiple edges. In a graph Γ, the neighborhood Γ(v) of a vertex v is the set of vertices adjacent to v. The degree of a vertex v is the size of Γ(v). A graph is regular if every vertex has the same degree. The i-neighborhood Γ i (v) is the set of vertices at geodetic distance i to v. The distance-i graph of Γ is the graph on the same vertex set where two vertices are adjacent if and only if the geodetic distance between them in Γ is i. A graph is distance-regular (DR) if for every pair or vertices ū and v at distance i apart the quantities

b i = |Γ i+1 (ū) ∩ Γ(v)|, c i = |Γ i-1 (ū) ∩ Γ(v)|,
which are referred to as the intersection numbers of the graph, solely depend on i and not on the special choice of the pair (ū, v). The automorphism group of a graph is the set of permutations of the vertices that preserve adjacency. Given a group G and an inversion-closed set G of non-identity elements of G, the Cayley graph Cal(G, G) is the graph on the vertex set G where two elements ū and v are adjacent if and only if ūv -1 ∈ G.

The Hamming graph H(n, q) is a distance-regular graph on the set Σ n of n-words over the alphabet Σ of size q, two vectors being adjacent if they are at Hamming distance one. We mostly consider the case when q is a prime power, Σ is the finite field F q of order q, and Σ n is an n-dimensional vector space over F q . The weight of a vertex of H(n, q) is the number of nonzero elements in the corresponding tuple. The weight distribution of a set C of vertices of H(n, q) is the sequence W = (W 0 , W 1 , . . . , W d ), where W i is the number of words of weight i in C. Sometimes, we will write the weight distribution in the form {i

W i 0 0 , i W i 1 1 , . . . , i W i k k }
, where W i 0 , . . . , W i k are the nonzero terms in W ; for example, (1, 0, 0, 7, 7, 0, 0, 1) = {0 1 , 3 7 , 4 7 , 7 1 }.

Equitable partitions, completely regular codes

A partition (P 0 , . . . , P r ) of the vertex set of a graph Γ is called an equitable partition (also known as regular partition, partition design, perfect coloring) if there are constants S ij such that

|Γ(v) ∩ P j | = S ij
for every v ∈ P i , i, j = 0, . . . , r.

The numbers S ij are referred to as the intersection numbers, and the matrix (S ij ) r i,j=0 as the intersection (or quotient) matrix. If the intersection matrix is tridiagonal, then P 0 is called a completely regular code of covering radius r and intersection array (S 01 , S 12 , . . . , S r-1 r ; S 10 , S 21 , . . . , S r r-1 ). That is to say, a set of vertices C is a completely regular code if the distance partition (C = C (0) , C [START_REF] Bassalygo | On uniformly packed codes[END_REF] , . . . , C (r) ) with respect to C is equitable. As was proven by Neumaier [?], for a distance-regular graph, this definition is equivalent to the original Delsarte definition [?]: a code C is completely regular if the outer distance distribution (|C ∩ Γ i (v)|) i=0,1,2,... of C with respect to a vertex v depends only on the distance from v to C.

Linear and additive codes

A linear code (that is, a linear subspace of F n q ) of dimension k and minimum distance d in H(n, q) is called an [n, k, d] q code. The duality is understood with respect to the standard inner product.

Linear codes are special cases of the additive codes, which are, by definition, the subgroups of the additive group of F n q . A coset of an additive code C is any translate of C by a constant vector. A coset leader is any coset element that minimizes the weight. The weight of a coset is the weight of any of its leaders. The coset graph Γ C of an additive code C is defined on the cosets of C, two cosets being adjacent if they differ in a coset of weight one.

Lemma 1 (see, e.g., [?]). An additive code with distance at least 3 is completely regular with intersection array {b 0 , . . . , b ρ-1 ; c 1 , . . . , c ρ } if and only if the coset graph is distance-regular with intersection numbers b 0 , . . . , b ρ-1 , c 1 , . . . , c ρ .

On the other hand, Lemma 2 (see [?, Theorem 11.1.10]). If a distance-regular graph of diameter at least 3 is a Cayley graph on an elementary abelian p-group, then it is the coset graph of a linear completely regular 1-code over F p .

So, the existence of a distance-regular graph of diameter at least 3 that is a Cayley graph on an elementary abelian p-group is equivalent to the existence of a linear completely regular 1-code in a Hamming graph over F p with the same intersection array.

Weight distributions Lemma 3 ([?],[?]

). If P is an equitable partition of the Hamming graph H(n, q) with quotient matrix S, then P is an equitable partition of the distancew graph of H(n, q) with quotient matrix S (w) := K w (K -1 1 (S)), where

K k (x) := K k (x; n, q) := k j=0 (-1) j (q -1) k-j x j n -x k -j
is the Krawtchouk polynomial. In particular, K w (K -1 1 (S)) is integer. (The same is true for any distance-regular graph, with the corresponding P-polynomial.)

The sequence of matrices (S (w) ) n w=0 can be called the generalized weight distribution of the equitable partition. Indeed, if the all-zero word belongs to P i , then the number of vertices of weight w in P j is precisely the (i, j)th element of the matrix S (w) . The following corollary shows that for some putative intersection arrays we can easily, without any calculations, see that S (w) is not integer for some w.

Corollary 1. Assume that we have a completely regular code C with intersection array {β 0 , ..., β ρ-1 ; γ 1 , ..., γ ρ } in a distance-regular graph with intersection array {b 0 , ..., b D-1 ; c 1 , ..., c D }. Then for any appropriate i ∈ {1, ..., ρ} and j ∈ {0, ..., ρ -i} the product β j β j+1 ...β j+i-1 of i consequent intersection numbers (similarly, for γ j+1 ...γ j+i ) is divisible by c 1 c 2 ...c i (for Hamming graphs, by i!).

Proof. We fix any vertex v from C (j) and denote by W l the set of vertices at distance l from v. We will prove by induction that

|W i ∩ C (i) | = i l=1 β j+l-1 c l . (1)
For i = 0, it is trivially 1 = 1 (by usual convention, the result of multiplying no factors is 1). For i > 0, it is sufficient to show that

c i |W i ∩ C (i) | = β j+i-1 |W i-1 ∩ C (i-1) |. (2)
By the definition of a completely regular code, the number of edges beginning in W i-1 ∩ C (i-1) and ending in

C (i-1) is β j+i-1 |W i-1 ∩ C (i-1) |. Since (C (t) ) t is a distance partition, all these edges end in W i .
By the definition of a distance-regular graph, the number of edges beginning in

W i ∩ C (i) and ending in W i-1 is c i |W i ∩ C (i) |. Since {C (t) } t is a distance partition, all these edges end in C (i-1) .
So, by double-counting the edges connecting W i-1 ∩ C (i-1) and W i ∩ C (i) we get (??). By induction, we have (??). Since the left part of (??) is integer, the claim of the corollary is now obvious. For {36, 28, 4; 1, 2, 24}, since 36 is the degree (q -1)n of the Hamming graph H(n, q), n ≥ 3, we have q ∈ {2, 3, 4, 5, 7, 10, 13}. In all the cases, direct calculations show that K 3 (K -1 1 (S)) is not integer, where K w (•) = K w (•; n, q) and S is the quotient matrix corresponding to the array {36, 28, 4; 1, 2, 24}.

For {63, 48, 10; 1, 8, 54}, the proof is similar.

5. {42,30,12;1,6,28}, q = 7

A putative distance-regular graph with intersection array {42, 30, 12; 1, 6, 28} has order 343 = 7 [START_REF] Brouwer | Parameters of strongly regular graphs[END_REF] . If such a graph can be realized as the coset graph of a linear completely regular code in a Hamming graph, then it is a 7ary code of length 7 = 42/(6 -1), dimension 4 = 7 -3, distance 3, and covering radius 3. The weight distributions of the code and of its dual are (1, 0, 0, 42, 42, 630, 840, 846) and (1, 0, 0, 0, 42, 0, 210, 90), respectively.

Proposition 1. If q is an odd prime power larger than 4, then there are no q-ary codes of length 7, dimension 3, and non-zero weights 4, 6, 7.

Proof. Assume that such a code C exists. W.l.o.g., it has a generator matrix   1 0 0 ? ? ? ? 0 1 0 ? ? ? ? 0 0 1 ? ? ? ?

 

We note that the weight of each row can only be 4, and two rows do not have two common zero positions (otherwise, for q > 4, they can be combined in a weight-5 codeword contradicting the hypothesis). So, w.l.o.g., we have Multiplying rows by coefficients, we get

  -1 0 0 0 -1 -1 -1 0 a -1 0 1 0 1 1 0 0 x -1 1 1 0 1   .
Now, the sum of the rows (-1, a -1 , x -1 , 2, 0, 0, 1) has weight 5, a contradiction.

Corollary 3. There are no linear completely regular codes with intersection array {42, 30, 12; 1, 6, 28} in H(7, 7).

We currently do not have a theoretical proof of the last proposition. The nonexistence of a code with the specified parameters was checked using the software [?] (1 hour of computation) and [?] (10 sec. of computation).

Corollary 5. There are no linear completely regular codes with intersection array {56, 42, 20; 1, 6, 28} in H(28, 3).

The following straightforward corollary of the fact above has an easy computer-free proof. Corollary 6. There are no linear completely regular codes with intersection array {56, 42, 20; 1, 6, 28} in H(7, 9).

Proof. If such a code exists, then its dual weight distribution is {0 1 , 4 56 , 6 392 , 7 280 }, which is impossible by Proposition ??.

8. {140,126,20,1;1,10,126,140}, q = 3

The only putative distance-regular graph in the [BCN] table with primepower order and more than 3 has order 2187 = 3

Conclusion

We have shown that no distance-regular graphs with an unsolved intersection array from the table [?, Ch. 14] can be constructed as a Cayley graph on an elementary abelian group, or, equivalently, as the coset graph of a linear completely regular code.

Finally, we observe that the existence of unrestricted (not necessarily linear) completely regular codes with the parameters considered in Sections ??-?? and the existence of few-distance codes with dual parameters remain unsolved.

  If b = c, then there is a weight-5 linear combination of the first two rows, a contradiction. So, b = c; similarly, y = z; and, if we assume w.l.o.g. that a = b, we also see x = y:

7 .Since C ( 4 )Corollary 7 .

 747 and intersection array {140, 126, 20, 1; 1, 10, 126, 140}. The corresponding completely regular code C in H(70, It is easy to see that C ∪ C (4) is a completely regular code with intersection matrix consists of 3 cosets of C, the code C ∪ C (4) is linear; its coset graph is a strongly regular graph whose parameters are also questionable, according to [?]. The dual code is a ternary two-weight code with weight distribution {0 1 , 45 588 , 54 140 }. However, a [70, 6, 45] 3 code does not exist [?]. There are no linear completely regular codes with intersection array {140, 126, 20, 1; 1, 10, 126, 140} or {140, 126; 1, 30} in H(70, 3).

Acknowledgments

The authors thank Sasha Kurz for helpful discussions and computer help.

The work of M. J. Shi is supported by the National Natural Science Foundation of China (12071001) and the Excellent Youth Foundation of Natural Science Foundation of Anhui Province (1808085J20); the work of D. S. Krotov is supported within the framework of the state contract of the Sobolev Institute of Mathematics (FWNF-2022-0017).

6. {52,42,16;1,6,28}, q = 5 A putative distance-regular graph with intersection array {52, 42, 16; 1, 6, 28} has order 625 = 5 [START_REF] Brouwer | Distance-Regular Graphs[END_REF] . If such a graph can be realized as the coset graph of a linear completely regular code in a Hamming graph, then it is a 5-ary code of length 13 = 52/(5 -1), dimension 9 = 13 -4, distance 3, and covering radius 3. The weight distributions of the code and of its dual are (1, 0, 0, 52, 260, 2028, 10660, 48620, 128076, 305240, 479596, 521040, 350480, 107072) and (1, 0, 0, 0, 0, 0, 0, 52, 0, 0, 364, 0, 208, 0) = {0 1 , 7 52 , 10 364 , 12 208 }, respectively. Proposition 2. There is no 5-ary code of length 13, dimension 4, and nonzero weights 7, 10, 12.

Proof. Let G be a 4 × 13 generator matrix of a putative code with considered parameters. The columns of G form a set S of 13 points of the vector space F 4 5 = GF(5) [START_REF] Brouwer | Distance-Regular Graphs[END_REF] . Any (3-dimensional) hyperplane H 0 in F 4 5 is determined by its dual vector h:

a codeword of weight 7, 10, or 12, any such hyperplane intersects with S in 6 = 13 -7, 3 = 13 -10, or 1 = 13 -12 points. Now, consider a (2-dimensional) plane P containing at least two points of S. It is included in exactly six 3-dimensional hyperplanes, call them H i , i = 0, 1, 2, 3, 4, 5. In the rest of the proof, we will show that any intersection numbers of S with P and H i , i = 0, 1, 2, 3, 4, 5, are contradictory. Denote a := |S ∩ P |, a > 1, and b

where k is the number of i in {0, 1, 2, 3, 4, 5} such that b i = 6. We derive k = 5(a -1)/3. Since k is integer, we have a = 4 and k = 5. However, a = 4 also implies that b i cannot be 3 and hence k = 6, a contradiction. Corollary 4. There are no linear completely regular codes with intersection array {52, 42, 16; 1, 6, 28} in H(13, 5).

7. {56,42,20;1,6,28}, q = 3 A putative distance-regular graph with intersection array {56, 42, 20; 1, 6, 28} has order 729 = 3 [START_REF] Krotov | On weight distributions of perfect colorings and completely regular codes[END_REF] . If such a graph can be realized as the coset graph of a linear completely regular code in a Hamming graph, then it is a 3-ary code of length 28 = 56/(3 -1), dimension 22 = 28 -6, distance 3, and covering radius 3. The weight distributions of the dual code is {0 1 , 12 56 , 18 392 , 21 280 }. Proposition 3. There is no 3-ary code of length 28, dimension 6, and nonzero weights 12, 18, 21.


come

Identification of deformable droplets from boundary measurements: the case of non-stationary Stokes problem

Christian Daveau, Stefan Bornhofen


Daniele Astolfi

Lorenzo Marconi

Andrew Teel

Low-power peaking-free high-gain observers for nonlinear systems

In this note the low-power observer structure presented in [1] is modified to avoid peaking of the observer state while preserving its good sensitivity properties to measurement noise and the implementation advantages in terms of powergain with respect to standard high-gain observers.

I. INTRODUCTION

High-gain observers have been extensively used in nonlinear control (see, for instance, the survey [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF] and the references therein) for their tunability property, namely the fact that the rate of convergence of the observer can be tuned. This important feature is motivated by the use of observers in output feedback control and it has been proved (see, among the others, the milestones [START_REF] Atassi | A separation principle for the stabilization of a class of nonlinear systems[END_REF], [START_REF] Teel | Global stabilizability and observability imply semi-global stabilizability by output feedback[END_REF]) that this tunability property plays a key role in establishing a nonlinear separation principle. However, the standard high-gain paradigm has many problems when used in real applications. Mainly:

1) the maximum gain to implement is increasing polynomially with the system dimension n, i.e. we need to implement a term n , where is the high-gain parameter. When nonlinear systems are considered, the value of the gain has to be taken large enough to dominate the nonlinear terms and, therefore, if the Lipschitz constant is very large, or the system dimension is high, the term n can be incredibly large; 2) high-gain observers are typically characterized by high sensitivity to high-frequency measurement noise by thus making their use practically impossible in a realistic noisy environment; 3) during the transient the variables present a peaking phenomenon which grows polynomially in , i.e. the value of the variables have an order of magnitude O( n-1 ). Many works have been published trying to overcome some of the above points. Among the others, the works [START_REF] Astolfi | Output feedback stabilization for SISO nonlinear systems with an observer in the original coordinates[END_REF], [START_REF] Maggiore | A separation principle for a class of non uniformly completely observable systems[END_REF], [START_REF] Prieur | Hybrid high-gain observers without peaking for planar nonlinear systems[END_REF] presented different solutions to deal with the peaking phenomenon. Nevertheless it is worth noting that these techniques are quite complicated to implement. More recently, the work presented in [START_REF] Teel | Further variants of the Astolfi/Marconi high-gain observer[END_REF] solves the peaking problem with a "nested-saturation" design.

A recent publication [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] presented a new "low-power" high-gain observer substantially overtaking the first two a Daniele Astolfi is with CASY -DEI, University of Bologna, Viale del Risorgimento 2, Bologna 40123, Italy, and MINES ParisTech, PSL Research University, CAS -Centre automatique et systémes, Paris 75006, France (email: daniele.astolfi@unibo.it).

b Lorenzo Marconi is with CASY -DEI, University of Bologna, Viale del Risorgimento 2, Bologna 40123, Italy (e-mail: lorenzo.marconi@unibo.it).

c Andrew Teel is with Electrical and Computer Engineering Department, University of California, Santa Barbara, CA 93106-9560 USA (e-mail: teel@ece.ucsb.edu).

problems. In particular the structure of the new observers is implemented by only using and 2 no matter the dimension n of the system is. Furthermore, the relative degree between the measure and the estimates of the new structure is higher than the ones of classical high-gain observers, resulting in a remarkable improvement of the sensitivity to high-frequency noise. The only price to pay is that the dimension of the new observer is 2n -2 and no more n. Furthermore, it can be easily proved that during the transient the values of the variables peak with an order of magnitude O( n-1 ).

The aim of this work is to present a modified version of the observer [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] which solves the three problems presented above and which is easy to implement and tune. The observer [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] is modified by adding saturations at various steps. The new structure preserves the benefits of the low-power observer presented in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF]. Furthermore, the values of the internal variables present a peaking whose order of magnitude is only O( ). This, from a practical point of view, is acceptable because the gains and 2 need to be implemented and therefore internal variables of order O( ) are numerically compatible. We observe that there are common aspects between the proposed observer design and the design in [START_REF] Teel | Further variants of the Astolfi/Marconi high-gain observer[END_REF], though the analysis in [START_REF] Teel | Further variants of the Astolfi/Marconi high-gain observer[END_REF] is complicated by the nested saturation structure. Furthermore, while the work in [START_REF] Teel | Further variants of the Astolfi/Marconi high-gain observer[END_REF] addresses the problem of peaking only indirectly (indeed, it was not the focus of [START_REF] Teel | Further variants of the Astolfi/Marconi high-gain observer[END_REF]), in this work we provide a rigorous peaking analysis.

A simulation result completes the presentation of the new observer, confirming the benefits of the proposed structure with respect to the observer proposed in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] and with respect to standard high-gain observers. Finally a new constructive Lemma (different from the one proposed in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF]) for the choice of the parameters is given in the Appendix.

Notation

We denote with R the set of real numbers. We denote with N, the set of positive integer i.e. N = {0, 1, . . . , ∞}.

Throughout the text, the variables i ∈ N, j ∈ N are used as positive integer indexes. The matrices

A ∈ R 2×2 , B ∈ R 2×1 , B i ∈ R i×1 , C ∈ R 1×2 and D 2 ( ) ∈ R 2×2 are defined as A := 0 1 0 0 , B := 0 1 , C := 1 0 , B i := 0 0 . . . 0 1 , D 2 ( ) := 0 0 2 , where ∈ R. Given the coefficients (k i1 , k i2 ) ∈ R 2 , i > 0 we define the following matrices K i ∈ R 2×1 , E i ∈ R 2×2 , Q i ∈ R 2×2 and N ∈ R 2×2 as K i := k i1 k i2 , E i := A -K i C = -k i1 1 -k i2 0 , Q i := K i B = 0 k i1 0 k i2 , N := B B = 0 0 0 1 .
We define in a recursive manner the matrices M i as

M 1 := E 1 , M i := M i-1 Ni Qi E i , i > 1,
with Ni := B 2(i-1) B , Qi := K i B 2(i-1) . The matrices M i are characterized by a block-tridiagonal structure

M i :=                E 1 N 0 . . . . . . 0 Q 2 E 2 N . . . . . . 0
. . . Q i-1 E i-1 N 0 . . . . . . . . . 0 Q i E i               
.

(1) We denote with M ∈ R (2n-2)×(2n-2) the block-tridiagonal matrix defined as M := M n-1 . We denote with Γ 1 ∈ R n×(2n-2) and Γ 2 ∈ R n×(2n-2) the following matrices

Γ 1 := blkdiag ( C, . . . , C (n-2) times , I 2 ) , Γ 2 := blkdiag (I 2 , B T , . . . , B T (n-2) times
) .

For any s ∈ R, the function sat is the function defined as sat(s) := min 1, 1 |s| s .

We denote with • a , respectively • ∞ , the asymptotic norm, respectively the infinity norm, of a time-varying signal:

s a := lim sup t→∞ |s(t)| , s ∞ := max t≥0 |s(t)| .

II. MAIN RESULTS

In this work we consider systems of the form (known also as phase-variable representation, see [START_REF] Gauthier | Deterministic observation theory and applications[END_REF])

ẋi = x i+1 ∀ i = 1, . . . , n -1 ẋn = ϕ(x, u) y = x 1 (2)
where x = (x 1 , . . . , x n ) ∈ X ⊆ R n is the state and y ∈ R is the output. The function ϕ may depend on some other known and bounded inputs u ∈ R ν . Without loss of generality and with a mild abuse of notation, we will write ϕ(x), omitting the dependence on u. The following assumption is made throughout the paper.

Assumption 1

The function ϕ is locally Lipschitz, i.e. for any compact set C ⊂ R n , there exists a L ϕ > 0 such that

|ϕ(x) -ϕ(x)| ≤ L ϕ |x -x| for all (x, x) ∈ C × C.
It is worth noting that this is the standard context in which high-gain observer are used (see, for instance [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF]). The class of system considered includes all the nonlinear systems which are uniformly observable and that can be put in the form (2) by a proper global diffeomorphism (see, for instance, [START_REF] Bornard | A high gain observer for a class of uniformly observable systems[END_REF], [START_REF] Gauthier | Observability for any u(t) of a class of nonlinear systems[END_REF], [START_REF] Gauthier | Deterministic observation theory and applications[END_REF]). We further assume that the set X where the state evolves is a compact set of R n .

A. The Low-Power High-Gain Observer

The structure of the observer proposed in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] has the following form

ξi = A ξ i + N ξ i+1 + D 2 ( ) K i e i i = 1, . . . , n-2 ξn-1 = A ξ n-1 + B ϕ s (x ) + D 2 ( ) K (n-1) e n-1 (3) with ξ = col(ξ 1 , . . . , ξ n-1 ) ∈ R 2n-2 , x = Γ 1 ξ e 1 = y-Cξ 1 , e i = B ξ i-1 -Cξ i i = 2, . . . , n-1 ,
and ϕ s (•) is an appropriate saturated version of ϕ(•) that agrees with ϕ(•) on X. The variable x represents an asymptotic estimate of the state x of (2). It is obtained by "extracting" n components from the state ξ according to the matrix Γ 1 defined above. The redundancy of the observer can be used to extract from ξ an extra state estimation that is x = Γ 2 ξ.

The following theorem, taken from [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF], shows that the observer (3) recovers the same asymptotic properties for the two estimates x and x of the "standard" high-gain observer. In the statement of the theorem we let

x = col(x , x ) , x = col(x, x) .
Theorem 1 Consider system (2) and the observer (3) with the coefficients (k i1 k i2 ) fixed so that the matrix M defined in (1) is Hurwitz (see Lemma 1 in the Appendix). Let ϕ s (•) be any bounded function that agrees with ϕ(•) on X. Then there exist c 1 > 0, c 2 > 0 and ≥ 1 such that for any ≥ and for any ξ(0) ∈ R 2n-2 the following bound holds

|x(t) -x(t)| ≤ c 1 n-1 exp(-c 2 t) |x(0) -x(0)| , ( 4
)
for all t ≥ 0 such that x(t) ∈ X.
The previous observer is proved to work in output feedback stabilization and principle of separations as shown for instance in [START_REF] Wang | Output Stabilization for a Class of Nonlinear Systems Via High-Gain Observer with Limited Gain Power[END_REF]. Furthermore, in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] the observer (3) has been shown also to improve the sensitivity to high-frequency noise with respect to standard implementation.

B. Saturated low-power high-gain observer

Having a look to the proof of Theorem 1 presented in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF], it is not so difficult to see that during the transient we have that the variables ξ i1 "peak" as O( i-1 ) whereas the variables ξ i2 "peak" as O( i ). This phenomenon can be prevented by modifying the low-power high-gain observer (3) as follows

ξi = A ξ i + B η i+1 + D 2 ( ) K i (η i-1 -Cξ i ) (5) for i = 1, . . . , n -1, with ξ = col(ξ 1 , . . . , ξ n-1 ) ∈ R 2n-2 , x = Γ 1 ξ, x = Γ 2 ξ, η 0 = y η i = L i+1 sat B ξ i L i+1 , i = 1, . . . , n -1, η n = L n+1 sat ϕ s (x ) L n+1 , (6)
and where the saturations level L i are defined as

L i := r i + ,
with a small positive real number and where we define the values r i > 0, i = 1, . . . , n + 1 as

r i := max x∈X |x i | , ∀ i = 1, . . . , n , r n+1 := max x∈X |ϕ(x)| . (7)
It can be proved that the observer (5) has the same asymptotic properties of the observer (3) and moreover that the peaking of the variables ξ i1 does not grow with , whereas the variables ξ i2 peak as O( ), with thus evident benefits in term of implementation. Moreover the variables x i for i = 1, . . . , n -1 present a peaking phenomenon which does not grow with , with only the n-th estimate of x affected by a peaking of order O( ).

The tuning of the parameters k i1 and k i2 , i = 1, . . . , n-1, must be done according to the forthcoming definition. With respect to the choice required in Theorem 1 in which the parameters were just asked to make the matrix M Hurwitz (see Lemma 1), in this case the presence of saturation functions asks the parameters also fulfill additional "common-Lyapunov" conditions.

Definition 1

The coefficients (k i1 , k i2 ), i = 1, . . , n-1, are said to be admissible if, for any i = 2, . . . , n -1, there exists a P i = P i > 0 satisfying

P i M i + M i P i < 0 , P i F i + F i P i < 0 , (8)
where

M i = M i-1 Ni Qi E i , F i := M i-1 0 Qi E i .
The set of admissible coefficients is denoted by G a .

According to Lemma 2 in the Appendix the set G a is nonempty. Therein, moreover, a constructive design procedure for the coefficients k i1 , k i2 , i = 1, . . . , n -1, is given. Theorem 2 Consider system (2) and the observer (5) and let the coefficients (k i1 k i2 ) be admissible according to Definition 1. Then, for any compact sets X ⊂ R n and Ξ ⊂ R 2n-2 and for any > 0, there exist µ 1 > 0, µ 2 > 0, δ 1 > 0, δ 2 > 0 and ≥ 1 such that, for any > and for any (x(0), ξ(0)) ∈ X × Ξ, the following bounds hold

|x(t) -x(t)| ≤ min µ 1 n-1 exp(-µ 2 t) |x(0) -x(0)| , δ 1 , ( 9
)
|x i (t)| ≤ δ 2 i = 1, . . . , n -1 , (10)
for all t ≥ 0 and as long as x(t) ∈ X.

By going throughout the proof of the theorem, it is easy to see that the observer (5) asymptotically coincides with the observer (3), namely the steady-state behaviour is not affected by the saturations. As a consequence the new design does not alter the sensitivity property of the observer (3).

III. PROOF OF THEOREM 2

The proof is divided in two parts. First we prove the asymptotic convergence of the observer. In this proof we get a conservative bound. Then, we make a more detailed analysis to show that this bound can be refined and that, indeed, the state variables are "peaking-free". In the following we define v i , i = 1, . . . , n -1 as

v i := η i+1 -x i+2 ∀ i = 1, . . . , n -2 v n-1 := η n -ϕ(x) .

Note that

|η i (t)| ∞ ≤ L i+1 , |v i (t)| ∞ ≤ r i+2 + L i+2 . ( 11
)
Furthermore, let the variables ξi , ζ i and ε i be defined as

ξi := ξ i -col(x i , x i+1 ) , ζ i := D 2 ( ) -1 ξi , ε i := 2-i D 2 ( ) -1 ξi . (12)
Finally, let R i > 0, i = 1, . . . , n be the numbers such that

| ξi (0)| ≤ R i , ∀ x ∈ X, ∀ ξ ∈ Ξ . (13)
A. Exit from saturation and convergence a) Case i = 1: Consider the first variable ξ 1 . Its dynamic is given by

ξ1 = Aξ 1 + Bη 2 + D 2 ( )K 1 (y -Cξ 1 ) .

By using the change of coordinates (12) the latter is trans

- formed into ε1 = E 1 ε 1 + -1 Bv 1
Let 1 ≥ 1 be chosen according to Lemma 4 with κ = 1 for some T 1 > 0 and 1 < . We obtain

| ε 1 (t)| ≤ 1 ∀ t ≥ T 1 ,
for all ≥ 1 Furthermore, by recalling (for any ≥ 1)

|ξ 12 | ≤ |ξ 12 -x 2 | + |x 2 | ≤ |ε 12 | + |x 2 | ≤ |ε 1 | + |x 2 |
we get

|ξ 12 (t)| ≤ |ε 1 | + |x 2 | ≤ 1 + r 2 ≤ L 2
for any ≥ 1 and for all t ≥ T 1 . With this in mind, we have that η 1 = ξ 12 for all t ≥ T 1 .

b) General Case i > 1: After time T i-1 > 0 we get η j = ξ j2 for j = 1, . . . , i-1. So consider the cascade system

ξ1 = A ξ 1 + N ξ j+1 + D 2 ( )K 1 (y -Cξ 1 ) . . . ξi-1 = A ξ i-1 + B η i + D 2 ( ) K i-1 (B ξ i-2 -Cξ i-1 ) ξi = A ξ i + Bη i+1 + D 2 ( ) K i-1 (B ξ i-1 -Cξ i )
By using the change of coordinates ( 12) the latter is transformed into

ε[1,i-1] εi = M i-1 0 Qi E i ε [1,i-1] ε i + 1 i B 2(i-1) 0 0 B v i-1 (t) v i (t) with the notation ε [1,i-1] = (ε 1 , . . . , ε i-1 ) . Note that v i-1 (t) = η i -x i+1 = L i+1 sat i ε i2 + x i+1 L i+1 -x i+1
By the Lipschitz mean-value theorem, for each t sufficiently large, there exists a continuous function

ρ(t) ∈ [0, 1] such that v i-1 (t) = ρ(t) i ε i2 .
As a consequence we get

ε[1,i-1] εi = Λ i (t) ε [1,i-1] ε i + 1 i 0 B v i (t)
where Λ i (t) is defined as in (16) in the Appendix. Let i ≥ i-1 be chosen according to Lemma 3 and 4 with κ = i and some T i > T i-1 and i < . We obtain

| i ε [1,i] (t)| ≤ i ∀ t ≥ T i ,
Furthermore, by recalling (for any ≥ 1)

|ξ i2 | ≤ |ξ i2 -x i+1 | + |x i+1 | ≤ i |ε i2 | + |x i+1 | ≤ i |ε [1,i] | + |x i+1 | we get |ξ i2 (t)| ≤ i |ε [1,i] | + |x i+1 | ≤ i + r i+1 ≤ L i+1
for any ≥ i and for all t ≥ T i . This prove that for t ≥ T i also η i = ξ i2 . c) Asymptotic Convergence: Finally, after the (n -1)th step we recover the observer [START_REF] Atassi | A separation principle for the stabilization of a class of nonlinear systems[END_REF]. By direct application of Theorem 1 we thus have that if is taken sufficiently large ( ≥ n ) then the following estimate holds

|x(t) -x(t)| ≤ µ 1 n-1 exp(-µ 2 t) |x(0) -x(0)| (14)
for some strictly positive real numbers µ 1 , µ 2 and for all t ≥ T n-1 . Note that the estimate bound n-1 is too conservative. Indeed, it can be proved that during the transient the peaking depends only on and not on n-1 , as shown in the next part of the proof.

B. Bound Estimate

Consider the ξ i dynamic. By using the change of coordinates [START_REF] Teel | Global stabilizability and observability imply semi-global stabilizability by output feedback[END_REF], it is transformed into

ζi = E i ζ i + -1 v i + K i-1 η i
(except for the ζ 1 for which the term η 1 is not present). Now consider the Lyapunov function V i defined as V i := η i S i η i with S i given by

S i E i + E i S i = -I, we get Vi ≤ -|ζ i | 2 + 2|ζ i ||S i |( -1 v i + |K i-1 |η i-1 ) .
Then, standard Lyapunov arguments show that there exist positive constants a ij , j = 1, . . . , 4, such that, if is taken sufficiently large, we get

|ζ i (t)| ≤ a i1 exp(-a i2 t)|ζ i (0)| + a i3 2 v i + a i4 |K i-1 |η i-1
for any i = 1, . . . , n -1. Furthermore, by recalling (for any ≥ 1)

| ξi1 | ≤ |ζ i | ≤ | ξi | , 1 | ξi2 | ≤ |ζ i | ≤ | ξi | , we get | ξi1 (t)| ≤ a i1 e -ai2 t | ξi (0)| + -2 a i3 v i + a i4 |K i-1 |η i-1 , | ξi2 (t)| ≤ a i1 e -ai2 t | ξi (0)| + -1 a i3 v i + a i4 |K i-1 |η i-1 , for all t ≥ 0. By using |x -x| ≤ n-1 i=1 | ξi1 | + | ξi2 | and by
defining δ 1 as

δ 1 := n-1 i=1 2a i1 R i + 2a i3 (r i+2 + L i+2 ) + a i4 |K i-1 |L i
with R i defined in [START_REF] Wang | Output Stabilization for a Class of Nonlinear Systems Via High-Gain Observer with Limited Gain Power[END_REF], we get |x(t) -x(t)| ≤ δ 1 for all t ≥ 0. By combining the latter with (14) we finally get [START_REF] Maggiore | A separation principle for a class of non uniformly completely observable systems[END_REF]. The bound [START_REF] Prieur | Hybrid high-gain observers without peaking for planar nonlinear systems[END_REF] follows by by noting

|x i | ≤ |ξ i1 -x i | + |x i | , i = 1, . . . , n -1,
and by setting

δ 2 := max i∈[1,n-1] r i +a i1 R i +a i3 (L i+2 +r i+2 )+a i4 L i |K i-1 | .

IV. SIMULATIONS

We consider as example a forced Van Der Pol oscillator, described by

ẅ + b(w 2 -1) ẇ + w = α cos(ωt + φ)
where b > 0. The system can be described in the phasevariable representation (2) with n = 4, By taking ω = 2 and α ≤ 1, it can be verified that, for any φ ∈ R, the solutions of the systems are ultimately bounded. Taking b = 0.01, a numerical simulation can be used to show that the solutions satisfy

ϕ(x) = -6bx 1 x 2 x 3 -2bx 3 2 -x 3 -b(x 2 1 -1)x 4 -ω 2 (x 3 + b(x 2 1 -1)x 2 + x 1 ) and initial conditions x 1 (0) = w(0) , x 2 (0) = ẇ(0) , x 3 (0) = -b(w(0) 2 -1) ẇ(0) -w(0) + α cos(φ), x 4 (0) = -2bw(0) ẇ(0) 2 -ẇ(0) -b(w(0) 2 -1)x 3 (0) + αω sin(φ).
x 1 a < 2.2 , x 2 a < 2.6 , x 3 a < 3 , x 4 a < 4.6 , ϕ (x
) a < 6.9 .

In the following we consider any compact set X of initial conditions satisfying (see also ( 7)) r 1 < 2.2 , r 2 < 2.6 , r 3 < 3 , r 4 < 4.6 , r 5 < 6.9 .

By following the procedure indicated by Lemma 1, we chose the coefficients of the observer as

K 1 = 3 6.4 , K 2 = 3 2.131 , K 3 = 3 0.7095 , such that the poles of M 3 are in (-2, -1.8, -1.6, -1.4, -1.2, -1).
It is easy to verify that they are also admissible, according to Definition 1. The saturations level are chosen as L 1 = 2.5 , L 2 = 3 , L 3 = 3.5 , L 4 = 5 , L 5 = 7.5 .

In the simulations we compared the observer (3) with the new observer [START_REF] Gauthier | Observability for any u(t) of a class of nonlinear systems[END_REF] for different choices of . The initial conditions of the forced Van Der Pol are chosen as w(0) = 1, ẇ(0) = -1, γ = 1, φ = 0, that is x(0) = (1, -1, 0, 0.98), whereas the two observers are initialized in the origin. The tables show the maximum peaking values in the error coordinates and the time needed to convergence to an error sufficiently small, i.e. the time T such that

|x(t) -x(t)| < ∀ t ≥ T .
In the tables we denote

x i = x i -x i , x i = x i -x i .
Note that x 1 = x 1 and x 4 = x 4 . The big numerical improvement during the transient, given by the saturations in the proposed design [START_REF] Gauthier | Observability for any u(t) of a class of nonlinear systems[END_REF], it is evident from the simulation, especially for the variables x 3 , x 3 and x 4 . The Table II confirms the theoretical result of Theorem 2. In particular the variables x 1 , x 2 and x 3 are always bounded with a bound independent of the high-gain parameter . Finally, it can be noticed from the simulations that the when the high-gain parameter is small, the speed of convergence of the novel observer ( 5) is worsen with respect to the standard lowpower high-gain observer (3), whereas when the high-gain parameter is high the rate of convergence is improved. This phenomenon is caused by the effect of the saturations which prevent the state of the observer to have very large transients.

V. CONCLUSION

By following the ideas in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF], we proposed a modified version of the low-power high-gain observer introduced therein. The new structure is designed by adding saturations at each block. The new design preserves the same benefits of the observer proposed in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] and overtakes the peaking phenomenon. The proposed observer has a dimension which is 2n -2 (where n is the system dimension) with a highgain parameter which grows only up to power 2 (whereas in standard high-gain observers the parameter is powered up to n) and the peaking phenomenon exhibits with a growth O( ) on the last component only (and not O( n-1 ) as in the observer of [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF]), giving evident benefits with respect to the standard implementations of high-gain observers. A new constructive Lemma for the choice of the parameter is given in the appendix. The new procedure may result easier with respect to the one proposed in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF]. Due to the asymptotic tunable convergence property [START_REF] Maggiore | A separation principle for a class of non uniformly completely observable systems[END_REF], the novel proposed observer can be used in any output feedback design where a standard high-gain observer applies (see, for instance, [START_REF] Teel | Global stabilizability and observability imply semi-global stabilizability by output feedback[END_REF] and [START_REF] Wang | Output Stabilization for a Class of Nonlinear Systems Via High-Gain Observer with Limited Gain Power[END_REF]).

APPENDIX

The eigenvalues of M i , i > 0, can be arbitrarily assigned by appropriately choosing the coefficients (k j1 , k j2 ), j = 1, . . . , i, as claimed in the next lemma (taken off-the-shelf from the Lemma 1 published in [START_REF] Astolfi | A High-Gain Nonlinear Observer with Limited Gain Power[END_REF] by setting i = n -1).

Lemma 1 Let

P Mi (λ) = λ 2i + m 1 λ 2i-1 + ... + m 2i-1 λ + m 2i
be an arbitrary Hurwitz polynomial. There exists a choice of k j1 , k j2 , j = 1, . . . , i, such that the characteristic polynomial of M i coincides with P Mi (λ).

Lemma 2

The set G a is non-empty.

Proof: We show constructively, though conservatively, a recursive procedure to assign the coefficients k i1 , k i2 .

Step i = 1) Let k 11 > 0 and k 12 > 0 be any (positive) real numbers.

Step i > 1) Let k i1 = k (i-1)1 , and let k i2 > 0 be chosen such that k i2 < k i1 /γ i-1 , with γ i-1 defined in the forthcoming equation (15).

To verify that this choice satisfies the Definition 1 let consider first the case i = 1. By choosing k 11 > 0 and k 12 > 0 the matrix E 1 is Hurwitz, and so is M 1 . Now consider the case i > 1 and the following two systems ẋi-1

y i-1 = = M i-1 x i-1 + B 2(i-1) u i-1 B 2(i-1) x i-1 ẋi y i = = E i x i + K i u i B x i
and let denote xi = (x i-1 , x i ). Let assume the matrix M i-1 is Hurwitz and note that the matrix E i is Hurwitz for any choice of positive real numbers k i1 , k i2 . Let H i-1 (s), H i (s) be the transfer function of the two systems H i-1 (s) := B 2(i-1) (sI 2(i-1) -M i-1 ) -1 B 2(i-1) , H i (s) := B (sI 2 -E i ) -1 K i . and let γ i-1 , γ i be the real numbers defined as

γ i-1 := max s∈R |H i-1 (s)| , γ i := max s∈R |H i (s)| . ( 15
)
It is easy to verify that

H i (s) = sk i2 s 2 + k i1 s + k i2 , γ i = k i2 k i1 .
By interconnecting the two subsystems with u i-1 = 0 and u i = y -1 we obtain a closed loop system described by ẋi = F i xi , whereas, by interconnecting the two subsystems with u i-1 = y i and u i = y -1 , we obtain a closed loop system described by ẋi = M i xi .

The choice of k i2 of the Step i satisfies the small-gain theorem (see, for instance, [START_REF] Isidori | Nonlinear control systems II[END_REF])

γ i-1 • γ i < 1
and therefore the interconnection (x i-1 , x i ) is stable in both cases. We conclude that F i and M i are Hurwitz and that there exists a common Lyapunov function satisfying [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF].

Note that, by using MATLAB, the values of γ i can be easily calculated with the command getPeakGain.

Lemma 3 Let the coefficients k j1 , k j2 , j = 1, . . . , i be admissible and let the matrix Λ i (t) be defined as

Λ i (t) := M i-1 ρ(t) Ni-1 Qi E i ( 16
)
with ρ(t) ∈ [0, 1] any continuous function. Then the origin of the system ẋ = Λ i (t) x (17

)
is globally exponentially stable.

Proof: When ρ(t) ≡ 0 we get Λ i = F i whereas when ρ(t) ≡ 1 we get Λ i = M i . As a consequence, by using the P i claimed in the Definition 1 we get the inequality P i Λ i (t) + Λ i (t) P i < 0 is satisfied for all ρ(t) ∈ [0, 1] and therefore the claim follows trivially.

Lemma 4 Consider the system ẋ = Λ i (t) x + -κ v i (t) where x ∈ R 2i , κ ≥ 0 and v i (t) satisfying |v i (t)| ∞ ≤ vi for some vi > 0. Suppose the origin of the system (17) is globally asymptotically stable and locally exponentially stable. Then, for any compact set X ⊂ R 2i , for any T > 0 and for any > 0, there exist a a ≥ 1 such that

| κ x(t)| ≤ , ∀ t ≥ T
for any > and for any initial condition x(0) ∈ X.

Proof: Trivial, by applying standard ISS-Lyapunov arguments.

TABLE I :

 I Features of the low-power high-gain observer (3).

		= 5	= 10	= 100	= 1000
	T 0.01	4.9898	1.8382	0.2544	0.0399
	x 1 ∞	1	1	1	1
	x 2 ∞	7.5	14.4	139	1.4 • 10 3
	x 2 ∞	7	13.6	133	1.3 • 10 3
	x 3 ∞	16.5	70	7.2 • 10 3	7.2 • 10 5
	x 3 ∞	14.1	60	6.4 • 10 3	6.4 • 10 5
	x 4 ∞	16.2	126	1.3 • 10 5	1.3 • 10 8

TABLE II :

 II Features of the low-power peaking-free high-gain observer[START_REF] Gauthier | Observability for any u(t) of a class of nonlinear systems[END_REF].

		= 5	= 10	= 100	= 1000
	T 0.01	5.6072	2.0222	0.2364	0.0312
	x 1 ∞	1	1	1	1
	x 2 ∞	7.2	13.4	127	1300
	x 2 ∞	4.3	4.4	4.4	4.4
	x 3 ∞	9.3	21	261	2600
	x 3 ∞	3.7	3.9	3.8	3.8
	x 4 ∞	4.7	9.4	92	900

ACKNOWLEDGMENT

The authors thanks Laurent Praly for a number of discussions and inputs that have influenced the outcome of this paper.


Soumaya Oueslati
email: soumaya.oueslati1@u-cergy.fr

Imen Balloumi
email: imen.balloumi@u-cergy.fr

Christian Daveau
email: christian.daveau@u-cergy.fr

Abdessatar Khelifi
email: abdessa-tar.khelifi@fsb.rnu.tn

Analytical method for the evaluation of singular integrals arising from boundary element method in electromagnetism

Keywords: Maxwell equations, boundary integral equations, singular integral, homogeneous functions AMS Subject Classification:

In the framework solving electromagnetic scattering problems with a Leontovich impedance boundary condition (LIBC), the numerical resolution requires the evaluation of singular integrals appearing in the discretization of the variational formulation. Our main interest is to pricisely evaluate these integrals. Thus, we propose an analytic method to approximate it. The performance of this method will be evaluated by calculating the radar cross section (RCS). Then, we compare RCS to the Mie series solution for the unit sphere in various configurations.

Introduction

Integral equation methods are an alternative for solving linear partial differential equations. It consists in transforming the partial differential equation posed in a domain, into an integral equation posed on the boundary. It is in particular a method adapted for diffraction problems that are posed in an unbounded domain (see [START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF] [START_REF] Bendali | Numerical analysis of the exterior boundary value problem for the time-harmonic Maxwell equations by a boundary finite-element method?Part 2: The discrete problem[END_REF]) and for which the consideration of the radiation condition is essential in order to insure a unique solution. We are interested in solving time-harmonic scattering problem by the integral method and we will take the constant impedance condition, known as standard or Leontovich impedence boundary condition (LIBC) [START_REF] Senior | Impedance boundary conditions for imperfectly conducting surfaces[END_REF]. This approximation does not depend on incident angle at all. It has been extensively used to simplify the formulation involved in the solution of complicated electromagnetic scattering problems. It was recognized that this type of boundary condition can be advantageously used to get a more tractable problem in numerous complex situations of electromagnetic scattering computations (see for example [START_REF] Artola | Diffraction d?une onde lectromagntique par un obstacle born permittivit et permabilit leves[END_REF] [7] [START_REF] Jingguo | Modified surface impedance boundary conditions fo 3-D eddy currents problems[END_REF]). Once the boundary problem has been converted into a boundary integral equation, it is necessary to discretize the boundary of the domain using plan polygons. Here we deal with plan triangles and Rao-Wilton-Glisson basis functions.

In this paper, we are particularly interested in the singular integrals resulting from the discretization of the variationnel problem of Maxwell 3-D equation. This singular character of the integrand makes the precise evaluation of these integrals difficult. It has been the subject of a large number of publications. M. Lenoire and N. Salles have developed in [START_REF] Salles | Evaluation of 3-D singular and nearly singular integrals in Galerkin BEM for thin layers[END_REF] [START_REF] Salles | Exact Evaluation of Singular and nearsingular integrals in Galerkin BEM[END_REF] a new approach to evaluate singular integrals. It allows to evaluate analytically and with great precision, singular integrals in the variational case. We adopt the same approach to evaluate a singular integral coming from a double layer operator.

To better understand this singular integrand problem, in the variational case, it should be noted that the domain of integration for these integrals with "problem", is the product of two elements of the mesh, precisely the product of two triangles. When the two triangles are distant, it is possible to evaluate the integral with standard numerical integration methods. But, when the intersection of the two triangles is not empty(adjacent triangles and triangles with a common vertex), it require separate treatment and this is a field of research still very active.

The approach developed in this work, makes it possible to reduce a 4-D integral to a linear combination of mono-dimensional integrals whose integrand is regular. It is then possible to numerically evaluate these 1-D integrals, but also explicitly. This method relies on several formulae that make it possible to reduce the dimension of the integration domain from n to n -1. So, repeating the same process, three times using these formulas, we obtain regular 1-D integrals. These reduction formulae are based on the homogeneity of the integrand; therefore, there are two conditions to check. The integrand must be homogeneous and we must use flat polygons to discretize the boundary of the domain in order to maintain the homogeneity of the integrand during the successive stages of the reduction process.

Scattering problem

We consider the scattering problem of electromagnetic waves (E, H) by a perfect conducting body with a complex coating. We refer by Ω -to the region of space embodying the scatterer, with a Lipschitz-continuous boundary denoted by Γ and we refer by n is the exterior unit normal to Γ, pointing to the exterior. Waves propagation medium is characterized by two physical quantities ϵ (electrical permittivity) and µ (magnetic permeability). The propagation of this wave is described in Ω + = R 3 \Ω -by the following Maxwell equations:

(1)

{ rotE + iωµH = 0, rotH -iωϵE = 0.
In order to ensure a unique solution of boundary value problem, it is necessary to apply the following boundary condition

(2)

Etg -Z(n × H) = 0 on Γ.
Where Z is the impedance operator that depends on incident angle, subscript tg denotes tangent component on the surface defined as:

Etg = n × (E × n).
Generally, we take the constant impedance operator, known as standard or Leontovich impedance boundary conditions Z = const [START_REF] Rahmat-Samii | Hoppe Impedance Boundary Conditions in Electromagnetics[END_REF][START_REF] Senior | Approximate boundary conditions in electromagnetics[END_REF][START_REF] Lange | Equations intégrales espace-temps pour les équations de Maxwell : calcul du champ diffracté par un obstacle dissipatif[END_REF]. We should notice that asymptotic behavior of the fields (E, H) is depicted by the Silver-Müller radiation condition:

(3) lim r→∞ r(E × nr + H) = 0.
where r = |x| and nr = x |x| , x ∈ R 3 . We introduce current densities J and M on the boundary Γ as follows

M = E × n , J = n × H.
We determine total electromagnetic fields (E, H) in Ω + as:

(4)

{ E = E inc + E sc in Ω + , H = H inc + H sc in Ω + .
Subscripts inc and sc characterize incident and scattered fields, respectively. Afterwards, using the Stratton-Chu formulation we can give the variational formulation of the problem. For all Ψ = (Ψ J , Ψ M ), find (J, M) such that:

(5)

< Z 0 (B -S)J, Ψ J > + < (P + Q)M, Ψ J >=< E inc , Ψ J >, ( 6
) -< (P + Q)J, Ψ M > + < 1 Z 0 (B -S)M, Ψ M >=< H inc , Ψ M > . Where (7) ⟨(B -S)ϕ, ψ⟩ = i ∫ ∫ Γ k G ϕ • ψ - 1 k G ∇y • ϕ ∇x • ψ dydx, (8) ⟨(P + Q)ϕ, ψ⟩ = 1 2 ∫ Γ ψ • (n × ϕ) dx + ∫ ∫ Γ (ψ × ϕ) • ∇xG dydx, (9) ⟨ϕ, ψ⟩ = ∫ Γ A • ψ dx,
with G(x, y) = e -ik∥x-y∥ 4π∥x -y∥ is the Green kernel in three dimensions.

For a numerical resolution of the problem we will adopt the finite element method of boundary. We assume that the initial boundary problem is defined in a 3-D domain. The boundary Γ of the domain is therefore a surface and can be discretized using triangles. We approach the surface of the obstacle by a surface Γ h composed of finite number of two dimensional elements. These elements are triangular facets denoted by T i for i = 1 to N T :

Γ h = N T ∪ i=1 T i .
We will call it an initial mesh (or original mesh).We denote by Ne the total number of edges of the mesh component Γ h . The discretization of unknowns J and M should verify a condition -flow conservation of these currents. One way to ensure this is to use the basis functions of Rao-Wilton-Glisson introduced below. We introduce local numbering of a triangle T . The vertices (a T j ) j=1,3 are arranged in clockwise order. Triangle edges are numbered so that the edge T ′ j connects vertices a T j and a T j+1 . Moreover, we give an orientation νn to each edge n. Consider the two triangles sharing this edge. We note T + n the triangle so that the direction of the edge n coincides with the forward direction (locally defined) of this triangle. For the other triangle, which we will denote T - n , the direction of the edge coincides with the indirect sense. Each basis function is associated with an edge and ensures the conservation of flux through this edge. If we denote |T | the area of a triangle T , the n th basis function is defined as follows:

Definition 2.1. If n is the i th local edge of triangle T + n and the j th of triangle T - n then (10) fn(x) =        ln 2|T + n | (x -a + i-1 ) if x ∈ T + n ln 2|T - n | (a - j-1 -x) if x ∈ T - n 0 if x / ∈ T + n ∪ T - n
The density is proportional to ∇ Γ • fn, where

(11) ∇ Γ • fn(x) =        + ln |T + n | on T + n -ln |T - n | on T - n 0 elsewhere.
We decompose the electric and magnetic currents:

J(y) = Ne ∑ i=1 J i f i (y), M(y) = Ne ∑ i=1 M i f i (y).
Considering to these assumptions we are looking for approximated solution of the problem [START_REF] Darrigrand | Operateurs intgraux, noyaux singuliers, mthodes rapides : Outils mathematiques pour la rsolution numerique de problmes denses et de grande dimension[END_REF][START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF]. We introduce the following system

(12) A h (U h , Ψ h ) = Ne ∑ i=1 < E inc , f i > + Ne ∑ i=1 < H inc , f i >,
where

A h (U h , Ψ h ) = Ne ∑ i,j=1 < Z 0 (B -S)f j , f i > J j + Z -1 0 Ne ∑ i,j=1 < (B -S)f j , f i > M j + Ne ∑ i,j=1 < Qf j , f i > M j - Ne ∑ i,j=1 < Qf j , f i > J j + a 0 2 Ne ∑ i,j=1 < f j , f i > J j + 1 2a 0 Ne ∑ i,j=1 < f j , f i > M j .
Since the basis function f i is defined on two triangles, namely adjacent to the i th edge , T + i and T - i . The same for the basis function f j , they are defined on two triangles, T + j and T - j . We define the matrix elements

(B -S) i,j = i ∫ T + i ∪T - i ∫ T + j ∪T - j k G(x, y)f j (y) • f i (x) - 1 k G(x, y)(∇ Γ • f i (x))(∇ Γ • f j (y)) dxdy.
Analogically

Q i,j = -i ∫ T + i ∪T - i ∫ T + j ∪T - j [f i (x) × f j (y)] • ∇ Γ • G(x, y) dydx. with G(x, y) = e -ik∥x-y∥ 4π∥x -y∥ = - 1 4π∥x -y∥ + H(∥x -y∥)
is the Green kernel in three dimensions and

∇xG(x, y) = - (1 + ik|x -y|) 4π|x -y| 3 e -ik|x-y| (x -y) = - (x -y) 4π∥x -y∥ 3 - k 2 (x -y) 8π∥x -y∥ -(x -y)R(∥x -y∥).
where H and R are regular analytic functions.

We then, notice that there are singular terms that appear. The regular integrals R i,j can be calculated numerically with the usual known methods such as the Gauss method.

We denote by S the triangle T + i and by T the triangle T + j . We evaluate the operators Q i,j and (B -S) i,j on S × T and it will be the same for other combinations of triangles.. On each triangle, the current is written as a linear combination of three RWG basis functions associated to three edges of a triangle. Now we come to the integrals developed in this paper. So we have to evaluate three type of singular integrals (13)

D i,j = ∫ S×T (x -a i ) • (x -y) ||x -y|| ξ+1 × (y -b j ) dxdy; i, j = 1, 2, 3; ξ = 0, 2 (14) S i,j = ∫ S×T (x -a i ) • (y -b j ) ||x -y|| dxdy; i, j = 1, 2, 3. ( 15
) M i,j = ∫ S×T 1 ||x -y|| dxdy; i, j = 1, 2, 3.
In this work, we adopt a purely analytic method for the evaluation of the singular integrals. It was developed in [START_REF] Salles | Evaluation of 3-D singular and nearly singular integrals in Galerkin BEM for thin layers[END_REF] [START_REF] Salles | Exact Evaluation of Singular and nearsingular integrals in Galerkin BEM[END_REF]. In these works, the term (15) was be evaluated and studied theoretically.

Definitions Tools and Notations

We present some notations and we give tools that we will need in our calculation.

Notations

Let the triangles S and T be defined (respectively) by its vertices a i and b j , i, j = 1, 2, 3 and its sides α i respectively β j , i, j = 1, 2, 3 such as α i and β j are opposed to a i respectively b j (see Figure 1). We introduce |α i | and |β j | the length of the side α i respectively β j . Let ⃗ α i and ⃗ β j be the vectors (a i+2 - 2). We denote by λ i and ϑ j the outside normal vector to α i respectively β j . Note p i (x) the orthogonal projection (see Figure 3) of a point x on the side α i of the triangle S and g i (x) = ||x -p i (x)|| the distance from x to the support of α i . When x = a i , we note that p i = p i (a i ) and that g i = g i (a i ) is the length of the height (see Figure 4). We will also use the distance signed from a point to the support of α i which involves the normal outgoing to the triangle on α i . This distance is defined by γ i = (p i (x) -x| ⃗ λ i ) (see Figure 5). We use notations for the triangle T similar to those introduced for the triangle S. The orthogonal projection of y on the support of β j , similar to p i (x), is denoted by q j (y) and the distance from y to the support of β j is denoted by d j (y) = ||y -q j (y)||. We define the signed distance δ j (y) from y to the support of β j by δ j (y) = (q j (y) -y| ⃗ ϑ j ). Finally, d j = d j (b j ) is the equivalent of g i , namely the length of the height from b j to β j .

a i+1 ) = (a + i -a - i ) respectively (b i+2 -b i+1 ) = (b + i -b - i ) (see Figure
On side α (resp β) the abscissa s (resp t) is defined with respect to an origin oα (resp o β ) and a unitary

direction vector ⃗ α |α| (resp ⃗ β |β|
). The abscissas of the ends a ± and b ± are respectively denoted by s ± and t ± .

s k = (a k -o| ⃗ α |α| ) t l = (b l -o| ⃗ β |β| ).

Reduction process

The reduction process is based on several simplification formulas presented in section . These formulas reduce the size of the domain integration. They can be used when the integrand is homogeneous, or positively homogeneous. Definition 3.1. A function is said to be homogeneous of degree q, or q-homogeneous, with respect to the variable z (16) f (λz) = λ q f (z). Definition 3.2. A function is said to be q-homogeneous with parameter if it is an homogeneous function with respect to a pair of variables (z, h), where z is the integration variable and h a real parameter [START_REF] Lange | Equations intégrales espace-temps pour les équations de Maxwell : calcul du champ diffracté par un obstacle dissipatif[END_REF] f (λz, λh) = λ q f (z, h). Theorem 1. Let f : Ω ⊂ R n → R a positively homogeneous function of degree q with q + n ̸ = 0 (18)

∫ Ω f (z)dz = 1 q + n ∫ ∂Ω ( -→ z | -→ ν )f (z) dsz,
where -→ ν is the exterior normal to ∂Ω and ( -→ z | -→ ν ) is the scalar inner product between z and -→ ν . Theorem 2. Let g : Ω ⊂ R n × R → R, be a homogeneous function of degree q, referred to as a homogeneous function with parameter. Then, provided that

(19) ( h h 0 ) q+n ∫ Ω g(z, h 0 )dsz -→ 0 quand h 0 -→ sgn(h)∞,
we get

(20) ∫ Ω g(z, h)dz = ∫ ∂Ω ( -→ z | -→ ν )F (z, h) dsz,
where

F (z, h) = h q+n ∫ sgn(h)∞ h g(z,u) u q+n+1 du. Theorem 3. Let f : Ω ⊂ R n → R
a positively homogeneous function of degree q with q = -n. We introduce an arbitrary homogeneous positive function θ : Ω → R of degree p > 0. Then we have

(21) ∫ Ω f (z)dz = 1 p ∫ ∂Ω ( -→ z | -→ ν )f (z)ln(θ(z)) dsz.
These formulae were proved in ( [START_REF] Salles | Evaluation of 3-D singular and nearly singular integrals in Galerkin BEM for thin layers[END_REF]).

In the following we present a method for the explicit evaluation of singular integrals. It is built on a recursive reduction of a 4-dimension integral into a linear combination of (m-1)-dimensional integrals. It leads to a linear combination of 1-dimensional regular integrals. The method we have used requires to distinguish between several geometric configurations and bases to the choice of the origin during the steps of reduction of the dimension of the integration domain. This technique is assured the homogeneity of the singular part of the Green function.

Application of the reduction formula

We propose to evaluate the integral

(22) I(x, y) = ∫ S×β j x -y ||x -y|| 3 dx dy,
where S is a triangle, x a point in the triangle S, β j is a segment in the same plan as S and y is a point in β j . Let b + j = b j+1 and b - j = b j+2 be the two extremities of β. The reduction process can be applied, provided that one takes some common point to β j and S as origin, so that the function f (z) = x-y ||x-y|| 3 is homogeneous on S × β j with q = 2, n = 3. The origin can be taken anywhere on the support of β j , but a wise choice is to choose one of the ends of β j as origin. Indeed, some of the distances being canceled, one will have a simplification. We choose for origin the point b -. Formula (18) provides ( 23)

I(x, y) = 1 5 ∫ ∂(S×β j ) ((x, y)| -→ ν ) x -y ||x -y|| 3 ∂(x, y),
where ν is the outward normal at the border of the edge of S × β j and ∂(x, y) is the surface element along

∂(S × β j ). We have ∂(S × β j ) = (∂S × β j ) ∪ (S × ∂β j ). ((x, y)| -→ ν ) b + j ×S = |β j |; ((x, y)| -→ ν ) |α i ×β j = γ i (b - j )
where |β j | denotes the length of the segment β j . Then, we get

(24) I(x, y) = |β j | 5 ∫ S x -b + j ||x -b + j || 3 dsx + 3 ∑ i=1 γ i (b - j ) 5 ∫ ∂(α i ×β j )
x -y ||x -y|| 3 dsx dsy.

So we get a linear combinaison of 2-D integrals.

Calculation of the singular integral in the forme double layer operator

Before explicitly evaluating the integral defined in (13), we will rewrite the RWG basis functions as a combination of the barycentric basis functions which allowed to transform the integral into integrals simpler to study. We introduce the basis functions ϕ i and ϕ j associated with the vertices of the triangle S respectively T . We use their properties (see appendix 6.1) above to determine θ, κ, ζ, ρ, η and µ verifying

(25) (x -a i ) = θϕ i (x) + κϕ i+1 (x) + ζϕ i+2 (x), (y -b j ) = ρψ j (y) + ηψ j+1 (y) + µψ j+2 (y).
We obtain

(26) (x -a i ) = (a i+1 -a i )ϕ i+1 (x) + (a i+2 -a i )ϕ i+2 (x), (y -b j ) = (a j+1 -a j )ψ j+1 (y) + (a j+2 -a j )ψ j+2 (y).
We are interested in calculating the singular integral [START_REF] Sauter | Boundary Element Methods[END_REF]. Using the basis functions as defined in (26), this integral becomes (27)

D i,j = ∫ S×T (κ × η) • (ϕ i+1 (x)ψ j+1 (y)) (x-y) ||x-y|| ξ+1 dxdy + ∫ S×T (κ × µ) • (ϕ i+1 (x)ψ j+2 (y)) (x-y) ||x-y|| ξ+1 dxdy + ∫ S×T (ζ × η) • (ϕ i+2 (x)ψ j+1 (y)) (x-y) ||x-y|| ξ+1 dxdy + ∫ S×T (ζ × µ) • (ϕ i+2 (x)ψ j+2 (y)) (x-y)
||x-y|| ξ+1 dxdy.

So we have to evaluate four integrals of the same type:

(28)

J i,j = ∫ S×T (ϕ i (x)ψ j (y)) (x -y) ||x -y|| ξ+1 dxdy; ξ = 0, 2; i, j = 1, 2, 3.
There is two different geometric configurations. Here, we deal with the case of adjacent triangles. We detail the formula for the evaluation of the expressions (28). In the case triangles have a common vertex, there no real additional difficulty. We can apply the same techniques.

Since the basis functions is affine, the integrant is no longer homogeneous, so we are led to decompose these into a linear part and a constant part by making a change of origin:

(29) ϕ i (x) = 1 + ( ---→ o -a i | -→ e i ) + ( ---→ x -o| -→ e i ) = ϕ i (o) + (x| -→ e i ); i = 1, 2, 3, ψ j (y) = 1 + ( ---→ o -b j | -→ l j ) + ( ---→ y -o| -→ l j ) = ψ j (o) + (y| -→ l j ); j = 1, 2, 3.
where o is the new origin, x = x -o and y = y -o. More generally, we note α and S the segment α and the triangle S after this change of origin. According to 29 we constate that the integrals J i,j must be decomposed into four integrals whose integrands are homogeneous J i,j

(30)

J i,j = ∫ S×T (x| - → e i )(y| - → l j )(x-y) ||x-y|| 1+ξ dxdy + ϕ i (o) ∫ S×T (x| - → e i )(y| - → l j )(x-y) ||x-y|| 1+ξ dxdy +ψ j (o) ∫ S×T (x| - → e i )(x-y) ||x-y|| 1+ξ dxdy + ψ j (o)ϕ i (o) ∫ S×T (x-y) ||x-y|| 1+ξ dxdy.
Once this decomposition has taken place, we obtain homogeneous integrands we can then apply the reduction method intermediate of the formulas (18), (20) or (21).

Three-dimensional reduction

We are dealing with the case of triangles with an edge in common, so there are two common vertices. For the first reduction, we take for origin one of the common vertices o = a i = b j . Let a i+1 = b j+1 be the second common vertex and a 1+2 ,b j+2 be the remaining vertices. We propose to calculate the integral (28). It will be wise to choose the vertex a i = b j , as origin o. In this case the basis functions are written (31)

ϕ i (x) = 1 + (x| -→ e i ), ψ j (y) = 1 + (y| -→ l j )
The four integrands are homogeneous of degrees q = 1 -ξ, q = 2 -ξ, q = 1 -ξ and q = -ξ with n = 4. We apply the formula (18), it gives the following result

J i,j = 3 ∑ i=1 γ i (a i ) 5 -ξ ∫ α i ×T (y| -→ l j )(x -y) ||x -y|| 1+ξ dsxdy + 3 ∑ j=1 δ j (b j ) 5 -ξ ∫ S×β j (y| -→ l j )(x -y) ||x -y|| 1+ξ dxdsy + 3 ∑ i=1 γ i (a i ) 6 -ξ ∫ α i ×T (x| -→ e i )(y| -→ l j )(x -y) ||x -y|| 1+ξ dsxdy + 3 ∑ j=1 δ j (b j ) 6 -ξ ∫ S×β j (x| -→ e i )(y| -→ l j )(x -y) ||x -y|| 1+ξ dxdsy + 3 ∑ i=1 γ i (a i ) 5 -ξ ∫ α i ×T (x| -→ e i )(x -y) ||x -y|| 1+ξ dsxdy + 3 ∑ j=1 δ j (b j ) 5 -ξ ∫ S×β j (x| -→ e i )(x -y) ||x -y|| 1+ξ dxdsy (32) + 3 ∑ i=1 γ i (a i ) 4 -ξ ∫ α i ×T x -y ||x -y|| 1+ξ dsxdy + 3 ∑ j=1 δ j (b j ) 4 -ξ ∫ S×β j x -y ||x -y|| 1+ξ dxdsy.
where γ i (x) and δ j (y) are the signed distances respectively of x to the support of α i and of y in support of β j (see Figure 4). In particular γ i (x) = 0 for all x belongs to the support of α i . γ i (a i ) = g i and δ j (b j ) = d j are the lengths of the height from the vertex a i or b j to the support of α i (respectively β j ). However the formula is simplified and we get (33)

J i,j = g i 5-ξ ∫ α i ×T (y| - → l j )(x-y) ||x-y|| 1+ξ dsxdy + d j 5-ξ ∫ S×β j (y| - → l j )(x-y)
||x-y|| 1+ξ dxdsy

+ g i 6-ξ ∫ α i ×T (x| - → e i )(y| - → l j )(x-y) ||x-y|| 1+ξ dsxdy + d j 6-ξ ∫ S×β j (x| - → e i )(y| - → l j )(x-y) ||x-y|| 1+ξ dxdsy + g i 5-ξ ∫ α i ×T (x| - → e i )(x-y) ||x-y|| 1+ξ dsxdy + d j 5-ξ ∫ S×β j (x| - → e i )(x-y)
||x-y|| 1+ξ dxdsy

+ g i 4-ξ ∫ α i ×T x-y ||x-y|| 1+ξ dsxdy + d j 4-ξ ∫ S×β j x-y ||x-y|| 1+ξ dxdsy.
We express the results in terms of ϕ i and ψ j in order to benefit the best of the properties (94) of the basic functions.

(x| -→ e i ) = ϕ i (x) -1, (y| -→ l j ) = ψ j (y) -1
we obtain

J i,j = g i 5-ξ ∫ α i ×T ψ j (y)(x-y) ||x-y|| 1+ξ dsxdy -g i 5-ξ ∫ α i ×T (x-y) ||x-y|| 1+ξ dsxdy + d j 5-ξ ∫ S×β j ψ j (y)(x-y) ||x-y|| 1+ξ dxdsy - d j 5-ξ ∫ S×β j (x-y) ||x-y|| 1+ξ dxdsy + g i 6-ξ ∫ α i ×T ϕ i (x)ψ(y)(x-y) ||x-y|| 1+ξ dsxdy -g i 6-ξ ∫ α i ×T ϕ i (x)(x-y) ||x-y|| 1+ξ dsxdy -g i 6-ξ ∫ α i ×T ψ j (y)(x-y) ||x-y|| 1+ξ dsxdy + g i 6-ξ ∫ α i ×T (x-y) ||x-y|| 1+ξ dsxdy + d j 6-ξ ∫ S×β j ϕ i (x)ψ(y)(x-y) ||x-y|| 1+ξ
dxdsy -

d j 6-ξ ∫ S×β j ϕ i (x)(x-y)
||x-y|| 1+ξ dxdsy -

d j 6-ξ ∫ S×β j ψ j (y)(x-y) ||x-y|| 1+ξ dxdsy + d j 6-ξ ∫ S×β j (x-y) ||x-y|| 1+ξ dxdsy + g i 5-ξ ∫ α i ×T ϕ i (x)(x-y) ||x-y|| 1+ξ dsxdy -g i 5-ξ ∫ α i ×T (x-y) ||x-y|| 1+ξ dsxdy + d j 5-ξ ∫ S×β j ϕ i (x)(x-y) ||x-y|| 1+ξ dxdsy - d j 5-ξ ∫ S×β j (x-y) ||x-y|| 1+ξ dxdsy + g i 4-ξ ∫ α i ×T x-y ||x-y|| 1+ξ dsxdy + d j 4-ξ ∫ S×β j x-y ||x-y|| 1+ξ dxdsy.
There are terms that vanish because of the fact that ϕ i|α i = ψ i|β j = 0. We have integrals in three dimensions we reapply the procedure to reduce dimention to two. The integrands are no more homogeneous. It is possible to continue the reduction by taking the second common vertex as the origin a 2 = b 2 = o. We reexpress the basis function as in (31) and we obtain

ϕ i (x) = ϕ i (a i+1 ) + (x| -→ e i ) = (x| -→ e i ); ψ j (x) = ψ j (b j+1 ) + (y| -→ l j ) = (y| -→ l j )
Finally, we get

J i,j = ( |S| 5-ξ - |S| 6-ξ ) 2 4-ξ P ψ j (a i+2 , T, ξ) + ( g i 5-ξ -g i 6-ξ ) d j+1 4-ξ Q(α, β j+1 , ξ) +( |S| 6-ξ -2 |S| 5-ξ + |S| 4-ξ ) 2 3-ξ P (a i+2 , T, ξ) + ( g i 6-ξ -2 g i 5-ξ + g i 4-ξ ) d j+1 3-ξ Q(α i , β j+1 , ξ) -( |T | 5-ξ - |T | 6-ξ ) 2 4-ξ P ϕ i (b j+2 , S, ξ) -( d 1 5-ξ - d j 6-ξ ) g i+1 4-ξ Q 1,ϕ i (β j , α i+1 , ξ) -( |S| 6-ξ - |S| 5-ξ ) 2 4-ξ P (b j+2 , S, ξ) + ( d j 4-ξ + d j 6-ξ - d j 5-ξ ) g i+1 4-ξ Q(α i+1 , β j , ξ)
, where

P (a, T, ξ) = ∫ T (a -y) ||a -y|| 1+ξ dsxdsy, Q(α i , β j , ξ) = ∫ α i ×β j (x -y) ||x -y|| 1+ξ dsxdsy, P ψq (a, T, ξ) = ∫ T ψq(y)(a -y) ||a -y|| 1+ξ dy, Q 1,ψq (α i , β j , ξ) = ∫ α i ×β j ψq(y)(x -y) ||x -y|| 1+ξ dsxdsy.
We give also the other terms calculated with the same procedure:

(34)

J i,j+1 = P (a i+2 , T, ξ)( 2|S| (5-ξ)(3-ξ) - 2|S| (5-ξ)(4-ξ) - 2|S| (6-ξ)(3-ξ) + 2|S| (6-ξ)(4-ξ) ) +P (b j+2 , S, ξ)(- 2|T | (5-ξ)(3-ξ) + 2|T | (5-ξ)(4-ξ) + 2|T | (6-ξ)(3-ξ) - 2|T | (6-ξ)(4-ξ) ) +Q(α i , β i+1 , ξ)( g i d j+1 (5-ξ)(3-ξ) - g i d j+1
(5-ξ)(4-ξ) -

g i d j+1 (6-ξ)(3-ξ) + g i d j+1 (6-ξ)(4-ξ) ) +Q(α i+1 , β j , ξ)( d j g i+1 (5-ξ)(3-ξ) - d j g i+1 (6-ξ)(3-ξ) + d j g i+1 (6-ξ)(4-ξ) - d j g i+1 (5-ξ)(4-ξ) ) +Q 1,ψ j+1 (α i+1 , β j , ξ)( d j g i+1 (5-ξ)(4-ξ) - d j g i+1 (6-ξ)(4-ξ) ) + P ψ j+1 (a i+2 , T, ξ)( 2|S| (5-ξ)(4-ξ) - 2|S| (6-ξ)(4-ξ) ) -Q 1,ϕ i (β j , α i+1 , ξ)( d j g i+1 (6-ξ)(4-ξ) - d j g i+1 (6-ξ)(5-ξ) ) + P ϕ i (b j+2 , S, ξ)(- 2|T | (6-ξ)(4-ξ) + 2|T | (6-ξ)(5-ξ) ) +Q ϕ i ,ψ j+1 (α i+1 , β j , ξ)( d j g i+1 (6-ξ)(5-ξ) ), where (35) Q ϕp,ψq (α i , β j , ξ) = ∫ α i ×β j ϕp(x)ψq(y)(y -x) ||x -y|| 1+ξ dsxdsy, ( 36
)
J i,j+2 = P (b j+2 , S, ξ)(- 2|T | (5-ξ)(4-ξ) + 2|T | (6-ξ)(4-ξ) ) -P ϕ i (b j+2 , S, ξ)( 2|T | (6-ξ)(5-ξ) ) +P ψ j+2 (a i+2 , T, ξ)( 2|S| (5-ξ)(4-ξ) - 2|S| (6-ξ)(4-ξ) ) + Q 1,ψ j+2 (α i , β i+1 , ξ)( d j+1 g i (5-ξ)(4-ξ) - d j+1 g i (6-ξ)(4-ξ) ) +Q 1,ψ j+2 (α i+1 , β j , ξ)( d j g i+1 (5-ξ)(4-ξ) - d j g i+1 (6-ξ)(4-ξ) ) + Q ϕ i ,ψ j+2 (α i+1 , β j , ξ)( d j g i+1 (6-ξ)(5-ξ) ), (37) J i+1,j+1 = P (a i+2 , T, ξ)( 2|S| (6-ξ)(3-ξ) - 4|S| (6-ξ)(4-ξ) + 2|S| (5-ξ)(6-ξ) ) + P ψ j+1 (a i+2 , T, ξ)( 2|S| (6-ξ)(4-ξ) - 2|S| (6-ξ)(5-ξ) ) +Q(α i , β i+1 , ξ)( g i d j+1 (6-ξ)(3-ξ) - 2g i d j+1 (6-ξ)(4-ξ) + g i d j+1 (6-ξ)(5-ξ) ) + Q 1,ϕ i+1 (β i+1 , α i , ξ)(- d j+1 g i (6-ξ)(4-ξ) + d j+1 g i (6-ξ)(5-ξ) ) +P (b j+2 , S, ξ)(- 2|T | (6-ξ)(3-ξ) + 4|T | (6-ξ)(4-ξ) - 2|T | (6-ξ)(5-ξ) ) + P ϕ i+1 (b j+2 , S, ξ)( 2|T | (6-ξ)(5-ξ) - 2|T | (6-ξ)(4-ξ) ) +Q(α i+1 , β j , ξ)( d j g i+1 (6-ξ)(3-ξ) - 2d j g i+1 (6-ξ)(4-ξ) + d j g i+1 (6-ξ)(5-ξ) ) + Q 1,ψ j+1 (α i+1 , β j , ξ)(- d j g i+1 (6-ξ)(5-ξ) + d j g i+1 (6-ξ)(4-ξ) ), ( 38
)
J i+1,j+2 = P (b j+2 , S, ξ)( 2|T | (6-ξ)(5-ξ) - 2|T | (6-ξ)(4-ξ) ) + P ψ j+2 (a i+2 , T, ξ)( 2|S| (6-ξ)(4-ξ) - 2|S| (6-ξ)(5-ξ) ) +P ϕ i+1 (b j+2 , S, ξ)(- 2|T | (6-ξ)(5-ξ) ) + Q 1,ψ j+2 (α i , β i+1 , ξ)( d j+1 g i (6-ξ)(4-ξ) - d j+1 g i (6-ξ)(5-ξ) ) +Q 1,ψ j+2 (α i+1 , β j , ξ)( d j g i+1 (6-ξ)(4-ξ) - d j g i+1 (6-ξ)(5-ξ) ) + Q ϕ i+1 ,ψ j+2 (α i , β i+1 , ξ)( d j+1 g i (6-ξ)(5-ξ) ),

and

(39)

J i+2,j+2 = P ψ j+2 (a i+2 , T, ξ)( 2|S| (6-ξ)(5-ξ) ) + P ϕ i+2 (b j+2 , S, ξ)(- 2|T | (6-ξ)(5-ξ) ) +Q ϕ i+2 ,ψ j+2 (α i , β i+1 , ξ)( d j+1 g i (6-ξ)(5-ξ) ) + Q ϕ i+2 ,ψ j+2 (α i+1 , β j , ξ)( d j g i+1 (6-ξ)(5-ξ) ).
We can also give the symmetrical terms:

(40)

J i+1,j = P (a i+2 , T, ξ)( 2|S| (5-ξ)(3-ξ) - 2|S| (5-ξ)(4-ξ) - 2|S| (6-ξ)(3-ξ) + 2|S| (6-ξ)(4-ξ) ) +P (b j+2 , S, ξ)(- 2|T | (5-ξ)(3-ξ) + 2|T | (5-ξ)(4-ξ) + 2|T | (6-ξ)(3-ξ) - 2|T | (6-ξ)(4-ξ) ) +Q(α i , β i+1 , ξ)( g i d j+1 (5-ξ)(3-ξ) - g i d j+1
(5-ξ)(4-ξ) -

g i d j+1 (6-ξ)(3-ξ) + g i d j+1 (6-ξ)(4-ξ) ) +Q(α i+1 , β j , ξ)( d j g i+1 (5-ξ)(3-ξ) - d j g i+1 (6-ξ)(3-ξ) + d j g i+1 (6-ξ)(4-ξ) - d j g i+1 (5-ξ)(4-ξ) ) +Q 1,ψ j (α i , β i+1 , ξ)( d j+1 g i (6-ξ)(4-ξ) - d j+1 g i (6-ξ)(5-ξ) ) -Q 1,ϕ i+1 (β i+1 , α i , ξ)( d j g i+1 (5-ξ)(4-ξ) - d j g i+1 (6-ξ)(4-ξ) ) +P ψ j (a i+2 , T, ξ)( 2|S| (6-ξ)(4-ξ) - 2|S| (6-ξ)(5-ξ) ) + P ϕ i+1 (b j+2 , S, ξ)(- 2|T | (5-ξ)(4-ξ) + 2|T | (6-ξ)(4-ξ) ) +Q ϕ i+1 ,ψ j (α i , β i+1 , ξ)( d j+1 g i (6-ξ)(5-ξ) ), (41) J i+2,j = P (a i+2 , T, ξ)( 2|S| (5-ξ)(4-ξ) - 2|S| (6-ξ)(4-ξ) ) + P ϕ i+2 (b j+2 , S, ξ)( 2|T | (6-ξ)(4-ξ) - 2|T | (4-ξ)(4-ξ) ) +P ψ j (a i+2 , T, ξ)( 2|S| (5-ξ)(6-ξ) ) + Q 1,ϕ i+2 (β j , α i+1 , ξ)( d j g i+1 (6-ξ)(4-ξ) - d j g i+1 (5-ξ)(4-ξ) ) +Q 1,ϕ i+2 (β i+1 , α i , ξ)( d j+1 g i (6-ξ)(4-ξ) - d j+1 g i (5-ξ)(4-ξ) ) + Q ϕ i+2 ,ψ j (α i , β i+1 , ξ)( d j+1 g i (6-ξ)(5-ξ) ),

and

(42)

J i+2,j+1 = P (a i+2 , T, ξ)( 2|S| (6-ξ)(4-ξ) - 2|S| (6-ξ)(5-ξ) ) + P ϕ i+2 (b j+2 , S, ξ)( 2|T | (6-ξ)(5-ξ) - 2|T | (6-ξ)(4-ξ) ) +P ψ j+1 (a i+2 , T, ξ)( 2|S| (6-ξ)(5-ξ) ) + Q 1,ϕ i+2 (β i+1 , α i , ξ)( d j+1 g i (6-ξ)(5-ξ) - d j+1 g i (6-ξ)(4-ξ) ) +Q 1,ϕ i+2 (β j , α i+1 , ξ)(- d j g i+1 (6-ξ)(4-ξ) + d j g i+1 (6-ξ)(5-ξ) ) + Q ϕ i+2 ,ψ j+1 (α i+1 , β j , ξ)( d j g i+1 (6-ξ)(5-ξ) ).

Reduction to mono-dimentional integrals

Calculation of P (a, T, ξ)

we will evaluate 43 • Case a don't belongs to the support of T In this case the reduction formulas are not usable. In order to operate the reduction, it is necessary to break P into two parts.

(44)

P (a, T, ξ) = ∫ T (a -â) (h 2 + ||â -y|| 2 ) 1+ξ 2 dsxdsy + ∫ T (â -y) (h 2 + ||â -y|| 2 ) 1+ξ 2
dsxdsy
where â is the projection of a into the support of T and h is the distance h = |a -â|. Taking â as a new origin we obtain homogeneous integrands of degrees -1 -ξ and -ξ with parameter h. We can then apply the formula (20) which gives (45) P (a, T, ξ) = (a -â)

3 ∑ j=1 δ j (â)S 1 (â, β j , h, ξ) + 3 ∑ j=1 δ j (â)R 2 (â, β j , h, ξ),
where S 1 and R 2 are defined by

(46) S 1 (a, β, h, ξ) = ∫ β h 1-ξ ∫ +∞ h du u 2-ξ (u 2 + ||a -y|| 2 ) 1+ξ 2 dsy = ∫ β f 1 (h, ||a -y||, ξ),
and

(47) R 2 (a, β, h, ξ) = ∫ β (a -y)h 2-ξ ∫ +∞ h du u 3-ξ (u 2 + ||a -y|| 2 ) 1+ξ 2 dsy = ∫ β (a -y)f 2 (h, ||a -y||, ξ), where (48) fn(h, s, ξ) = h n-ξ ∫ +∞ h du u n+1-ξ (u 2 + s 2 ) 1+ξ 2
. S 1 (a, β, h, ξ) and R 2 (a, β, h, ξ) are calculated in (6.3) and (6.5).

• Case a belongs to the support of T Since a belongs to the support of T we can consider a as an origin and get homogeneous integrant with q = -ξ and n = 2 but we must distingue two different cases in order to apply the suitable reduction formula.

* If n + q = 0, we apply (21) and we get See appendix (6.2).

Calculation of Q(α, β, ξ)

Here we deal with

(53) Q(α i , β j , ξ) = ∫ α i ×β j (x -y) ||x -y|| 1+ξ dsxdsy.
We distinguish two cases.

• Case α and β are neither parallel nor secants

First of all, we must find a new origin to ensure homogeneity of the integrand. Let α be the projection of α on the plane parallel to α going through β and let o be the intersection of α with β. We will consider o as a new origin. Define z = ǎ -a and h = ||z||. We express Q as a sum of two integrals whose integrands are homogeneous of degrees q = -1 -ξ and -ξ.

(54) Q(α, β, ξ) = ∫ α×β (x -x) (h 2 + ||x -y|| 2 ) 1+ξ 2 dsxdsy + ∫ α×β (x -y) (h 2 + ||x -y|| 2 ) 1+ξ 2 dsxdsy,
where x is a point α. This decomposition allows us to apply the formula (20) and gives

Q(α, β, ξ) = (x -x)[ ∑ k=± ks k S 1 (ǎ k , β, h, ξ) + ∑ l=± lt l S 1 (b l , α, h, ξ)] (55) + [ ∑ k=± ks k R 2 (ǎ k , β, h, ξ) - ∑ l=± lt l R 2 (a, β, h, ξ)],
where the abscissae are evaluated by

s k = (ȃ k -o| ⃗ α |α| ) and t l = (b l -o| ⃗ β |β| ).
• Case α and β are secants Let o = α ∩ β be the new origin. Note that can be o = a i+1 = b j+1 or o = a i+2 = b j+2 . In this case the integrand is homogeneous of degree q = -ξ and n = 2. We must distinguish two cases: * If q + n ̸ = 0, the formula (18) gives

(56) Q(α, β, 0) = ∑ k=± ks k 2 S(a k , β, 0) - ∑ l=± lt l 2 S(b l , α, 0).
* If q + n = 0, we apply (21) and we get

(57) Q(α, β, 2) = ∑ k=± ks k L(a k , β) - ∑ l=± lt l L(b l , α).
• Case α and β are parallel When the supports of α and β are parallel, we must use the second formula (20), called with parameter. This is possible by orthogonally projecting one of the segments (we choose β) on the support of the other (α). Thus, the distance d between the supports of the two segments, appears and plays the role of the parameter. Let β and ỹ, be the respective orthogonal projections of β and y on the support of α, ⃗ z = ỹ -y and d = ||ỹ -y|| the distance from α to the support of β. Then, the integral Q(α, β, ξ) becomes ( 58)

Q(α, β, ξ) = ⃗ z ∫ α× β 1 (d 2 + ||x -ỹ|| 2 ) 1+ξ 2 dsxds ỹ + ∫ α× β (x -ỹ) (d 2 + ||x -ỹ|| 2 ) 1+ξ 2 dsxds ỹ .
Taking any point of α as a new origine, for exemple one of the two vertices of the triangle a -and a + , we get two homogeneous integrands with parameter d and degrees q = -1 -ξ respectively q = -ξ.

(59)

Q(α, β, ξ) = ⃗ z[s + ∫ β d 1-ξ ∫ +∞ d du u 2-ξ (u 2 + ||a + -ỹ|| 2 ) 1 2 ds ỹ + ∑ l=± lt l ∫ α d 1-ξ ∫ +∞ d du u 2-ξ (u 2 + ||x -bl || 2 ) 1 2 dsx] +s + ∫ β d 2-ξ (a + -y) ∫ +∞ d du u 3-ξ (u 2 + ||a + -y|| 2 ) 1 2 ds ỹ + ∑ l=± lt l ∫ α d 2-ξ (x -bl ) ∫ +∞ d du u 3-ξ (u 2 + ||x -bl || 2 ) 1 2 dsx, which is equivalent to (60) Q(α, β, ξ) = ∑ l=± lt l [⃗ zS 1 ( bl , β, d, ξ) -R 2 ( bl , β, h, ξ)] + s + [S 1 (a + , β, d, ξ) + R 2 (a + , β, d, ξ)].

Calculation of P ψq (a, T, ξ)

To calculate the singular integral we are inspired by the section (4.2.1). There are two different geometrical cases.

• Case a don't belongs to the support of T Note that a = a i+2 is a vetex of S. To maintain the homogeneity of the integrands the origin must belongs to the support of the triangle T . So we consider â the projection of a on the support of T (see Figure ??) and we break P ψq (a, T, ξ) into two parts (62)

P ψq (a, T, ξ) = (a -â) ∫ T ψq(y) (h 2 + ||â -y|| 2 ) 1+ξ 2 dy + ∫ T ψq(y)(â -y) (h 2 + ||â -y|| 2 ) 1+ξ ,
where h = |a -â| and ξ = 0, 2. Next, we express ψq(y) with the new origin o = â, ψq(y) = ψq(â) + (y|lq), which gives

P ψq (a, T, ξ) = (a -â)[ψq(â) ∫ T 1 (h 2 + ||â -y|| 2 ) 1+ξ 2 dy + ∫ T (y|lq) (h 2 + ||â -y|| 2 ) 1+ξ 2 dy] (63) + ψq(â) ∫ T (â -y) (h 2 + ||â -y|| 2 ) 1+ξ 2 dy + ∫ T (y|lq)(â -y) (h 2 + ||â -y|| 2 ) 1+ξ 2
dy.

We notice that the four intergrands are homogeneous with parameter h and degrees (respectively) q = -1 -ξ, q = -ξ, q = -ξ and q = 1 -ξ. We can then apply (20) and we obtain (64)

P ψq (a, T, ξ) = (a -â)ψq(â) ∑ 3 j=1 δ j (â) ∫ β j h 1-ξ ∫ +∞ h du u 2-ξ (u 2 +||â-y|| 2 ) 1+ξ 2 dy +(a -â) ∑ 3 j=1 δ j (â) ∫ β j h 2-ξ (y|lq) du u 3-ξ (u 2 +||â-y|| 2 ) 1+ξ 2 dy +ψq(â) ∑ 3 j=1 δ j (â) ∫ β j h 2-ξ (â -y) ∫ +∞ h du u 3-ξ (u 2 +||â-y|| 2 ) 1+ξ 2 dy + ∑ 3 j=1 δ j (â) ∫ β j h 3-ξ (y|lq)(â -y) ∫ +∞ h du u 4-ξ (u 2 +||â-y|| 2 ) 1+ξ 2
dy.

To simplify the formula, we express the basic functions in this form (y|lq) = ψq(y) -ψq(â) and we delete the terms that vanish, thanks to (94) (65)

P ψq (a, T, ξ) = (a -â)ψq(â) ∑ 3 j=1 δ j (â) ∫ β j f 1 (h, ||â -y||, ξ) -f 2 (h, â -y, ξ)dsy +(a -â)[δ q+1 (â) ∫ β q+1 ψq(y)f 2 (h, â -y, ξ)dsy + δ q+2 (â) ∫ β q+2 ψq(y)f 2 (h, â -y, ξ)dsy] +ψq(â) ∑ 3 j=1 δ j (â) ∫ β j (â -y)(f 2 (h, ||â -y||, ξ) -f 3 (h, â -y, ξ))dsy +δ q+1 (â) ∫ β q+1 ψq(y)(â -y)f 3 (h, ||â -y||, ξ)dsy + δ q+2 (â) ∫ β q+2 ψq(y)(â -y)f 3 (h, â -y, ξ)dsy.
otherwise we can write (66) Integrands are homogeneous of degrees q = -ξ and q = 1 -ξ. We apply (18) and we get (69) P ψq (a, T, ξ) = ψq(a)

P ψq (a, T, ξ) = (a -â)ψq(â) ∑ 3 j=1 δ j (â)(S1(â, β, h, ξ) -S2(â, β, h, ξ)) +(a -â)[δ q+1 (â)S ψq 2 (â, β q+1 , h, ξ) + δ q+2 (â)S ψq 2 (â, β q+2 , h, ξ)] +ψq(â) ∑ 3 j=1 δ j (â)R 2 (â, β j , h, ξ) -R 3 (â, β j , h, ξ) +δ q+1 (â)D ψq 3 (â, β q+1 , h, ξ) + δ q+2 (â)D
3 ∑ j=1 δ j (a) 2 -ξ ∫ β j (a -y) ||a -y|| 1+ξ dsy + 3 ∑ j=1 δ j (a) 3 -ξ ∫ β j (y|lq)(a -y) ||a -y|| 1+ξ dy.
We replace (y|lq) by (y|lq) = ψq(y) -ψq(a), the we get

P ψq (a, T, ξ) = ψq(a) 3 ∑ j=1 δ j (a) 2 -ξ ∫ β j (a -y) ||a -y|| 1+ξ dsy. (70) + 3 ∑ j=1 δ j (a) 3 -ξ ∫ β j ψq(y)(a -y) ||a -y|| 1+ξ dsy -ψq(a) 3 ∑ j=1 δ j (a) 3 -ξ ∫ β j (a -y) ||a -y|| 1+ξ dsy, which is equivalent to (71) P ψq (a, T, ξ) = 3 ∑ j=1 δ j (a) 3 -ξ S ψq (a, β j , ξ) + ψq(a) 3 ∑ j=1 [ δ j (a) 2 -ξ S(a, β j , ξ) - δ j (a) 3 -ξ S(a, β j , ξ)].

Calculation of

Q 1,ψq (α, β, ξ)
We want to evaluate the singular integral

(72) Q 1,ψq (α i , β j , ξ) = ∫ α i ×β j ψq(y)(x -y) ||x -y|| 1+ξ dsxdsy, where β ̸ = α ̸ = α i+2 , β ̸ = β i+2 ̸ = βq.
We are going to study two different cases.

• Case α and β are neither parallel nor secants Let α be the projection of α on the plane containing β and parallel to α, x a point in α, h = |a -ǎ| and

⃗ z = x -x. We decompose Q 1,ψq in two parts (73) Q 1,ψq (α, β, ξ) = ⃗ z ∫ α×β ψq(y) (h 2 + ||x -y|| 2 ) 1+ξ 2 ds xdsy + ∫ α×β ψq(y)(x -y) (h 2 + ||x -y|| 2 ) 1+ξ 2
ds xdsy .

We will follow the same approach as previously. We consider o = α ∩ β the new origin and we replace in the equation (73) the basic function ψq(y) by the equivalent expression ψq(y) = ψq(o) + (y|lq).

Once we have homogeneous integrands, the reduction formula (20) gives

Q 1,ψq (α, β, ξ) = ψq(o)[ ∑ k=± ks k R 2 (â k , β, h, ||â k -y||, ξ) -R 2 (b l , α, h, ||x -b l ||, ξ)] + ∑ k=± ks k [D ψq 3 (â k , β j , ||â k -y||, ξ) -ψq(o)R 3 (â k , β, ||â k -y||, ξ)] + ∑ l=± lt l [(-ψq(b l ) + ψq(o))R 3 (b l , α, ||x -b l ||, ξ)] +⃗ zψq(o)[ ∑ k=± ks k S 1 (â k , β, h, ξ) + ∑ l=± lt l S 1 ( bl , α, h, ξ)] +⃗ z ∑ k=± ks k [S ψq 2 (â k , β, h, ||â k -y||, ξ) -ψq(o)S 2 (â k , β, h, ||â k -y||, ξ)] +⃗ z ∑ l=± lt l [(ψq(b l ) -ψq(o))S 2 (b l , α, h, ||x -b l ||, ξ)].
• Case α and β are secants Let o = α ∩ β be the new origin. The basis function ψq(y) can be expressed as following

ψq(y) = ψq(o) + (y|lq), which gives (74) Q 1,ψq (α, β, ξ) = ψq(o) ∫ α×β (x -y) ||x -y|| 1+ξ dsxdsy + ∫ α×β (y|lq)(x -y) ||x -y|| 1+ξ dsxdsy
The integrants are homogeneous of degrees q = -ξ and q = 1 -ξ with n = 2.

* If q + n ̸ = 0 the formula (18) gives Q 1,ψq (α, β, ξ) = ψq(o)[ ∑ k=± ks k 2 S(a k , β, 0) - ∑ l=± lt l 2 S(b l , α, 0)] + ∑ k=± ks k 3 S ψq (y) (a k , β, 0) (75) -ψq(o) ∑ k=± ks k 3 S(a k , β, 0) - ∑ l=± lt l 3 (ψq(b l ) -ψq(o))S(b l , α, 0), where (76) S ψq (y) (a k , β, 0) = ∫ β ψq(y)(a k -y) ||a k -y|| dsy.
* If q + n = 0 it will be appropriate to apply the formula (21), which leads us to this expression

Q 1,ψq (α, β, ξ) = ψq(o)[ ∑ k=± ks k L(a k , β) - ∑ l=± lt l L(b l , α)] + ∑ k=± ks k L ψq (a k , β) (77) - ∑ k=± ks k ψq(o)L(a k , β) - ∑ l=± lt l (ψq(b l ) -ψq(o))L(b l , α), where (78)
L ψq (a k , β) = ∫ β ψq(y)(a k -y) ||a k -y|| 3 dsy.

Calculation of Q ϕp,ψq (α, β, ξ)

Here we deal with the singular integral

(79) Q ϕp,ψq (α, β, ξ) = ∫ α×β ϕp(x)ψq(y)(y -x) ||x -y|| 1+ξ
dsxdsy.

• Intersecting triangles

We project α on the plan parallel to α and going through β. Let α be the projection of α and x the projection of x. Then, we provide an equivalent expression of the basis functions: 2 ,

(80) ϕp(x) = 1 + µ (x -ap|⃗ α) |α| 2 = 1 + µ (x -âp| ⃗ α) |α| 2 , (81) ψq(y) = 1 + ν (y -bq)| ⃗ β |β|
where µ = ±1(µ = 1 if ap = a + and µ = -1 if ap = a -). Likewise, ν = ±1(ν = 1 if bp = b + and ν = -1 if bp = b -).
We choose the intersection of the supports of α and β, as origin and we note it o. we obtain :

ϕp(x) = 1 + µ (o -âp|⃗ α) |α| 2 + µ (x -o|⃗ α) |α| 2 (82) = Cp + µ (x -o) ⃗ α , ψq(y) = 1 + ν (o -bq| ⃗ β) |β| 2 + ν (y -o| ⃗ β) |β| 2 (83) = Cq + µ (y -o| ⃗ β) |β| 2 where Cp = 1 + µ (o-âp)|⃗ α |α| 2 and Cq = 1 + ν (o-bq )| ⃗ β |β| 2
Then the expression of Q ϕp,ψq becomes

Q ϕp,ψq (α, β, ξ) = CqCq ∫ α×β (x -y) (h 2 + ||x -y|| 2 ) 1+ξ 2 ds xdsy + Cpν |β| 2 ∫ α×β (y -o| ⃗ β)(x -y) (h 2 + ||x -y|| 2 ) 1+ξ 2 ds xdsy (84) + Cqµ |α| 2 ∫ α×β (x -o|⃗ α)(x -y) (h 2 + ||x -y|| 2 ) 1+ξ 2 ds xdsy + µµ | α| 2 |β| 2 ∫ α×β (y -o| ⃗ β)(x -o|⃗ α)(x -y) (h 2 + ||x -y|| 2 ) 1+ξ 2
ds xdsy .

In the formula below the integrands are homogeneous with variable h and of degrees q = -ξ, q = 1 -ξ, 1 -ξ and 2 -ξ. Formula (20) gives (85)

Q ϕp,ψq (α, β, ξ) = CpCq[ ∑ k=± kŝ k ∫ β (â k -y)h 2-ξ ∫ +∞ h du u 3-ξ (u 2 +||â k -y|| 2 ) 1+ξ 2 dsy + ∑ l=± l tl ∫ α(x -b l )h 2-ξ ∫ +∞ h du u 3-ξ (u 2 +||x-b l || 2 ) 1+ξ 2 dsx] + Cpν |β| 2 [ ∑ k=± kŝ k ∫ β (y -o| ⃗ β)(â k -y)h 3-ξ ∫ +∞ h du u 4-ξ (u 2 +||â k -y|| 2 ) 1+ξ 2 dsy + ∑ l=± l tl ∫ α(b l -o|⃗ α)(x -b l )h 3-ξ ∫ +∞ h du u 4-ξ (u 2 +||x-b l || 2 ) 1+ξ 2 dsx] + Cq µ |α| 2 [ ∑ k=± kŝ k ∫ β (â k -o|⃗ α)(â k -y)h 3-ξ ∫ +∞ h du u 4-ξ (u 2 +||â k -y|| 2 ) 1+ξ 2 dsy + ∑ l=± l tl ∫ α(x -o|⃗ α)(x -b l )h 3-ξ ∫ +∞ h du u 4-ξ (u 2 +||x-b l || 2 ) 1+ξ 2 dsx] + νµ |α| 2 |β| 2 [ ∑ k=± kŝ k ∫ β (â k -o|⃗ α)(y -o| ⃗ β)(â k -y)h 4-ξ ∫ +∞ h du u 5-ξ (u 2 +||â k -y|| 2 ) 1+ξ 2 dsy + ∑ l=± l tl ∫ α(x -o|⃗ α)( bl -o|⃗ α)(x -b l )h 4-ξ ∫ +∞ h du u 5-ξ (u 2 +||x-b l || 2 ) 1+ξ 2 dsx],
which is equivalent to (86)

Q ϕp,ψq (α, β, ξ) = CpCq[ ∑ k=± kŝ k R 2 (â k , β, h, ξ)ξ - ∑ l=± l tl R 2 (b l , α, h, ξ)] + ∑ k=± kŝ k [ Cq µ |α| 2 (â k -o|⃗ α)R 3 (â k , β, h, ξ) + Cpν |β| 2 M 3 (â k , o, β, h, ξ)] + ∑ l=± l tl [- Cpν |β| 2 (b l -o| ⃗ β)R 3 (b l , α, h, ξ) - Cq µ |α| 2 M 3 (b l , o, α, h, ξ)] + νµ |α| 2 |β| 2 [ ∑ k=± kŝ k (â k -o|⃗ α)M 4 (â k , o, β, h, ξ) - ∑ l=± l tl (b l -o| ⃗ β)M 4 (b l , o, α, h, ξ)].
where

(87) M 3 (â k , o, β, h, ξ) = ∫ β (â k -o|⃗ α)(â k -y)f 3 (h, ||â k -y||, ξ)dsy. (88) M 4 (â k , o, β, h, ξ) = ∫ β (a -o|⃗ α)(a -y)f 4 (h, ||â k -y||, ξ)dsy.
• Coplanar triangles As we saw in section (4.2.2), we need to distinguish two geometrical configurations for the supports of the sides α and β.

Two intersecting segments

When the supports of α and β are intersecting, we choose the intersection of the supports, denoted o, as the origin. Then we use the same decompositions as (82) and ( 83).

(89)

ϕp(y) = Cp + µ (x -o|⃗ α) |α| 2 ψq(y) = Cq + µ (y -o| ⃗ β) |β| 2 .
In case ξ -2 ̸ = 0: If we combine (89) with (18) we obtain (90)

Q ϕp,ψq (α, β, 0) = CpCq 2 [ ∑ k=± ks k S(a k , β, 0) - ∑ l=± lt l S(b l , α, 0] + Cpµ 3|β| 2 [ ∑ k=± ks k N (a k , o, β, 0) - ∑ l=± lt l (b l -o| ⃗ β)S(b l , α, 0)] + Cq ν 3|α| 2 [ ∑ k=± ks k (a k -o|⃗ α)S(a k , β, 0) - ∑ l=± lt l N (b k , o, α, 0)] + νµ 4|α| 2 |β| 2 ∑ k=± ks k [(a k -o|⃗ α)N (a k , o, β, 0) - ∑ l=± lt l (b l -o| ⃗ β)N (b k , o, α, 0)].
where

(91) N (a k , o, β, ξ) = ∫ β (y -o| ⃗ β)(a k -y) ||a k -y|| 1+ξ dsy.
5 Numerical application

Calculation of the Radar Cross Section

In radiation theory it's well known that the energy intercepted by an object can be reflected, transmitted or absorbed through the target. We can assume that most of the energy is reflected. The spatial distribution of this energy depends on the size, shape and composition of the target, and on the frequency and nature of the incident wave. This distribution of energy is called scattering, and the target itself is often referred to as a scatterer. The radar cross section (RCS) of the body is a measure of the energy scattered in a particular direction for a given illumination [START_REF] Rahmat-Samii | Hoppe Impedance Boundary Conditions in Electromagnetics[END_REF]. Bistatic scattering is the name given to the situation when the scattering direction is not back toward the source of the radiation. If E and H represent fields scattered by an object illuminated by incident plane wave E inc traveling in the direction of the unit vector k, the bistatic radar cross section in the observation direction r is

σ(r, k) = lim r→∞ 4πr 2 |E| 2 |E inc | 2 .
This cross section is defined as the area through which an incident plane wave carries sufficient power to produce, by omnidirectional radiation, the same scattered power density as that observed in a given far field direction. The monostatic radar cross section is defined as the radar cross section observed in the back scattering direction, σ(-k, k). Radar cross section is the measure of a target's ability to reflect radar signals in the direction of the radar receiver, i.e. it is a measure of the ratio of backscatter power per steradian (unit solid angle) in the direction of the radar (from the target) to the power density that is intercepted by the target. The radar cross section of a target can be viewed as a comparison of the strength of the reflected signal from a target to the reflected signal from a perfectly smooth sphere. The units for RCS are square meters. As RCS can span a wide range of values, a logarithmic decibel scale is also used with a typical reference value σ ref equal to 1m 2 : (92) σ dBm 2 = 10 log 10 ( σ σ ref )

Numerical results

To illustrate our approach, we compute the RCS of a unit sphere coated with a thin dielectric layer. The radius of the inner conductor is denoted by r and the thickness of the coating by d. The material properties are described by a Leontovich Impedance Z = const.

Z = z 0 √ µr ϵr tan ( √ µrϵrk 0 d)
where µr and ϵr denote the magnetic permeability and electical permittivity of the coating. To validate the code, we test it on a sphere with different sphere-meshes and we compare it with Mie series. We take the sphere's radius is 1m and incident waves pulsation ω = 0.9GHz.

In figure [START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF], we plot the exact RCS calculated by Mie Series(continuous black line) and we compare it with our code results in different sphere mesh with inner radius r = 0.95 and coating thikness d = 0.05. In 6, the azure continuous line is the ucp code calculated on the mesh λ/10, red continuous line in our code calculated on the mesh λ/15 and the blue discontinue line is our code on the mesh λ/20.

In the second test, we decrease the coating thikness to d = 0.125 and we increase the inner sphere radius, for the same frequence as the previous test. In figure [START_REF] Huddleston | An impedance boundary condition approach to radiation by uniformly coated antennas[END_REF], we plot the exact RCS calculated by Mie Series(continuous black line) the red continuous line is our code calculated on the mesh λ/10, azure continuous line represent our code calculated on the mesh λ/15 and the blue discontinuous line represent our code on the mesh λ/20 In figure 8, we consider a coating thikness d = 0.125, a frequence f = 014313GHz, the same constant for the magnetic permiability µr = 1.0 and electical permittivity ϵr = 2073.42 . In figure [START_REF] Jingguo | Modified surface impedance boundary conditions fo 3-D eddy currents problems[END_REF], we plot the exact RCS calculated by Mie Series(continuous red line) and the RCS calculated with our code 6 Appendix

Barycentric basis function

We begin by presenting the basic functions used and some of their properties that will simplify calculations in the section 4. Let S be the triangle of vertices a i and side α i . We introduce the basic functions ϕp associated with

(93) ϕ i (x) = 1 + ( ----→ x -a i | -→ e i ); i = 1, 2, 3
where

⃗ e i = µ 1 --→ α i+1 + µ 2 --→ α i+2 with µ 1 = || --→ α i+2 || 2 + ( --→ α i+1 | --→ α i+2 ) || --→ α i+2 || 2 || --→ α i+1 || 2 -( --→ α i+1 | --→ α i+2 ) 2 µ 2 = - || --→ α i+1 || 2 + ( --→ α i+1 | --→ α i+2 ) || --→ α i+2 || 2 || --→ α i+1 || 2 -( --→ α i+1 | --→ α i+2 ) 2
The basis function

ϕ i verify ϕ i (a i ) = 1, ϕ i (a i+1 ) = ϕ i (a i+2 ) = 0 besides (94) ϕ i|α i = 0, ϕ i (x) |α i+1 = 1 - x -a i |α i+1 | ; ϕ i (x) |α i+2 = 1 - x -a i |α i+2 | . 6.2 Calcul of S(a, β, ξ) (95) S(a, β, ξ) = ∫ β a -y ||a -y|| 1+ξ dsy
Let q(a) be the projection of a on the support of β and d = ||q(a) -a||. Then (96) S(a, β, ξ) = (a -q(a))

∫ t + t - 1 (d 2 + s 2 ) 1+ξ 2 ds + ∫ t + t - s (d 2 + s 2 ) 1+ξ 2
ds
where s is the abscissa on β which will be calculated as follows s = y -q(a) and t ± where t ± are abscisses of the extremities of beta with q(a) origin and d = q(a) -a.

6.3 Calcul of S 1 , S 2 and S 3

• S 1 (a, β, h, ξ) S 1 (a, β, h, ξ) = ∫ β f 1 (h, ||a -y||, ξ)dsy (97) = ∫ β h 1-ξ ∫ +∞ h du u 2-ξ (u 2 + ||a -y|| 2 ) 1+ξ 2
ds y where fn(h, s, ξ) were be defined in 48. Let q(a) be the projection of a into the support of β, then

(98) S 1 (a, β, h, ξ) = ∫ t + t - h 1-ξ ∫ +∞ h du u 2-ξ (u 2 + d 2 + s 2 ) 1+ξ 2
dsy where t ± are abscisses of the extremities of beta with q(a) origin and d = q(a) -a.

if ξ = 0 S 1 (a, β, h, 0) = ∫ t + t - h 1 ∫ +∞ h du u 2 (u 2 + d 2 + s 2 ) 1 2 dsy = ∫ t + t - √ h 2 + d 2 + s 2 -h d 2 + s 2 ds (99) = ∫ t + t - 1 √ h 2 + d 2 + s 2 -h ds (100) S 1 (a, β, h, 0) = [ h s - √ h 2 + s 2 s + lg(s + √ h 2 + s 2 )] t + t - for ξ = 2 S 1 (a, β, h, 2) = ∫ t + t - h -1 ∫ +∞ h du (u 2 + d 2 + s 2 ) 3 2 dsy = ∫ t + t - 1 h(d 2 + s 2 ) (1 - h √ h 2 + d 2 + s 2 )ds (101) = [ 1 dh arctan( s d ) - 1 dh arctan( hs d √ h 2 + d 2 + s 2 )] t + t - -----------------

------------------------------------------• S

ds -----------------------------------------

• T 2 (a, β, h, ξ)

(

-----------------------------------------

• T 3 (a, β, h, ξ) = (a -q(a))

(
∫ t + t - h 2-ξ ∫ +∞ h du u

----------------------------------------------------------

• R 2 (a, β, h, ξ)

Figure 1 :Figure 2 :

 12 Figure 1: Definition of Triangle S Figure 2: Orientation of the triangle S

Figure 3 : Projection Figure 4 :

 34 Figure 3: Projection Figure 4: Projections and distances

Figure 5 :

 5 Figure 5: Distances signees

  y) ||a -y|| 1+ξ dsxdsy, with a = a i+2

*

  ||a -y|| 3 ln(||a -y||)dsy. If n + q ̸ = 0, we apply (18) and we obtain (51) P (a, T, ξ) ||a -y|| 1+ξ dsy.

  (61) P ψq (a, T, ξ) = ∫ T ψq(y)(a -y) ||a -y|| 1+ξ dy.

ψq 3 (

 3 â, β q+2 , h, ξ), where (67) D ψq n (a, βp, h, ξ) = δ q+1 (a) ∫ βp ψκ(y)(a -y)fn(h, a -y, ξ)dsy. • Case a belongs to the support of T It suffices to consider a as origin to have homogeneity . Since o = a is the new origin, ψq(y) can be written ψq(y) = ψq(a) + (y|lq). We obtain (68) P ψq (a, T, ξ) = ψq(a) ∫ T (a -y) ||a -y|| 1+ξ dy + ∫ T (y|lq)(a -y) ||a -y|| 1+ξ dy.

Figure 6 :Figure 7 :

 67 Figure 6: Bistattic RCS for a coated unit sphere, when d = 0.05, r = 0.95, µ r = 1.0, ϵ r = 2 and ω = 0.9GHz

Figure 8 :

 8 Figure 8: Bistattis RCS for a coated unit sphere

  -------------------------• S 2 (a, β, h, ξ) S 2 (a, β, h, ξ) = ∫ β f 2 (h, ||a -y||, ξ)dsy

  -------------------------------------------------------• R 2 (a, β, h, ξ) R 2 (a, β, h, ξ) = ∫ β (a -y)f 2 (h, ||a -y||, ξ)dsy = (a -q(a)) ∫ β f 2 (h, ||a -y||, ξ)dsy + ∫ β (q(a) -y)f 2 (h, ||a -y||, ξ)dsy (117)

R 3 ((u 2 + d 2 + s 2 )(u 2 + d 2 + s 2 )R 3

 32222223 a, β, h, ξ) = ∫ β (a -y)f 3 (h, ||a -y||, ξ)dsy = (a -q(a)) ∫ β f 3 (h, ||a -y||, ξ)dsy + ∫ β (q(a) -y)f 3 (h, ||a -y||, ξ)dsy (119) = (a -q(a)) (a, β, h, ξ) = (a -q(a))S 3 (a, β, h, ξ) -⃗ β |β|T 3 (a, β, h, ξ)

  3-ξ (u 2 + ||a -y|| 2 ) (u 2 + d 2 + s 2 ) (u 2 + d 2 + s 2 )

													1+ξ	dsy
													2
	(103)			S 2 (a, β, h, ξ) =	∫ t + t -	h 2-ξ	∫ +∞ h	du u 3-ξ 1+ξ 2	ds
	if ξ = 0
					S 2 (a, β, h, 0) =	∫ t + t -	h 2	∫ +∞ h	du u 3 1 2	ds
	(104)	[	h 2 + d 2 2d 2 arcsinh(	√	s h 2 + d 2	) -	2d 2	sh 2 √ s 2 + d 2	arcsinh(	√	d 2 + s 2 h	)] t + t -
	for ξ = 2
							S 2 (a, β, h, 2) =	∫ t + t -	∫ +∞ h	du u(u 2 + d 2 + s 2 )	2 3	ds
		=	∫ t + t -	-	(d 2 + s 2 )	1 h 2 + d 2 + s 2 √	+	1 (d 2 + s 2 )	3 2	arcsinh	(	√ d 2 + s 2 ) h	ds
	(105)	= [	d 2 ( √	s d 2 + s 2 )	arcsinh(	√	d 2 + s 2 h	) -	1 d 2 log(s +	√	d

2 + h 2 + s 2 )] t + t -

  3 (a, β, h, ξ) S 3 (a, β, h, ξ) = (u 2 + d 2 + s 2 ) + h 2 + s 2 )(d 2 + s 2 -2h 2 ) 3(d 2 + s 2 ) 2ds Calculation of T 1 , T 2 and T 3• T 1 (a, β, h, ξ)

	6.4 (109)		T 1 (a, β, h, ξ) =	∫ t + t -	sf 1 (h,	√	d 2 + s 2 , ξ)ds
	if ξ = 0
		T 1 (a, β, h, 0) =	∫ t + t -	sh	∫ +∞ h	du u 2 (u 2 + ||a -y|| 2 )	2 1	ds
	if ξ = 2
		T 1 (a, β, h, 2) =	∫ t + t -	sh -1	∫ +∞ h	du (u 2 + ||a -y|| 2 )	3 2
								∫
									f 3 (h, ||a -y||, ξ)dsy
									β
	(106)	S 3 (a, β, h, ξ) =	∫ t + t -	h 3-ξ	∫ +∞ h	du u 4-ξ 1+ξ 2	ds
	if ξ = 0
		S 3 (a, β, h, 0) = d 2 = ∫ t + t -2h 3 + ( √ 1 3 [ h 3 s -h 2 s √ d 2 + h 2 s 2 d 2 (d 2 + s 2 )	+	h 3 arctan( s d ) d 3
	(107)	-	h 3 d 3 arctan(	d	hs d 2 + h 2 + s 2 √	) + lg(s +	√	d 2 + h 2 + s 2 )]
	if ξ = 2
	S (108) = [ s(-h + s d 2 (d 2 + s 2 ) √ d 2 + h 2 + s 2 )	-	h d 3 (arctan(	s d	) -arctan(	d √	hs d 2 + h 2 + s 2	))]

3 (a, β, h, 2) = ∫ t + t - d 2 + s 2 + 2h 2 (h -√ d 2 + h 2 + s 2 ) (d 2 + s 2 ) 2 √ d 2 + h 2 + s 2 ds

  110) T 2 (a, β, h, ξ) =

												∫ t + t -	sf 2 (h,	√	d 2 + s 2 , ξ)ds
	if ξ = 0
				T 1 (a, β, h, 0) =	∫ t + t -	sh 2	∫ +∞ h	du u 3 (u 2 + ||a -y|| 2 )	2 1	dsy
		=	∫ t + t -	s	√ 2(d 2 + s 2 ) d 2 + h 2 + s 2	-	sh 2 2(d 2 + s 2 )	3 2	arcsinh(	√	d 2 + s 2 h	)ds
	(111)			= [	√	d 2 + h 2 + s 2 2	+	2 √	h 2 d 2 + s 2	arcsinh(	√	d 2 + s 2 h	)] t + t -
	if ξ = 2
				T 1 (a, β, h, 2) =	∫ t + t -	s	∫ +∞ h	du u(u 2 + ||a -y|| 2 )	2 3	ds
	=	∫ t + t -	-	(d 2 + s 2 )	s d 2 + h 2 + s 2 √	+	s (d 2 + s 2 )	3 2	arcsinh(	√ d 2 + s 2 h	)ds
	(112)								[-	√	1 d 2 + s 2	arcsinh	(	√ d 2 + s 2 h	)] t + t -

  113) T 3 (a, β, h, ξ) = + s 2 ) 2 ( √ d 2 + h 2 + s 2 (d 2 + s 2 -2h 2 ) + 2h 3 )ds + h 2 + s 2 ) √ d 2 + h 2 + s 2 (d 2 + s 2 ) 2 Calcul of R 1 , R 2 and R 3 • R 1 (a, β, h, ξ) R 1 (a, β, h, ξ) = ∫ β (a -y)f 1 (h, ||a -y||, ξ) = ∫ β (a -y)h 1-ξ ∫ +∞ h du u 2-ξ (u 2 +||a-y|| 2 )

	6.5 (116)
	∫ t + t -	sf 3 (h,	√	d 2 + s 2 , ξ)ds
	if ξ = 0
	∫ t + t -3(d 2 (114) T 3 (a, β, h, 0) = [ -h 3 + (d 2 + h 2 + s 2 ) 3 s 2 3(d 2 + s 2 )	] t + t -
	if ξ = 2
	T 3 (a, β, h, 2) = d 2 (115) ∫ t + t -sd 2 + s 3 + 2sh(h -√ [ h -√ d 2 + h 2 + s 2 d 2 + s 2 ] t + t -

  3-ξ (u 2 + d 2 + s 2 ) (u 2 + d 2 + s 2 )

		1+ξ 2	dsy -	⃗ β |β|	∫ t + t -	sh 2-ξ	∫ +∞ h	du u 3-ξ 1+ξ 2	dsy
	wich is equivalent to
	(118)	R

2 (a, β, h, ξ) = (a -q(a))S 2 (a, β, h, ξ) -⃗ β |β| T 2 (a, β, h, ξ)