Yuanrui Zhang

Hengyang Wu

Yixiang Chen

Frédéric Mallet

Embedding CCSL into Dynamic Logic: A Logical Approach for the Verification of CCSL Specifications

The Clock Constraint Specification Language (CCSL) is a clock-based specification language for capturing causal and chronometric constraints between events in Real-Time Embedded Systems (RTESs). Due to the limitations of the existing verification approaches, CCSL lacks a full verification support for 'unsafe CCSL specifications' and a unified proof framework. In this paper, we propose a novel verification approach based on theorem proving and SMT-checking. We firstly build a logic called CCSL Dynamic Logic (CDL), which extends the traditional dynamic logic with 'signals' and 'clock relations' as primitives, and with synchronous execution mechanism for modelling RTESs. Then we propose a sound and relatively complete proof system for CDL to provide the verification support. We show how CDL can be used to capture RTES and verify CCSL specifications by analyzing a simple case study.

Introduction

UML/MARTE [START_REF]OMG: UML profile for MARTE: Modeling and analysis of real-time embedded systems[END_REF] is an extension of UML dedicated to the modelling and analysis of Real-Time Embedded Systems (RTESs). Its time model relies on so-called clocks to identify control and observation points in the UML model. These clocks can be used to specify how the system behaves. The Clock Constraint Specification Language (CCSL) [START_REF] Mallet | Clock constraint specification language: specifying clock constraints with UML/MARTE[END_REF][START_REF] André | Syntax and Semantics of the Clock Constraint Specification Language (CCSL)[END_REF] is a formal declarative language defined in an annex of MARTE to specify the expected behaviour of the model. Given a system model (or a concrete implementation) and a CCSL specification, the question to answer is whether the system can only perform behaviors that are accepted by the CCSL specification [START_REF] Mallet | Correctness issues on MARTE/CCSL constraints[END_REF]. When a CCSL specification can be encoded as a finite transition system, it is called 'safe' [START_REF] Mallet | Safe CCSL specifications and marked graphs[END_REF], then the verification task mainly consists in making reachability analysis on the product of the system and the CCSL specification. Most recently SMT encoding of CCSL [START_REF] Zhang | Towards SMT-based LTL model checking of clock constraint specification language for real-time and embedded systems[END_REF] proved to be a promising way to verify unsafe CCSL specification, however, there is no proof environment available so far for reasoning on general specifications.

In this paper, we propose a novel approach for the verification of CCSL, which is based on the combination of theorem proving and SMT-checking. To capture both the system model and the CCSL specification, we choose dynamic logic [START_REF] Harel | Dynamic logic[END_REF], since it contains both dynamic program and static logic as its primitives. We propose a variation of dynamic logic, called 'CCSL Dynamic Logic' (CDL), which extends the traditional First-Order Dynamic Logic (FODL) [START_REF] Harel | First-Order Dynamic Logic[END_REF] with 'signal' and 'CCSL clock relations' as primitives in its syntax. CDL also supports synchronous events in order to capture synchronous system models [START_REF] Halbwachs | Synchronous programming of reactive systems[END_REF]. We propose a sound and relatively complete proof system for CDL in order to verify CDL formulas in a modular way. Our approach for verification of CCSL specifications can be illustrated in a verification framework given in Fig. 1. The verification task can be captured as a CDL formula of the form [p]ξ, where part of the CCSL specification, called 'clock relations' (will be introduced in Sect. 2), are expressed by a formula ξ, and the product of the system model and 'clock definitions' (the other part of the CCSL specification) can be captured by a program of CDL p. In CDL, a formula [p]ξ can be transformed into Quantifier-Free, Arithmetical First-Order Logic (QF-AFOL) formulas through a deduction procedure in the proof system of CDL. Then the validity of these formulas can be handled by an SMT-checking procedure in an efficient way [START_REF] Barrett | The SMT-LIB Standard: Version 2[END_REF], and according to which the verification result is obtained. With CDL, CCSL specifications can be verified in a unified proof framework, provided with strong tool support, e.g. Isabelle [START_REF] Nipkow | Isabelle/HOL -A Proof Assistant for Higher-Order Logic[END_REF] and Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development -Coq'Art: The Calculus of Inductive Constructions[END_REF].

The rest of this paper is organized as follows: Sect. 2 gives a general introduction to CCSL and FODL. Sect. 3 introduces the syntax and semantics of CDL. In Sect. 4, we propose the proof system for CDL. In Sect.5, we give a simple case study to show how CDL can express and verify CCSL verification problems. Sect. 6 introduces the related works, and Sect. 7 concludes this paper and discusses about future work.

Preliminaries of CCSL and FODL

We present the syntax and semantics of CCSL based on [START_REF] Mallet | Correctness issues on MARTE/CCSL constraints[END_REF][START_REF] Zhang | An SMT-based approach to the formal analysis of MARTE/CCSL[END_REF]. In CCSL, a logical clock actually models a sequence of occurrences of a signal in synchronous models [START_REF] Berry | The Esterel synchronous programming language: design, semantics, implementation[END_REF]. A logical clock c is defined as an infinite sequence of instants (c i ) i∈N + , where each c i can be 'tick' or 'idle', representing that the signal associated to c occurs or not at a discrete time i. N + is the set of natural numbers. Clock relations describe binary relationships between clocks. The syntax of clock relations is defined by:

Rel ::= c 1 ⊆ c 2 | c 1 #c 2 | c 1 ≺ c 2 | c 1 c 2 ,
where c 1 , c 2 are arbitrary clocks. We use C to denote a finite set of clocks. A schedule σ : N → P(C) is a finite or infinite sequence of clock ticks, N = N + ∪{0}.

It gives a global view of how each clock ticks at each instant. For any 'Intersection' defines the clock that ticks whenever both c 1 and c 2 tick; '(Strict) Sample' defines the clock that (strictly) samples c 1 based on c 2 ; 'Interruption' defines the clock that ticks as c 1 until c 2 ticks; 'Periodicity' defines the clock that ticks every n ticks of clock c ; 'Delay' defines the clock that ticks when c ticks but is delayed for n ticks of c . 'Infimum' ('Supremum') defines the slowest (fastest) clock that is faster (slower) than both c 1 and c 2 . e.g., Fig. 2 shows a possible schedule of clock definitions c c ∝ n and c c $n (when n = 2), which are used in the case study we give in Sect. 5.

i ∈ N + , σ(i) = {c | c ∈ C ∧ c i = tick}. σ(0) = ∅ indicates
σ = ∅{c 1 }{c 2 }{c 1 }{c 1 , c 2 }∅∅{c 1 , c 2 }{c 2 }∅{c 1 }{c 2 }∅.... X σ (c 1 , 1) = 1, X σ (c 1 , 2) = 1, X σ (c 1 , 3) = 2. X σ (c 2 , 1) = 0, X σ (c 2 , 2) = 1. b c 1 c 2 c 1 ≺ c 2 b c ′ c c c ′ ∝ n (n = 2) b c ′ c c c ′ $n (n = 2)
Table 1. Semantics of CCSL 1. σ ccsl c1 ⊆ c2 iff ∀i ∈ N + .c1 ∈ σ(n) → c2 ∈ σ(n) (Subclock) 2. σ ccsl c1#c2 iff ∀i ∈ N + .c1 / ∈ σ(i) ∨ c2 / ∈ σ(i) (Exclusion) 3. σ ccsl c1 ≺ c2 iff ∀i ∈ N + .(Xσ(c1, i) = 0 ∧ Xσ(c2, i) = 0) ∨ Xσ(c1, i) > Xσ(c2, i) (Precedence) 4. σ ccsl c1 c2 iff ∀i ∈ N + .Xσ(c1, i) ≥ Xσ(c2, i) (Causality) 5. σ ccsl c c1 + c2 iff ∀i ∈ N + .c ∈ σ(i) ↔ (c1 ∈ σ(i) ∨ c2 ∈ σ(i)) (Union) 6. σ ccsl c c1 * c2 iff ∀i ∈ N + .c ∈ σ(i) ↔ (c1 ∈ σ(i) ∧ c2 ∈ σ(i)) (Intersection) 7. σ ccsl c c1 c2 iff ∀i ∈ N + .c ∈ σ(i) ↔ (c2 ∈ σ(i)∧ (∃0 < j < i)(∀j ≤ k < i).c1 ∈ σ(j) ∧ c2 / ∈ σ(k)) (Strict Sample) 8. σ ccsl c c1 c2 iff ∀i ∈ N + .c ∈ σ(i) ↔ (c2 ∈ σ(i)∧ (∃0 < j ≤ i)(∀j ≤ k < i).c1 ∈ σ(j) ∧ c2 / ∈ σ(k)) (Sample) 9. σ ccsl c c1 c2 iff ∀i ∈ N + .c ∈ σ(i) ↔ (c1 ∈ σ(i) ∧ (∀0 < j ≤ i).c2 / ∈ σ(j)) (Interruption) 10. σ ccsl c c ∝ n iff ∀i ∈ N + .c ∈ σ(i) ↔ (c ∈ σ(i) ∧ ∃m ∈ N + .Xσ(c , i) = m • (n + 1)) (Periodicity) 11. σ ccsl c c $n iff ∀i ∈ N + .Xσ(c, i) = max(Xσ(c , i) -n, 0) (Delay) 12. σ ccsl c c1 ∧ c2 iff ∀i ∈ N + .Xσ(c, i) = max(Xσ(c1, i), Xσ(c2, i)) (Infimum) 13. σ ccsl c c1 ∨ c2 iff ∀i ∈ N + .Xσ(c, i) = min(Xσ(c1, i), Xσ(c2, i)) (Supremum)
A CCSL specification is a conjunction of clock relations and clock definitions, denoted as a triple SP ::= C, fi Cdf , Rel , where C is the set of clocks. fi Cdf is a set of clock definitions and Rel is a set of clock relations. σ ccsl C, fi Cdf , Rel is defined s.t. σ ccsl Rel and σ ccsl Cdf hold for all Rel ∈ Rel and Cdf ∈ fi Cdf . FODL is an extension of propositional dynamic logic with assignment x := e and testing P ? in its program model. The FODL we present here is based on [START_REF] Harel | Dynamic logic[END_REF]. The program of FODL is a regular program, defined as follows:

