Simon Archipoff
email: simon.archipoff@labri.fr

David Janin
email: david.janin@labri.fr

Unified Media Programming: An Algebraic Approach

Keywords: temporal animation, 3D programming, model based programming language

In this paper, we aim at defining a simple and sound mathematical framework for describing temporal media programming language semantics. It occurs that semigroup theory offers various concepts that are especially well suited for this purpose.

As a result, a fairly general programming scheme can be defined in order to specify, compose and render both spatial media objects (e.g. 3D drawings) and timed media objects (e.g. Animation or Music). Each of these constructs is specified in Haskell via an adequate type class definition and an associated uniform data type construct.

A simple monoid based semantics model of the turtle command language of Logo is detailed and extended throughout the paper. This allows for providing step by step introductions and usage examples of the algebraic concepts and constructs our proposal is based on.

Introduction

In this paper, we develop a semantics model for the synthesis of temporal media as defined by Hudak [START_REF] Hudak | A Sound and Complete Axiomatization of Polymorphic Temporal Media[END_REF] such as animations, as illustrated in this paper, but also applicable to music, as illustrated in former work [START_REF] Archipoff | Structured Reactive Programming with Polymorphic Temporal Tiles[END_REF] or in Euterpea [START_REF] Hudak | The Haskell School of Music : From signals to Symphonies[END_REF].

Spatio-temporal modeling. Focusing our attention on the synthesis of temporal media, programs can handle the time dimension quite in the same way drawing programs handle space dimension. It is only at the rendering stage, that is, when executing a program, that the time dimension specifically differs from space dimensions.

More precisely, when executing a temporal media program, in other words, when rendering a defined temporal media, time flow is irreversible, i.e. a figure that has been displayed on the screen cannot be "undisplayed". However, at the programming stage, that is, when specifying a temporal media program, nothing prevents a programmer/designer from going back and forth along the time dimension. A temporal media program is out-of-time while its execution is in-time [START_REF] Hudak | From out-of-time design to in-time production of temporal media[END_REF]. This means that, when defining temporal media, both space and time dimensions can be treated similarly.

It follows that any programming construct that allows for drawing a picture, say on the two dimensional space defined by the screen, can also be used for defining the transformation of that picture along the time dimension when defining an animation. This simple observation and a careful study of the underlying semantics of a typical drawing programming language: Logo's turtle command, leads us to the definition of a kind of unified model for temporal media synthesis programming.

Semantics driven DSL definition. Aiming at defining a domain specific language (DSL) for temporal media programming, our approach is governed by semantics.

More precisely, we develop a semantic model for temporal media from which syntax can be derived. Our syntax will essentially reflect the way temporal media can be transformed and combined. It follows that, at the present stage, it remains mostly implicit. What does matter throughout our presentation is semantics.

Then one may ask what is the mathematical framework we use for our semantic model. Our answer is, somehow paradoxically, Haskell programs, that is, pieces of Haskell syntax! There is (essentially) no need for distinguishing program's syntax from program's semantics: every (pure) Haskell function (implicitly) denotes a unique (mathematical) function.

Incidentally, purity also forces all semantical dependencies to be made explicit; there is no possible "magic semantic trick" hidden in a side effect. However, despite such a complete and explicit semantic model definition, a simple application programming interface (API) can still be achieved. All the semantical dependencies that can uniformly (and automatically) be realized can be hidden into higher order programming constructs: mostly monoids, but also monads.

The program constructs we propose throughout this paper are essentially based on combining, by means of a single uniform associative binary operator, various pieces of programs that both act on space (or time) dimensions and produce media values.

Main contribution. We first show, as announced above, that spatial and temporal dimensions semantics can indeed be treated (mostly) in the same way. For this purpose, we develop an original algebraic semantics for temporal media programming languages.

Our semantic model is based on semigroup theoretical concepts such as semigroup actions over a set, semi-direct products, inverse semigroups and their extension to semigroups with local units (called resettable semigroup in the present paper).

Though well-known in semigroup theory, these concepts are yet quite unknown in programming language theory1 despite the fact they perfectly fit, as shown below, our needs.

Structure of the paper. We start our presentation by examining the semantics of Logo's turtle programs in which commands for movement and drawing produce line graphics either on screen or with a small robot called a turtle. The language [..] enable [..] "bodysyntonic reasoning", where students can understand, predict and reason about the turtle's motion by imagining what they would do if they were the turtle. 2 .

After reviewing a basic (equivalent) syntax for turtle programs, we show that their semantics can be described by means of a monoid. This semantic monoid is shown to be an inverse monoid [START_REF] Lawson | Inverse Semigroups : The theory of partial symmetries[END_REF] with elements reflecting both the moves and the drawings performed by the turtle.

The algebra theoretic concepts that support our semantic models are detailed throughout. Then we show, in the remainder of the text, how these concepts are applicable for lifting Logo's turtle program to time and additional space dimensions, i.e. for animation and 3D drawings. Little incursions into monads also provide a rather comfortable syntax with do notation and timely rendering of media programs within the IO monad.

Notations.

Most examples, formulas and programs are given using Haskell syntax. In particular, as in most functional programming language, we write f x for the application of the function f to the argument x. Also, our paper makes intensive use of the notion of semigroups: sets equipped with a associative binary operator. The semigroup binary operator is denoted by . We also use the notion of monoids: semigroups with a neutral element, denoted by mempty, with the binary operator possibly denoted by mappend when used as a function, and the concatenation of monoid elements in a list denoted by mconcat.

The Turtle Example

Our approach can be illustrated in a fairly simple way by analyzing the semantics of basic logo programs that allow for defining pictures by commanding turtle movements. These programs are, from now on, called turtle programs.

It is an easy observation that turtle programs act both on the space of possible pen positions and the figures that are drawn. Together, positions and figures form the states that are eventually transformed by programs.

Actions on positions are called moves, actions on figures are called drawings. Together they form the turtle program semantic model we aim at defining.

It occurs that this simple drawing language offers a deep insight into the mathematically well structured semantics we propose for temporal media programming language. Example. With angles defined in degree and provided we start with the pen down, the series of actions Walk 1, Turn 90, Walk 1, Turn 90, Walk 1, Turn 90, Walk 1, Turn 90, draws a triangle with unitary length borders as depicted in Figure 1 where a final TogglePen action has also been performed. In this figure, each move has been grouped with its following turn.

Remark. For the sake of simplicity, we use the same type d as both the type distance in a walk or the type angle in a turn. Later in the text, the same type d is even used as the type of time duration. Of course, distinguishing all these types would certainly lead to a more constrained hence a safer type system.

Moves, Drawings and Semantics Types

Defined by lists of basic actions, turtle's programs form a monoid where sequential composition of program is modeled by list concatenation. Aiming at defining a compositional semantic model for these programs, we thus look for a semantic monoid. The turtle performs, both at the same time, moves when changing positions, and drawings when creating new figures. They both constitute the elements of our semantics.

Formally, a turtle position is defined as a triple composed of a boolean (true when the pen is in drawing position, false otherwise), a 2D vector (the pen coordinates in the plane), and an angle (the turtle orientation). Moves are defined as position transformations. Clearly, the type Semantics d that is given above is not a functional type as one could expect. However, turning a program semantics into a program that draw figures, or into a function over states that compute a new state can be done as follows. with figures forming an idempotent and commutative monoid with union as product and empty set as neutral element.

Producing a figure from a given position, as done by function runS, just amounts to applying the drawing function specified by the second component of the semantic model.

Applying a semantics to a program state, as done by function applyS, is a little more tricky. The output position is computed by applying the move part of the semantics m to the input position p. The output figure is computed by merging the input figures f with the new figure computed by applying the drawing part of the semantics d to the input position p.

Remark. We aim at making the control structures, the functional dependencies and the algebraic properties of the underlying semantics explicit: a much less obvious task. The generalization of these constructions, presented below, may look uselessly (or even pedantically) theoretical. However, it has a merit: as we shall see, it generalizes to a fairly wide range of temporal media programs.

Moves, Drawings and Semantic Monoids

As already mentioned, we are looking for a semantic monoid. It occurs that the types above lead to the definition of three monoids instead of one. The third one, the program semantic monoid, is obtained by combining uniformly the first two, the moves and the drawings monoids.

First, we observe that both moves and drawings form monoids With the help of the above monoids, program semantics, that is, pairs of moves and drawings, should be sequentially combined as follows:

combine ::

Semantics d → Semantics d → Semantics d combine (m 1 , d 1 ) (m 2 , d 2 ) = (m 1 m 2 , d 1 (d 2 • m 1 ))
On the move side, the sequential application of two turtle programs combines the moves defined by these programs. On the drawing side, the above expression says that, from a given input position p, a first figure This follows from the similarities one can observe between our definition of the applyS and the combine functions.

d

Basic Moves and Drawings Implementation

So far, we have defined types and combinators for composing semantics of turtle programs. These types are based on the decomposition of our semantic values into moves and drawings. It thus remains to explicit the moves and drawings associated to basic actions in order to recover a complete turtle semantic model. Basic moves semantics is defined as follows.

moveA ::

BasicAction d → Move d moveA TogglePen (b, v, a) = (¬ b, v, a) moveA (Turn d) (b, v, a) = (b, v, a + d) moveA (Walk d) (b, v, a) = (b, v + w, a) where w = (V 2 (d * cos ra) (d * sin ra)) ra = 2 * pi * a/360
The move semantics associated to every program, that is, list of basic actions, is then defined by:

move :: Program d → Move d move prog = mconcat (map moveA prog)
As expected, the semantics of a sequential composition of programs is mapped to the (flipped) composition of their semantics. Basic drawing semantics is defined similarly. Observe that drawing semantics is slightly more complex to define than move semantics because it depends on positions which may be altered by moves. Nevertheless, this corresponds to a fairly well known construction in algebra called semi-direct product that is detailed below. This semantics is compositional in the sense that, for all program p1 and p2, we have sem mempty == mempty sem (p1 p2) == sem p1 sem p2 Indeed, the function sem is the unique monoid morphism from the (free) monoid of turtle programs into turtle semantics that is induced by (moveA, drawA) on elementary actions.

Remarks on Basic Turtle Moves

At first sight, one may think that a turtle move essentially traverses the 2D space defined by the (x,y)-coordinates of every position. However, the capacity to turn says that the turtle can also act on the space itself: namely its underlying 2D basis.

Indeed, by combining basic turtle walks and turns, one can define any isometric and orientation preserving transformation of the underlying 2D space, that is, the transformations generated by arbitrary translations and rotations.

Extended with additional basic moves, turtle semantics captures arbitrary isometries (adding reflexion across, say, the x-axis), arbitrary affine bijections (adding reflexion and non zero contraction/walking speed along the x-axis) or even arbitrary affine transform (also adding projection of the entire 2D space into a line by means of zero scaling along the x-axis).

This observation plays a key role when extending turtle moves to additional space or time dimensions.

Monoid Semantics Revisited

We have seen above that monoids (or even just semigroups) play a central role in the definition of the turtle semantics. Moreover, most semantics definitions are generic. Functions move, draw and sem are uniformly defined in terms of basic actions (BasicAction d) and their move and drawing semantics (functions moveA and drawA).

This suggests that the turtle semantics relies on some general structures that are worth being made explicit. In this section, we study more in depth the generic constructions that have been (yet implicitly) applied and the rather rich mathematical properties these constructions satisfy.

Semi-direct Product

We first observe that the turtle semantic monoid is a fairly standard construction known in semigroup theory as the semi-direct product of two monoids.

This notion is defined via what algebraists call an action of a monoid over a set. class (Monoid m) ⇒ MonoidAction m a where act :: m → a → a with the requirement that properties should hold for all a, m 1 and m 2 of adequate types.

