Marjolein De Jong-Hagelstein

Mieke Van De Sandt-Koenderman

MD Niels D Prins

MD Diederik W J Dippel

MD Peter J Koudstaal

MA Evy G Visch-Brink

De Jong-Hagelstein

MA W M E Van De Sandt-Koenderman

The Efficacy of Early Cognitive-Linguistic Treatment and Communicative Treatment in Aphasia after Stroke -A Randomized Controlled Trial (RATS-2)

Keywords:

come

INTRODUCTION

Aphasia is present in about 30% of all acute stroke patients and affects their daily communication and social participation. There are two main approaches in aphasia treatment: cognitive-linguistic treatment (CLT) and communicative treatment. The ultimate goal of both approaches is to improve patients' everyday communication. They differ fundamentally in how they achieve this: CLT focuses on the impairment and aims at improving the underlying linguistic processing at the linguistic levels affected, e.g. semantics (word meaning) [START_REF] Kiran | Semantic complexity in treatment of naming deficits in aphasia: evidence from well-defined categories[END_REF], phonology (word sound) [START_REF] Kendall | Phoneme-based rehabilitation of anomia in aphasia[END_REF] or syntax (understanding and building sentences). [START_REF] Thompson | Complexity in treatment of syntactic deficits[END_REF] Communicative treatment focuses on the disability: patients are trained to use their residual language skills combined with compensatory strategies in order to optimize information transfer. [START_REF] Elman | The efficacy of group communication treatment in adults with chronic aphasia[END_REF] It is unclear which of both approaches is best for which patients in which period of recovery. In the extensive literature about the efficacy of aphasia treatment there is more evidence for the efficacy of CLT -recommended as Practice Standard in 2005 [START_REF] Cicerone | Evidence-Based Cognitive Rehabilitation: Updated Review of the Literature From 1998 Through 2002[END_REF] -than of communicative treatment, which has been evaluated less frequently. However, a metaanalysis and recent reviews on cognitive rehabilitation have emphasized the need for welldesigned trials on aphasia treatment in general and on specific treatments [START_REF] Cicerone | Evidence-Based Cognitive Rehabilitation: Updated Review of the Literature From 1998 Through 2002[END_REF][START_REF] Robey | A meta-analysis of clinical outcomes in the treatment of aphasia[END_REF][START_REF] Greener | Speech and language treatment for aphasia following stroke[END_REF][START_REF] Cappa | EFNS guidelines on cognitive rehabilitation: report of an EFNS task force[END_REF] with a sufficient sample size, a functional outcome measure and well-defined methods of intervention. [START_REF] Cappa | EFNS guidelines on cognitive rehabilitation: report of an EFNS task force[END_REF] In the last three decades 11 randomized controlled trails (RCTs) were conducted on a specific treatment method for aphasic stroke patients delivered by a speechlanguage therapist (SLT). To date, no RCTs have been conducted in which the benefits of both approaches, in the form of individual treatment, are compared.

A factor that may influence the efficacy of treatment is timing. Meta-analyses of uncontrolled studies and RCTs suggested that the largest improvements after language treatment occur within one year post-injury [START_REF] Robey | A meta-analysis of clinical outcomes in the treatment of aphasia[END_REF][START_REF] Rohling | Effectiveness of Cognitive Rehabilitation Following Acquired Brain Injury: A Meta-Analytic Re-Examination of Cicerone et al.'s[END_REF] and mainly when treatment was started within the first three months. [START_REF] Robey | A meta-analysis of clinical outcomes in the treatment of aphasia[END_REF] In these analyses the type of treatment was not controlled. It is very well possible that CLT and communicative treatment differentially interact with time post-onset.

Code [START_REF] Code | Multifactorial Processes in Recovery from Aphasia: Developing the Foundations for a Multileveled Framework[END_REF] poses that treatment aimed at restoration of impaired cognitive processes is probably more appropriate in acute stages when natural recovery occurs. Addressing specific neural networks, involved in semantics and phonology, by specific treatment activities (CLT) might facilitate or speed up neural recovery processes. Hence, starting early may be crucial for the efficacy of CLT, but less important for communicative treatment. This view is reflected in the current preference in many centers to give CLT in the acute stage followed by communicative treatment when a plateau in improvement has been reached.

In our previous RCT, RATS-1 [START_REF] Doesborgh | Effects of Semantic Treatment on Verbal Communication and Linguistic Processing in Aphasia After Stroke: A Randomized Controlled Trial[END_REF], semantic treatment (BOX) [START_REF] Visch-Brink | een semantisch therapieprogramma[END_REF] was compared with phonological treatment (FIKS) [START_REF] Van Rijn | FIKS, een fonologisch therapieprogramma[END_REF], applied 4-12 months post-stroke. The two treatments appeared to be equally effective in improving verbal communication (Amsterdam-Nijmegen Everyday Language Test, ANELT-A). [START_REF] Blomert | Amsterdam-Nijmegen Test voor Alledaagse Taalvaardigheden[END_REF] In the current study, RATS-2, we therefore combined BOX and FIKS and compared this CLT with communicative treatment to evaluate their efficacy in an earlier stage of aphasia.

Our objective was to measure the efficacy of CLT, applied in the first six months starting within three weeks post-stroke, on everyday verbal communication and on semantic and phonological processing. We hypothesized that CLT would be more effective than communicative treatment and that its effect would be the largest in the first three months.

METHODS

Participants

All patients with aphasia after intracerebral hemorrhage or ischemic stroke of less than three weeks duration were screened for eligibility by the local speech-language therapist (SLT) of 15 hospitals in the Netherlands and Belgium. We included patients aged 18-85 with a life expectancy of more than six months.

Apart from a disorder in verbal communication as measured with the ANELT-A (score <44/50), a semantic and/or phonological disorder had to be present. A semantic disorder implied a score on Semantic Association Test-verbal [START_REF] Visch-Brink | De Semantische Associatie Test[END_REF] of less than 26/30 and/or a score on Semantic Association (PALPA) [START_REF] Bastiaanse | Psycholinguïstische testbatterij voor de taalverwerking van afasiepatiënten[END_REF] of less than 12/15. A phonological disorder implied a score on Nonword Repetition Task [START_REF] Bastiaanse | Psycholinguïstische testbatterij voor de taalverwerking van afasiepatiënten[END_REF] of less than 20/24 and/or on Auditory Lexical Decision [START_REF] Bastiaanse | Psycholinguïstische testbatterij voor de taalverwerking van afasiepatiënten[END_REF] of less than 76/80.

Exclusion criteria were severe dysarthria, developmental dyslexia or visual perceptual disorder; pre-existing aphasia, premorbid dementia and recent psychiatric disorder.

Interventions

Experimental treatment

Cognitive-linguistic treatment (CLT) consisted of BOX, a semantic treatment program, and/or FIKS, a phonological treatment program (paper and computer versions). BOX contains many semantic decision tasks using written words, sentences and texts that may also be presented orally. BOX aims to enhance semantic processing. FIKS has a similar structure but is directed at the phonological input and output routes. The SLT determined which treatment program(s) and which subparts the patient needed.

Control treatment