p ::= x := e | P ? | p; p | p ∪ p | p * ,
where e is an arithmetical expression. P ? means at current state, P is true. p; q means the program first executes p, and after p terminates, it executes q. p ∪ q means the program either executes p, or executes q, it is a non-deterministic choice. p * means the program executes p for a finite number of times. An FODL formula is defined as follows:

φ ::= tt | e ≤ e | [p]φ | ¬φ | φ ∧ φ | ∀x.φ,
where tt is the boolean true, ≤ represents the 'less than' relation in number theory. [p]φ is the dynamic formula, meaning after all executions of program p, formula φ holds.

The semantics of FODL is based on Kripke structure [START_REF] Harel | Dynamic logic[END_REF]. A Kripke structure is a pair (S, val) where S is a set of states, val is a function that interprets a logic to data structures on S.

In FODL, val interprets a regular program as a set of state pairs (s, s ) and interprets a formula as a set of states. Intuitively, each pair (s, s ) ∈ val(p) means that starting from state s, after execution of p, the program may terminate at state s . Each state s ∈ val([p]φ) means that for all pairs (s, s ) ∈ val(p), s satisfies φ. For a formal definition of the semantics of FODL, refer to [START_REF] Harel | Dynamic logic[END_REF].

The deductive system for FODL is sound and relatively complete. Except for the rule for atomic program 'x := e', all rules can be found in Table 3, 4 below, as a part of CDL proof system. Refer to [START_REF] Harel | Dynamic logic[END_REF] for more details.

Syntax and Semantics of CDL

CDL enriches the traditional FODL with a synchronous program model that contains 'signal' as a primitive, and 'clock relation' as an ingredient of logic formulas. We first give the syntax of the CDL program model and the CDL formula, and then define their semantics.

The Syntax of CDL

Syntax of Synchronous Event Programs. CCSL essentially describes the logical and chronometrical constraints between signals in synchronous models, where the time model is discrete and at each time, several signals can be triggered simultaneously. To capture CCSL constraints in dynamic logic, we need to introduce the synchronous execution mechanism in the regular program of FODL. Synchronous systems often involve infinite executions, thus to support it we also import 'infinite loop'. The program after enriched turns out to be an 'omega program', with the support of synchronous mechanism. We call it 'Synchronous Event Program' (SEP).

Definition 1 (Syntax of SEP). The syntax of SEP is based on the regular program of FODL, defined as follows:

p ::= ε | α | P ?α | p; p | p ∪ p | p * | p ω ,
where α is a combinational event, defined as:

α ::= | Cmb, Cmb ::= c | x := e | (Cmb|Cmb).
Arithmetical expression e, testing condition for signals and testing condition P are defined as follows: P ?α is a testing event, it means that if condition P is true, event α proceeds, otherwise the program causes a deadlock. In SEP, testing P ? must combines with an event α, because P ? does not consume time. P can be expressed with a QF-AFOL formula, where tt represents the boolean true, ≤ represents the 'less than' relation between two integers. Operator ; , ∪, * are defined just as in FODL [START_REF] Harel | Dynamic logic[END_REF]. ω represents the infinite loop. p ω means that program p executes for infinite number of times and never terminates. e.g., program f = 1?α 2 ; p * where p ::

= n = 0 ∧ f = 0?α 3 ∪ n > 0 ∧ f = 0?α 4 firstly executes α 2 if f = 1 holds, then it executes program p for finite number of times. In p, it either executes α 3 (if n = 0 ∧ f = 0 holds), or executes α 4 (if n > 0 ∧ f = 0 holds).
The precedence of operators are listed as follows from the highest to the lowest: ω, * , ;, ∪. We stipulate that ; is right-associative, ∪ is left-associative. e.g., program α 1 ∪p 1 ; p

2 ; p • 3 ∪P 1 ?α • 2 ∪P 2 ?α 3 means (((α 1 ∪p 1 ; (p 2 ; p • 3 ))∪P 1 ?α • 2 )∪ P 2 ?α 3 ).
As in synchronous models (e.g., Esterel [START_REF] Berry | The Esterel synchronous programming language: design, semantics, implementation[END_REF]), we do not allow two signals with the same name triggered at the same time. e.g. event (c|c). For simplification, we also do not allow two assignments with the same target variable executing simultaneously, e.g. event (x := 5|x := y + 1).

Syntax of CDL.

In CDL formula, we need to introduce a special kind of variable which is related to clock. These variables help record the 'information' of each clock at current time, just as the roles the schedule σ and the configuration X σ play in CCSL.

Definition 2 (Clock Related Variables). For each clock c ∈ C, we define two variables related to it: c n , c s . Variable c n is of type N, it records the number of times the clock has ticked at current time. Variable c s is of type {0, 1}, and it records the status of the clock (1 for present and 0 for absent) at current time.

Given a clock set C, we denote the set of variables related to C as V ar(C).

Definition 3 (Syntax of CDL Formula). The CDL formula φ is defined as:

φ ::= tt | E ≤ E | [p]ξ | [p]φ | ¬φ | φ ∧ φ | ∀x.φ where ξ ::= Rel | (Rel 1 , ..., Rel n ), E ::= x | c n | c s | n | E + E | E • E. E ≤ E
σ ccsl ∼ cr iff σ ccsl cr, (ii) σ ccsl (cr 1 , ..., cr n ) iff σ ccsl ∼ (∼ cr 1 , ..., ∼ cr n ), where cr, cr i ∈ {Rel i , ∼ Rel i }(1 ≤ i ≤ n).
[p]φ is the dynamic formula in FODL, meaning that after all executions of p, formula φ is satisfied.

We often call ξ or ∼ ξ "path formulas", denoted by π. Other arithmetic expressions, relations, and logic expressions, e.g., E -E, E/E, E = E, E < E, ff , p ∼ ξ, p φ, φ ∨ φ, φ → φ, ∃x.φ, etc, can be expressed using the formulas given above. e.g., p ∼ ξ can be expressed as

¬[p]ξ, E 1 -E 2 , E 1 /E 2 can be expressed as ∃x.(E 2 + x = E 1 ), ∃x.(x • E 2 = E 1 ) respectively.
In FODL, given a formula φ, a variable whose value changes with the execution of a program is called a 'dynamic variable' [START_REF] Harel | Dynamic logic[END_REF] of the formula φ. Here in CDL, for convenience sake, any clock-related variable c s , c n is defined as a dynamic variable. As we will see in Def. 5(ii), they can be seen 'changed' after the execution of any event at current time. Any general variable that appears on the left side of an assignment is defined as a dynamic variable as well. Variables which are not dynamic variables are called 'static variables'. e.g., the set of dynamic variables of formula z = 5 → [(c

1 |x := y+1); c 2 ]c 1 c 3 is {x, c n 1 , c s 1 , c n 2 , c s 2 , c n 3 , c s 3 }
, where c n 3 , c s 3 can be seen as 'changed' after the the set of static variables is {y, z}. Like in FODL, we say a variable x is 'bound' in φ iff: 1. x is in the scope of the effect of some quantifier ∀x, or 2. x is in the scope of the effect of some event α which has x on the left side of an assignment of the form x := e. A variable is not bound in φ is called 'free'. e.g., in formula φ ::= (x = 1 ∧ z = 2 ∧ ∃z.x = z) → [(x := z + 1|c|y := 1); x := y + 1]x > z, the first and second variable x is free, while the third one (in expression 'x > z') is bounded by the assignment 'x := y + 1'.

Given a formula φ, a substitution φ[E/x] in CDL replaces all the free occurrences of variable x with expression E (of the same type). Given a formula multi-set Γ , Γ [E/x] means to carry out the substitution φ[E/x] for each formula φ in Γ . Given two vectors (E 1 , ..., E n ), (x 1 , ..

., x n ), φ[E 1 , ..., E n /x 1 , ..., x n ] is the shorthand of φ[E 1 /x 1 ][E 2 /x 2 ]...[E n /x n ].
A substitution is admissible with respect to a formula φ if there are no variables x, y such that y is in E, and after the replacement φ[E/x], y is bound in φ. e.g., in the formula φ given above, φ[z + 1/z] = (x = 1 ∧ z + 1 = 2 ∧ ∃z.x = z) → [(x := (z + 1) + 1|c|y := 1); x := y + 1]x > z + 1 is admissible, while φ[x + 1/z] = (x = 1 ∧ x + 1 = 2 ∧ ∃z.x = z) → [(x := (x + 1) + 1|c|y := 1); x := y + 1]x > x + 1 is not admissible. Intuitively, in φ[x + 1/z], it is about to prove x > x + 1 which is generally not true. In the rest of paper, unless we specially point out, all substitutions we discuss are admissible.

The Semantics of CDL

In the Kripke structure (introduced in Sect.2) of CDL, val interprets a program as a set of traces on S and a logic formula as a set of states. A trace tr is a finite or infinite sequence of states. Given a finite trace tr 1 = s 1 s 2 ...s n and a (possibly infinite) trace tr 2 = u 1 u 2 ...u n ..., we define:

tr 1 •tr 2 ::= s 1 s 2 ...s n u 2 u 3 ... if s n = u 1 .
Given any tr 1 , tr 2 , we define tr 1 • tr 2 ::=