act mempty a == a (1)

act m 1 (act m 2 a) == act (m 1 m 2 ) a (2)
The associated notion of semi-direct product of monoids is defined by:

instance (MonoidAction m a, Monoid a) ⇒ Monoid (m, a) where mempty = (mempty, mempty) mappend (m 1 , a 1 ) (m 2 , a 2 ) = (m 1 m 2 , a 1 act m 1 a 2 )
Back to the turtle semantics, we simply put:

instance MonoidAction (Move d) (Drawing d) where act m d = d • m
The generic monoid instance given above applies to Semantics d, turning it into a monoid and we have combine == ( ).

Inverse Monoid

So far, we haven't used the fact that (turtle) moves are reversible. It turns out that those properties open the way to a transformation of turtle program, called inverseP that have, up to semantic equivalence, fairly interesting properties.

which states, following inverse semigroup theory [START_REF] Lawson | Inverse Semigroups : The theory of partial symmetries[END_REF], that inverse a is the semigroup inverse of a, and

a inverse a b inverse b == b inverse b a inverse a (5)
which states that elements of the form a inverse a commute. Then, inverse semigroup theory ensures that idempotents 5 , necessarily of the form a inverse a, commute and that inverse are uniquely defined, i.e. inverse a is the unique element b such that both properties a b a == a and b a b == b hold.

As another derivative of inverse semigroup theory, the relation norder is a partial order relation known as the natural order [START_REF] Nambooripad | The natural partial order on a regular semigroup[END_REF]. The word natural refers in semigroup theory to the fact it is uniformly defined over all inverse semigroups.

Then we can show that turtle programs form, up to semantic equivalence, an inverse monoid: for some turtle program p, one can check that we have

m ′ • m == id == m • m ′ and d ′ == d • m ′
This confirms that, as already mentioned, the moves induced by Turtle programs are in fact bijections, or, stated differently, the moves induced by turtle programs form a group. Then, it is a known fact in semigroup theory (see [START_REF] Lawson | Inverse Semigroups : The theory of partial symmetries[END_REF], Section 7.1) that a sufficient condition for the semi-direct product (m, a) of a monoid m by a monoid a to produce an inverse semigroup is that m is in fact a group and a is a lattice with product as meet. This is the case for (the subtype of bijective) moves and drawings induced by turtle programs. This fact concludes the proof that the program monoid is, up to semantic equivalence, an inverse monoid.

Resettable Monoid

The inverse semigroup detailed above also induces an additional program transformation, called the reset operation, that is also of interest. It is definable even in the case of non reversible moves, and, as we shall see, it behaves like a fork operator hence allowing concurrent evaluation.

More precisely, we define the function resetP over turtle programs by In other words, the program resetP p can also be read as: from any given initial position, send a copy of the current turtle to perform the expected drawing and disappear upon termination while the initial turtle stays in that position waiting for further instructions.

This says at least two things. First, the reset function acts over a program p as a sort of fork operation that allows parallel evaluation. Second, implementing the reset directly at the semantic level is more efficient and can be defined even in the case of non reversible moves.

In semigroup theory, resettable semigroups are known as left semiadequate semigroups (see [START_REF] Fountain | A Munn type representation for a class of E-semiadequate semigroups[END_REF] and the related works mentioned at the end). Property [START_REF] Hines | The algebra of Self-Similarity and its Applications[END_REF] says that reset a is a left local unit of a. Properties ( 7) and ( 8) say that these left local units are idempotent and commute henceforth form a semi-lattice with product as meet. Property [START_REF] Hudak | Tiled Polymorphic Temporal Media[END_REF] says that reset a is the least left local unit of a in this semi-lattice. From these properties, one can also prove that the property reset • reset == reset is satisfied. In other words, the reset function is a projection. Indeed, let a be some monoid element, let e = reset a and let e ′ = reset e. It remains to prove that e == e ′ . By property [START_REF] Hines | The algebra of Self-Similarity and its Applications[END_REF] we have e ′ e == e. This implies that e ⩽ e ′ in the semi-lattice of local units. But since e is idempotent, we have e e = e hence, by property (9) e ′ ⩽ e. Together, this implies that e == e ′ .

Of course, every monoid can be turned into an instance of ResettableSemigroup by taking reset = const mempty. However, this instance is of very little interest. Incidentally, we do not require m to be a monoid as shown by the premises Semigroup m of the class definition. In the absence of a neutral element, existence of a reset is by no mean a trivial property.

Inverse semigroups are particular cases of resettable semigroups thanks to the following (non trivial) instance:

instance InverseSemigroup m ⇒ ResettableSemigroup m
where reset m = m inverse m and, as already observed above, over semi-direct product, one may prefer the following and more efficient instance. Of course, there is no syntactic way to reset efficiently a turtle program. But, as already noticed, thanks to our uniform semantics, this is really just a matter of syntax. We can add a basic action Reset (Program d) that takes an entire program as argument and whose semantics is defined by moveA (Reset ) = mempty and drawA (Reset p) = draw p. Also, we do not use here the fact that the set of moves induced by turtle programs are reversible. This means that we may even add irreversible basic move commands such as WalkToPosition v, a command that is irreversible, without losing the possibility to perform a reset. So far, our turtle programming language is essentially based on lists over some parametrized alphabets (basic actions) with mathematically rich but not (yet) clearly usable properties. For instance, conditionals, loops and variable bindings are missing for drawing complex figures. Also, effectively drawing figures necessitates IO actions. These leads us to embed the turtle semantics into monads.

Drawing Monad Semantics

Embedded in Haskell, our proposal inherits all the programming framework proposed by Haskell. However, combining Haskell functions and turtle programs may not be that convenient. The monoid product syntax ( ) does not sound much like programming at all. Also, the above discussion about the reset function and possible extension to non reversible moves suggests that our programs should be defined directly at the turtle semantics level as pairs of moves and drawings with potential unnecessarily heavy turtle state transfer.

The monad approach (and the do notation) offers a clean and efficient approach for this purpose. It occurs that Haskell GHC Prelude even offers an easy way to embed a monoid into a monad.

The following instance, recalled below for the sake of reasonable self completeness, can relevantly be used: It follows that turtle semantics can instead be seen as an IO action that takes a position as input and actually produces a drawing while also sending back the new current position. The type of moves over a space type s can thus be described by: Various equalities could be written down in order to validate all these definitions. For instance, one can check that we have liftM (load p) == liftS (sem p) for all turtle program p. The important point is that the monoid, the monad or the IO monad approaches essentially define equivalent semantics.

Remark. Command by command rendering may be useful in programming context such as live coding. However, for proper 2D or even 3D rendering of complex pictures, it may be worth accumulating in the underlying state all drawing commands so that, ultimately, after evaluating an entire program, features that are context dependent, such as face covering in 3D drawings, can be computed and rendered correctly.

Time Extension

So far, our turtle is capable of performing figures in a two dimensional space, henceforth performing 2D drawings. We aim now at extending our turtle with time in order to define 2D animation. All the algebraic tools developed so far can be used for this purpose.

Animation

At the semantics level, extending our turtle to the time dimension is rather easy. We just define the type:

type Animation d = d → Semantics d
The meaning of an element a :: Animation d is the following: the move performed and the image drawn at some time stamp t is given defined by a t. In other words, we can define the function Just given as a reminder as this is a particular instance of the generic monoid instance of the monoids of functions into a monoid.

Example. Every figure f can simply be embedded into an animation by defining liftA f = const (id, const f ), that is, the turtle programming displaying the figure f at every instant and from any position. Then we have

liftA (f 1 f 2 ) == liftA f 1 liftA f 2
Observe that, with arbitrary animations and turtle moves, the monoid composition of two animation defined above is non trivial. It is not a mere superposition of two animations.

Remark. The figure produced by (a 1 a 2 ) t at every instant t, is not in general the superposition of the two figures produced by a 1 t and a 2 t. The space position from which the drawing induced by a 2 t is performed is indeed defined by the position where the turtle stops when drawing a 1 t. It turns out that composing two animations from the same position necessitates a reset that can be defined by instance ResettableSemigroup (Animation d) where reset a = \t → reset (a t) that inherits from the ResetableSemigroup instance of the type Semantics d. The timed superposition of two animations a 1 and a 2 at the same instant and position can then be defined by reset a 1 a 2 .

Temporal Moves

In terms of semantics, the animation type defined above already allows for defining any 2D animation. However, in terms of temporal programming, one may be in need of uniformly defined functions that allows to start, stop or run an animation from, to or between two specified timestamps. One may also wish to modify the speed of a given animation. Even more importantly, one may feel the need to combine animations over time. This observation leads us to the definition of the temporal moves type and its associated monoid instance: Remark. The usage of the type Maybe d as codomain may look a bit awkward at first sight. However, this turns out to be the most elegant way to "block" passing time while preserving the property that an instant has zero duration. Indeed, when blocking time, no instant value will be given.

Examples Many other examples of timed moves can be defined such as, for instance, the time stretch \t → Just (d * t) that amounts to compressing (when d < 1) or stretching (when d > 1) the time dimension by some (generally strictly positive) constant.

However, in the absence of "visible" effect of these timed moves over animations makes these examples a bit abstract. The real power of timed moves appears when associated with animations as shown right below.

Temporal Semantics Extension

The temporal extension of the turtle semantics type is obtained by pairing timed moves and animations.

(d 1 , a 1 ) (d 2 , a 2 ) = (d 1 d 2 , a 1 act d 1 a 2 )
all these instances given just as a reminder. Reset is extended to timed semantics (and the underlying animations) by the following instance:

instance ResetableSemigroup (TExtension d) where reset ( , a) = (Just, reset a)

which amounts to superpose two (embedded) animations at the same time and position. One can check that we have

reset (d 1 , a 1 ) (d 2 , a 2 ) == (d 2 , reset a 1 a 2 )
for all time moves d 1 and d 2 and all animations a 1 and a 2 . Running a timed semantics from a given instant t and a given space position p can be defined by the function Examples. As a first application example, delaying an animation a by d units of time can be defined by delay d liftTS a with delay d = (\t → Just (d + t), mempty).

As another application example, stretching an animation by a non-zero factor f can be defined by stretch f liftTS a with stretch f = (\t → Just (f * t), mempty).

Gathering these timed semantics elements that only act on the time dimension eventually leads us to the definition of the following class type.

Timed Monoid Class Type

The class type TimedMonoid is defined by:

class (Monoid m) ⇒ TimedMonoid m d | m → d where liftTM :: (d → Just d) → m delay :: d → m delay d = liftTM (Just • (d+)) stretch :: d → m stretch d = liftTM (Just • (d * )) start :: m → m stop :: m → m play :: d → d → m → m play t 1 t 2 m = start (delay (t 1 -t 2 ) stop (delay (t 2 -t 1 ) m))
In the class type defined above, the method liftTM method shall lift any time transformation into the timed monoid m in such a way that the derived methods delay and stretch behave as already illustrated in the examples above.

Methods start and stop shall start or stop a timed process at a given (fixed) time origin, say, the zero instant. The derived method play then plays a timed element from the instant t 1 to the instant t 2 . Timed moves form an instance of timed monoid:

instance TimedMonoid (TMove d) d where liftTM = id start d t = if (t < 0) then Nothing else d t stop d t = if (t > 0) then Nothing else d t
The timed monoid instance of timed semantic extension follows from the above as it uniformly derives from the timed monoid instance of timed moves 6 .

instance Example. Running the timed semantics play 0 d (liftTS f ) will display the figure f from the instant 0 to the instant d, inclusive, performing no drawings when outside this interval.