Communicative treatment aimed at improving communicative ability using all verbal and nonverbal strategies available to the patient, e.g. written choice communication and communication books. By definition, exercises are personally relevant and embedded in a communicative setting. Examples of methods used are PACE [START_REF] Davis | Adult aphasia rehabilitation: Applied pragmatics[END_REF] (Promoting Aphasics' Communicative Effectiveness), role playing and conversational coaching.

Assessment

Baseline measures, including patient demographics and date and type of stroke, were recorded before randomization. The assessment at baseline and at three and six months post-stroke consisted of various linguistic measures and a measure of disability (Table 1).

Table 1 Reported assessments

Linguistics

Semantic measures

• Semantic Association Test (SAT) [START_REF] Visch-Brink | De Semantische Associatie Test[END_REF], verbal version. The SAT is based on the principles of the Pyramids and Palm Trees Test. [START_REF] Howard | The Pyramids and Palm Trees Test[END_REF] The patient chooses from four written words (three semantically related words and an unrelated word) the word that is semantically closest to the target word. • Semantic Association with low-imageability words (PALPA). [START_REF] Bastiaanse | Psycholinguïstische testbatterij voor de taalverwerking van afasiepatiënten[END_REF] • Semantic Word Fluency: animals and professions. Phonological measures

• Nonword Repetition Task (PALPA).

• Auditory Lexical Decision (PALPA). The patient decides if words are existing or nonexisting. • Letter Fluency: D, A and T. Other linguistic measures

• Amsterdam-Nijmegen Everyday Language Test (ANELT). [START_REF] Blomert | Amsterdam-Nijmegen Test voor Alledaagse Taalvaardigheden[END_REF] Verbal responses to ten everyday language scenarios are scored on a 5-point scale for informational content (scale A). • Aachen Aphasia Test [START_REF] Graetz | Akense Afasie Test. Nederlandstalige versie[END_REF] (only at 6-8 weeks post-stroke).

Disability

• Modified Rankin Scale. [START_REF] Van Swieten | Interobserver agreement for the assessment of handicap in stroke patients[END_REF] The Amsterdam-Nijmegen Everyday Language Test (ANELT) consists of ten scenarios to which the patient has to respond verbally. For example: You are in a store and you want to buy a television. I am the salesperson here. "Can I help you?". Patients' verbal responses are rated for informational content on scale A "understandability", and for articulation on scale B (intelligibility"). The scales are from 1 (bad) to 5 (good) so the total score on both scales ranges from 10 to 50. The ANELT is both a valid test (ecological validity is strong, criterion-related validity is .81, construct validity is good) and a reliable test (inter-rater reliability is .92, test-retest reliability is .92, goodness of fit is >0.91). [START_REF] Blomert | Amsterdam-Nijmegen Everyday Language Test: Construction, reliability and validity[END_REF] Procedure This trial was approved by the Ethical Committee of Erasmus MC and is registered (ISRCTN67723958). Written informed consent was obtained from all patients and their proxy before enrollment.

Patients were assessed as soon as possible as from day three. Subsequently, they were included in the study and the allocated treatment was started three weeks post-stroke at the latest. Treatment was provided for six months or shorter if the patient had completely recovered. Patients were retested at three and six months. The assessment and treatment took place in patients' subsequent treatment settings or at home.

Treatment was applied with a minimum of two and preferably for five hours per week, partly individual and partly as homework. The SLTs wrote down the content and amount of treatment their patients received on registration forms that were returned to us and discussed this with us every two to three weeks.

Outcomes

The primary outcome was the ANELT scale A (understandability) at six months. This scale measures functional verbal communication. The ANELT-A was scored from audiotapes by two independent, experienced SLTs, blinded to test moment and treatment allocation. The means of both raters' scores were used in the analyses. In case of a difference between two scores of ≥7 points, the raters were asked, without giving further information, to score the particular test again. In the few cases in which the difference was still ≥7 points, the scores were averaged with that of a third rater.

The secondary outcome measures were the ANELT-A at three months, three semantic tests: Semantic Association Test-verbal, Semantic Association with lowimageability words and Semantic Word Fluency; and three phonological tests: Letter Fluency, Auditory Lexical Decision and Nonword Repetition Task.

The primary effect measure was the difference between the two treatment groups in mean score on the ANELT-A at six months. This difference at three months was a secondary effect measure. The other secondary effect measures were the following: the difference in proportion of patients who improved ≥7 points, the clinically relevant change for individual patients on the ANELT (critical difference) according to the test manual [START_REF] Blomert | Amsterdam-Nijmegen Test voor Alledaagse Taalvaardigheden[END_REF]; the difference in proportion of patients who, after treatment, fell in the ANELT category "moderate" or "mild/no communication disorder", score 30-50; and finally, the difference in score on the three semantic and three phonological tests, all at three and six months. The assessments were done by members of the research group, of a supporting aphasia team, trainees, or involved or uninvolved colleagues of the treating SLT. In the minority of cases it was unavoidable that the tests were administered by the treating SLT. One hundred of 158 follow up assessments were carried out by a person who was blind for treatment allocation. The assessors were instructed both orally and through a manual on how to administer the tests. In addition, they were trained specifically in administering the ANELT.

Statistical analyses

Analyses were performed on the basis of intention to treat. We also performed an on treatment analysis by limiting the analysis to patients who had completed treatment. We used ANCOVA and adjusted for baseline severity [START_REF] Vickers | Analysing controlled trials with baseline and follow up measurements[END_REF] to test group differences in score on the ANELT-A and on the semantic and phonological tests at three and six months, with 95% confidence interval (CI). Beside baseline severity, we planned to adjust for age, gender, level of education, aphasia type, and intensity of treatment. The proportion of patients in each group who improved ≥7 points after three and six months, and the proportion of patients in each group who fell in the ANELT category "moderate" or "mild/no communication disorder" were compared by Odds Ratio with 95% CI by means of multiple logistic regression. All analyses were performed in SPSS 15.0 for Windows.

Sample size

We calculated that a sample of 70 patients would provide a power of 0.87 to detect a difference of seven points, the clinically relevant change for individual patients, on the ANELT between the two treatment groups at a 5% 2-sided significance level. To compensate for non-evaluable patients we randomized 80 patients.

Randomization and blinding

Treatment allocation was stratified by center. An independent statistician provided the computer-generated random allocation sequence per center. An uninvolved secretary put the assignments per center in sequentially numbered, opaque envelopes that were sealed and stored in a locked drawer. The research coordinator wrote the research number of the patient on the next appropriate envelop and then opened it to assign the intervention.

The patients and SLTs could not be blinded to individual treatment allocations, but the assessment of the primary outcome was blinded. Patients' responses on the ANELT were tape-recorded and scored by two independent raters, blinded to test moment and treatment allocation.

RESULTS

From September 2006 to April 2008, 85 patients were enrolled in 27 treatment centers in The Netherlands and Belgium. In 3 of 41 patients assigned CLT and in 2 of 44 patients assigned communicative treatment only baseline assessment was obtained and no follow-up due to serious concomitant illness, death or refusal to further participate. Because no outcomes could be determined in these five patients, we do not report on them. The intention-to-treat group therefore consisted of 80 patients who had received at least one follow-up assessment and in whom outcomes could be determined. For the on-treatment analyses we excluded 5 of 80 patients who prematurely aborted treatment (see details in Figure 1).

Baseline characteristics were well balanced (Table 2), except for gender: there were more men in the control group (57% versus 37% in the CLT group).

Primary outcome

In both treatment groups the average ANELT-A scores improved. There was no significant difference in the mean ANELT-A scores of the two treatment groups, neither at three nor at six months post-stroke (Table 3). Almost all improvement occurred in the first three months. In total, four ANELT scores were lacking, two in the CLT group and two in the control group, due to loss of the audiotape, death, residence abroad and refusal. We replaced these lacking scores by the sum of the patient's ANELT-A score on the previous assessment and the mean improvement of the whole group in the previous period.

There was also no difference in improvement of ≥7 points on the ANELT-A. At three months 22/38 patients (58%) in the CLT group improved ≥7 points, compared with 26/42 (62%) in the communicative treatment group (OR=0.85, 95% CI=0.35 to 2.07). At six months 27/38 patients (71%) in the CLT group improved ≥7 points, compared with 31/42 (74%) in the communicative treatment group (OR=0.87, 95% CI=0.33 to 2.33).

At three months there was a trend regarding the proportion of patients in each group who fell in the ANELT category "moderate" or "mild/no communication disorder": 27/38 patients (71%) from the CLT group, versus 23/42 (55%) in the communicative treatment group (OR=2.0, 95% CI=0.80 to 5.13). But this trend was not present anymore at six months (29/38 [76%] in the CLT group versus 30/42 [71%] in the communicative treatment group, OR=1.3, 95% CI=0.47 to 3.52).

In the on-treatment analyses, with five patients less than in the intention to treat group, the treatment effects were much the same. Adjustment for neither the baseline characteristics age, gender, and level of education, nor for the variables aphasia type and intensity of treatment changed the results of the main outcomes.

The intraclass correlation coefficients between the two independent raters of the ANELT indicated excellent agreement (at baseline 0.95, at three months 0.97 and at six months 0.96).

Secondary outcome measures

Both treatment groups improved on all secondary tasks. There was a statistically significant difference between the groups on the fluency tasks, in favor of CLT: on Semantic Word Fluency at three months (adjusted difference=3.2, 95% CI=0.4 to 6.0) and on Letter Fluency at six months (adjusted difference=3.1, 95% CI=0.3 to 6.0). On the remaining secondary tasks, Semantic Association Test-verbal, Semantic Association with low imageability words, Auditory Lexical Decision and Nonword Repetition Task, there was no significant difference in improvement between the groups (Figure 2).

Treatment intensity

The mean intensity of treatment was 2.1 hours. The mean number of hours of therapy was 45.4 hours of which 33.8 hours were face-to-face with an SLT and 11.6 hours were spent on homework. More details are given in Table 4. Intensity did not differ significantly between the two treatment groups (0-3 months: p=0.2; 3-6 months: p=0.5). A smaller percentage of patients in the communicative treatment group received the minimum treatment intensity of two hours compared with patients in the CLT group. However, this difference was not statistically significant (0-3 months: χ 2 =1.7, p=0.2; 3-6 months: χ 2 =2.8, p=0.1).

DISCUSSION

We compared the efficacy of CLT (semantic and phonological treatment) and communicative treatment, applied in the first six months post-stroke. The two groups showed an equal improvement on the ANELT-A at three and six months with the largest increase in the first three months. The proportion of patients who fell in the ANELT category "moderate" or "mild/no communication disorder" showed a trend in favor of CLT at three months, but not at six months. This difference at three months might be due to chance: more patients in the CLT group than in the control group scored near the cut-off and thus were more likely to shift from the severe to the moderate category. The scores on nearly all specific semantic and phonological tests were higher after CLT than after communicative treatment, but the difference was only significant for the semantic and letter fluency tasks.

This treatment effect in favour of CLT could be meaningful. Both fluency tasks are explicitly related to the aim of CLT, i.e. to improve semantic and phonological processing, which has a positive influence on word finding. The fluency tasks are productive tasks that require self generation of words, a stage in the pursuit of adequate verbal communication in general, measured by the ANELT. Although there was no significant overall treatment effect on our primary outcome measure, nearly all differences between groups were in favour of the CLT group and therefore the efficacy of CLT remains to be evaluated in future studies.

Some methodological aspects of our study should be discussed. To our knowledge, RATS-2 is the first RCT that has evaluated aphasia treatment started in the acute stage, with everyday language use as primary outcome. Other strengths were the relatively large sample size (n=80), very good compliance and minimal loss to follow-up.

In line with recommendations for efficacy research on cognitive rehabilitation [START_REF] Cicerone | Evidence-Based Cognitive Rehabilitation: Updated Review of the Literature From 1998 Through 2002[END_REF], we used a functional communication measure as primary outcome. In the Netherlands, the ANELT, in origin Dutch, was the only adequate test available. It has a high ecological validity [START_REF] Blomert | Amsterdam-Nijmegen Everyday Language Test: Construction, reliability and validity[END_REF] and is increasingly applied as primary outcome in treatment studies, both in the Netherlands [START_REF] Doesborgh | Effects of Semantic Treatment on Verbal Communication and Linguistic Processing in Aphasia After Stroke: A Randomized Controlled Trial[END_REF][START_REF] Doesborgh | Cues on request: the efficacy of Multicue, a computer program for wordfinding therapy[END_REF][START_REF] Bastiaanse | The training of verb production in Broca's aphasia: A multiple-baseline across-behaviours study[END_REF] and elsewhere.[25] It is considered a weakness if, as in most studies, the outcome measure is trained material. Fewer studies considered generalisation to untrained material. The ultimate result is generalisation of the intervention to everyday communication, which we aimed to show. Therapy-specific findings in our previous trial [START_REF] Doesborgh | Effects of Semantic Treatment on Verbal Communication and Linguistic Processing in Aphasia After Stroke: A Randomized Controlled Trial[END_REF] support the view that improving the underlying linguistic processing results in improved verbal communication.

Although we could not avoid that about one third of the outcome assessments was non-blinded, the final judgement of all ANELT samples, the primary outcome measure, was blinded.

A limitation is that we did not include a control group without language treatment. Therefore, we are not able to specify the potential effect of treatment over natural recovery. Because our aim was to measure the efficacy of a 6-months treatment period, we considered it both unethical and impracticable to withhold treatment from patients with a recent stroke for such a long period. Also, one may question whether the contrast between both treatments was large enough. Although only one patient in our study received less than 75% of treatment according to protocol, it is obvious that in any communicative exercise, semantic and phonological processes are implicitly addressed. Therefore, we cannot exclude that overlap between the two treatments may have played a role. Finally, treatment intensity is currently an important issue in efficacy research and is assumed to be vital for the efficacy of treatment. [START_REF] Cicerone | Evidence-Based Cognitive Rehabilitation: Updated Review of the Literature From 1998 Through 2002[END_REF][START_REF] Robey | A meta-analysis of clinical outcomes in the treatment of aphasia[END_REF][START_REF] Breitenstein | Intensive language training for aphasia. Contribution of cognitive factors[END_REF] Possibly, the intensity in our study did not reach the threshold necessary to exceed natural recovery and find potential treatment effects [START_REF] Bakheit | A prospective, randomized, parallel group, controlled study of the effect of intensity of speech and language therapy on early recovery from poststroke aphasia[END_REF] as we did not succeed in achieving the preferred intensity of five hours per week (the mean was 2.1 hours). A metaanalysis [START_REF] Bhogal | Intensity of aphasia treatment, impact on recovery[END_REF] suggested that 8.8 hours of treatment per week is needed to obtain a treatment effect and that two hours per week is insufficient. Recommendations in the remaining literature on treatment intensity range from 1.5 to two hours per week as being too little [START_REF] Robey | A meta-analysis of clinical outcomes in the treatment of aphasia[END_REF][START_REF] Lincoln | Effectiveness of speech therapy for aphasic stroke patients. A randomised controlled trial[END_REF], to two or three hours as the minimum to obtain positive results. [START_REF] Robey | A meta-analysis of clinical outcomes in the treatment of aphasia[END_REF][START_REF] Salter | Evidence-based review of stroke rehabilitation: Aphasia. Version 11[END_REF][START_REF] Basso | Aphasia and its therapy[END_REF] Of the few well-designed RCTs on the efficacy of aphasia treatment, the one of Wertz et al. [START_REF] Wertz | Veterans administration cooperative study on aphasia: a comparison of individual and group treatment[END_REF] is most comparable to ours. These authors compared treatment of specific language deficits with communicative treatment, started at four weeks post-onset, and found that the two were equally effective. Communicative treatment, however, was provided in a group instead of individually. In most other studies, conventional treatment was used, so approaches were mixed. Constraint-induced Aphasia Therapy [START_REF] Pulvermueller | Constraint-induced therapy of chronic aphasia after stroke[END_REF], a treatment that is to date evaluated best, also combines a cognitive-linguistic and communicative approach. Elman&Bernstein-Ellis [START_REF] Elman | The efficacy of group communication treatment in adults with chronic aphasia[END_REF] examined the efficacy of group communication treatment separately, and found higher scores on communicative and linguistic measures compared to no treatment. CLT directed to semantic and phonological processing, which are crucial to word finding, has not been contrasted with no treatment in an RCT.

Our study does not support the widespread notion that CLT is more appropriate in an early stage and that communicative training is more suitable at a later stage. [START_REF] Code | Multifactorial Processes in Recovery from Aphasia: Developing the Foundations for a Multileveled Framework[END_REF] Nor do the results support the recommendation that, in cognitive rehabilitation, clinicians should focus on training cognitive skills directly rather than broad interventions with the expectation of subsequent generalization to broader use in daily life. [START_REF] Rohling | Effectiveness of Cognitive Rehabilitation Following Acquired Brain Injury: A Meta-Analytic Re-Examination of Cicerone et al.'s[END_REF] The results also do not support the hypothesis that treatment of communication via the activation of the underlying processes, i.e. semantics and phonology, would be more effective in early aphasia, when natural recovery takes place, than a direct training of the communication itself. Because of the possible overlap between the two treatments and the low treatment intensity, the question whether CLT is efficacious particularly in the acute stage remains open. Therefore, in our next study, RATS-3, we aim to compare the effect of intensive CLT in aphasia patients very early post-stroke, with deferred treatment.

COPYRIGHT LICENCE STATEMENT

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non-exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd and its Licensees to permit this article (if accepted) to be published in Journal of Neurology, Neurosurgery & Psychiatry and any other BMJPGL products to exploit all subsidiary rights, as set out in our licence.

Figure 1 Figure 2 4 *

 124 Figure 1 Flow chart of enrolled patients Figure 2 Difference between the mean improvement on the secondary outcome measures of the CLT (n=38) and the communicative treatment group (n=42)

Table 2

 2 Baseline characteristics of enrolled patients

	CLT (n=38)	Communicative treatment (n=42)

Table 3

 3 Primary outcome measure: mean ANELT-A scores for the CLT and the communicative treatment group

		CLT	Communicative	Difference	Adjusted	P
		(n=38)	treatment	(95% CI)	difference	value
			(n=42)		(95% CI)*
	3 months	33.4	31.6	1.8	1.5	0.48
	post-			(-3.8 to 7.4)	(-2.6 to 5.6)
	stroke
	6 months	35.2	33.2	1.9	1.6	0.42
	post-			(-3.4 to 7.3)	(-2.3 to 5.6)
	stroke †

* Adjusted for baseline score; † Primary effect measure

Table 4

 4 Data on treatment intensity for both treatment groups

		0 -3 months	3 -6 months
		CLT	Communicative	CLT	Communicative
		(n=38)	(n=41)	(n=38)	(n=41)
	Mean intensity (SD)	2.4 (0.9)	2.0 (1.0)	2.1 (1.2)	2.0 (1.2)
	Number of patients	25 (66%)	21 (51%)	21 (55%)	15 (37%)

Funding

This study was funded by the Stichting Nuts Ohra (T-07-71). Stichting Nuts Ohra had no involvement in the study design, the data, the report or publication.

ACKNOWLEDGEMENTS, COMPETING INTERESTS, FUNDING Acknowledgements

We wish to thank the patients and the speech-language therapists who provided and treated study patients of the following centers: Hospitals: Albert Schweitzer, AMC, Atrium MC, academisch ziekenhuis Maastricht, Erasmus MC, Haven, Ikazia, Maasstad, Maria Middelares, Middelheim, Reinier de Graaf Gasthuis, Sint Franciscus Gasthuis, Sint Lucas Andreas, Vlietland, IJsselland, VUmc. Rehabilitation centers: Hoensbroeck, RC Amsterdam, Rijndam, Sophia Den Haag and Delft and Gouda, De Waarden Dordrecht and Gorinchem. Nursing homes: Antonius Binnenweg and IJsselmonde, Berkenstede, Bernardus, Bieslandhof, De Burcht, Cicero, Crabbehoff, Elf Ranken, Invia, Klevarie, Meerweide, Naaldhorst, Pniël, De Poort, Rheuma, Schiehoven-Wilgenplas, Sevagram, Sint Jacob, Sterrelanden, Twee Bruggen, Vijf Havens, Zonnehuis. Private practices and other institutions: Daphne Salek, Gert Vrancken, Hanneke Sistermans and Peggy Penris, Ingrid Veltman, Lieve de Witte, Maddy Meijer, Mantelers-Nijssen, Petra Roukens, Suzanne Cornelussen; Afasiecentrum Rotterdam, Afasie trainingscentrum, Stichting Afasietherapie Amsterdam, Zorghotel De Stromen SFG. We are very grateful to Dineke Blom, Sandra Wielaert and Jane van Gelder-Houthuizen, speech-language therapists, for scoring the ANELT. We wish to thank Cynthia Thomson, Albyn Davis, David Howard and Malcolm McNeil for their scientific advises.

Competing interests

None.


Editorial

Software products are increasingly dominating our lives. Their pervasiveness makes it crucial for both economical and social progress that producing software of high quality and at low cost becomes routine. In order to achieve this goal, software engineers need to have at their disposal theories, languages, methods, and tools that derive from both the systematic research of the academic community and the distilled experience of practitioners.

This special issue originates from the 2008 edition of the International Conference on Fundamental Approaches to Software Engineering (FASE), one of the European Joint Conferences on Theory and Practice of Software (ETAPS), held in Budapest, Hungary, from March 31st to April 2nd. The conference proceedings were published by Springer in the Lecture Notes in Computer Science Series (volume 4961). After the conference, the authors of the highest-ranked papers were invited to submit significantly extended papers for inclusion in this special issue; ten of these were accepted and have been included here.

"Deriving Non-Zeno Behaviour Models from Goal Models using ILP", by D. Alrajeh, J. Kramer, A. Russo and S. Uchitel, addresses the problem of constructing behaviour models from declarative goal specifications using a combination of model checking and machine learning techniques.

"From StPowla Processes to SRML Models", by Laura Bocchi, Stephen Gorton and Stephan Reiff-Marganiec, discusses the way business goals captured in the policy language StPowla can be aligned with high-level models of service provision expressed in the language SRML.

"An Algebraic Semantics for MOF", by Artur Boronat and José Meseguer, defines a reflective, algebraic, executable framework for precise metamodeling based on membership equational logic (mel) that supports the Meta-Object Facility (MOF) standard.

"Automating the Transformation-Based Analysis of Visual Languages", by Juan de Lara and Hans Vangheluwe, shows how, given a description of the operational semantics of a domain-specific visual language in the form of graph transformation rules, model-to-model transformations can be automatically generated for its translation into a semantic domain for analysis.

"Consistent Integration of Models Based on Views of Meta Models" by Hartmut Ehrig, Karsten Ehrig, Claudia Ermel and Ulrike Prange, extends the framework of typed attributed graph morphisms to allow modellers to formulate type hierarchies and views of visual languages.

"From a Domain Analysis to the Specification and Detection of Code and Design Smells", by Naouel Moha, Yann-Gäl Guéhéneuc, Anne-Franoise Le Meur, Laurence Duchien, and Alban Tiberghien, introduces a new way of automating the generation of detection algorithms for code and design smells from specifications written using a domain-specific language.

"A Generic Complete Dynamic Logic for Reasoning about Purity and Effects", by Till Mossakowski, Lutz Schröder and Sergey Goncharov, introduces a generic framework for reasoning about purity and effects in which effects are modelled abstractly and axiomatically (using Moggis idea of encapsulation of effects as monads) and a dynamic logic is defined whose logical formulae are pure programs in a strong sense.

Correspondence and offprint requests to: "Reasoning with Graph Constraints", by Fernando Orejas, Hartmut Ehrig and Ulrike Prange, concentrates on the problem of proving the consistency of specifications based on graph constraints as they arise in the areas of visual software modelling and of semi-structured documents or websites.

"Automatic Verification of Java Programs with Dynamic Frames", by Jan Smans, Bart Jacobs, Frank Piessens and Wolfram Schulte, demonstrates that the dynamic frames approach for the verification of objectoriented programs can be integrated into an automatic verifier based on verification condition generation and automated theorem proving.

"Preserving Correctness During Business Process Model Configuration", by Wil M.P. van der Aalst, Marlon Dumas, Florian Gottschalk, Arthur H.M. ter Hofstede, Marcello La Rosa and Jan Mendling, proposes a formal foundation for individualizing configurable business process models incrementally, while preserving correctness both with respect to syntax and behavioral semantics.

We are grateful to 'Formal Aspects of Computing' for allowing us to publish this collection of papers, and John Cooke in particular for his invitation to organise this special issue and help throughout the editorial process. We are also grateful to all Program Committee members of FASE 2008 and reviewers involved in selecting the papers included in this special issue.

José Luiz Fiadeiro1 (FASE 08 Programme Committee co-chair) Paola Inverardi (FASE 08 Programme Committee co-chair) Tom Maibaum (Member of the Editorial Board of FACJ)

Department of Computer Science, University of Leicester, Leicester LE1 7RH, UK , jose@mcs.le.ac.uk


S Tzirti

K Tsiganis

H Varvoglis

Effect of 3 rd -degree gravity harmonics and Earth perturbations on lunar artificial satellite orbits

Keywords: Lunar artificial satellites, 3 rd -degree gravity coefficients, 3 rd -body effect

In a previous work we studied the effects of (I) the J 2 and C 22 terms of the lunar potential and (II) the rotation of the primary on the critical inclination orbits of artificial satellites. Here, we show that, when 3 rd -degree gravity harmonics are taken into account, the long-term orbital behavior and stability are strongly affected, especially for a non-rotating central body, where chaotic or collision orbits dominate the phase space. In the rotating case these phenomena are strongly weakened and the motion is mostly regular. When the averaged effect of the Earth's perturbation is added, chaotic regions appear again for some inclination ranges. These are more important for higher values of semi-major axes. We compute the main families of periodic orbits (POs), which are shown to emanate from the 'frozen eccentricity' and 'critical inclination' solutions of the axisymmetric problem ('J 2 + J 3 '). Although the geometrical properties of the orbits are not preserved, we find that the variations in e, I and g can be quite small, so that they can be of practical importance to mission planning.

Introduction

During the past decades there has been extensive study of satellite orbits around an oblate primary. The consideration of the J 2 term of the spherical harmonic expansion of the potential leads to a critical value of the inclination, equal to 63 o .43, for which the orbit does not precess secularly [START_REF] Garfinkel | The global solution in the problem of the critical inclination[END_REF][START_REF] Hughes | The 'Critical Inclination': Another Look[END_REF]; Coffey, [START_REF] Coffey | The critical inclination in artificial satellite theory[END_REF][START_REF] Jupp | The Critical Inclination problem -30 years of progress[END_REF]). For the case of the Earth, [START_REF] Cutting | Orbit Analysis for SEASAT-A[END_REF] examined the conditions under which the long-periodic argument of pericentre and eccentricity remain both constant with the additional effect of the J 3 term (frozen orbits). [START_REF] Coffey | Frozen Orbits close to an Earth-like planet[END_REF] searched for frozen orbits close to an Earth-like planet, including up to 6 th -degree zonal coefficients, while [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF] took into consideration all the available zonal coefficients for the case of a polar lunar orbiter. The tesseral and sectorial coefficients were, in all cases, neglected, as their values were considered too small to affect the long term orbital behavior or the critical inclination value.