ß tr 1 • tr 2 , if tr 1 is finite tr 1 , otherwise . Given two sets of traces S 1 , S 2 , S 1 • S 2 is defined as {tr 1 • tr 2 | tr 1 ∈ S 1 , tr 2 ∈ S 2 }. Let tr(i)
denotes the i th element of trace tr, i ≥ 0; tr b denotes the first element of trace tr, tr b = tr(0). Let tr e denotes the last element of trace tr, provided that tr is a finite trace.

In CDL, we assume an interpretation which interprets arithmetical operators '+, -, •, /' and relation '≤' as their usual meanings in the traditional number theory, and interprets relations '⊆, ≺, , #' as their corresponding clock relations in CCSL. Next we first define the concept of 'state' in CDL.

Definition 4 (State and Evaluation in CDL).

A state s in CDL is a total function defined as follows: (i) s maps each variable c n in V ar(C) to a value in domain N. (ii) s maps each variable c s in V ar(C) to a value in domain {0, 1}. (iii) s maps each variable x in V ar to a value in domain Z.

Given an expression E and a state s, an evaluation Eval s (E) is defined as:

(i) If E = a, where a ∈ {x, c n , c s }, then Eval s (a) ::= s(a). (ii) If E = n, then Eval s (n) ::= n. (iii) If E = f (E 1 , E 2 ), where f ∈ {+, •}, then Eval s (E) ::= f (Eval s (E 1 ), Eval s (E 2 )).
e.g., given a state s ::

= {x → 9, c n → 2, (c ) s → 0, ...}, there is Eval s (2) = 2, Eval s (x) = 9, Eval s (x + c n ) = Eval s (x) + Eval s (c n ) = 11.
Semantics of SEP. Different from traditional FODL, the semantics of SEP is based on traces, since our CDL contains a path formula π which is satisfied by a program trace.

Definition 5 (Semantics of SEP). Given a Kripke structure (S, val), for any SEP p, let C be a finite set of clocks, the semantics of SEP is given as follows:

(i) val(ε) := S, S is the set of all traces of length 1. Note that ε defines a set of traces of length 1, so val(p; ε) = val(ε; p) = val(p), which means that ε can be taken as a unit element of operator ;. Event α defines a transition from a state s to a state s . In s , for each clock c in α, the variable c n that records the number of ticks is added by 1 and the variable c s is set to 1, indicating at current time, clock c is emitted. For each clock d not in α, its variable d n in s is kept the same while d s is set to 0. For any assignment x := e in α, the value of x in s is set to the value of expression e in state s, while other variables in both s and s are kept the same. Traces satisfying P ?α are exactly those traces satisfying p adding that their beginning states must satisfy P . e.g., let α = (c|x := x + 1),

P = x > 1, C = {c, c }, V ar = {x, y}, if s = {x → 0, y → 0, c n → 0, c s → 0, c n → 0, c s → 0}, s = {x → 1, y → 0, c n → 1, c s → 1, c n → 0, c s → 0}, then trace ss ∈ val(α). If u = {x → 2, y → 0, c n → 1, c s → 1, c n → 0, c s → 1}, u = {x → 3, y → 0, c n → 2, c s → 1, c n → 0, c s → 0}, then trace uu ∈ val(P ?α).
The semantics of p; q, p ∪ q, p * are directly inherited from the traditional FODL [START_REF] Harel | Dynamic logic[END_REF]. The traces of program p ω consists of all infinite traces of the form tr 1 • tr 2 ... where each tr i ∈ val(p) is finite (i ∈ N + ), or of the form

tr 1 • tr 2 • ... • tr n , where n ≥ 1, tr 1 , ..., tr n-1 ∈ val(p) is finite, but tr n ∈ val(p) is infinite. e.g., suppose val(p) = {s 1 s 2 }, val(q) = {u 1 u 2 , t 1 t 2 } where s 2 = u 1 , s 2 = t 1 , then val(p; q) = {s 1 s 2 u 2 }, val(p ∪ q) = {s 1 s 2 , u 1 u 2 , t 1 t 2 }, val(p * ) = {ε, s 1 s 2 , s 1 s 2 s 1 s 2 , ..., s 1 s 2 s 1 s 2 ...s 1 s 2 2n , ...} (n ≥ 1), val(p ω ) = val(p * )∪{s 1 s 2 s 1 s 2 ...s 1 s 2 ... ∞ }.
Semantics of CDL. For each trace tr, we can actually build a corresponding schedule σ tr s.t. for all clock c ∈ C and i ∈ N + , there is:

1. tr(i)(c n ) = X σ (c, i). 2. tr(i)(c s ) = 1 iff c ∈ σ tr (i).
In this way, we can actually define tr ccsl X given a clock relation or definition X: tr ccsl X iff σ tr ccsl X. Note that we do not require any relationships between tr(0) and σ(0). e.g., consider the trace ss discussed above, we have a schedule σ ss defined as:

σ ss ::= ∅{c}. So X σ ss (c, 0) = X σ ss (c , 0) = 0, X σ ss (c, 1) = 1, X σ ss (c , 1) = 0.
Definition 6 (Semantics of CDL Formula). Given a Kripke structure (S, val), the semantics of CDL formula is given as follows:

v 0 ∈ Z, s ∈ val(φ[v 0 /x])}.
The semantics of CDL formula is based on states. In (iii), a trace satisfying a clock relation is from the second state of the trace due to the definition of ccsl in Sect. 2. So state s itself is unrelated to ξ. (iv)-(vii) are similar to the definition in FODL [START_REF] Harel | Dynamic logic[END_REF], except that the semantics of SEP is based on traces. (iv) requires the trace must be finite, indicating that it only matters whether φ holds on those states on which program p terminates.

The first two figures in Fig. 3 give an illustration of [p]φ and [p]ξ, where the 'snake arrow' indicates an execution path (could be infinite) of program. Some states are tagged with a formula aside that they satisfy. States and paths are colored red to stress that they satisfy the corresponding formulas (φ, ξ).

At last we define the satisfaction relation of the CDL logic. Given a state s and any CDL formula φ, the satisfaction relation s |= cdl φ is defined as: s |= cdl φ iff s ∈ val(φ). If for all state s, s |= cdl φ holds, then we say φ is valid, denoted as |= cdl φ.

Proof System of CDL

In this section we propose a proof system, which forms the foundation of the verification of CDL. The proof system provides a modular way of transforming a CDL formula into a QF-AFOL formula. Our proof system is based on that of FODL, which is only for regular program model [START_REF] Harel | Dynamic logic[END_REF].

ϕ φ to represent Γ ⇒ ϕ ⇒ ∆ Γ ⇒ φ ⇒ ∆ if
Γ, ∆ can be neglected. We call Γ, ∆ the context of formula φ in sequent Γ ⇒ φ, ∆ or Γ, φ ⇒ ∆.

Proof Rules for CDL

The proof rules of CDL we present are divided into three categories: rules for path formulas π (in Table 2), rules for non-path formulas (in Table 3) and rules of First-Order Logic (FOL) (in Table 4).

In Table 2, rule (π) is for a single event, where we set α = (c|x := e) as an example of combinational events.

Γ [V /V ], c n = (c n ) + 1, c s = 1, x = e[V /V ], (d n 1 , ..., d n n ) = ((d n 1 ) , ..., (d n n ) ), (d s 1 , ..., d s n ) = (0, ..., 0) n ⇒ (ξ) ⇒ ∆[V /V ] Γ ⇒ [α]ξ ⇒ ∆ (π)
where α = (c|x := e), {d1, ..., dn} = C -C(α), V = V(α),

V is the set of new variables (w.r.t. Γ, [α]ξ, ∆, ) corresponding to V.

P → [α]A [P ?α]A (P ?)

where A ∈ {ξ, φ}

tt [ε]ξ (πε) [p * ]ξ [p ω ]ξ (π[ω]) [p]ξ ∧ [p][q]ξ [p; q]ξ (π[;]) [p]ξ ∧ [q]ξ [p ∪ q]ξ (π[∪]) [p; p * ]ξ [p * ]ξ (π[ * ]u) [p * ][p]ξ [p * ]ξ (π[ * ]i)
The rule says that for any state s, the conclusion holds at state s, iff there exists a state s with ss ∈ val((c|x := e)), s.t. the premise holds at s . The vector equation (x 1 , ..., x n ) = (e 1 , ..., e n ) is the shorthand of x 1 = e 1 , ..., x n = e n . d 1 , ..., d n are all clocks not appeared in α. Given a CDL formula φ (or an SEP p), let C(φ) (C(p)) returns all clocks appeared in φ (p), V(φ) (V(p)) returns all dynamic variables appeared in φ (p). V is the set of new variables corresponding to V , for each variable x ∈ V , there is a new variable x with respect to Γ, [α]ξ, ∆ corresponding to it. Function (ξ) maps each relations to an AFOL formula which should hold at state s . It is defined as follows: for any c 1 , c 2 , (i

) (c 1 ⊆ c 2 ) ::= c s 1 = 1 → c s 2 = 1. (ii) (c 1 #c 2 ) ::= c s 1 = 0 ∨ c s 2 = 0. (iii) (c 1 ≺ c 2 ) ::= (c n 1 = 0 ∧ c n 2 = 0) ∨ (c n 1 > c n 2 ). (iv) (c 1 c 2 ) ::= c n 1 ≥ c n 2 . (v) ( (Rel 1 , ..., Rel n ) ::= 1≤i≤n (Rel i ).
(P ?) is a rule for both path-formulas and non-path formulas. Rule (P ?) says that the conclusion at a state is true, iff if P is true, then [α]A is true. In rule (πε), tr ccsl ξ always holds for trace tr of length 1. Rule (πω) is based on two ) means every trace of p; q satisfies ξ iff every trace of p satisfies ξ, and after p every trace of q satisfies ξ. (π[∪]) says every trace of p∪q satisfies ξ iff every trace of p and q satisfies ξ. Rule (π[ * ]u) unwinds the star operator * . It is due to the fact that every trace (whose length ≥ 2) of p * are the trace of p; p * . (π[ * ]i) states that ξ holds along all paths of any times of repetitions of p, iff after any times of repetitions p , ξ holds along all paths of p. Fig. 3 gives a graphical illustration of rule (π), (π[; ]), (π[∪]), (π[ * ]u), (π[ * ]i).