Remark. It can be shown that the monoid of timed semantics generated by delays, non-zero stretches, starts and stops, and arbitrary lifts of animation is not an inverse monoid.

It could be turned into one by propagating time cuts induced by start and stop. However, such a propagation of time cuts does not seem to be desirable. On the contrary, the product of timed semantics computes instead some kind timed envelope of the combined animation much like in the way spatial envelopes are computed in Diagrams [START_REF] Yorgey | Monoids: Theme and Variations (Functional Pearl)[END_REF].

Example. There is the notion if partial identities over d that can be defined as function p :: d → Just d such that if p t = Just t ′ for some t :: d then t = t ′ , i.e. seen as a partial function, a partial identity is essentially defined by its domain. Then, as an example of implicit computation of temporal envelopes, one can check that for any partial identities d 1 , d 2 , for all animation a 1 , a 2 , we have

(d 1 , a 1 ) (d 2 , a 2 ) == (d, a 1 a 2 )
where d is the partial identity (uniquely) defined by the union of the domain of d 1 and d 2 .

Space Extension

So far, our turtle is spatio-temporal. It performs timed 2D figures made of curves. Now we aim at extending its space to perform 3D figures made of surfaces. Again, most of the algebraic tools developed above are reusable for this purpose.

Extruding 3D Turtle

There are many possibilities for defining a 3D drawing language. Following the idea that drawings should be achieved by means of combining turtle pen position and moves one possibility is to define 3D surface by mean of extrusion: a well known 3D technique that aims to create 3D surfaces by moving in the space some varying 2D curve.

More precisely, in the 2D case, we observe that figures are drawn by "extruding" turtle positions that defines points (when pen is down) into segments. Extending turtle positions to 3D space with basic 2D curves as possible turtle pen shapes eventually leads us to the desired programming language.

As preliminary code, we define 3D affine transform by

type Affine3D d = (M33 d, V3 d) instance Monoid (Affine3D d) where mempty = (V3 (V3 1 0 0) (V3 0 1 0) (V3 0 0 1), 0) mappend (m 1 , v 1 ) (m 2 , v 2 ) = (m 1 ! * ! m 2 , v 1 + m 1 ! * v 2 )
Then, limiting our example to circle shapes with a given radius, turtle pen shape, basic 3D actions and 3D turtle programs are defined by

data Penshape d = Empty | Circle d data BasicAction3D d = Pen (Penshape d) | Transform (Affine3D d) type Program3D d = [BasicAction3D d ]
3D turtle moves are defined by arbitrary affine transformation of the underlying 3D space (command Tranform). Extrusion is performed between two successive non empty Pen command, separated by one or more Transform command that are implicitly combined. This allows for defining complex movements by composition of simpler ones, realizing the extrusion only when the resulting complex movement is defined. Then, as it shall become clear with the extrusion semantics detailed below, the expression cylinder r l indeed specifies a cylinder with radius r and length l along the z-axis. where a position describes the last pen shape, its associated affine transform and the current affine transform, and a scene describes a set of triangles. Again, scenes inherit from the monoid structure of sets with union as product and empty set as neutral element. Technically, our proposal is a lot more evolved than what is described above. Octopus programs, implemented in Haskell and executed on the central processing unit (CPU), define 3D scene specifications that consists in trees of (2D shaped) pen positions. Extrusion, tessellation, normal computation and lighting are executed in parallel on the graphics processing unit (GPU) via a modern OpenGL graphic pipeline with vertex, geometry and fragment shaders.

As example, the image in Figure 3 is extracted from a 3D animated scene example, described by means of 5000 timed elements, extruded into more than 500,000 triangles per image, displayed 7 hence inducing only a segment that is omitted 8 see http://poset.labri.fr/octopus/ at 50 frames per second. Its associated triangulation is given in Figure 4 and the normal colors9 are given in Figure 5.

It must be noted that the efficiency of the rendering is also due to Haskell itself that automatically exploit the multicore structure of the CPU, easily reaching more than 180% CPU usage. Although we aim in the long term at defining a 3D animation programming language dedicated to a broad class of users, be they trained with basic geometry notions or not, this goal is yet not achieve. The Octopus language necessitates further experiments and extensions before being considered as a mature front end.

A significant difference with the 2D turtle case is that, in our 3D version, a pen "action" creates now a drawing by extruding a surface from the pen shape and position defined by the previous pen "action". The moves between these two actions are thus implicitly combined to be treated as a single one. More precisely, extrusion allows for combining several simple moves into a complex one which facilitates the task of drawing regular shapes.

For instance, a regular polygon with n vertices and (outer) radius n is defined in Logo by n repetitions of walks of 2 × r × sin(π /n) distance units followed by turns of 2×π /n. Its definition necessitates some basic knowledge in trigonometry. With (2D) extrusion, it suffices instead to define n + 1 successive pen positions (with point shape) separated by the combined movements Turn 90, Walk r, Turn (360 / n), Walk (-r) and Turn (-90), two successive pen positions being extruded into a segment. The resulting drawing is the same without any explicit trigonometric function.

Related Works

Between groups and monoids, inverse semigroups have already been studied and developed quite a while in algebra [START_REF] Lawson | Inverse Semigroups : The theory of partial symmetries[END_REF][START_REF] Meakin | Groups and semigroups: connections and contrasts[END_REF]. In computer science, they already appears in connection with Girard's Linear Logic [START_REF] Hines | The algebra of Self-Similarity and its Applications[END_REF] and related studies of reversible computations [START_REF] Abramsky | A structural approach to reversible computation[END_REF][START_REF] Danos | Reversible, Irreversible and Optimal lambda-Machines[END_REF]. More recently, inverse semigroup have been applied to music modeling [START_REF] Archipoff | Structured Reactive Programming with Polymorphic Temporal Tiles[END_REF][START_REF] Hudak | Tiled Polymorphic Temporal Media[END_REF][START_REF] Janin | A robust algebraic framework for high-level music programming[END_REF].

Resettable monoids, first studied in [START_REF] Lawson | Semigroups and ordered categories. I. The reduced case[END_REF], are known in semigroup theory as semi-adequate monoids [START_REF] Fountain | A Munn type representation for a class of E-semiadequate semigroups[END_REF][START_REF] Lawson | Semigroups and ordered categories. I. The reduced case[END_REF]. Although not much popular even in theoretical computer science, it is known that their free algebras capture deterministic unranked trees [START_REF] Kambites | Free adequate semigroups[END_REF]. They have recently been studied in formal language theory by the second author for developing a new algebraic characterization of regular languages of finite trees [START_REF] Janin | On Labeled Birooted Trees Languages: Algebras, Automata and Logic[END_REF]. It is conjectured that dropping the commutation hypothesis (property (8)) will eventually lead to languages of ranked trees.

To the best of our knowledge, resettable monoids have never been considered in programming language semantics even though interleaving semantics for concurrency relies for long on the notion of commuting processes : a property explicitly axiomatized in resettable monoids.

In the context of functional programming languages, drawing and animation languages have already led to many proposals. As already suggested above, the novelty of our approach does not lay our yet fairly implicit EDSL language proposal. Our front-end is still a toy language example essentially designed for studying the applicability of inverse semigroup theory to the definition of a generic and efficient back-end for media programming languages.

Nevertheless, our proposal can still be related with existing approaches among which Euterpea [START_REF] Hudak | The Haskell School of Music : From signals to Symphonies[END_REF] in computer music and Diagrams [START_REF] Yorgey | Monoids: Theme and Variations (Functional Pearl)[END_REF] in computer graphics.

Indeed, both proposals provide semigroup based uniform and robust program structures that are applicable to many distinct purposes. Paul Hudak's notion of polymorphic temporal media, based on a 2D structuration of time with its (horizontal) sequential and its (vertical) parallel product is an example. The automatic computation of object envelopes of figures in Diagrams is another example. In Diagrams, semigroup actions are even made explicit (see [START_REF] Yorgey | Monoids: Theme and Variations (Functional Pearl)[END_REF] for many other examples).

Our approach clearly follows these paths already open by previous authors.

Conclusion

We thus propose a formal approach for defining (and implementing) a temporal media program semantics model based on semigroup theoretical concepts.

The main novelty of our approach is to make explicit the usage of semi-direct products. This allows, in particular, for combining, within semigroup structures, both (spacial or temporal) moves and (musical or graphical) renderings.

To the best of our knowledge, the notion of semi-direct product, though deriving from the notion of monoid action that is already used in Diagrams [START_REF] Yorgey | Monoids: Theme and Variations (Functional Pearl)[END_REF], have not yet been considered in temporal media programming.

Octopus, an extension of Logo's turtle for programming 3D animation is based on this approach. Efficient rendering is achieved by combining simple (but possibly large) tree shaped temporal 3D scene specifications defined in Haskell (on the CPU) with efficient rendering on a modern OpenGL graphic pipeline (on the GPU).

Though applicable to the definition of non trivial 3D animation, the resulting programming language is stable enough to be viewed as a programming language proposal. In particular, we believe that the uniform semantics modeling of both space and time dimensions that we have proposed throughout these pages is neither fully exploited not fully understood.

Our belief is that inverse (and resettable) semigroups could be much more developed towards programming semantics. For instance, the approach presented here could be generalized for combining in a fairly abstract way reversible computations, e.g. moves, with non reversible ones, e.g. drawings.

Aside the Octopus experiment, one possible development of our proposal may consist in rephrasing the FRP definitions [START_REF] Elliott | Functional Reactive Animation[END_REF] in terms of monoid constructions. Doing so, the inderlyng semigroup syntax will certainly restrict the way behaviors and events can be defined and combined, and this may eventually lead to a more efficient operational semantics. As a matter of fact, recent works on FRP (see [START_REF] Van Der Ploeg | Practical principled FRP: forget the past, change the future, FRPNow![END_REF]) tend to restrict FRP programs to somewhat automata (henceforth semigroup) like specifications for better efficiency in a reactive context.

2. 1

 1 Program Syntax Turtle programs are defined by means of series of basic actions. Up to a minor technical change from the original definition 3 , basic actions are defined by data BasicAction d = TogglePen | Turn d | Walk d which respectively allows for flipping up or down the pen's head (TogglePen), pivoting by the specified angles (Turn) or moving forward by the specified distance (Walk). Turtle programs are then defined as lists of basic actions. type Program d = [BasicAction d ]

type

  Position d = (Bool, V 2 d, d) type Move d = Position d → Position d A figure is defined by a set of drawn segments, that is, a set of pairs of points. Since the turtle never reads a picture, but, instead, produces a picture from a given position, drawings also take positions as input. type Figure d = Set (V 2 d, V 2 d) type Drawing d = Position d → Figure d Program states are defined as pairs of positions and figures, and, program semantics as pairs of moves and drawings. type State d = (Position d, Figure d) type Semantics d = (Move d, Drawing d)

runS : :

 : Semantics d → Position p → Figure d runS ( , d) p = d p applyS :: Semantics d → State d → State d applyS (m, d) (p, f ) = (m p, f d p)

  4 . instance Monoid (Move d) where mempty = id mappend = flip (•) instance Monoid (Drawing d) where mempty = const mempty mappend d 1 d 2 = \p → d 1 p d 2 p