The situation for the Moon is different than for the Earth, as the C 22 coefficient is only nine times smaller than J 2 . This problem was first examined by De Saedeleer and Henrard (2006), who showed that the value of the critical inclination depends on the initial phase of the ascending node. In their derivation the rotation of the Moon and the time variation of the averaged inclination were neglected, most likely both considered to be insignificant on the assumed time scale. In a previous work [START_REF] Tzirti | Quasi-critical orbits for artificial lunar satellites[END_REF]) we showed that, under the combined effect of J 2 and C 22 harmonics, the averaged inclination does not remain constant but performs long periodic oscillations. The time evolution of the inclination causes variations in the argument of pericentre, of similar amplitude. Under those circumstances, the term 'critical inclination' becomes meaningless. Instead, we used the term 'quasi-critical' for the solutions in which the argument of pericentre librates. On the other hand we showed that the strong dependence of the 'critical inclination' on the longitude of the ascending node, found by De Saedeleer and Henrard (2006), is considerably weakened when the rotation of the Moon is taken into account.

Besides C 22 , even the higher-degree harmonics are important for the case of Moon. In this work we take into account perturbations up to 3 rd -degree. So, we include the 'J 3 ' effect and the non-axisymmetric part as well. The averaged system has now two degrees of freedom and the long-term stability of the orbits is modified. As we show below, for a non-rotating Moon, extended regions of chaotic motion or regions of collision orbits (due to the large amplitude oscillations of the eccentricity) occupy a large portion of phase space. However, the rotation of the Moon smooths out these phenomena. The motion appears to be mostly regular for all semi-major axes and orientations studied. We compute the main families of stable periodic orbits, which are shown to emanate from the POs of the 'J 2 + J 3 ' axisymmetric case. There are two families, one corresponding to the low-e 'frozen' orbits (∀I) and the other to the large-e 'critical inclination' solution (at I ≈ 63 o .5).