α [p; q]ξ [q]ξ p ξ q p; q [p]φ [p]ξ (π) (π[; ]) p q ξ p ∪ q [p * ]ξ [p; p * ]ξ p p * ξ p * [p * ]ξ [p]ξ p * ξ p p * [p; q]φ [q]φ φ p q p; q (π[∪]) (π[ * ]u) (π[ * ]i) ([; ]) φ φ p q p ∪ q φ [p * ]φ [p; p * ]φ φ p p * p * φ φ φ φ → [p]φ φ p * p * p p * φ(n) φ(n -1) ∃x ≥ 0.φ ∃x ≤ 0.φ(x) p * p * p p * ([∪]) ([ * ]u) ([ * ]ind) ( *
All non-path formula rules in Table 3 except for (φ), (ε), (ω) are based on the corresponding structure rules of FODL in [START_REF] Harel | Dynamic logic[END_REF]. (φ) is similar to (π), except that φ is kept unchanged in the premise. (ε) is obvious because the traces of ε all have length 1. Rule ([; ]) describes that φ holds after p; q iff [q]φ holds after p. Rule ([∪]) says φ holds after p ∪ q iff φ holds after p, and also holds after q. ([ * ]u) means that φ holds after any times of repetitions of p, iff φ holds at current state, and φ holds after p; p * .

([]gen), ( gen) strengthen the conclusions by extending the proposition φ → ϕ into dynamic situations. ([]gen) (( gen)) expresses that if φ → ϕ holds under all context of Γ, ∆, after any (some) executions of p, φ implies ϕ. ([ * ]ind) is the mathematical induction by the number of repetitions of program p: to prove φ holds after any repetitions (including 0), we need to prove that under any context of Γ, ∆, if φ holds, then it also holds after p. ( * con) is from the Harel's convergence rule in [START_REF] Harel | Dynamic logic[END_REF] where integer x indicates the existing number of

Γ ⇒ φ, ∆ Γ, ϕ ⇒ ∆ Γ, φ → ϕ ⇒ ∆ (→l) Γ ⇒ φ[tn/x], ∆ Γ ⇒ ∃x.φ, ∆ (∃r) Γ, ∃x.φ, φ[x /x] ⇒ ∆ Γ, ∃x.φ ⇒ ∆ (∃l)
where Γ, ∆ are multi-sets of QF-AFOL formulas. x is a new variable w.r.t. Γ, φ, ∆, φ[tn/x] is admissible.

A Case Study

In this section, we illustrate how our proposed CDL can be used to capture RTES models and verify CCSL specifications, by analyzing a simple RTESthe Digital Filter (DF) system. The DF system we analyze here is based on [START_REF] André | Specification and verification of time requirements with CCSL and Esterel[END_REF].

As Fig. 4 shows, the DF is used in a video system, it reads image pixels from a memory, filters them and sends the result out to a video device. The explicit structure of DF is shown in the right figure of Fig. 4. The DF consists of two modules: a Feeder and a Filter. They interact with each other and with their environment through ports p, r and o. Ports are the only way for different modules to communicate in synchronous models. They can be modelled as signal c p , c r and c o in SEP. The behaviour of the DF, as a whole system of two modules, is as follows: the Filter sends a 'Ready' message to the Feeder through port r, to tell it 'I am ready for the pixels'. The Feeder receives this message and the next time it begins feeding pixels towards the Filter, one pixel per unit of time. SP 1 expresses the property that 'the result can be obtained only after the "Ready" message is sent', i.e., clock c r ticks strictly before c o . SP 2 says that 'only after the last result is computed, the new pixels can be received'. SP 2 contains two clock definitions. c p , c p are generated clocks not appeared in program DF . Two specifications are indicated by red and blue arrows in Fig. 5 respectively. The ticks of clock c p , c p are indicated by green and yellow circles respectively.

For SP 1 , the verification problem can be captured by a CDL formula:

I → [DF ]c r ≺ c o ,

Fig. 1 .

 1 Fig. 1. Verification framework of CDL

Fig. 2 .

 2 Fig. 2. A possible schedule for selected clock constraints

e

  ::= x | n | e + e | e -e | n • e | e/n, P ::= tt | e ≤ e | ¬P | P ∧ P. ε represent an 'empty program', it does nothing nor consumes time. A combinational event α consumes a unit of time, it consists of an 'idle event' , or several signals or assignments that occur simultaneously. An idle event does nothing but waits for a unit of time. Several signals and assignments can be composed by operator '|'. A signal 3 c in an SEP represents that its corresponding clock (with the same name c) ticks at current time. Since CCSL constraints only captures the logical relationships between signals which are not related to the value of signals, we only consider 'pure signals' (signals without values) in SEP. e is a Presburger arithmetic expression. In e, n ∈ Z is an integer number, +, -, •, / are the addition, subtraction, multiplication and division signs respectively.

  (ii) val(α) := {ss | s, s ∈ S; for each clock c ∈ α, s (c s ) = 1 ∧ s (c n ) = s(c n ) + 1; for other clock d ∈ C, s (d n ) = s(d n ) ∧ s (d s ) = 0; for each x := e in α, s (x) = Eval s (e); for other x ∈ V ar, s (x) = s(x)}. (iii) val(P ?α) ::= {ss | s ∈ val(P ), ss ∈ val(α)}. (iv) val(p; q) ::= val(p) • val(q). (v) val(p ∪ q) ::= val(p) ∪ val(q). (vi) val(p * ) ::= n≥0 val n (p), where val n (p) = val(p) • ... • val(p) n , val 0 (p) = S. (vii) val(p ω ) ::= val(p) • val(p) • ... ∞ .

  (i) val(tt) ::= S. (ii) val(E ≤ E ) ::= {s | Eval s (E) ≤ Eval s (E )}. (iii) val([p]ξ) ::= {s | for all tr s.t. s = tr b and tr ∈ val(p), tr ccsl ξ}. (iv) val([p]φ) ::= {s | for all finite tr ∈ val(p) s.t. tr b = s, tr e ∈ val(φ)}. (v) val(¬φ) ::= {s | s / ∈ val(φ)}. (vi) val(φ ∧ ϕ) ::= val(φ) ∩ val(ϕ). (vii) val(∀x.φ) ::= {s | for any

Fig. 4 .Fig. 5 .

 45 Fig. 4. The Digital Filter system

  the beginning of the sequence where no clock ticks. X σ : C × N + → N keeps track of the number of ticks for each clock. X σ (c, i) = |{j | j ∈ N + , j ≤ i, c ∈ σ(j)}| is called a configuration of clock c at time i. The semantics of clock relations is defined as items 1-4 in Table1. 'Subclock' says that c 1 can only tick if c 2 ticks; 'Exclusion' means that c 1 , c 2 can not tick at the same instant; 'Precedence' means that c 1 always ticks faster than c 2 ; 'Causality' expresses that c 1 ticks not slower than c 2 .For example, the leftmost figure of Fig.2shows a possible schedule σ for clock relation c 1 ≺ c 2 , where clock b = tick tick tick tick tick tick tick tick tick tick tick tick ..., c 1 = tick idle tick tick idle idle tick idle idle tick idle idle ..., c 2 = idle tick idle tick idle idle tick tick idle idle tick idle ....

b is a based clock representing the minimal granularity of time. Schedule

  is an atomic AFOL formula. E is an integer arithmetic expression. Different from e, it also includes clock-related variable c n , c s , and multiplication between variables. x ∈ V ar. [p]ξ is a dynamic formula, where p is an SEP. [p]ξ is the dynamic formula special in CDL, it means that all execution paths of program p satisfies ξ. (Rel 1 , ..., Rel n ) represents the conjunction of clock relations Rel 1 , ..., Rel n , we define σ ccsl (Rel 1 , ..., Rel n ) iff σ ccsl Rel 1 ,...,σ ccsl Rel n . In order to express the negation of [p]ξ in CDL, we also import the negation ∼ and the disjunction of clock relations: (i)

Table 2 .

 2 Rules for path formulas

Table 1 .

 1 With them not hard to see the premise and conclusion of rule (πω) are logical equivalent. With rule (πω) we can reduce the proof case of [p ω ]ξ to the proof case of [p

con) Fig. 3. Graphical illustrations of [p]φ, [p]ξ and some proof rules facts about clock relation ξ and SEP traces: (i) For any infinite trace tr ∈ val(p ω ) and any state s in tr, there exists a finite trace tr ∈ val(p * ) that contains s. (ii) For any relation ξ and trace tr, tr ccsl ξ iff tr(i) |= cdl (ξ) for any i ∈ N + . These two facts can be easily obtained according to Def. 5 and the definition of ccsl in * ]ξ. (π[; ]), (π[∪]), (π[ * ]u), (π[ * ]i) are structure rules for path formulas. (π[; ]

In SEP, for convenience, we use the same name 'c' to represent the signal corresponding to clock c, which should not cause any ambiguities. Sometimes we also say a signal c in p 'a clock c in p'.

A sequent[START_REF] Gentzen | Untersuchungen über das logische Schließen[END_REF] is defined as follows: Γ ⇒ ∆ ::= φ∈Γ φ → ϕ∈∆ ϕ, where Γ, ∆ are two finite multi-sets of logic formulas. It means that every formula in Γ holds can conclude that at least one of formulas in ∆ holds. The conditions when either (both) Σ or (and) ∆ is (are) empty set(s) is (are) expressed as follows: 1.• ⇒ ∆ ::= tt → ϕ∈∆ ϕ, 2. Γ ⇒ • ::= φ∈Γ φ → ff , 3. • ⇒ • ::= tt → ff , where we use • to indicate Γ or ∆ is empty. A rule in sequent calculus is of the form: Γ 1 ⇒ ∆ 1 ... Γ n ⇒ ∆ n Γ ⇒ ∆ , which means that if Γ 1 ⇒ ∆ 1 , ..., Γ n ⇒ ∆ n are all valid, then Γ ⇒ ∆ is valid. Each Γ i ⇒ ∆ i in the upper part is called a 'premise', while Γ ⇒ ∆ in the lower part is called 'conclusion'. We use Γ ⇒ ϕ ⇒ ∆ Γ ⇒ φ ⇒ ∆ to represent a pair of sequent rules: Γ, ϕ ⇒ ∆ Γ, φ ⇒ ∆ and Γ ⇒ ϕ, ∆ Γ ⇒ φ, ∆, i.e., φ, ϕ can be on both side of the sequent. Sometimes we write

Table 4. Rules of first order logic|=cdl φ∈Γ φ → ϕ∈∆ ϕ Γ ⇒ ∆ (o) Γ, φ ⇒ φ, ∆ (ax) Γ ⇒ φ, ∆ Γ, φ ⇒ ∆ Γ ⇒ ∆ (cut) Γ, ¬φ ⇒ ∆ Γ ⇒ φ, ∆(¬r)Γ ⇒ ¬φ, ∆ Γ, φ ⇒ ∆ (¬l) Γ ⇒ φ, ∆ Γ ⇒ ϕ, ∆ Γ ⇒ φ ∧ ϕ, ∆ (∧r) Γ, φ, ϕ ⇒ ∆ Γ, φ ∧ ϕ ⇒ ∆ (∧l) Γ ⇒ φ[x /x], ∆ Γ ⇒ ∀x.φ, ∆(∀r)Γ, ∀x.φ, φ[tn/x] ⇒ ∆ Γ, ∀x.φ ⇒ ∆ (∀l) Γ ⇒ φ, ϕ, ∆ Γ ⇒ φ ∨ ϕ, ∆(∨r)Γ, φ ⇒ ∆ Γ, ϕ ⇒ ∆ Γ, φ ∨ ϕ ⇒ ∆ (∨l) Γ, φ ⇒ ϕ, ∆ Γ ⇒ φ → ϕ, ∆(→r)

This work was partly funded by the French Government, through program #ANR-11-LABX-0031-01 allowing new clocks to be defined using different clock expressions.

repetitions of p. ([ * ]i) and ( * i) are rules for eliminating the star operator * in practical verification. They can be derived by ([ * ]ind), ( * con) with generalisation ([]gen), ( gen) (see [START_REF] Harel | Dynamic logic[END_REF][START_REF] Platzer | Logical analysis of hybrid systems -proving theorems for complex dynamics[END_REF]). ϕ is the loop invariant of p. ([ * ]i) says that to prove φ holds after any repetitions of p, we need to prove that there exists an invariant ϕ such that: (i) ϕ holds at the beginning. (ii) Under any context of Γ, ∆ if ϕ holds, then ϕ holds after p as well. (iii) Under any context of Γ, ∆, ϕ implies φ. Fig. 3 gives a graphical illustration of rule ([∪]), ([; ]), ([ * ]u), ([ * ]ind), ( * con).

Table 3. Rules for non-path formulas

where α = (c|x := e), {d1, ..., dn}

Other FOL rules are listed in Table 4. As indicated in Sect.1, after a QF-AFOL formula is obtained we can adopt SMT-checking procedure to check the validation of it. Since the SMT-checking procedure is independent from the CDL proof system, we propose an 'oracle' rule (o) in our proof system to indicate the termination of the proof. We assume that the validity of this QF-AFOL formula can be SMT-checked in a 'black box', through this oracle rule. Other rules comes from the traditional FOL and we omit the details of them. Now we define the deduction relation of CDL. For any CDL formula φ and a formula multi-sets Φ, Φ cdl φ iff the sequent Φ ⇒ φ can be derived according to rules in Table 2, 3, 4. If Φ is empty, we also write cdl φ. As a variation of dynamic logic, the soundness and relative completeness of proof system cdl can be analyzed in a similar way as those of FODL in [START_REF] Harel | Dynamic logic[END_REF][START_REF] Harel | First-Order Dynamic Logic[END_REF]. For the soundness, above we have explained the intuitive meaning of each rule and their relations to the corresponding rules in FODL. For the relative completeness, intuitively, to prove it we show that each formula of form [p]ξ can be transformed into an AFOL formula by applying the rules of the CDL proof system, which is similar for the formula [p]φ in FODL. Due to space limit, we omit the complete proof.

where I = c∈{cr,cp,co} (c n = 0 ∧ c s = 0) represents the initial environment of the program. The deduction procedure of this formula is illustrated in Fig. 6. Starting from the root node (node 1 ), the procedure answers 'yes' iff every leave node of the proof tree returns a valid QF-AFOL formula, which can be checked through an SMT-checking procedure. Inference rules with a 'double line' indicate that there are more than one deductions between the premises and the conclusion, and rules they are based on are listed on the right side. e.g., from node 2 , by applying rule (π[; ]), (∧r) in sequence, we obtain two premises: node 3 and 4 . The meanings of other symbols and contexts are explained in the table of Fig. 6. Note that variables of the form like '(c n r ) 1 ' is a new general variable (corresponding to c n r ), not a clock-related variable.

ϕ2 ⇒ [P1?α3]ϕ2

(P ?,→r,φ)

ϕ2 ⇒ [P2?α4]ϕ2

ϕ1

(P ?,→r,φ)

, where x = 1, 2, 3, 4. v x , x1, x2, y1, z1, z2, u1 are new variables corresponding to their counterparts in substitutions.

Fig. 6. The deduction procedure of I → [DF ]cr ≺ co Due to the limit of space, we omit the details of the branch from node 7 . At node 5 , 8 , we apply rule ([ * ]i) to eliminate the loop operator * . Here we need to manually decide the loop invariants ϕ 1 , ϕ 2 . The selecting of a suitable invariants is according to the loop body (here p 1 , p 2 ) and the formulas we want to verify after the loop program (here [p 1 ]c 1 ≺ c 2 , ϕ 1 ). e.g., in ϕ 1 , we have to guarantee 'c n r > c n o ' always holds during each execution of p 1 , because if not so, [p 1 ]c 1 ≺ c 2 would not hold for some state during the execution of p * 1 . 'C 1 ∨ C 2 ' is to make sure that n, f can only be 'reasonable values' during the execution of p 1 . At last, easy to see that each leave node is a valid QF-AFOL formula. e.g., at node 9 , clearly from Γ 4 , we have C 2 holds. In Γ 4 , since z 1 = 0 ∧ z 2 = 0(in

For SP 2 , just like in previous approaches [START_REF] Mallet | Correctness issues on MARTE/CCSL constraints[END_REF][START_REF] Suryadevara | Verifying MARTE/CCSL mode behaviors using UPPAAL[END_REF][START_REF] Zhang | Timed automata semantics of spatio-temporal consistency language STeC[END_REF], we firstly make the product of the system model DF and the clock definitions c p c p $1, c p c p ∝ 3. As indicated in Fig. 1, this product can then be captured by an SEP program. A similar verification procedure as above can be carried out.

Related Work

Previous approaches [START_REF] Suryadevara | Verifying MARTE/CCSL mode behaviors using UPPAAL[END_REF][START_REF] Zhang | Timed automata semantics of spatio-temporal consistency language STeC[END_REF] for the verification of CCSL specifications are mainly based on model checking, where the reachability analysis is made for the product of the system model and the CCSL specification. When the CCSL specification is unsafe, a bound needs to be set to avoid the enumeration of infinite number of states. Our approach is based on theorem proving and SMT-checking, which provides a unified framework under which both safe and unsafe CCSL specifications can be analyzed.

Another subject of analysis for CCSL is to find a schedule of a given CCSL specification [START_REF] André | Syntax and Semantics of the Clock Constraint Specification Language (CCSL)[END_REF][START_REF] Zhang | An SMT-based approach to the formal analysis of MARTE/CCSL[END_REF][START_REF] Zhang | Periodic scheduling for MARTE/CCSL: Theory and practice[END_REF], where no system models were involved. The earliest approach [START_REF] André | Syntax and Semantics of the Clock Constraint Specification Language (CCSL)[END_REF] combined BDD-based boolean solving and the rewriting on clock expressions, while the method in [START_REF] Zhang | Periodic scheduling for MARTE/CCSL: Theory and practice[END_REF] was based on the rewriting logic in Maude. In [START_REF] Zhang | An SMT-based approach to the formal analysis of MARTE/CCSL[END_REF], the schedule was found by solving an UFLIA formula that encodes the CCSL specification through an SMT-checking procedure. Comparing with [START_REF] Zhang | An SMT-based approach to the formal analysis of MARTE/CCSL[END_REF], we propose a proof system to transform the CDL formula into QF-AFOL formulas, which are more efficient for an SMT-checking procedure to solve.

CDL is largely based on the traditional FODL [START_REF] Harel | First-Order Dynamic Logic[END_REF] and its rules (π[; ]), (π[∪]), (π[ * ]u), (π[ * ]i) are inspired from the Differential Dynamic Temporal Logic (DDTL), a dynamic logic for verification of hybrid systems [START_REF] Platzer | Logical analysis of hybrid systems -proving theorems for complex dynamics[END_REF][START_REF] Platzer | A temporal dynamic logic for verifying hybrid system invariants[END_REF]. In DDTL, the program supports a continuous time model with differential equations embedded into it. Our SEP supports a discrete time model with synchronous mechanism which we think would be more friendly for modelling RTESs.

Conclusion and Future Work

In this paper, we propose a logical approach for verification of CCSL specifications. We build a variation of dynamic logic called CDL to capture the verification problem, and a proof system to provide the verification support. We give a case study to illustrate how CDL can be used for verifying CCSL specifications.

Unlike traditional synchronous programming languages, SEP only supports sequential models. We shall present a concurrent extension in a future work by adding a ' ' operator. We also consider mechanizing CDL with the popular theorem prover Coq in order to see more practical potentials for this method.


A Labibian
email: a.labibian@gmail.com

S H Pourtakdoust

A Alikhani
email: aalikhani@ari.ac.ir

H Fourati
email: hassen.fourati@gipsa-lab.grenoble-inp.fr

Development of a radiation based heat model for satellite attitude determination

Keywords:

This paper is focused on the development and verification of a heat attitude model (HAM) for satellite attitude determination. Within this context, the Sun and the Earth are considered as the main external sources of radiation that could effect the satellite surface temperature changes. Assuming that the satellite orbital position (navigational data) is known, the proposed HAM provides the satellite surface temperature with acceptable accuracy and also relates the net heat flux (NHF) of three orthogonal satellite surfaces to its attitude via the inertial to satellite transformation matrix. The proposed HAM simulation results are verified through comparison with commercial thermal analysis tools. The proposed HAM has been successfully utilized in some researches for attitude estimation, and further studies for practical implementations are still ongoing.

Introduction

Attitude determination (AD) of satellites has been an important issue since early developments of the first space mission, and it is still the subject of many scientific and applied researches. Proper knowledge of attitude is required for pointing capabilities via the satellite Attitude Control System (ACS). The history of research in AD goes back to 1965, when G. Wahba initially proposed her method of attitude estimation [START_REF] Wahba | A least-square estimation of satellite attitude[END_REF]. In general, when the AD of satellites is studied, the following three main parts need to be addressed: attitude representation, attitude sensors and estimation algorithms. There exists many researches on each of the above mentioned parts, under the so called AD problem. For the attitude representation, Euler angles, quaternion and Gibbs Vector are widely used. Conventional attitude sensors include magnetometer, gyroscope, Sun sensor, Earth sensor, as well as star tracker that can either be used separately [START_REF] Psiaki | Global magnetometer-based spacecraft attitude and rate estimation[END_REF][START_REF] Chen | High precision attitude estimation algorithm using three star trackers[END_REF]o r in combination [START_REF] Psiaki | Autonomous low-earth-orbit determination from magnetometer and sun sensor data[END_REF][START_REF] Son-Goo | Kalman filtering for relative spacecraft attitude and position estimation[END_REF][START_REF] Springmann | The attitude determination system of the RAX satellite[END_REF]. Sensors choice is usually dictated by mission statements and/or is dependent on the desired autonomy and accuracy. Estimation algorithms are also widespread and still subject of many current researches. Extended Kalman Filter (EKF) is known as one of the most popular techniques for attitude estimation in many applications. EKF uses nonlinear attitude equations in a forward estimation process and updates the estimated states in a recursive way [START_REF] Farrell | Attitude determination by Kalman filter[END_REF][START_REF] Murrel | Precision attitude determination for multimission spacecraft[END_REF][START_REF] Lefferts | Kalman filtering for spacecraft attitude estimation[END_REF][START_REF] Markley | Attitude error representation for Kalman filtering[END_REF]. Sometimes, some constraints are a critical issue, then the structure of EKF can be modified in order to satisfy the constraints. Multiplicative EKF is developed in that case, with a fairly acceptable performance [START_REF] Shuster | Constraint in attitude estimation, part I: constrained estimation[END_REF][START_REF] Shuster | Constraint in attitude estimation, part II: unscented estimation[END_REF]. Unscented Kalman Filter (UKF) and particle filter (PF) are among other alternatives for EKF. Although UKF increases the computation burden, it is suitable for applications involving complex measurement equations, and its performance is better than those of EKF [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF][START_REF] Wan | The unscented Kalman filter, in: Kalman Filtering and Neural Networks[END_REF][START_REF] Norgaard | New development in state estimation for nonlinear systems[END_REF][START_REF] Crassidis | Unscented filtering for spacecraft attitude estimation[END_REF][START_REF] Kiani | Concurrent orbit and attitude estimation using minimum sigma point unscented Kalman filter[END_REF]. PF has emerged as a new filter structure that is based on sequential Monte Carlo Simulation (MCS) in which distributions are approximated by weighted particles [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF][START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Cheng | Particle filtering for sequential spacecraft attitude estimation[END_REF][START_REF] Oshman | Estimating attitude from vector observations using a genetic algorithm-embedded quaternion particle filter[END_REF].

During last years, the thermal modelling and testing tools and techniques used in satellite programs have become a major area of research. Multiple thermal nodes are considered and numerical methods are used to solve the resulting thermal coupled nonlinear partial differential equations [START_REF] Gilmor | Spacecraft Thermal Control Handbook[END_REF][START_REF] Karam | Satellite Thermal Control for System Engineers[END_REF]. Since the computational cost of these methods is high, there are analytical schemes that are less costly to use depending on the application [START_REF] Tsai | Overview of satellite thermal analytic model[END_REF][START_REF] Gaite | Nonlinear analysis of a simple model of temperature evolution in a satellite[END_REF][START_REF] Gaite | Nonlinear analysis of spacecraft thermal models[END_REF]. Further, analytical single node thermal models can also produce acceptable results [START_REF] Anh | Thermal radiation analysis for small satellites with single-node model using techniques of equivalent linearization[END_REF]. Within this framework, there are some recent studies involving utility of temperature data for flux computation as primary step in AD [START_REF] Khaniki | Determining the heat flux absorbed by satellite surfaces with temperature data[END_REF], where limitations, resolution and accuracies of thermal data for heat flux determination are analyzed. Heat flux for AD is proposed in [START_REF] Khaniki | Satellite attitude determination using absorbed heat fluxes[END_REF], where a basic Levenberg-Marquardt (LM) deterministic least square algorithm is used. A major drawback of the LM algorithm is its inability for noise filtering that affects the attitude accuracy. In addition, the performance of the proposed HAM for attitude estimation is also demonstrated via EKF and UKF using MCS [START_REF] Labibian | Performance of a novel heat based model for spacecraft attitude estimation[END_REF]. Moreover, experimental validation of the proposed HAM has been investigated [START_REF] Labibian | Experimental validation of a novel radiation based model for spacecraft attitude estimation[END_REF], where a cubic laboratory satellite is used for testing within a vacuum chamber lit by a Sun simulator. There is also a commercial sensor system called Advanced Coarse Earth Sun Sensor (ACESS) [START_REF]products/satellite-equipment/cess-a-cess[END_REF]. It can provide an omnidirectional coarse estimation of the Earth and Sun position vector in the satellite reference frame. A combination of six ACESS packs (a total of 36 thermistors) is used to give the satellite navigation data. However, the proposed HAM requires only three non-contact thermopiles (NCT) for AD. A basic part of most estimation algorithms is the measurement model that is often a nonlinear function relating the measured data to variables whose proper estimates are needed. Therefore, the key contribution of this paper is focused on development of a novel measurement model for AD, namely the HAM. This model is based on the net heat flux (NHF) absorbed by three orthogonal, properly insulated, surfaces of a satellite. For a typical Earth-orbiting satellite, Sun and Earth are the main heat sources to be considered. As a result, the direct Sun radiation, the Earth flux and its albedo are the major means of heat exchange with the satellite surfaces. In order to establish a suitable HAM, later for estimation process, we need to establish a relation between the NHF of each satellite surface (measurements) and the satellite attitude (states).

In this paper, we represent each surface's NHF in terms of the components of the normal to surface unit vectors using the concept of view factors. Subsequently, the knowledge of satellite relative position vectors with respect to Sun and Earth is used to relate the satellite body system to the Earth Centred Inertial (ECI) system, which leads to establish the inertial satellite attitude. The proposed HAM is verified using thermal analysis software. The satellite attitude dynamics is useful in this study to simulate the satellite attitude behaviour in orbit. Therefore, we demonstrate effectiveness of the surface temperatures resulting from the HAM for an Earth-orbiting satellite whose attitude knowledge is simulated (or known).

This paper is organized as follows: Section 2 covers the necessary thermal relations for the development of the HAM. Section 3 integrates different parts of the HAM for attitude estimation. Simulation, evaluation and verification of the HAM for an Earth-orbiting cubic satellite is discussed in Section 4. Finally, Section 5 concludes the paper, and addresses future researches currently under development.

Main relations of radiation in space for the design of HAM

One of the major characteristics of space environment is its high vacuum. Hence, for most high orbiting satellites, the residual atmospheric pressure and the resulting drag will be negligible for the problem of attitude estimation. Therefore, the aerodynamic heating and convective interaction between the satellite and the space environment can be ignored [START_REF] Fortescue | Spacecraft Systems Engineering[END_REF].

Since the radiation is the key mean of heat exchange in space between a satellite and its surrounding, it can be characterized by the following factors (see Fig. 1).

• Radiation from the satellite to deep space.

• Direct solar radiation.

• Albedo radiation.

• Planetary (Earth) radiation.

Radiation emission from a satellite

Radiation is usually related to heat exchange through electromagnetic energy (with wavelength ranging from 0.1t o 100 µm) between a surface and the surrounding environment [START_REF] Karam | Satellite Thermal Control for System Engineers[END_REF]. All objects with a temperature more than absolute zero emit radiation rates proportional to the fourth power of their temperature. The most efficient emitter surface is called a black surface b. Radiation emission, denoted by q e b (W/m 2 ≡ Watt per unit Area), from a black surface at temperature T (K ≡ Kelvin) is expressed by the Stefan-Boltzmann law:

q e b = σ T 4 (1)
or

Q e b = A r σ T 4 (2)
where σ = 5.67 × 10 -8 W/m 2 K 4 is the Stefan-Boltzmann constant and A r (m 2 ) is the radiating area.

Radiation emission from a nonblack surface nb with similar area and temperature is expressed by:

Q e nb = ε A r σ T 4 (3)
where ε < 1i s known as the total hemispherical (or average) emissivity, which is the ratio of energy emission by the nonblack surface nb to that of the black surface b at the same temperature [START_REF] Karam | Satellite Thermal Control for System Engineers[END_REF].

Considering now two black surfaces, the amount of radiation heat exchange between them depends on their temperatures and how they view one another. The visibility of one surface as seen from the other is referred to a view (shape) factor. It defines the fraction of radiation leaving one surface and intercepted by another. The net radiation from a black surface element dA 1 at temperature T dA 1 to another interacting black surface element dA 2 at temperature T dA 2 is given by [START_REF] Karam | Satellite Thermal Control for System Engineers[END_REF]:

d 2 Q r dA 1 -dA 2 = σ T 4 dA 1 -σ T 4 dA 2 cos θ ′ 1 cos θ ′ 2 dA 1 dA 2 /π s 2 (4)
This results in a net radiation from A 1 to A 2 , which is given by:

Q r A 1 -A 2 = Ä 1 A 2 σ T 4 dA 1 -σ T 4 dA 2 cos θ ′ 1 cos θ ′ 2 dA 1 dA 2 π s 2 (5)
where θ ′ 1 and θ ′ 2 are the angles between surface normal vectors (of elements dA 1 and dA 2 ) and the line, of length s, that connects them (see Fig. 2). The superscript r indicates net radiation exchange, while the superscript e denotes the emission alone.

Assuming isothermal conditions for the two surfaces, the net radiation Q r 12 from surface 1 can be written in the following form:

Q r 12 = F 12 A 1 σ T 4 1 -σ T 4 2 ( 6
)
where F 12 is the shape factor expressed by [START_REF] Karam | Satellite Thermal Control for System Engineers[END_REF]:

F 12 = 1 A 1 Ä 1 A 2 cos θ ′ 1 cos θ ′ 2 π s 2 dA 1 dA 2 (7)
If dA 1 and dA 2 are considered as the elements on the satellite and the Earth, respectively, the integration of Eq. ( 7)y i e l d s the required shape factors which are shown in Fig. 3.

According to Fig. 3, shape factors are determined using a parameter, ν = a/R in which a and R are the satellite height and Earth radius respectively.

Solar radiation flux

Most measurements above the Earth's atmosphere have expressed a solar radiation flux value between 1365 and 1373 W/m 2 for LEO satellites. Solar radiation exists within a broad range of electromagnetic spectrum that includes about 7% in the ultraviolet (UV) wavelength (0.31-0.40 µm), 46% visible (0.40-0.69 µm), and 47% infrared (IR) radiation (above 0.70 µm). Since solar IR has shorter wavelength against the satellite emitted IR at normal satellite temperatures, it can be utilized for conditioning of a surface to have high reflectivity in the solar spectrum with high emissivity in IR at the same time. This property is reflected in solar absorptivity, α S , which is the fraction of straight solar energy absorbed by a surface. Therefore, the absorbed solar energy S α for a surface, whose normal vector forms an angle β with Sun direction, can be computed using Eq. ( 8):

S α = α S S cos β ( 8
)
where S is the solar flux.

Albedo radiation flux

The reflected heating sunlight from the Earth surface is called albedo. Usually, it is considered as a fraction of the solar constant and can be written as:

Q ′ Albedo = f a S (9)
where f a is known as the albedo factor. Albedo appears more significant at the Earth's polar ice caps, and can be estimated in those regions with some accuracy as a function of Sun's elevation and the satellite orbital parameters. However, predictions for overland and above oceans are distorted by variable effects of cloud formations and water distribution in the atmosphere.

Earth radiation flux

Emitted radiation from the Earth is considered equivalent in intensity and wavelength to the heat transferred from a black surface at -20 • C, which gives a nominal value of Earth flux approximately equal to 236 W/m 2 . For this value, a tolerance of ±38 W/m 2 (16%) is usually imposed in the satellite thermal analysis [START_REF] Gilmor | Spacecraft Thermal Control Handbook[END_REF][START_REF] Karam | Satellite Thermal Control for System Engineers[END_REF]. Since the Earth radiation is in the same band as that normally emitted by satellites, the fraction of impinging Earth flux absorbed by the satellite radiator is given by its emissivity ε. Hence, a surface treatment intended to reflect the Earth radiation will also reduce the space radiation through the same proportion as the surface's emission ability.

Design of heat-attitude model (HAM)

In order to establish a direct relation between the received heat fluxes by the surface of a typical satellite and its corresponding attitude, the orientation of a satellite in space thermal environment is presented in Fig. 4, where r Sat-Sun and r Sat-Earth denote the known satellite relative position vectors with respect to Sun and Earth in the ECI respectively.

The ultimate purpose is to find the attitude of the satellite with respect to ECI. Since the normal vectors of satellite surfaces that receive the heat fluxes are defined in body coordinate system, the relation between these vectors and satellite position vectors ( r Sat-Sun , r Sat-Earth ) can be established using shape factors computed along the parallel and the perpendicular directions with respect to the satellite position vectors. This concept is further clarified through Fig. 5.

As it is shown in Fig. 5, the unit vector k could be either in the direction of the satellite-Earth or satellite-Sun position vectors depending on the heat source. Therefore, k can be expressed as:

k = r Sat-Earth r Sat-Earth or k = r Sat-Sun
r Sat-Sun [START_REF] Markley | Attitude error representation for Kalman filtering[END_REF] While the distance between an Earth-orbiting satellite and the Sun is considerable, the satellite-Sun position vector is assumed to contain only the parallel component and thus the perpendicular component is neglected.

On the other hand, the relation between the ECI and the satellite coordinate system is also established through a rotation matrix, for which various attitude representations are possible. Considering the Euler's angles representation, the inertial to body rotation matrix is expressed as [START_REF] Sidi | Spacecraft Dynamics and Control[END_REF]:

C B I = cos θ cos ψ
cos θ sin ψ sin θ sin φ sin θ cos ψcos φ sin ψ sin φ sin θ sin ψ + cos φ cos ψ sin φ cos θ cos φ sin θ cos ψ + sin φ sin ψ cos φ sin θ sin ψsin φ cos ψ cos φ cos θ [START_REF] Shuster | Constraint in attitude estimation, part I: constrained estimation[END_REF] where φ, θ, ψ are the Euler's angles. C B I is only a function of the Euler's angles to be determined.

NHF decomposition

The NHF for each surface of a satellite can be represented as follows: [START_REF] Shuster | Constraint in attitude estimation, part II: unscented estimation[END_REF] where Q is the total heat received from three thermal sources, i.e. [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF] and Q lost is the emitted heat from satellite surfaces to the space environment that can be computed via Eq. ( 3).

Q net = Q -Q lost
Q = Q Sun + Q Earth + Q Albedo
Using the satellite-Earth and satellite-Sun position vectors, the heat can be represented as follows:

Q Sun = Aα S S cos β = Aα S S n• r Sat-Sun r Sat-Sun ( 14
)
Q Earth = AGε F h P E + F v 1 -P 2 E ( 15
)
where G is the Earth flux and F h and F v are the shape factors parallel and perpendicular to the satellite-Earth position vector, respectively to be determined from Fig. 3. P E is the projection of the satellite surface normal unit vector along the direction of the satellite-Earth position vector, and is defined as:

P E = n• r Sat-Earth r Sat-Earth ( 16
)
For the reflective albedo, that is a fraction of direct solar radiation in the direction of the satellite-Earth position vector, the heat can be written as [START_REF] Moghanipour | Modeling and Identification of Satellite Thermal Radiation Model for Attitude Estimation Using Temperature Sensors[END_REF]:

Q Albedo = Af a Sα S F h P E + F v 1 -P 2 E cos θ A ( 17
)
Fig. 6. Reflection of solar ray from the Earth surface [START_REF] Moghanipour | Modeling and Identification of Satellite Thermal Radiation Model for Attitude Estimation Using Temperature Sensors[END_REF].

in which θ A is the angle between the Sun-Earth vector and albedo reflection, see Fig. 6.

HAM for satellite attitude estimation

In order to establish the proposed HAM in a form of measurement equation required for estimation problem, we consider the relation between heat and temperature for each satellite surface:

Q i = m i c i dT i dt ( 18
)
where i = 1, 2, 3 denotes the number of orthogonal satellite surfaces, m i and c i are the mass and specific heat capacity of each surface respectively. Based on Eq. ( 12) and Eq. ( 18), the temperature rate for each satellite surface will be:

dT i dt = Q net m i c i ( 19
)
On the other hand, substituting Eq. ( 12)t h r o u g h Eq. ( 17)i n t o Eq. ( 19), yields the temperature rate of each orthogonal satellite surfaces:

dT i dt = 1 m i c i A i α S i S n i .C B I r I Sat-Sun r I Sat-Sun + A i Gε i F h P E i + F v 1 -P 2 E i + A i f a Sα S i F h P E i + F v 1 -P 2 E i cos θ A -A i ε i σ T 4 i ( 20
)
Finally, using the inertial to body transformation matrix, C B I from Eq. ( 11)i n t o Eq. ( 20) complete the HAM for three satellite orthogonal surfaces,

dT 1 dt = A 1 m 1 c 1 -ε 1 σ T 4 1 + Sα S 1 r S x cos(φ) cos(θ ) -r S z sin(θ ) + r S y cos(θ ) sin(ψ ) + a 1 Gε 1 + Sf a α S 1 cos(θ A ) a 1 = F v 1 -r E x cos(ψ) cos(θ ) -r E z sin(θ ) + r E y cos(θ ) sin(ψ) 2 + F h r E x cos(ψ ) cos(θ ) -r E z sin(θ ) + r E y cos(θ ) sin(ψ ) (21)
dT 2 dt = -A 2 m 2 c 2 ε 2 σ T 4 2 + b 1 Gε 2 + Sf a α S 2 cos(θ A ) + Sα S 2 -r S x b 3 + r S y b 4 + r S z cos(θ ) sin(φ) b 1 = F h b 2 + F v 1 -b 2 2 b 2 =-r E x b 3 + r E y b 4 + r E z cos(θ ) sin(φ) (22) b 3 = cos(φ) sin(ψ ) -cos(ψ ) sin(φ) sin(θ ) b 4 = cos(φ) cos(ψ ) + sin(φ) sin(ψ ) sin(θ ) dT 3 dt = A 3 m 3 c 3 ε 3 σ T 4 3 + c ′ 1 Gε 3 + Sf a α S 3 cos(θ A ) + Sα S 3 r S x c ′ 4 -r S y c ′ 3 + r S z cos(φ) cos(θ ) c ′ 1 = F h c 2 + F v 1 -c ′ 2 2 c ′ 2 = r E x c 4 -r E y c 3 + r E z cos(φ) cos(θ ) (23)
c ′ 3 = cos(ψ ) sin(φ) -cos(φ) sin(ψ ) sin(θ ) c ′ 4 = sin(φ) sin(ψ ) + cos(φ) cos(ψ ) sin(θ )
where, r S x , r S y , r S z are the components of r I Sat-Sun vector and, similarly, r E x , r E y , r E z are the components of r I Sat-Earth vector, described previously.

In general, the measurement equation takes usually the following form:

z = h(x) + v (24)
where h (Eq. ( 21)t h r o u g h Eq. ( 23)) relates the measurements z (temperature rates) to the states x (attitude) and v is the sensor measurement noise, which is assumed to be a white Gaussian noise.

Therefore, by considering Eq. ( 21)t h r o u g h Eq. ( 23) and Eq. ( 24), the measurement vector is composed of the temperature rates (to be measured via temperature sensors in appropriate time intervals):

z = dT 1 dt dT 2 dt dT 3 dt T + v (25)

Results and discussion

To resume, the following summary establishes clearly the relation between the NHF and the satellite's Attitude.

• Based on Eqs. ( 14), [START_REF] Norgaard | New development in state estimation for nonlinear systems[END_REF] and [START_REF] Kiani | Concurrent orbit and attitude estimation using minimum sigma point unscented Kalman filter[END_REF], the received heat is related to the shape factors and the components of n for each surface of the satellite. Note that the satellite relative position vectors toward the major heat sources (Sun and Earth) are known in ECI.

• The body components of each surface normal unit vector is considered as a row in the transformation matrix, namely C B I to be determined (see Eq. ( 11)).

• Since the nine elements of the rotation matrix can be written in terms of three independent variables (namely the Euler's angles) that are explicitly present in Eq. ( 21)t h r o u g h Eq. ( 23), there are three equations and three unknowns to be determined.

• Solution of this system of equations will deterministically establish the satellite attitude using the proposed HAM.

It is important to emphasize that the key contribution of this paper is to develop a novel HAM and verify its reliability for an Earth-orbiting satellite. The verification process, which requires the temporal variation of the rotation matrix for continuous attitude determination, is accomplished using the satellite attitude dynamics and will be discussed in the next section.

Verification of heat-attitude model

For the purpose of this paper, the Euler's angles representation is chosen to express the rotational kinematics of the satellite. Through the use of the Euler's equations ( 26) and the kinematic equations ( 27), the satellite rotational dynamics and attitude can be written as [START_REF] Wie | Spacecraft Vehicle Dynamics and Control[END_REF]: where ω =[ω x ω y ω z ] T is the angular velocity of the satellite expressed in the body coordinate system. By considering Eq. ( 26) and Eq. ( 27) and using Eq. ( 21)t h r o u g h Eq. ( 23), the temperature, as well as the NHF that each surface of satellite receives, will be attainable. Table 1 presents the orbital characteristics used for the simulation of HAM. In addition, Table 2 provides some system features required for the simulations that are performed using MATLAB [START_REF][END_REF], as well as the initial states of the nadir pointing satellite. Specific surface features and preparations required for proper applicability of the developed HAM are beyond the scope of this work. However, some of the thermal data which are used in thermal part of the model are presented in Table 3. The results of simulations using the data of Tables 1, 2 and 3 are presented in Figs. 78910.

In order to partially verify the proposed HAM, commercial thermal analysis tools are used for comparison of results. Thermal Desktop [38] and SINDA [39]a r e used to obtain the NHFs and the satellite surface temperatures, respectively. For comparative purposes, the temperature of three orthogonal surfaces resulted out of HAM are compared with those of the Thermal Desktop in Fig. 11 through Fig. 13, whose trends and accuracy are well compatible. As shown in Figs. 11 through 13, the HAM surface temperature results closely follow those obtained via the thermal analysis tools with good accuracy. This observation verifies the HAM appropriateness to represent a measurement model for AD stochastic formulation. Based on the above results and those presented in [START_REF] Labibian | Performance of a novel heat based model for spacecraft attitude estimation[END_REF] and [START_REF] Labibian | Experimental validation of a novel radiation based model for spacecraft attitude estimation[END_REF], a performance comparison is presented between the proposed HAM package against the other existing AD systems in Table 4.

Conclusion

A novel radiation based heat model for attitude estimation of satellites is proposed in this paper. The heat-attitude model (HAM) uses external radiating heat sources affecting the satellite external surface temperatures to produce the net heat flux (NHF) that each satellite surface receives. Decomposition of the NHFs in parallel and perpendicular components with respect to the heat sources allows for generation of the inertial to satellite rotation matrix that, in turn, leads to its attitude knowledge. Simulation results of the proposed HAM are verified via commercially available thermal analysis tools for a typical LEO satellite. The proposed HAM satellite temperature results are in good agreement with those of the available commercial thermal softwares. Thus, the proposed HAM can be considered as a viable measurement model for stochastic AD via a large class of filtering techniques, some of which have been already utilized by researchers.

Fig. 1 .

 1 Fig. 1. Thermal environment of a typical spacecraft [33].

Fig. 2 .

 2 Fig. 2. Radiation heat exchange between two surfaces.

Fig. 3 .

 3 Fig. 3. Shape factors from a surface element to a sphere [23].

Fig. 4 .

 4 Fig. 4. Orientation of a typical satellite with respect to Sun and Earth.

Fig. 5 .

 5 Fig. 5. Decomposition of the satellite surface normal vector n.

  ωx = M x -ω y ω z (I z -I y ) I x ωy = M y -ω x ω z (I x -I z ) I y (26) ωz = M z -ω x ω y (I y -I x ) I z φ = ω x + ω y sin(φ) + ω z cos(φ) tan(θ ) θ = ω y cos(φ)ω z sin(φ) (27) ψ = ω y sin(φ) + ω z cos(φ) sec(θ )

Fig. 7 .

 7 Fig. 7. Satellite angular velocities.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Satellite Euler angles.

Fig. 10 .

 10 Fig. 10. Satellite surface temperatures.

Fig. 11 .

 11 Fig. 11. Temperature variations of surface 1.

Fig. 12 .

 12 Fig. 12. Temperature variations of surface 2.

Fig. 13 .

 13 Fig. 13. Temperature variations of surface 3.

  = 0.005 deg to 3 deg 0.1 to 2 0 to 3 Star sensors Accuracy = 1a r c s e c o n dt o1a r c m i n u t e 2t o5 5t o2 0 (Scanners and cameras) (= 0.0003 deg to 0.01 deg) Horizon sensors Accuracy: 1 to 4 5 to 10 Scanner/pipper 0.05 deg to 1 deg (0.1 deg is best for LEO) 0.5 to 3.5 0.3 to 5 F i x e dh e a d( s t a t i c ) <0.1 deg to 0.25 deg Magnetometer Accuracy = 0.5 deg to 3 deg 0.3 to 1.2

Table 1

 1 Satellite orbital parameters.

	Eccentricity (e)	0
	Right ascension ( )6 0 •
	Inclination (i)	40.5 •
	Argument of perigee (ω)	270 •
	Orbital height	460 km

Table 2

 2 Satellite system features.

	Satellite mass	49.86 kg
	Satellite volume	1 m × 1m× 1m
	Initial angular velocity (ω 0 )	[0 -0.063 0] T deg/sec
	Initial Euler angles	[-49.506 0 ] T deg

Table 3

 3 Satellite thermal model properties.

	Solar radiation heat flux	1414 W/m 2
	Earth radiation heat flux	257 W/m 2
	Solar absorptivity	0.33
	Albedo factor	0.26
	Surface material	Aluminium
	Heat specific capacity	961 J/kg K
	Surface emissivity factor	0.509

Table 4

 4 Existing AD systems.

	Sensors	Typical performance range	Mass (kg)	Power (W)
	Gyroscopes	Drift rate = 0.003 deg/h to 1 deg/h	<0.1 to 15	<1t o2 0 0
		Drift rate stability varies widely
	Sun sensors	Accuracy

Powered by TCPDF (www.tcpdf.org)

Acknowledgement

Authors thank Dr. Mohsen Abedi for his consultation on thermal modelling utilized in parts of this research.

Conflict of interest statement

The authors declare that there is no conflict of interest over this article.