drawA : :

 : BasicAction d → Drawing d drawA (Walk d) (False, , ) = mempty drawA (Walk d) p@(True, v, ) = insert (v, nv) mempty where ( , nv, ) = moveA (Walk d) p drawA = mempty The drawing semantics of program is defined as the second projection of the complete semantics of programs. draw :: Program d → Drawing d draw prog = snd (sem prog) semA :: BasicAction d → Semantics d semA prog = (moveA prog, drawA prog) sem :: Program d → Semantics d sem prog = mconcat (map semA prog)

  inverseP :: Program d → Program d inverseP p = map (invA) (reverse p) where invA TogglePen = TogglePen invA (Turn d) = Turn (-d) invA (Walk d) = Walk (-d) Then, one can check that, for every program p, programs p and p inverseP p p have the same semantics, as well as programs inverseP p and inverseP p p inverseP p. In semigroup theory, inverseP p is called the semigroup inverse of p. This leads us to the definition of the class of inverse semigroups. class Semigroup m ⇒ InverseSemigroup m where inverse :: m → m norder :: (Eq m) ⇒ m → m → Bool norder m n = m == (m inverse m n) with the requirement that the following properties should hold for arbitrary elements a and b of the semigroup m. a inverse a a == a (3) inverse a a inverse a == inverse a

instance

  InverseSemigroup (Program d) where inverse = inverseP with properties (3)-(5) valid up-to semantic equivalence. From the semantics side, given (m, d) = sem p and (m ′ , d ′ ) = sem (inverse p)

  resetP :: Program d → Program d resetP p = p inverse p When executing the program p inverse p, the turtle performs the actions specified by p, producing some figure, and then it executes the inverse of these actions in the reverse order, reproducing the same exact figure while moving back to the starting position. In other words, the program resetP p describes a double drawing of the figure defined by p; a forward drawing with p and a backward drawing with inverse p. At the semantic level, such a redundancy can be avoided as follows. Given the semantics (m, d) = sem p of the program p, that says that the turtle should move following m while drawing following d, one can check that we have (mempty, d) == sem (resetP p)

  Generalizing these observations leads us to the definition of resettable semigroups. class (Semigroup m) ⇒ ResettableSemigroup m where reset :: m → m with the requirement that the properties should hold for arbitrary elements a and b of the semigroup m. reset a a == a (6) reset a reset a == reset a (7) reset a reset b == reset b reset a (8) reset b a == a ⇒ reset b reset a == reset a

  instance (SemigroupAction m a, Semigroup a) ⇒ ResettableSemigroup (m, a) where reset ( , a) = (mempty, a) Applied to turtle semantics, this can be made explicit by: instance ResettableSemigroup (Semantics d) where reset ( , d) = (mempty, d)

instance

  Functor ((, ) a) where fmap f (x, y) = (x, f y) instance Monoid a ⇒ Applicative ((, ) a) where pure x = (mempty, x) (u, f ) < * > (v, x) = (u 'mappend' v, f x) instance Monoid a ⇒ Monad ((, ) a) where (u, a) > > = k = case k a of (v, b) → (u 'mappend' v, b) with > > = the bind operator of monads. Our turtle semantics is then embedded into a monad by defining: type MSemantics d a = (Semantics d, a) This allows for (re)defining basic moves as monadic actions: nop = load mempty toggle = load [TogglePen] walk d = load [Walk d ] turn d = load [Turn d ] with the loading and running function defined by: load :: Program d → MSemantics d () load p = (sem p, ()) runM :: MSemantics d a → Figure a runM (( , d), ) = d (True, 0, 0) Then, turtle programs inherit from monad do notation. Example. The following program draws a triangle. prog = do { walk 1; turn 120; walk 1; turn 120; walk 1 } triangle = runM prog 4.2 Drawing IO-monad Semantics Viewing figures on the screen necessitates an additional function render :: Figure d → IO () which code, depending on the chosen graphical library, is not detailed here.

type

  IOSemantics s a = s → IO (s, a) with the Monad instance instance Monad (IOSemantics s) where return a s = return (s, a) (> > =) m f s = do (ns, a) ← m s (f a) ns Then, thanks to the render function, one can lift turtle monoid semantics into IO semantics by: liftS :: Semantics d → IOSemantics (Position d) () liftS (m, d) p = do render (d p) return (m p, ()) or turtle monad semantics into IO semantics by: liftM :: MSemantics d a → IOSemantics (Position d) a liftM ((m, d), a) p = do render (d p) return (m p, a)

runA : :

 : Animation d → d → Position d → Figure d runA a t = runS (a t) p that simply draws the figure specified by the animation at every instant.Since Semantics d is known to be a monoid, animations can also be turned into a monoid by letting instance Monoid (Animation d) where mempty = const mempty mappend a 1 a 2 = \t → (a 1 t) (a 2 t)

type

  TMoves d = d → Maybe d instance Monoid (TMoves d) where mempty = Just mappend f g = \t → case (f t) of Just t 1 → g t 1 Nothing → Nothing that generalizes the flipped function composition.

type

  TExtension d = (TMoves d, Animation d) The monoid of timed moves acts over animations by instance MonoidAction (TMoves d) (Animation d) where act d a = \t → case (d t) of Just t 1 → a t 1 Nothing → mempty where, in the case time is stopped by a temporal move, the animation no longer produces any figure but the empty one. Then, following the definition of semi-direct product construction, we have: instance Monoid (TExtension d) where mempty = (mempty, mempty) mappend

  runTS :: TExtension d → d → Position d → Figure d runTS ( , a) t p = runA a t p which explicits the fact that the temporal move d in the timed semantics (d, a) only affects the animations that will be combined to the right of the animation (d, a). This generalizes to the time dimension the fact that spatial move performed by a given 2D turtle program p affects the drawings performed by programs combined to the right of p as illustrated by the following examples. Every animation a can simply be embedded into a timed semantics by defining liftTS :: Animation d → TExtension d liftTS a = (Just, a) Then we have liftTS (a 1 a 2 ) == liftTS a 1 liftTS a 2 for all animation a 1 and a 2 .

  TimedMonoid (TExtension d) d where liftTM d = (d, mempty) start (d, a) = (d, act (start mempty) a) stop (d, a) = (d, act (stop mempty) a)

Example.

  As program examples, we define circle :: d → Program3D d circle d = [Pen (Circle d) ] shift :: (d, d, d) → Program3D d shift (x 1 , x 2 , x 3 ) = [Transform (0, V3 x 1 x 2 x 3 ) ] cylinder :: d → Program3D d cylinder r l = circle 0 circle r shift (0, 0, l) circle r circle 0

6. 2 3D

 2 Position and Scene Three dimensional positions and figures, called Scene, are defined by type Position3D d = (Penshape d, Affine3D d, Affine3D d) initialPosition3D = (Empty, mempty, mempty) type Scene d = Set (V3 d, V3 d, V3 d)

6. 3

 3 3D Moves, Drawings and Semantics Types Three dimensional moves, drawings are defined much like in the 2D case: type Move3D d = Position3D d → Position3D d instance Monoid (Move3D d) where mempty = id mappend = flip (•) type Drawing3D d = Position3D d → Scene d instance Monoid (Drawing3D d) where mempty = const mempty mappend d 1 d 2 = \p → (d 1 p d 2 p) with a similar semigroup action and derived semi-direct product for semantics, defined by: instance SemigroupAction (Move3D d) (Drawing3D d) where act m d = d • m type Semantics3D d = (Move3D d, Drawing3D d) instance Monoid (Semantics3D d) where mempty = (mempty, mempty) mappend(m 1 , d 1 ) (m 2 , d 2 ) = (m 1 m 2 , d 1 act m 1 d 2 )6.4 Running 3D TurtleAs for running is concerned, we define the move and drawing semantics for 3D turtle basic actions and programs by: move3DA :: BasicAction3D d → Move3D d move3DA (Pen ps) = \ ( , , aff ) → (ps, aff , aff ) move3DA (Transform aff ) = \ (ps, p, aff0) → (ps, p, aff0 aff ) move3D :: Program3D d → Move3D d move3D p = mconcat (map move3DA p) draw3DA :: BasicAction3D d → Drawing3D d draw3DA (Pen Empty) = (const mempty) draw3DA (Pen (Circle d)) = (extrude d) draw3DA = mempty draw3D :: Program3D d → Drawing3D d draw3D p = let ( , d) = sem3D p in d sem3D :: Program3D d → Semantics3D d sem3D p = mconcat (map (\a → (move3DA a, draw3DA a)) p) with the extrude function defined by extrude :: d → Position3D d → Scene d extrude (Empty, , ) = mempty in the case the last turtle pen position is empty, extrude 0 (Circle 0, , ) = mempty in the case the last and current turtle pen position have zero radius hence both define a point 7 extrude d (Circle 0, ( , v), aff ) = span3D v (d, aff ) in the case the previous turtle pen shape has zero radius hence define a single point, extrude 0 (Circle d, aff , ( , v)) = cospan3D (d, aff ) v in the case the current turtle pen shape has zero radius hence define a single point, and extrude d 2 (Circle d 1 , aff1, aff2) = tube3D (d 1 , aff1) (d 2 , aff2) in the case both last and current turtle pen shapes are ellipse with non zero radius.Functions span3D, tube3D and cospan3D allow for triangulating the surface defined by extrusion from one position to another. Their semantics is informally depicted in Figure2.

Figure 2 .

 2 Figure 2. Triangulation examples respectively performed by functions span3D, tube3D and cospan3D.

Figure 3 .

 3 Figure 3. A scene image example.

Figure 4 .

 4 Figure 4. Associated triangulation.

Figure 5 .

 5 Figure 5. Associated surface normals.

  1 p is drawn by the first program, moving the turtle to a position m 1 p from which a second figure d 2 (m 1 p) is drawn.

	This is, indeed, the expected sequential composition of semantics
	we seek for.
	The function combine turns the type Semantics d into a monoid
	by
	instance Monoid (Semantics d) where
	mempty = (id, const mempty)
	mappend = combine
	Then, one can check that, for every turtle semantics s 1 and s 2 , we
	have
	applyS (s 1 s

2 ) == (applyS s 2 ) • (applyS s 1 )

  . First examples of timed moves are delay functions of the form \t → Just (t + d) for some constant value d. In other words, this function allows for modeling the fact that a turtle move may last d unit of time.The (monoid) composition of delays \t → Just (t + d 1 ) and \t → Just (t + d 2 ) equals delay \t → Just (t + d 1 + d 2 ). As expected, this reflects the fact that the two delay values have been summed.

with noticeable exceptions as detailed at the end of the paper

We replace both PenDown and PenUp (irreversible) commands of Logo by a single (reversible) TogglePen command. This is just a technical point since PenUp and Pen-Down command can also be turned into reversible commands provided their history is recorded in Turtle states.

Many instances defined in this pages clash with Haskell Prelude instances or even one on another. In our real implementation newtype embeddings are used. However, for the sake of clarity, they have been omitted throughout, as well as have been omitted most type constraints.

The ScopedTypeVariables option, necessary to remove any ambiguity on the type of mempty has been omitted in this instance.

The surface vertex colors are set to the vertex normal coordinates.

Acknowledgments

Thanks to Jean-Michaël Celerier, Pierre Bénard and Mikhail Ruskin for their help with modern 3D technology.


Christophe Bourlier

Gildas Kubické
email: gildas.kubicke@intradef.gouv.fr

Philippe Pouliguen
email: philippe.pouliguen@intradef.gouv.fr

come

Accelerated computation of the Physical Optics approximation for near-field single-and double-bounces backscattering

Accelerated computation of the Physical Optics approximation for near-field single-and doublebounces backscattering

Christophe Bourlier, Gildas Kubické and Philippe Pouliguen