For distances greater than 3R (where R is the mean equatorial radius of the Moon), it is necessary to consider the effect of the Earth (as e.g. in De Saedeleer 2006b), which changes drastically the orbital behavior, especially at high inclinations. We find that in some regions of phase space the motion seems to be chaotic or to lead to collision with the surface of the Moon. For relatively large semi-major axis values, the effect of the Earth becomes dominant, while for orbits close to the Moon, the 3 rd -degree gravitational harmonics seem to play the decisive role. In both cases we examine the distribution of the main families of POs (in the averaged problem).

Hamiltonian of the problem and equations of motion

Gravity harmonics of the Moon

We consider an artificial satellite in orbit around the Moon. We use a rotating frame whose origin is at the center of the Moon, the x axis passes through the longest lunar

V = - µ r ∞ n=0 R r n n m=0 Pnm(sin φ) [Cnm cos mλ + Snm sin mλ] (1)
where µ = GM, R is the mean equatorial radius of the Moon, (λ, φ) are the selenographic longitude and latitude, respectively, r the selenocentric distance of the satellite, Cnm and Snm the non-normalized gravity coefficients and Pnm the Legendre polynomials of degree n and order m. For m = 0 we get the zonal harmonic coefficients, for which we will use the notation Jn = -C n0 . The numerical values used for the constants are G = 6.6739 × 10 -11 m 3 kg -1 s -2 , M = 7.349 × 10 22 kg and R = 1737530 m. To simplify the calculations, we use dimensionless units, by taking µ = 1 and ares 1-1 = 1, where ares 1-1 (= 88459747.12676102 km) represents the semi-major axis value for which the period T of the satellite is in 1 -1 resonance with the rotational period of the Moon in the unperturbed problem. In this system of units T = 2π. The numerical results presented in this paper follow the above unit conventions.