Abstract-This paper deals with the acceleration of the Physical Optics (PO) approximation for the backscattering computation in the near and far field, for both single-(SB) and doublebounces (DB) PO contributions. It is based on physical arguments, which allow us to apply either a closed-form expression (rapid calculation) or a numerical double integration to evaluate the scattered field from sub-surfaces of a given mesh. The use of the closed-form expression is governed by two criteria, which must be extended to the bistatic case and also to the double bounce: the conventional Fraunhofer (related to the far-field zone) criterion and the chord error (related to the surface curvature) criterion. The proposed method is tested on a dihedral geometry.

Index Terms-Physical optics (PO), Radar cross section, highfrequency methods, fast algorithm

I. INTRODUCTION

F Or more clarity, the following acronyms are defined:

• DB: Double Bounce.

• IPO: Iterative Physical Optics.

• PO: Physical Optics.

• PO1: Physical Optics at the first order (contribution of the SB). • PO2: Physical Optics at the second order (contribution of the DB). • RCS: Radar Cross Section.

• SB: Single Bounce.

• SC: Scattering Coefficient.

• SD: Sub-surfaces Decomposition. In the low-and intermediate-frequency regimes, exact numerical methods, such as the method of moments (MoM) [START_REF] Harrington | Field Computation by Moment Methods[END_REF], can be used for solving the problem of scattering from a target. At very high frequencies, such methods become unattractive due to their high computational complexity. Fortunately, for such frequencies, asymptotic techniques, such as the ray-based [START_REF] Bhalla | 3D scattering center representation of complex targets using the shooting and bouncing ray technique: A review[END_REF] and PO-based methods, become applicable for structures with radii of curvature large compared to the wavelength. The ray methods provide a phenomenological solution, cast in terms of reflected and diffracted ray contributions, but suffer from high sensitivity to geometrical details and occasional failures.

Using the PO approximation for general (non convex) geometries, in order to capture the relevant scattering mechanisms, one may describe the solution as a series of contributions, each generated by its respective term in a series of source distributions on the scatterer's surface, named "bounces". The first term in the series, the "single bounce" contribution, provides a sufficient description of the scattering mechanism, asymptotically associated with the specular reflection, which can be significant for convex geometries. For concave geometries, accurate computation of the scattered field requires, in addition, to take into account at least the second-order term in the series, the "double-bounce" contribution, which asymptotically corresponds to two reflections on surfaces of the scatterer. For particular geometries, like a dihedral, the DB scattering can become the main contribution.

For some applications, the target can be close to the transmitter/receiver and then the Fraunhofer criterion is not satisfied, which means that the scattered field must be computed in the near field. The computation of the multi-bounce contributions involves multiple surface integrations. For instance, for the calculation of the SB contribution, two foldintegrations are required, whereas for the DB contribution, four fold-integrations are required. These integrations dominate the computational complexity of the entire PO solution and become a computational bottleneck for large scatterers.

For canonical geometries, the calculation of the SB integral can be done analytically [START_REF] Pouliguen | Analytical formulae for radar cross section of flat plates in near field and normal incidence[END_REF], [START_REF] Vallecchi | Physical Optics Curved-Boundary Dielectric Plate Scattering Formulas for an Accurate and Efficient Electromagnetic Characterization of a Class of Natural Targets[END_REF], [START_REF] Bourlier | Useful analytical formulae for nearfield monostatic radar cross section under the Physical Optics: far-field criterion[END_REF], [START_REF] Corbel | Rough Surface RCS Measurements and Simulations Using the Physical Optics Approximation[END_REF], [START_REF] Kubické | Bridging the Gap Between the Babinet Principle and the Physical Optics Approximation : Vectorial Problem[END_REF], [START_REF] Karaca | Closed-Form Time Domain PO Expressions of the Electric Field Scattered From PEC Objects Illuminated by an Electric Dipole[END_REF] by introducing special functions, especially in the near field and for complex shapes. For complex geometries, numerical techniques have also been developed [START_REF] Ludwig | Computation of radiation patterns involving numerical double integration[END_REF], [START_REF] Pouliguen | Radar reflectivity of a target illuminated by a spherical wave[END_REF], [START_REF] Boag | A Fast Physical Optics (FPO) Algorithm for High Frequency Scattering[END_REF], [START_REF] Zhang | Efficient Evaluation of the Physical-Optics Integrals for Conducting Surfaces Using the Uniform Stationary Phase Method[END_REF], [START_REF] Corucci | Near Field Physical Optics modelling for Concealed Weapon Detection[END_REF], [START_REF] An | Improved multilevel physical optics algorithm for fast computation of monostatic radar cross section[END_REF], [START_REF] Zhang | An alternative treatment of saddle stationary phase points in Physical optics for smooth surfaces[END_REF] to evaluate the SB integral, especially for different excitations [START_REF] Boag | A Fast Physical Optics (FPO) Algorithm for High Frequency Scattering[END_REF], [START_REF] Zhang | Efficient Evaluation of the Physical-Optics Integrals for Conducting Surfaces Using the Uniform Stationary Phase Method[END_REF] (see also [START_REF] Roudstein | Multilevel Physical Optics Algorithm for Near-Field Double-Bounce Scattering[END_REF] for a complete review of this method). For the calculation of the DB integral, it is impossible to derive closed-form expressions even for canonical geometries, without the introduction of simplifying assumptions [START_REF] De Adana | Calculation of the RCS From the Double Reflection Between Planar Facets and Curved Surfaces[END_REF], [START_REF] Cátedra | Physical Optics analysis of multiple interactions in large scatters using current modes[END_REF], [START_REF] Delgado | Analytical Field Calculation Involving Current Modes and Quadratic Phase Expression[END_REF], [START_REF] Kubické | Polarimetric bistatic signature of a faceted octahedron in high-frequency domain[END_REF], [START_REF] Kubické | High-frequency bistatic scattering by depolarizing, nearly omni-directional reflectors: higher order polyhedral reflectors[END_REF], [START_REF] Jackson | Analytic Physical Optics Solution for Bistatic, 3D Scattering From a Dihedral Corner Reflector[END_REF] (like the geometric optics approximation or by expanding the kernel into a Taylor series expansion). For any geometry, to overcome this issue, the multilevel physical optics (MLPO) [START_REF] Roudstein | Multilevel Physical Optics Algorithm for Near-Field Double-Bounce Scattering[END_REF] has been developed both in the near and far fields, for which the time saving is mainly obtained from the treatment of multiple source positions/angles/frequencies. Other accelerations of PO were proposed in the context of open-ended cavities with the Iterative Physical Optics (IPO)

A c c e p t e d M a n u s c r i p t method [START_REF] Obelleiro | An Iterative Physical Optics approach for analyzing the electromagnetic scattering by large open-ended cavities[END_REF], [START_REF] Obelleiro | A segmented approach for computing the electromagnetic scattering of large and deep cavities[END_REF], [START_REF] Boag | A fast iterative physical optics (FIPO) algorithm based on non-uniform polar grid interpolation[END_REF], [START_REF] Burkholder | Iterative Physical Optics for radar scattering predictions[END_REF], [START_REF] Thomet | Optimizing Iterative Physical Optics by using an ACA compression on interaction matrices[END_REF], [START_REF] Thomet | Low Computational cost method for scattering of large cavities based on ACA compression of Iterative Physical Optics[END_REF].

This paper addresses the general issue of the acceleration of the calculation of the SB and DB contributions, both in the near and far field, in the backscattering direction by developing the sub-surfaces decomposition physical optics "PO+SD" method. It consists in applying either a closedform expression (rapid calculation) or a numerical double integration to evaluate the scattered field from sub-surfaces of a given mesh. This idea has been tested in [START_REF] Pouliguen | Radar reflectivity of a target illuminated by a spherical wave[END_REF]. The starting point of this algorithm is similar to that of MLPO, which also decomposes the structure into sub-surfaces, but it is different for the calculations of the fields scattered by the sub-surfaces. To apply the closed-form expression, two criteria must be satisfied. The first one corresponds to the Fraunhofer criterion, which gives the distance r 0 , from which, both the transmitter (which can be a sub-surface) and the receiver (which can be a sub-surface) can be considered to be in the far field. For the monostatic case, this distance is often defined by r 0 ≥ 4d 2 /λ 0 , in which d is the largest dimension of the sub-surface and λ 0 , the electromagnetic wavelength. In this paper, the distance r 0 is derived rigorously for the bistatic case and extended to the monostatic case for the DB. The second criterion is related to the curvature of the sub-surface. The use of a closedform expression of the scattered field, without introducing special functions, implies that the curvature of the sub-surface can be neglected. This leads to the derivation of the "chord error" criterion for the bistatic case and its extension to the monostatic case for the DB.

This paper is organized as follows. Using the PO method, Section II derives the bistatic scattering coefficient (the RCS is obtained from the SC by taking its square modulus) for the SB and DB and simplifies them for the monostatic case. In Section III, method "PO+SD" is presented, in which the two criteria are derived for the SB and DB and the resulting monostatic SC is derived. Section IV presents numerical results of the monostatic SC and the time saving. The last section gives concluding remarks.

The time convention e -jωt is used throughout this paper.

II. SCATTERING COEFFICIENT COMPUTED BY THE PHYSICAL OPTICS IN NEAR FIELD

A. Bistatic case

In the near field, Pouliguen et al. [START_REF] Pouliguen | Analytical formulae for radar cross section of flat plates in near field and normal incidence[END_REF] showed that the POscattered magnetic field of the SB is expressed as (for a perfectly-conducting scatterer)

H s = - j λ 0 S I G(r s ) [ n ∧ H i ] ∧ rs e jk0rs r s dS, (1)
where

G(r) = 1 + j k 0 r . (2)
Here (see Fig. 1)

• λ 0 (k 0 = 2π/λ 0 ) is the electromagnetic wavelength in free space. • r s = r -R s , where -r s is the observation point defined from the point M on the surface of the scatterer. • The vector R s = (x s , y s , z s ) stands for the location of the receiver.

• n = n(r) is the normal to the surface S at the point M of coordinates r = (x, y, z).

• H i = H i (r)
is the incident magnetic field on the surface. • S I is the portion of the surface illuminated by the incident wave. The points of the surface are illuminated if ri • n < 0, where ri stands for the direction of the incident wave. • The boldface stands for a vector and the hat ˆindicates that the vector is unitary (û = u/ u ). As shown in Fig. 2, the scattered magnetic field at the second order of the DB, H s,12 , resulting from the scattering from scatterer 2 illuminated by scatterer 1, is obtained from Eq. (1), in which the incident magnetic field H i is substituted for H s,1 , the magnetic field scattered by object 1. This leads to

H s,12 = - 1 λ 2 0 S I1 ,S I2 G(r 12 )G(r s,2 )e jk0(r12+rs,2) r 12 r s,2 × n2 ∧ [( n1 ∧ H i,1 ) ∧ r12 ] ∧ rs,2 dS 1 dS 2 , (3)
where n1 = n(r 1 ), n2 = n(r 2 ), H i,1 = H i (r 1 ), r 12 = r 2 -r 1 and r s,2 the observation point on the surface defined from scatterer 2. The subscripts "1" and "2" are added to distinguish the scatterers 1 and 2. The points of the surface S 2 are illuminated (corresponding to the surface S I2 ) from S 1 if r12 • n2 < 0.

The incident magnetic field can be expressed as

H i (r) = V i Z 0 e jk0r r pi , (4)
where V i is the voltage source, Z 0 the wave impedance in free space and pi the polarization state of the transmitter, which can be either vertical (v i ) or horizontal ( ĥi ). In addition, in the near field, the bistatic SC (the RCS is obtained from the SC by taking its square modulus) can be defined as [START_REF] Pouliguen | Analytical formulae for radar cross section of flat plates in near field and normal incidence[END_REF]

ρ pips = 2 √ πR i R s Z 0 H s • ps V i , (5)
where {R i,s } are the distances from the transmitter and the receiver, respectively, to the origin O, and ps is the polarization state of the receiver, which can be either vertical (v s ) or horizontal ( ĥs ). The vectors are defined by ĥi,s = ri,s ∧ ẑ/ ri,s ∧ ẑ , where ẑ = (0, 0, 1) and vi,s = ĥi,s ∧ ri,s . Then, the substitution of Eqs. ( 4) and (1) into Eq. ( 5) leads to

ρ pips = 2 √ π jλ 0 S I [( n ∧ pi ) ∧ rs ] • ps × G(r s )R i R s r i r s e jk0(ri+rs) dS, (6)
where r i = r -R i , in which R i = (x i , y i , z i ) stands for the location of the transmitter. The bistatic SC of the DB is derived from the substitution of Eqs. ( 4) and (3) into Eq. ( 5) (in which

H s = H s,12 , R i = R i,1 and R s = R s,2 ), leading to ρ pips,12 = - 2 √ π λ 2 0 S 1I ,S 2I R i,1 R s,2 G(r 12 )G(r s,2 ) r i,1 r 12 r s,2 × e jk0(ri,1+r12+rs,2) × n2 ∧ n1 ∧ pi,1 ∧ r12 ∧ rs,2 • ps,2 dS 1 dS 2 . (7)
Eqs. ( 6) and ( 7) are valid if the following conditions are satisfied:

• To apply the PO approximation, the radii of curvature of the two surfaces are much larger than the wavelength λ 0 . • For the definition of the SC, the electromagnetic field radiated by the antenna has a spherical wave structure. • For the definition of the SC, the scattered electromagnetic wave field has a spherical wave structure, locally plane on the receiving antenna. In conclusion, the total SC up to the second order is obtained from four sub-SCs

ρ pips,tot = ρ pips,1 + ρ pips,2 + ρ pips,12 + ρ pips,21 , (8) where
• ρ pips,1 stands for the SC of scatterer 1 assumed to be alone (obtained from Eq. ( 6) by adding the subscript 1). • ρ pips,2 stands for the SC of scatterer 2 assumed to be alone (obtained from Eq. ( 6) by adding the subscript 2). • ρ pips,12 stands for the SC of scatterer 2 illuminated by scatterer 1. • ρ pips,21 stands for the SC of scatterer 1 illuminated by scatterer 2.

B. Monostatic case

For a monostatic configuration, r s = r i , Eq. ( 6) leads to

ρ pips = 2 √ π jλ 0 S I (r i • n) pi G(r i )e j2k0ri R 2 i r 2 i dS, (9)
where

( n ∧ pi ) ∧ ri = (r i • n) pi -(r i • pi ) n = (r i • n) pi .
In addition, (r i,1 , r s,2 ) involved in Eq. ( 7) becomes (r i,2 , r s,1 ) for ρ pips,21 , in which r s,2 = r i,2 and r s,1 = r i,1 . Then

ρ pips,12 + ρ pips,21 = - 2 √ π λ 2 0 S 1I ,S 2I dS 1 dS 2 R i,1 R i,2 G(r 12 )G pips r i,1 r 12 r i,2 e jk0(ri,1+r12+ri,2) , (10)
where

G pips = G(r i,2 ) n2 ∧ n1 ∧ pi,1 ∧ r12 ∧ ri,2 • pi,2 -G(r i,1 ) n1 ∧ n2 ∧ pi,2 ∧ r12 ∧ ri,1 • pi,1 . (11)
To accelerate the computations of integrals ( 9) and ( 10), the surface decomposition method is presented in the next section.

III. EVALUATION OF SUB-SURFACES CONTRIBUTIONS

This method consists in applying either a closed-form expression (rapid calculation) or a numerical double integration (done from a conventional trapezoidal rule) to evaluate the scattered field from the sub-surfaces {S n } of a given mesh. The closed-form expression is valid if the bistatic Fraunhofer and chord error criteria are both satisfied. They are derived thereafter for the single and double bounces.

A. Bistatic Fraunhofer and chord error criteria

From Fig. 1, one has r i = r i + δ. If δ = δ r i = r i , then a Taylor series expansion up to the second order over δ and around zero leads to

r i ≈ r i + δ cos φ i + δ 2 2r i sin 2 φ i , (12)
where φ i = ( δ, r i ). The δ term is related to the local behavior of a plane wave, whereas δ 2 is related to the local behavior of a spherical wave. The Fraunhofer criterion is obtained from Eq. ( 12) by neglecting the δ 2 term. In other words, this approximation is satisfied if δ 2 sin 2 φ i /(2r i ) does not exceed A c c e p t e d M a n u s c r i p t λ 0 /n 0 (typically n 0 is an integer ranging from 8 to 16) where λ 0 is the radar wavelength. This leads for k 0 (r i + r s ), to

δ 2 2 sin 2 φ i r i + sin 2 φ s r s ≤ λ 0 n 0 . ( 13
)
The maximum value of δ, named ∆, equals

∆ = max (x -x ) 2 + (y -y ) 2 + (z -z ) 2 ≈ d 2 + d 2 + d 2 (γ 2 x + γ 2 y ) ≈ |d | 2 + γ 2 x + γ 2 y , ( 14
)
where d is the radius of the circle circumscribing the subsurface and z -z is expressed from Eq. ( 17). Then, the subsurface can be considered to be in far field if

r is ≥ n 0 ∆ 2 4λ 0 , 1 r is = 1 2 sin 2 φ i r i + sin 2 φ s r s . (15)
In other words, condition ( 15) is satisfied if the arithmetic mean of the curvature radii weighted by sin 2 φ i,s of the incident and scattered wave is smaller than 4λ 0 /(n 0 ∆ 2 ).

For a monostatic configuration, r s = r i and φ s = φ i , which leads to r i ≥ n 0 ∆ 2 sin 2 φ i /(4λ 0 ). For n 0 = 16 and sin φ i = 1, the conventional Fraunhofer criterion (for RCS application) is retrieved,

r i ≥ 4∆ 2 /λ 0 ≈ 8d 2 /λ 0 .
If criterion ( 15) is satisfied, from Eq. ( 12), r i ≈ r i +δ cos φ i and

r i ≈ r i + (x -x )r ix + (y -y )r iy + (z -z )r iz , (16)
where r i = r i (r ix , r iy , r iz ) and the point O has coordinates (x , y , z ).

To obtain a closed-form expression of the double integral over x and y, the term z -z is expanded as

z -z ≈ (x -x )γ x + (y -y )γ y , (17)
where γ x = ∂z/∂x| x=x ,y=y and γ y = ∂z/∂y| x=x ,y=y . Using the same way as for criterion [START_REF] Corucci | Near Field Physical Optics modelling for Concealed Weapon Detection[END_REF], expansion ( 17) is valid if

|r iz + r sz | F is < λ 0 n 0 , ( 18
)
where

F is = γ xx (x -x ) 2 2 + γ yy (y -y ) 2 2 + γ xy (x -x )(y -y ) , ( 19
) and γ xx = ∂ 2 z/∂x 2 | x=x ,y=y , γ yy = ∂ 2 z/∂y 2 | x=x ,y=y and γ xy = ∂ 2 z/∂x∂y| x=x ,y=y .
Introducing d , Eq. ( 18) becomes

|r iz + r sz | d 2 γ xx 2 + γ yy 2 + γ xy < λ 0 n 0 . ( 20
)
Criterion [START_REF] Kubické | Polarimetric bistatic signature of a faceted octahedron in high-frequency domain[END_REF] is commonly named the "chord error" criterion.

As a conclusion, if criteria ( 13) and ( 20) are satisfied, from Eqs. ( 16) and [START_REF] De Adana | Calculation of the RCS From the Double Reflection Between Planar Facets and Curved Surfaces[END_REF], the distance

r i + r s -(r i + r s ) = δ is is approximated by δ is = (x -x ) [r ix + r sx + γ x (r iz + r sz )] + (y -y ) r iy + r sy + γ y (r iz + r sz ) . ( 21
)
B. Resulting scattering coefficient and complexity of PO1

In Eq. ( 6), the kernel omitted of the phase term e jk0(ri+rs) is approximated at the point O (x , y , z ). Then, this term becomes independent of the integration variables x and y. This leads to

ρ pips = 2 √ π jλ 0 N n=1 n ∧ p i ∧ r s • p s G(r s )R i R s r i r s × e jk0(r i +r s ) S I e jk0δ is dS , (22)
where the sum is done over all the sub-surfaces, which satisfy criteria [START_REF] Corucci | Near Field Physical Optics modelling for Concealed Weapon Detection[END_REF] and [START_REF] Kubické | Polarimetric bistatic signature of a faceted octahedron in high-frequency domain[END_REF]. The symbol prime indicates that the variables are constant for a given sub-surface but depend on the sub-surface n.

If the sub-surface has a rectangular shape of lengths l x and l y with respect to the directions x and y, respectively, then

S I e jk0δ is dS = l x l y sinc a x l x 2 sinc a y l y 2 , (23)
where

a x = r ix + r sx + γ x (r iz + r sz ) a y = r iy + r sy + γ y (r iz + r sz ) , sinc(x) = sin(x)/x. (24)
The double integration is done analytically instead of numerically and the resulting complexity is O(M x M y )η PO1+SD , where M x and M y are the numbers of the surface samples with respect to the x and y directions, respectively and η PO1+SD ∈ [0; 1] is defined as

η PO1+SD = 1 M x M y N n=1 M n , (25)
where

M n = 1 if Eqs. (13) and (20) are satisfied else M x,n M y,n , (26)
where M x,n and M y,n are the numbers of the sub-surface samples n with respect to the x and y directions, respectively. If all the sub-surfaces do not satisfied the two criteria, then

η PO1+SD = ( n M x,n M y,n )/(M x M y ) = 1
, corresponding to the complexity of the conventional PO1.

C. Resulting scattering coefficient and complexity of PO2

In this subsection, the formulation is extended to the double bounce.

In Eq. [START_REF] Pouliguen | Radar reflectivity of a target illuminated by a spherical wave[END_REF], the distance r i,1 +r 12 +r i,2 can be approximated by

r i,1 + r 12 + r i,2 ≈ r i,1 + r 12 + r i,2 + δ 1 (cos φ i,1 -cos φ 1 ) + δ 2 (cos φ i,2 + cos φ 2 ), (27)
where

δ 1 = ---→ M 1 O 1 , δ 2 = ---→ M 2 O 2 , φ i,1 = ( δ 1 , r i,1 ), φ i,2 = ( δ 2 , r i,2 ), φ 1 = ( δ 1 , r 12 ) and φ 2 = ( δ 2 , r 12 ), in which r 12 = ---→ O 2 O 1 and {O 1 , O 2 }
are the origins of the subsurfaces of the scatterers 1 and 2, respectively. In addition,

r 12 = r 12 , r i,1 = ---→ O 1 P i , r i,2 = ---→ O 2 P i .
A c c e p t e d M a n u s c r i p t

For the phase k 0 (r i,1 + r 12 + r i,2 ), Eq. ( 27) is valid if Using the same approach as for PO1, from Eq. ( 27), the chord error criterion becomes

∆ 2 1 sin 2 φ i,1 2r i,1 + ∆ 2 2 sin 2 φ i,2 2r i,2 + (∆ 1 sin φ 1 -∆ 2 sin φ 2 ) 2 2r 12 ≤ λ 0 n 0 , (28)
(r iz,1 -r 12z )d 2 1 F 1 + (r iz,2 + r 12z )d 2 2 F 2 ≤ λ 0 n 0 , (29)
where 29), the component rz = r • ẑ/ r .