In this paper we consider the effect of terms up to 3 1.

The Hamiltonian of the problem has the following form:

H = H 0 + H J2 + H C22 + H J3 + 3 j=1 H C3j + 3 j=1 H S3j + Hn M (2)
where H 0 = u 2 /2 -µ/r is the Keplerian part, H J 3 , H C 3j and H S 3j correspond to the 3 rd -degree gravity harmonics and Hn M = -n M p λ (p λ is the momentum conjugate to λ) describes the rotation of the Moon. In selenographic coordinates, H will contain powers of the distance r, powers of sin φ, coming from the Legendre polynomials Pnm, and cos mλ, sin mλ, which are also easily expressed as powers of cos λ and sin λ. The angles φ, λ and the distance r can be written as functions of the semi-major axis, a, the eccentricity, e, the inclination, I, of the orbital plane relative to the equatorial one, the longitude of the ascending node in the rotating system, i.e. h = Ω -n M t, the argument of pericentre, g, and the true anomaly, f , using the relations (7.10) from De Saedeleer (2006a):

sin φ = sY (3)
X = cos φ cos (λ -h) (4)
cos φ = cY sin (λ -h) + X cos (λ -h) (5)
where s = sin I, c = cos I, Y = sin (f + g), X = cos (f + g).

After performing the trigonometric expansions in (4) and ( 5), we can solve the first of them for sin λ and replace it in the second one, which can now be solved for cos λ. Using the relations

X 2 = 1 -Y 2 , s 2 = 1 -c 2 and cos φ = √ 1 -s 2 Y 2
where it is essential, we can derive the simplified expressions for cos λ and sin λ as functions of the orbital elements Table 2 3 rd degree terms that appear in the Hamiltonian. K = 3µeR 3 / 16a 4 1 -e 2 5/2 , c = cos I, s = sin

I parameter K ± 3j K ′± 3j C 31 K(1 ± 11c -5c 2 ∓ 15c 3 ) 0 C 32 0 10K(-s ± 2cs + 3c 2 s) C 33 30K(1 ± c -c 2 ∓ c 3 ) 0 S 31 0 ±K(1 ± 11c -5c 2 ∓ 15c 3 ) S 32 ∓10K(-s ± 2cs + 3c 2 s) 0 S 33 0 ±30K(1 ± c -c 2 ∓ c 3 ) cos λ = X cos h -cY sin h √ 1 -s 2 Y 2 (6) sin λ = X sin h + cY cos h √ 1 -s 2 Y 2 (7)
Substituting cos λ and sin λ in H, the latter becomes a function of [a, e, I, f, h, g].

Following Roy (1982) we can derive the averaged Hamiltonian to first order in (Jn, Cnm, Snm), by integrating H over the mean anomaly, using standard trigonometric formulas and applying the following change of variables:

dM = r 2 a 2 (1 -e 2 ) 1/2 df (8)
The steps described above can be applied to obtain the averaged Hamiltonian for up to every term of degree n and order m in the perturbing potential. In our case, the averaged terms of ( 2) are given below:

H 0 = - µ 2a , Hn M = -n M µa(1 -e 2 )c (9)
H J2 = µJ 2 R 2 1 -3c 2 4a 3 (1 -e 2 ) 3/2 , H C22 = - 3µC 22 R 2 s 2 2a 3 (1 -e 2 ) 3/2 cos 2h (10)
H J3 = 2J 3 K -4 + 5s 2 s sin g (11)
H C3j = C 3j K + 3j cos (g + jh) + K - 3j cos (g -jh)+ K ′+ 3j sin (g + jh) + K ′- 3j sin (g -jh) (12)
H S3j = S 3j K + 3j cos (g + jh) + K - 3j cos (g -jh)+ K ′+ 3j sin (g + jh) + K ′- 3j sin (g -jh) (13)
where K ± 3j and K ′± 3j are given in Table 2.

Before proceeding to the inclusion of a 3 rd -body perturbation and analyzing the methods of work, it would be of interest to have a look at the solutions of the axisymmetric problem (the averaged Hamiltonian contains only the terms J 2 and J 3 ). In that way, we will be able to examine how the asymmetry introduced by the additional sectorial and tesseral terms affects the phase space and POs. Following [START_REF] Vallado | Fundamentals of Astrodynamics and Applications[END_REF] ġ and ė will be:

ġ = - 3n M J 2 (1 -5c 2 ) 4(1 -e 2 ) 2 R a 2 N (e, I, g), (14)
where

N (e, I, g) = 1 + J 3 2J 2 R a 1 1 -e 2 s 2 -e 2 c 2 s sin g e (15)
and

ė = 3J 3 8 n M (1 -e 2 ) 2 R a 3 s 1 -5c 2 cos g (16)
One can reach the same results using Eqs. [START_REF] Murray | Solar System Dynamics[END_REF] in the next page. A periodic solution ( ġ = 0, ė = 0) is obtained in two cases: (a) for I = 63 • .43 = const. ('critical inclination' orbits), so that Eqs. ( 14) and ( 16) become both equal to zero. This solution occurs for every value of e 0 and g 0 and it already exists in the J 2 problem (b) for g = ±π/2, which satisfies ė = 0 and appropriate values of e 0 , I 0 , so that N (e, I, ±π/2) = 0 or equally ġ = 0 ('frozen eccentricity' orbits). For I = 0 or e = 0 a singularity appears in Eq. ( 15), which does not allow us to give initial conditions exactly equal to these values.

Frozen satellite orbits have been also numerically computed by [START_REF] Lara | Numerical continuation of families of frozen orbits in the zonal problem of artificial satellite theory[END_REF], who included in their model zonal coefficients up to order 9, while Abad et al. ( 2009) developed an analytic model for locating frozen orbits keeping zonal harmonics J 2 and J 7 . Other work in the same direction is that of [START_REF] Elipe | Frozen orbits in the tesseral artificial satellite theory[END_REF], in which, starting from results on frozen orbits in the zonal problem, they examine their continuation under the inclusion of the C 22 term. Moreover, [START_REF] Lara | Preliminary design of low lunar orbits[END_REF] computed frozen orbits in the zonal problem and they compared them to results from analytic solutions.

We now turn to the computation of the averaged perturbation exerted by the Earth.

The effect of the Earth as a 3 rd -body perturbation

As we go further away from the surface of the Moon, the effect of the Earth, as a 3 rdbody perturbation, cannot be neglected. To include the interaction with the Earth, we work in a rotating frame with the Moon at its origin (Fig. 2). The perturbing potential arising from Earth has the form (see e.g. Murray and Dermott 1999) where µ ′ = GM Earth (M Earth = 5.9736 x 10 24 kg) and r, r ′ are the position vectors of the satellite and the Earth relative to the Moon. V Earth can be expressed as a series of Legendre Polynomials P k