F p = γ p,xx 2 + γ p,yy 2 + γ p,xy , (30)
In Eq. ( 10), the kernel omitted of the phase term equals e jk0(ri,1+r12+ri,2) , is approximated at the points O 1 (x 1 , y 1 , z 1 ) and O 2 (x 2 , y 2 , z 2 ). Then, this term becomes independent of the integration variables (x 1 , y 1 , x 2 , y 2 ). This leads to Eq. ( 10

) = - 2 √ π λ 2 0 N1 n1=1 N2 n2=1 R i,1 R i,2 G pips G(r 12 )
r i,1 r 12 r i,2 λ 0 × e jk0(r i,1 +r 12 +r i,2 ) S 1I ,S 2I e jk0δ 12 dS 1 dS 2 , (31)

where

δ 12 = p=2 p=1 (x p -x p )a x,p + (y p -y p )a y,p , (32)
in which

       a x,1 = r ix,1 -r 12x + γ x,1 (r iz,1 -r 12z ) a y,1 = r iy,1 -r 12y + γ y,1 (r iz,1 -r 12z ) a x,2 = r ix,2 + r 12x + γ x,2 (r iz,2 + r 12z ) a y,2 = r iy,2 + r 12y + γ y,2 (r iz,2 + r 12z ) , (33)
where γ x,p = ∂z p /∂x p | xp=x p ,yp=y p and γ y,p = ∂z p /∂y p | xp=x p ,yp=y p .

In Eq. ( 31), the sums are done over all the sub-surfaces of the scatterers p = {1, 2}, which satisfy criteria [START_REF] Thomet | Low Computational cost method for scattering of large cavities based on ACA compression of Iterative Physical Optics[END_REF] and [START_REF] Thomet | Improvment of the Iterative Physical Optics using the Physical Optics shadow radiation[END_REF]. The symbol prime indicates that the variables are constant for a given sub-surface but depend on the sub-surfaces n 1 and n 2 .

If the sub-surface of the scatterer p has a rectangular shape of area l x,p × l y,p , then the quadruple integrations can be done analytically, leading to

The resulting complexity is O(M x,1 M y,1 M x,2 M y,2 )η PO2+SD , where M x,p and M y,p are the numbers of samples on the surface p with respect to the x and y directions, respectively, and η PO2+SD ∈ [0; 1] is defined by

η PO2+SD = 1 M x,1 M y,1 M x,2 M y,2 N1 n1=1 N2 n2=1 M n1,n2
where M x,ni and M y,ni are the numbers of samples on the sub-surface n i of the scatterer i with respect to the directions x and y, respectively. If all the sub-surfaces dot not satisfied the two criteria, then n1 n2 = M x,n1 M y,n1 M x,n2 M y,n2 and η PO2+SD = 1, corresponding to the complexity of the conventional PO2.

For the DB and from Eq. ( 34), the integral is expressed as the product of two SB identical functions with different arguments related to sub-surfaces of the scatterers 1 and 2.

D. Numerical implementation

The main steps of the algorithm are : 1) Mesh the geometry.

2) Decompose the surface into sub-surfaces of rectangular shapes. For a more general polygon shape (with the condition that two adjacent sub-surfaces must remain connected), the formula of Gordon [START_REF] Gordon | Far-field approximations to the Kirchhoff-Helmholtz representations of scattered fields[END_REF] can be applied. 3) One loop on the number of sub-surfaces for PO1 and two loops on the number of sub-surfaces of each scatterer. 4) Calculate the center coordinates and the circumscribed radii of the sub-surfaces . 5) If the Fraunhofer and chord criteria are satisfied, then the integration is analytical, else numerical. 6) Cumulated sum of the resulting contribution. 7) Stop the loop(s).

IV. NUMERICAL RESULTS

A. Single bounce

First, the monostatic SC of the SB is computed from PO1. The geometry is shown in Fig. 3: a corrugated surface defined as

z(x, y) = a cos 2πx Λ x cos 2πy Λ y , ( 37
)
where a is the amplitude and {Λ x , Λ y } the periods with respect to the x and y directions, respectively. The surface area is L x × L y and the wavelength equals λ 0 = 3 cm (the Radar frequency is f = 10 GHz). In Fig. 3, a = λ 0 , L x = L y = 400λ 0 , Λ x = Λ y = L x /4 = 100λ 0 . The maximum values of the surface slopes are γ x,max = max(|∂z/∂x|) = 2πa/Λ x = 0.0628 (a > 0) and γ y,max = max(|∂z/∂y|) = 2πa/Λ y = 0.0628 with respect to the x and y directions, respectively. These values are chosen such as the multiple reflections on the surface can be neglected. The sampling step with respect to the x and y directions, respectively, is ∆x = ∆y = λ 0 /10 and M x = M y = 4001. In addition, the positions of the receiver/transmitter are x i = y i = 0. z i ranges from 0.06 to 40, 000 m > 4(L 2

x + L 2 y )/λ 0 = 38, 400 m, which is the conventional Fraunhofer distance.

At the top, Fig. 4 plots the square modulus |ρ| 2 of the monostatic SC versus z i . In addition, at the bottom, to better see the differences, the ratio |ρ PO1+SD /ρ PO1 | 2 is plotted in dB scale (becoming a difference). In the legend, the labels

• "PO1" means that the SC is computed from Eq. ( 9).

• "PO1+SD (l x × l y )" means that the SC is computed from Eqs. ( 9) and ( 22) and n 0 = 10. In addition, the lengths l x = l y = {1, 2, 4}λ 0 (M x,n = M y,n = {10, 20, 40} ∀p and N = {400 2 , 200 2 , 100 2 }) are constant versus z i . As z i increases, the difference between SD+PO1 and PO1 decreases and does no exceed ±2 dB. Computing time ratio t PO1 /t PO1+SD versus z i . l x = l y = (1, 2, 4)λ 0 . The parameters are the same as in Fig. 4. Fig. 6. Percentages of the sub-surfaces, which satisfy the Fraunhofer (C 1 , Eq. ( 13)), the chord error (C 2 , Eq. ( 20)) and (C 1 and C 2 ) criteria, versus z i . (lx, ly) = (2, 2)λ 0 . Fig. 5 plots the computing time ratio t PO1 /t PO1+SD versus z i . To explain the behaviours of the curves, for l x = l y = 2λ 0 , Fig. 6 plots the percentages of the sub-surfaces {p C1 , p C2 , p C1,2 }, which satisfies the Fraunhofer (C 1 , Eq. ( 13)) and the chord error (C 2 , Eq. ( 20)) criteria, respectively, versus z i . Fig. 7 plots the same variations as in Fig. 6, but for l x = l y = 4λ 0 and the y scale differs. As we can see in Fig. 5, for l x = l y = {1, 2}λ 0 , the time saving increases with z i since the Fraunhofer C 1 criterion, p C1 , increases as z i increases, whereas the chord error C 2 criterion is always satisfied (p C2 = 100% ∀z i ). The same remark holds for l x = l y = 1λ 0 (not shown).

In addition, for small z i and l x = l y = 4λ 0 , the agreement in Fig. 4 is better because most contributions of the subsurfaces are calculated without acceleration (comparison of Fig. 6 with Fig. 7, where p C1,2 = {99.7, 95}% for z i near zero, respectively). The computing time ratio t PO1 /t PO1+SD is directly related to the complexity. From Eq. ( 25), it is expressed as 1/η PO1+SD . For all the N sub-surfaces, if the two criteria are satisfied then the time ratio is M x M y /N . For l x = l y = (1, 2)λ 0 , we obtain {4001 2 /400 2 , 4001 2 /200 2 } ≈ {100, 400}. From Fig. 5, the asymptotic values are nearly these values divided by 2.

Fig. 5 also shows that the choice of the lengths (l x , l y ) must be optimized with respect to r i , but cannot exceed an upper limit obtained from the chord error criterion [START_REF] Kubické | Polarimetric bistatic signature of a faceted octahedron in high-frequency domain[END_REF]. Typically, for a plane surface and for r i → ∞, the numbers (N, l x , l y ) tend to (1, L x , L y ), but if the surface is strongly curved, these asymptotic values cannot be reached.

Indeed, for l x = l y = 4λ 0 , Fig. 5 shows that the time saving increases and next decreases. This behaviour is directly related to Fig. 7. For z i close to 0, the criterion C 1 is not fully satisfied p C1 ≈ 95%, whereas the criterion C 2 is fully satisfied p C2 ≈ 100%. As z i increases, p C1 increases to tend toward 100%, whereas p C2 decreases and then p C1,2 decreases to tend toward 55%. Then, there is an antagonist behaviour between C 1 and C 2 , which produces an optimized value, for which the method is the most efficient.

From Eq. ( 20), for r i → ∞, r iz → 1, the criterion is then satisfied for the monostatic case if d = d max = λ 0 /(2n 0 ) max(γ xx /2 + γ yy /2 + γ xy ) ≈ 3.559λ 0 . This value is consistent with Fig. 5.

At the top of Fig. 8, the percentages of the sub-surfaces, which satisfy the Fraunhofer (C 1 , Eq. ( 13)) and the chord error (C 2 , Eq. ( 20)) criteria, are plotted versus z i . In addition, in the legend, the label "C 10 " means that in Eq. ( 13) sin φ i = sin φ s = 1, which corresponds to the conventional Fraunhofer criterion. At the bottom, the corresponding computing time ratio t PO1 /t PO1+SD is plotted versus z i .

If the conventional criterion C 10 is applied instead of C 1 , then for small values of z i , the percentage of sub-surfaces satisfying C 10 is smaller than that satisfying C 1 because C 10 is an upper limit of C 1 . Then, the resulting computing time is approximately multiplied by two. 13)) and the chord error (C 2 , Eq. ( 20)) criteria, versus z i . (lx, ly) = (2, 2)λ 0 . Bottom: Corresponding computing time ratio t PO1 /t PO1+SD versus z i . Fig. 9. Geometry for the computation of the monostatic SC of the DB. The scatterer 1 is defined by y < 0 and the scatterer 2 by y ≥ 0.

B. Double bounce

In this subsection, the monostatic SC ρ of the DB is computed from PO2 1) Case a = 0: To produce a DB, a dihedral composed of two connected rectangular plates of same lengths (L x , L y ) and shown in Fig. 9, is considered. The scatterer 1 is defined by y < 0 and the scatterer 2 by y ≥ 0. The lengths with respect to the x and y directions are L x = 20λ 0 and L y = 10 √ 2λ 0 , respectively, in which λ 0 = 3 cm (the Radar frequency is f = 10 GHz). The angles between the y direction and the plates are π/4 + π/2 and π/4, respectively. The sampling step with respect to the x and y directions is ∆x = ∆y = λ 0 /10 and (M x,1 , M y,1 , M x,2 , M y,2 ) = (201, 101, 201, 101). In addition, the positions of the receiver/transmitter are x i = y i = 0 and z i ranges from 0.1 to 100 m > 4 max(2L y cos(π/4), L x ) 2 /λ 0 = 48 m, which is the conventional Fraunhofer distance. A the top, Fig. 10 plots the square modulus (RCS) of the monostatic SC, ρ, versus z i . In addition, at the bottom, to better see the differences, the ratio |ρ PO2+SD /ρ PO2 | 2 is plotted in dB scale (becoming a difference). In the legend, the labels

• "PO2" means that the SC is computed from Eq. [START_REF] Pouliguen | Radar reflectivity of a target illuminated by a spherical wave[END_REF].

• "PO2+SD (l x,1 × l y,1 )" (l x,1 = l x,2 , l y,1 = l y,2 ) means that the SC is computed from Eqs. ( 9) and [START_REF] Jackson | Analytic Physical Optics Solution for Bistatic, 3D Scattering From a Dihedral Corner Reflector[END_REF] and n 0 = 10.