V Earth = - µ ′ |r -r ′ | + µ ′ r • r ′ r ′3 (17
V Earth = - µ ′ d ∞ k=2 r d k P k (cos φ) (18)
from which we will only keep the k = 2 term. We assume circular motion of the Earth around the Moon at the mean distance d = r ′ = a 0 1 + e 2 0 /2 , where a 0 = 384, 400 km and e 0 = 0.0549, and that the orbital plane is the Moon's equatorial plane. In this way the selenographic latitude, φ, is written as a function of the angles f, g, h by using the first term (in our case the other two terms are considered equal to zero because of the previous assumptions) of Eq. (8.10) from De Saedeleer (2006a).

cos φ = X cos h -cY sin h (19)
After expanding the trigonometric sums in X (= cos (f + g)) and Y (= sin (f + g)), using standard formulas for expressing cos f , sin f and r as functions of the eccentric anomaly (E) and changing consistently the variable from M to E, we obtain the first order averaged Hamiltonian by integrating over E. The choice of working with E simplifies the calculations. The final result is

H Earth = µ ′ a 2 16d 3 1 -3c 2 2 + 3e 2 -3s 2 5e 2 cos (2g) + 2 + 3e 2 cos (2h) - 15e 2 2 (-1 + c) 2 cos [2 (g -h)] -(1 + c) 2 cos [2 (g + h)] (20)
and it formally agrees with that of De Saedeleer (2006b).

In order to find the equations of motion, either with or without the inclusion of Earth, we use the well-known equations

ġ = ∂ H/∂G, Ġ = -∂ H/∂g (21) ḣ = ∂ H/∂H, Ḣ = -∂ H/∂h (22)
where [L, G, H] are the Delaunay actions.

L = √ µa, G = µa (1 -e 2 ), H = µa (1 -e 2 ) cos I (23)
In the following, we focus our study on three models: the '3G' model, which contains only up to 3 rd -degree gravity harmonics, the '3G+R' model, which contains the additional effect of the rotation of the Moon and, finally, the more complete '3G+R+E' model that also introduces the effect of the Earth as a 3 rd -body perturbation. In all models, the dynamical system has two degrees of freedom. To study the global behavior of orbits, we will use Poincaré maps, as well as Fast Lyapunov Indicator (FLI) maps [START_REF] Froeschlé | The fast Lyapunov indicator: a simple tool to detect weak chaos. Application to the structure of the main asteroidal belt[END_REF]. Finally, we will compute the main families of POs and analyze their geometrical and dynamical properties.

Poincaré sections and FLI maps

The Poincaré maps presented here are taken on the G -g plane, defined from the section h = π and ḣ < 0. Note that G is dimensionless in the system of units described above. The choice h = π is just indicative. A different choice of h for the section plane (i.e. h = 0) would give similar results to the h = π. The section is defined for ḣ < 0, because h decreases as time increases. The choice of the opposite sign would give no points on the section. For given initial conditions (a, g 0 , G 0 , h 0 ) and 'energy' values, we use a Newton-Raphson method to determine H 0 . The equations of motion are solved using a 4 th -order Runge-Kutta scheme with step equal to T /10, where T is the period of the satellite in the unperturbed problem (the order of magnitude of T is several hours for the semi-major axis values considered here. For example, T = 2.23 hours for a = 2000 km and T = 6.3 hours for a = 4000 km). To achieve better accuracy on the crossing point of the section surface, we first find a point immediately after the crossing, we change the independent variable to h (the equations of motion now are t

′ (h) = dt/dh = 1/ ḣ, g ′ (h) = ġ/ ḣ, H ′ (h) = Ḣ/ ḣ, G ′ (h) = Ġ/ ḣ)
and integrate backwards until h = π with accuracy 10 -12 . As the energy level and h are constant on each section, each curve on the map corresponds to a different inclination range.

The computation of FLI allows us to discriminate ordered and chaotic motion or even resonant and non resonant regular orbits on a relatively short time interval. Here we use the revisited definition given in Lega and Froéschle (2001). FLI maps are not useful for integrable systems, since the motion is regular in the whole phase space. For given initial a, g 0 and h 0 , we chose a 50 × 90 grid in (e 0 , I 0 ) in the range 0 < e < ec and 0 < I < 180 • , where ec is the collision eccentricity. For each set of initial conditions, FLI is equal to the logarithm of the maximum value that the norm of the tangent vector reaches up to time t. As the evolution of the system is slowed when we increase the semi-major axis value, we selected to integrate for t = 1500 lunar months. This time scale is quite large and it exceeds the typical duration of a lunar mission, but it is essential for semi-major axis values until a = 12R, for which the characteristics of the orbits (for example chaotic behavior) are exhibited after some hundreds of lunar months. Regular orbits close to resonance (F LI ≤ 5) are represented by black color on the FLI maps, while the white colored regions indicate chaotic or collision regions. Below we discuss the results for each model separately.

'3G' model: As it was shown in [START_REF] Tzirti | Quasi-critical orbits for artificial lunar satellites[END_REF], under the combined effect of J 2 and C 22 terms of the lunar potential the orbits remain regular. This is because the system is effectively reduced to one degree of freedom, since in this case G is an integral of motion. The inclusion of the 3 rd -degree terms, however, seems to change significantly the secular orbital behavior. Poincaré sections and FLI maps contain extended chaotic regions, encompassing regions where regular orbits have secular oscillation amplitude in e large enough to lead to collision on the surface of the Moon (collision regions -see Figs. 3,4a). The collision limits are indicated on the Poincaré sections by horizontal, solid lines. The semi-major axis values used in this model are a = (2 × 10 3 , 3 × 10 3 , 4×10 3 ) km. For a = 2×10 3 km, chaotic and collision regions dominate the phase space to such an extent, that regular regions are hardly detectable (that is why this is not shown). It should be noted that points in Fig. 3b, that seem to form an unusual series crossing the other lines, belong to different curves. The 3 rd -degree terms introduce an asymmetry relative to the I = 90 • plane (Fig. 4a). This practically means that two orbits with initial inclinations I and 180 • -I do not exhibit the same behavior (see Fig. 4a).

'3G+R' model: In this case the motion seems to be mostly regular for all semi-major axes studied, as the effect of the lunar rotation smooths out the chaotic behavior. Terms depending on h are now dominated by rotation, so the asymmetry with respect to I = 90 • is strongly weakened. Poincaré maps show the existence of stable POs at small eccentricities (< 0.01) for any value of inclination (Fig. 5). These correspond to the 'frozen-e' solution of the axisymmetric problem. For inclinations close to 63 • .5 or 116 • .5, another family of eccentric POs seems to bifurcate from the 'critical inclination' solution of the 'J 2 + J 3 ' problem.

'3G+R+E' model: We studied orbits with semi-major axes in the range R + 100 km ≤ a ≤ 20000 km under the combined effects of 3 rd -degree gravity harmonics and the perturbation of Earth. The characteristics of the motion are different close to and far from the surface of the Moon, where the lunar gravity harmonics or the 3 rd -body effects are dominant, respectively. In Fig. 6 we present the Poincaré sections for a = 3000 km and a = 10000 km, for three different inclination ranges (again it is chosen ho = π for the surfaces of section. This value is just indicative. The same analysis performed for other values of ho gives similar results). For a = 3000 km and I < 58 • , the motion appears to be regular and there are two main islands where motion does not lead to collision, centered around g 0 ≈ ±π/2 (Fig. 6a). For 58 • < I < 78 • , the orbits are either chaotic or lead to collision, due to regular secular oscillations (Fig. 6c). The main islands do not anymore correspond to orbits with acceptable eccentricity values, so there are no POs corresponding to 'collision-free' motion (Fig. 4b). For inclinations close to 90 • (near-polar orbits), the phase space mainly consists of regular orbits and two main islands of 'collision-free' motion appear again, but now at g 0 ≈ 0, π (Fig. 6e). For a > 3000 km the structure of phase space is almost the same as that for a < 3000 km only for I < 60 • (Fig. 6b). As the inclination increases, the main islands disappear and the phase space is dominated by chaotic regions that become more extensive as the inclination approaches 90 • (Figs. 6d,6f). The distribution of POs in phase space is essential for understanding the dynamics of a nearly integrable system, as they form the 'backbone' of phase space. From a practical point of view, POs are essential for efficient design of lunar orbiter missions, since, as a rule, minimal spacecraft control is required in the vicinity of that kind of orbits (see [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF].

In order to determine the initial conditions of the POs, we use differential corrections [START_REF] Deprit | Natural Families of Periodic Orbits[END_REF]. For a PO, we have that g(t, T ) = g(t), G(t + T ) = G(t) and h(t + T ) = h(t), where t is the time, T is the period of the orbit and g, h ∈ [0, 2π]. Starting from initial conditions (g 0 , G 0 ), h = h Section , H = f (g 0 , G 0 , h 0 , H), we find the next section point in time t * with coordinates (g(t * ), G(t * )). This method is based on the idea of consecutive corrections of the initial conditions, in order to ensure that at time t * the orbit returns at the same point on the section. This process is repeated, until succeeding the setting accuracy. In our calculations tolerance was set equal to 10 -12 . The method converges if we start close enough to the PO (i.e. in the region of an island).

We calculated the main, stable POs for a variety of semi-major axis values, with emphasis on the '3G+R' and '3G+R+E' models, which are more realistic. The description and the diagrams that follow, refer to a = 3000 km and h 0 = π. '3G' model: Stable POs appear in regions of regular motion and are not symmetrically distributed with respect to I = 90 • (see Fig. 7). The eccentricity of these POs is quite large, almost for all inclinations. Since the model has 2-degrees of freedom, the quantities e, g and I of POs are not expected to be constant. Instead they perform librations with amplitudes ∆e, ∆g and ∆I (defined here as the difference between the largest and the smallest value) 0.08 < ∆e < 0.28 and 5 • < ∆I < 35 • , while g of the POs is circulating for this model almost for all inclinations. The libration amplitudes decrease as the semi-major axis of the orbit increases. The distribution of POs is 'broken' close to I = 60 • , 90 • and 125 • , where the ratio fg/f h (fg and f h represent the frequencies of the angles g and h, respectively) changes resonance.

'3G+R' model: The eccentricity of the POs remains close to zero for all inclinations, except for I ≈ 63 • .5 and 116 • .5, where two branches of POs emanate from the main families of the axisymmetric model (Fig. 7). By varying the energy, one of them moves to lower eccentricities and the other to higher ones, both keeping I almost constant. The nearly-circular POs correspond to the island appearing at large G (small e) in the Poincaré sections (see Fig. 5). The high eccentricity POs correspond to the island that appears at I ≈ 63 • .5, i.e. correspond to the 'critical' inclination orbits of the J 2 problem. The effect of the rotation is to reduce the libration amplitudes, so that ∆e < 0.003, ∆I < 0 • .4. Also, ∆g is smaller than 20 • almost for all inclinations (libration), but it increases as I tends to 0 • or 180 • , where g starts to perform rotations (see Fig. 8). We studied the dependence of the POs on a for R + 100 km ≤ a ≤ 4000 km with ∆a = 100 km, or even smaller in some intervals. Some of these results are presented on an e -I diagram (Fig. 11a). The results for all the cases studied are not included for better readability of the diagram. '3G+R+E' model: In this case, the most complete model in our study, there appear two 'gaps' with no POs (see Fig. 7 and4b). This is because near these gaps the eccentricity of a PO exceeds the collision value, ec (see Fig. 6d) and thus we do not mark it in the diagram. In the other regions, the eccentricity remains very close to zero (e < 0.001). As far as libration amplitudes are concerned, we find ∆g < 10 • , ∆e < 0.01 and ∆I < 0 • .5 (Fig. 8), slightly increasing for larger values of a. It is worth mentioning that for a = R + 100 km, ∆e is such that the minimum distance between the satellite and the surface of the Moon for polar orbits equals to 50 km, (I ≈ 90 • ). Again, we find the POs for different values of a, but now we include the cases a = 5000, 10000, 15000, 20000 km. Some of these results are presented in Fig. 11b, while others are not included for better readability of the diagram. Nevertheless, cases for a > 5000 km have similar behavior as those for a = 5000 km. For small values of a (for example a = R + 500 km), the 'gap' of POs around I ≈ 63 • .5 is narrow and it widens as a increases, an effect attributed to perturbations by Earth. Moreover, for a ≤ 2900 km we find polar POs at small eccentricities, while for larger semi-major axes the eccentricity of the POs around I = 90 • has higher values. The time evolution of some POs for a = 3, 000 km and a = R + 100 km is presented in Figs. 9. and 10, respectively (the orbits on the right are under the additional effect of the Earth, while these on the left are not). It is quite clear that for a values around 3, 000 km the effect of Earth modifies the shape of the orbits, something that is not true for low semi-major axis values. In all cases, the libration amplitude of e is smaller than 0.02.

As mentioned earlier in this paper, [START_REF] Abad | Analytical model to find frozen orbits for a lunar orbiter[END_REF] have done similar work, for different semi-major axis values, taking into account coefficients J 2 and J 7 . The problem they studied is of 1-degree of freedom, so they searched for frozen orbits (i.e. e and g remain constant with time). In our case there is an additional degree of freedom because of the tesseral and sectorial coefficients, so we study POs instead of frozen orbits, the orbital elements of which perform oscillations. Fig. 9 Time evolution of some POs, for a = 3000 km and h 0 = π, for the '3G+R' (left) and the '3G+R+E' models.

Conclusions

We have shown that when 3 rd -degree terms are included in the lunar potential, the secular behavior of the orbits changes significantly, as compared with the 'J 2 +C 22 (+n M )' models, for which the motion, as a rule is regular. When the rotation of the Moon is not taken into account, the phase space contains extensive chaotic or collision regions. The inclusion of the lunar rotation weakens these phenomena, so that the motion appears to be mostly regular. When the effect of Earth is added, chaotic or collision regions appear again in phase space. There is an obviously different behavior between a < 3000 km, where the 3 rd -degree lunar gravity terms are dominant, and a > 3000 km, where the effect of Earth is more important. Note that we have selected the value 3000 km after extensive numerical computations, as an easy to remember approximate value. The behavior of the orbits changes close to that semi-major axis value and not exactly at that value.

The main families of POs emanate from the families of POs of the axisymmetric problem ('J 2 + J 3 '). These are the equivalent of nearly-circular 'frozen-eccentricity' orbits (∀I) and 'critical inclination' orbits (I ≈ 63 • .5 and 116 • .5). In the '3G' model, POs usually appear at large eccentricities, while in the '3G+R' model there exist elliptic POs for every I. For I close to 63 • .5 new POs appear for higher eccentricities. When Earth is included, there are (main, stable) POs for every inclination, except for a gap of chaotic/collision orbits. G and I are no longer integrals of motion, but their variations (∆I, ∆e, ∆g) are small, when n M = 0. These POs can be used as the basis for an efficient design of artificial satellite orbits, requiring minimal control, as already discussed in [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF].

In the future, we intend to extend this work to the numerical computation of POs in the full gravitational potential of the Moon. To this end, it would be desirable to have a closed form relation (in e and I) for the 1 st -order averaged Hamiltonian (De Saedeleer has done this for the zonal problem, 2005). As the computation of the POs will become quite cumbersome, it would be desirable to have an analytic (perturbative) method that would be able to provide an accurate determination of the POs location. This work is under way and we hope to report on these issues in a forthcoming paper.

Fig. 1

 1 Fig. 1 Orientation of the orbit in space

Fig. 2

 2 Fig.2The position vectors of the satellite and the Earth with respect to the Moon

Fig. 3 Fig. 4 Fig. 5 Fig. 6

 3456 Fig. 3 Poincaré sections of the '3G' model. The horizontal, solid lines indicate the border under which the satellite collides with the surface of the Moon

Fig. 7 e

 7 Fig. 7 e -I diagram of the main, stable POs for the three models studied (a = 3000 km)

Fig. 10

 10 Fig.10Time evolution of some POs, for a = R + 100 km and h 0 = π, for the '3G+R' (left) and the '3G+R+E' models.

Fig. 11 e

 11 Fig. 11 e -I diagrams of the main POs without (a) and with (b) the effect of the Earth for various semi-major axis values

  rd -degree, omitting C 21 , S 21 and S 22 , whose values are at least three orders of magnitude smaller than C 22 and have considerable error bars (see Table1or De Saedeleer 2006a -Tab. 1.5. Note, however, that the values come from different sources: Bills 1980 for 2 nd -degree and Konopliv 1998 for 3 rd -degree). The inclusion of these terms would complicate the equations, while it would give negligible improvement to the results. It is already obvious that the rd -degree terms are J 3 and C 31 . In the following, we will use the parameter values given in Table

	Table 1 Gravitational parameters values and uncertainties (De Saedeleer 2006a)
	parameter	(value ± uncertainty) x 10 -6
	J 2	202.43 ±1.14
	C 21	-0.0035 ±0.0058
	C 22	22.26 ±0.13
	S 21	-0.00098 ±0.0045
	S 22	0.0108 ±0.0020
	J 3	8.476 ±0.017
	C 31	28.437 ±0.0081
	C 32	4.846 ±0.002
	C 33	1.7132 ±0.0006
	S 31	5.9022 ±0.0048
	S 32	1.6709 ±0.0021
	S 33	-0.249 ±0.0006
	leading 3