• "PO1+PO2" means that the SC is computed from Eqs. ( 9) and (10). • "MoM" means that the SC is computed from the fullwave Method of Moments, for which the sampling step is λ 0 /8.

The sub-surface lengths are l x,1 = l y,1 = l x,2 = l y,2 = {0.5, 1}λ 0 and are constant according to z i . In addition, the numbers of the sub-surfaces samples of scatterers i = {1, 2}

are M x,pi = M y,pi = {5, 10} ∀(p, i) and the numbers of the sub-surfaces of scatterers 1 and 2 are P 1 = P 2 = {800, 200}. For small values of z i , Fig. 10 shows a good agreement between PO2 and PO2+SD because, as shown in Fig. 11 (and also Fig. 12 for (l x,p , l y,p ) = (0.5, 0.5)λ 0 (p = {1, 2})), most of the sub-surface pair contributions are computed without acceleration. Next, as z i increases, the difference increases and in the Fraunhofer region, the opposite effect occurs and the difference does not exceed 0.5 dB.

Fig. 10 also reveals that for small values of z i , the main contribution is given from PO1, whereas the opposite effect occurs for large values of z i . In addition, a satisfactory agreement is obtained between the MoM and "PO1+PO2" results, especially in far field. In near field, the difference can be attributed to the fact that PO does not account for the edge diffraction.

Fig. 11 plots the computing time ratio t PO2 /t PO2+SD versus z i . As expected, the time saving is significant for hight values of z i and it is larger for l x,1 = l y,1 = l x,2 = l y,2 = l x = λ 0 , but the range over z i is larger for l x = 0.5λ 0 .

Like PO1, to give an explanation, Fig. 12 plots the percentages {p C1 , p C2 , p C1,2 } of the pairs of sub-surfaces (between scatterers 1 and 2), which satisfy the Fraunhofer (C 1 , Eq. ( 28)), the chord error (C 2 , Eq. ( 29)) and (C 1 and C 2 ) criteria, versus z i . As expected, for a = 0, the chord error criterion is satisfied for any z i (p C2 = 100%), whereas p C1 differs from 100% for z i = 0, next decreases (it is not a monotonic function of z i ) and tends toward 100% for large values of z i . For l x = λ 0 , simulations not reported here, showed that p C1,2 has the same behaviour as that computed for l x = 0.5λ 0 with the main difference that the time saving is significant for a larger value of z i . This explains the behaviours in Fig. 11.

The computing time ratio t PO2 /t PO2+SD is directly related to the complexity. From Eq. (35), it is expressed as 1/η PO2+SD . For all the N = N 1 N 2 sub-surfaces, if the two criteria are satisfied then the time ratio is M x,1 M y,1 M x,2 M y,2 /N . For l x = (0.5, 1)λ 0 , we obtain {(201 × 101) 2 /640, 000, (201 × 101) 2 /40, 000} ≈ 10 4 ×{0.64, 10.3}. From Fig. 11, we obtain {143, 290}, which is 4.5 and 35.5 times smaller that the theoretical time savings. This difference can be attributed by the fact that the codes are written in MatLab (not optimized for the loops) and that many additional tests and auxiliary variables are introduced for SD+PO2 method.

2) Case a = 0: Fig 13 plots the same geometry as in Fig. 9 but a = 0.5λ 0 and the periods Λ x = L x and Λ y = L y /2. Fig. 14 plots the same variations as in Fig. 10 but the object is defined from Fig. 13 (a = 0). Fig. 15 plots the corresponding computing time ratio t PO2 /t PO2+SD versus z i . Fig. 16 plots the same variations as in Fig. 12, but the object is defined from Fig. 13 (a = 0).

The comparison of Fig. 14 with Fig. 10 shows that the RCS computed from PO1+PO2 strongly differs and the difference increases slightly in comparison to a = 0. Since a = 0, the percentage of the sub-surface pairs p C2 associated to the chord error criterion C 2 can differ from 100%. For l x = 0.5λ 0 , simulations not reported here, showed that p C1,2 is nearly the same as that plotted in Fig. 12 and that p C2 = 100% for any z i . This explains why the time saving is nearly the same in For l x = λ 0 , Fig. 16 shows that {p C1 , p C2 , p C1,2 } strongly differ from those plotted in Fig. 12. Firstly, p C2 has a similar behaviour as that obtained from PO1 and depicted in Fig. 7: it is a decreasing function of z i . In addition, since a = 0, for small values of z i , p C2 is smaller than 100%, unlike in Fig. 12. Secondly, p C1 is not a monotonic function and does not tend toward 100% because sub-surface pairs are in nearfield independently of z i . The resulting percentage p C1,2 ≤ (p C1 , p C2 ) (by definition) and since p C1,2 = p C1 or p C1,2 = p C2 , the sub-surface pairs that satisfy criteria C 1 and C 2 are not the same. Indeed, C 1 and C 2 are independent. In far field, from Fig. 16 a percentage of p C1,2 ≈ 1% is not enough to have a significant time saving (t PO2 /t PO2+SD ≈ 1.16 instead of 147 for l x = 0.5λ 0 ) in Fig. 15.

V. CONCLUSION

To accelerate the computation of the SB and DB contributions of PO backscattering (in the near and far fields), the SD is developed. It is based on two physical criteria, the Fraunhofer and chord criteria, for which their derivations are generalized in this paper. For the computation of the SB contribution, the resulting algorithm PO1+SD provides significant time savings in comparison to the direct computation of the double integral. For the computation of the DB contribution, PO2+SD method is also very efficient in terms of time saving especially in far field, for which the computing time can be reduced by a factor ranging from 100 to 300 in comparison to the conventional quadruple numerical integrations.

For objects with more complex shapes, e.g, object with self-shadowing, the visible function can be computed from ray-based techniques or similar algorithms. For any scattering problem solved from the PO, it is important to underline that the computation of the visible function is inherent to the problem whatever the method used to accelerate the PO. In addition, recently, A. Thomet et al. [START_REF] Thomet | Improvment of the Iterative Physical Optics using the Physical Optics shadow radiation[END_REF] showed that the shadowing effect can be modeled by using the Physical Optics shadow radiation. This avoids to calculate the visible function from ray-tracing based algorithms which can be time consuming.

The SD algorithm is highly parallelizable and can be further accelerated if it is implemented on multi-cores or GPU architectures. For PO1, all the integrals of the subsurfaces can be evaluated independently and then can be parallelized. The RAM requirement depends on the choice of the integration technique. From a simple trapezoidal rule, this requirement is near zero. For PO2, all the integrals between the sub-surfaces of the scatterers 1 and 2 can also be evaluated independently and they can be parallelized. From a simple trapezoidal rule, this memory requirement is near zero. To sum up, only the fields on the sub-surfaces must be stored and for a parallelized version of the code, the memory requirement is proportional to the number of sub-surfaces.

The most restrictive criterion (chord error) is related to the curvature of the surface. To overcome this issue, the prospect of this paper is to incorporate the spherical structure up to the second order (in Eq. ( 12), the term in δ 2 would be kept) of the incident wave into the calculation of the surface integral. Then, as shown in [START_REF] Bourlier | Useful analytical formulae for nearfield monostatic radar cross section under the Physical Optics: far-field criterion[END_REF], in Eqs. ( 23) and (34), the sinc function is changed by the modified Fresnel integrals, but simplifying assumptions must be introduced. Another means is to extend the algorithm to multi levels, which consists in dividing sub-surfaces into sub-sub-surfaces of smaller areas, for which the two criteria become valid.

Fig. 1 .

 1 Fig. 1. Geometry for the single bounce case (PO1).

A c c e p t e d M a n u s c r i p t

Fig. 2 .

 2 Fig. 2. Geometry to calculate the scattered field from PO2 (double bounce).

  where ∆ p = |d i | 2 + γ 2 p,x + γ 2 p,y . In addition, γ p,x = ∂z p /∂x p | xp=x p ,yp=y p , γ p,y = ∂z p /∂y p | xp=x p ,yp=y p and d p is the radius of the circle circumscribing the prime sub-surface of the scatterer p = {1, 2}. This corresponds to the monostatic Fraunhofer criterion extended to the DB.

  p=2 p=1 l x,p l y,p sinc a x,p l x,p 2 sinc a y,p l y,p 2 .

.

Fig. 3 .

 3 Fig. 3. Geometry for the computation of the monostatic SC of the SB. At the bottom, the full line is for x = 0 (left) and y = 0 (right). a = λ 0 = 3 cm, Lx = Ly = 400λ 0 and Λx = Λy = Lx/4 = 100λ 0 .

Fig. 5 .

 5 Fig. 5.Computing time ratio t PO1 /t PO1+SD versus z i . l x = l y = (1, 2, 4)λ 0 . The parameters are the same as in Fig.4.

Fig. 7 .

 7 Fig.7. Same variations as in Fig.6, but for l x = l y = 4λ 0 and the y scale differs.

Fig. 8 .

 8 Fig. 8. Top: Percentages of the sub-surfaces, which satisfied the Fraunhofer (C 1 , Eq. (13)) and the chord error (C 2 , Eq. (20)) criteria, versus z i . (lx, ly) = (2, 2)λ 0 . Bottom: Corresponding computing time ratio t PO1 /t PO1+SD versus z i .

  A c c e p t e d M a n u s c r i p t

Fig. 10 .

 10 Fig. 10. Top: Square modulus |ρ| 2 in dBm 2 scale of the monostatic SC computed from PO2 and PO1+PO2 (without acceleration) versus z i . Bottom: Difference in dB scale |ρ PO2+SD | 2 (dB) -|ρ PO2 | 2 (dB) (against PO2 without acceleration). l x,1 = l y,1 = l x,2 = l y,2 = {0.5, 1}λ 0 and n 0 = 10.

Fig. 11 .

 11 Fig. 11. Computing time ratio t PO2 /t PO2+SD versus z i .

Fig. 12 .

 12 Fig. 12. Percentage {p C 1 , p C 2 , p C 1,2 } of the pairs of the sub-surfaces (between the scatterers 1 and 2), which satisfy the Fraunhofer (C 1 , Eq. (28)), the chord error (C 2 , Eq. (29)) and (C 1 and C 2 ) criteria, versus z i . (l x,p , l y,p ) = (0.5, 0.5)λ 0 (p = {1, 2}).

Fig. 13 .

 13 Fig. 13. Same geometry as in Fig. 9 but a = 0.5λ 0 and the periods Λx = Lx and Λy = Ly/2.

Figs. 11

 11 Figs. 11 and 15.

Fig. 14 .

 14 Fig. 14. Same variations as in Fig. 10 but the object is defined from Fig. 13 (a = 0).

Fig. 15 .

 15 Fig. 15. Computing time ratio t PO2 /t PO2+SD versus z i . The object is defined from Fig. 13 (a = 0).

  A c c e p t e d M a n u s c r i p t

Fig. 16 .

 16 Fig.[START_REF] Roudstein | Multilevel Physical Optics Algorithm for Near-Field Double-Bounce Scattering[END_REF]. Same variations as in Fig.12, but the object is defined from Fig.13(a = 0) and (l x,p , l y,p ) = (1, 1)λ 0 (p = {1, 2}).

  and γ p,xx = ∂ 2 z p /∂x 2 p | xp=x p ,yp=y p , γ p,yy = ∂ 2 z p /∂y 2 p | xp=x p ,yp=y p and γ p,xy = ∂ 2 z p /∂x p ∂y p | xp=x p ,yp=y p . In Eq. (