The analysis system for astrophysical plasmas (ASAP) of the Osservatorio Astronomico di Palermo

Computer Physics Communications ELSEVIER

Computer Physics Communications 81 (1994) 105-119

The analysis system for astrophysical plasmas (ASAP) of the Osservatorio Astronomico di Palermo A. Maggio *, F. Reale, G. Peres, A. Ciaravella lstituto ed Osservatorio Astronomico, Piazza del Parlamento 1, 90134 Palermo, Italy Received 20 January 1994

Abstract

We present the software package ASAP, mainly devoted to the presentation and analysis of fluid models of astrophysical plasmas. ASAP allows to generate optically thin X-ray and UV spectra emitted by plasma of given density, temperature and velocity, by means of well established spectral synthesis codes; the emission can then be easily folded with the instrumental response of interest and compared directly with measurements, or can be used to predict the performance of planned instruments. As examples of typical applications of the ASAP package, we show high resolution emission spectra derived from 1D static models of coronal loops, detailed fitting with coronal loop models of observations made with the ROSAT/PSPC instrument, hydrodynamic simulations of a solar flaring loop as seen by the Yohkoh satellite, and the evolution of line emission from a dense blob of interstellar plasma subject to the passage of a shock. The package, designed to be highly flexible, modular, and easily expandable, is expected to evolve and grow rapidly in response to new needs and required tasks. It operates within the IDL environment, and it is under consideration for being included in the SOHO Archive System.

1. Introduction

The development of the analysis system described in this p a p e r stems from the active inv o l v e m e n t of the group at the Osservatorio Astronomico di Palermo in the modeling, and comparison with observations of astrophysical environments as diverse as solar and stellar coronae, supernova remnants, and accretion flows in galactic halos. To this aim, our group has developed, in collaboration with other researchers, various

* Corresponding author.

numerical models of astrophysical plasmas, such as a one-dimensional hydrostatic model of coronal plasma confined in a loop [37], the P a l e r m o Harvard hydrodynamic plasma model [22], a hydro-code suited to study winds and cooling flows [25], the Palermo two-dimensional hydrodynamic plasma model based on the F C T technique [30], and - more recently - a siphon flow model of coronal loops [19]. The approach is to derive with these models the evolution of distributions of plasma temperature, density, and velocity, which in turn allow to compute the source emission according to an appropriate radiation model, and to compare the model results with observations, usually at X-ray or U V wavelengths.

0010-4655/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0 0 1 0 - 4 6 5 5 ( 9 4 ) 0 0 0 1 4 - S

106

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

This approach, and in particular the need of performing detailed spectral synthesis, has already been called for by the spatially-resolved wide-band or spectroscopic observations of the solar corona obtained in the seventies with the S-054 and S-055 experiments on-board of the Skylab [42,33], and later with SMM [1,20]. A fortiori, this occurs today because of the spectroscopic and imaging capabilities of the newly operative detectors on-board the solar observatories Yohkoh [18,7], NIXT [10,11], and the near at hand SOHO [9], as well as on space missions dedicated to observe galactic and extragalactic sources in EUV and X-rays, such as EUVE [4], ROSAT [40] and ASTRO-D [39]. The analysis system for astrophysical plasmas (ASAP) has been developed mostly to achieve three goals: presentation of model results, spectral synthesis, and direct comparison with observational data. The first task is intended to favor physical insight and to help the researcher to disentangle the interplay of the several physical effects present in the plasma model (e.g. radiative losses, hydrodynamics, heat conduction, energy and momentum inputs, interaction of non-thermal particles with thermal plasma). Some of such effects are sometimes equally important, at work simultaneously, and their functional form can be strongly non-linear. Therefore, the models are complex and the interpretation of their results non-obvious. The difficulties increase when the model is multi-dimensional, or time-dependent. The second task of ASAP, namely the spectral synthesis, is simple in principle, because in most cases the plasma is optically thin to EUV and X-ray radiation, so that source emission can be computed by integrating all contributions along the line of sight. Nevertheless, the task of computing the spectral properties of the plasma emission may need a detailed and sophisticated emissivity model and, sometimes, also taking into account effects such as non-equilibrium of ionization. The last step, i.e. the comparison with observational data, requires a detailed knowledge of the instrument characteristics, which are often subject to several revisions during the lifetime of the experiment. Our approach is to compute the

model emission with the proper spatial and spectral resolution for a direct comparison with the available measurements. We have preferred this approach, instead of trying a deconvolution of the observations with respect to the instrument response, because in the latter case we are often faced with an ill-defined problem from a mathematical point of view, which may admit nonunique solutions [8,14]. Our interest in the above fields has brought, over the years, the development of an extensive library of programs and routines for presentation and analysis of model results, which reflect in their diversity the development of the software resources, the model sophistication, and the growing information content of the observational data. We have recently realized that an integrated system with a rapid and flexible interface with model results and observational data, would allow a qualitative leap of our efficiency to address current research problems. Therefore, we started the effort to develop ASAP, reorganizing the existing library of routines, eliminating redundancy and incompatibility, and building up new facilities required by typical modern scientific projects. Presently, we plan to manage the distribution of ASAP to other selected sites. In particular, ASAP is under consideration for being available as a subpackage of the SOHO archive system, whose strategy and general features are currently being addressed by the SOHO project scientists office (Martens 1993, private communication). In this paper we present an overview of this system; more specifically, in Section 2 we describe the rationale and the strategy adopted, we present in Section 3 its structure and capabilities, and in Section 4 some application examples; we discuss future prospects in Section 5.

2. Rationale and development strategy

ASAP has been designed with the aim of being flexible, simple and effective. These characteristics stem naturally from two requirements: (a) the need to interface with minimal effort different existing or feasible plasma models, usually corn-

A. Maggio et al. /Computer Physics Communications 81 (1994) 105-119

puted by Fortran codes, and (b) the choice to build a system easy to expand and update, i.e. able to grow in time, where any user could easily develop new customized routines, analysis procedures, or presentation schemes, which may be eventually included in ASAP, if useful. Following the above guidelines, ASAP has been developed as a library of specialized routines, within a high-level programming environment, namely IDL (Interactive Data language). Although we have tried other environments, IDL has resulted as the most suitable base system for our purposes, since it is a commercially available, widely diffused software environment, and allows numerical computations, interactive analysis and presentation of scientific data through graphs,

107

images, or combinations of the two. We expect that IDL will follow the hardware development for several years to come, and this guarantees that ASAP can be easily moved to the next generation of hardware platforms. An additional advantage of using the package within a tested analysis environment is the ability of the user to develop by himself complex procedures, which might, eventually, become themselves part of the package. The present version of ASAP is based on a core of routines widely used within our group. These routines are highly modular, and allow different combinations of analysis and display of results, as shall be evident in the examples in Section 4.

PRESENTATION A s t r o p h y s iMcoadle l s

- Model -

Temperature, velocity vs time

Spectral

density, space and

->

Codes

-

Raymond-Smith

-

Landini&Monslgnorl

-

High resolution synthesis of specific lines

-

source

spectra

-

emission distributions

-

-

Synthesized

spectra and emission

Folding

spatial distributions

synthetic

spectra

Fossl

- Plasma /

-

results

Instrument

Synthesis

~Isslvlty

Instrument

of

Plasma

emissivities

responses

EINSTEIN/IPC ROSAT/PSPC

- ROSAT/HRI - ROSAT/WFC - EUVE spectrometers - ASTROD/SIS - SMM/FCS

-

/

Instrument characteristics

- YOHKOH/SXT - NIXT

ICampmrlson withJ ~]

experiments

MEASUREMENTS

Fig. 1. Block diagram of ASAP in its present configuration. The arrows indicate the information flow needed to perform any of the ASAP tasks. Intermediate results can be scrutinized and displayed at any stage of the analysis.

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

108

3. Structure and present capabilities The ASAP concept is schematically explained by the block diagram in Fig. 1. The arrows indicate the information flow needed to perform any of the ASAP tasks. Starting points for the analysis are data arrays of distributions of plasma quantities computed with models, emissivity functions, and instrument characteristics (see Table 1 for more details). These data are employed by the ASAP routines in several ways, to produce other data arrays as a result, and to proceed further with the analysis. Partial results can be graphically displayed at any stage. The final goal is to produce model data which can be directly compared with real observations from space or even laboratory experiments. Note that we can consider also oversimplified plasma models (e.g. isothermal, isodensity plasma distributions), and

simplified (flat) instrument responses, in order to get different end products: for example, it is possible to obtain an instrument-folded emissivity table (cf. Section 3.6) assuming a trivial plasma model, or we may be interested to study the model emission distribution incident on the telescope, considering a boxcar instrument response. In the following we illustrate in some detail the main routines available, and how they achieve their tasks.

3.1. The standard ASAP data format We have developed a flexible I / O interface between ASAP and various disk-resident data arrays, adopting a single file format. The most general case is represented by the plasma models, which provide data arrays with the distribution of several plasma physical quantities, either time-de-

Table 1 ASAP structure and content Category

Description

Type

Plasma models

1D hydrostatic coronal loop models either with Spitzer or with limited free-streaming thermal conduction [37,6] Palermo-Harvard 1D hydrodynamic loop model [22] Palermo 2D FCT hydrodynamic model [30] 1D siphon flow loop model [19]

external programs

Emissivity models

Raymond-Smith model [28,29] Landini and Monsignori Fossi [16] Model

external program external tables

Instrument responses

Einstein/IPC effective area [15] ROSAT/PSPC effective area and response matrix a ROSAT/HRI effective area a ROSAT/WFC filter effective areas a ASTRO-D/SIS effective area [39] EUVE spectrometers effective areas a SMM/FCS [1] Yohkoh/SXT [41] NIXT [13] AXAF-S [38]

external tables

General data I / O

Reading and writing data files in ASAP standard format

ASAP routines

Spectral synthesis

Spectra emitted by the model plasma distributions ISM absorption

ASAP routines

Instrument folding

Folding with effective area and energy response matrix

ASAP routines

Graphics

Presentation of model results and comparison with measurements

ASAP routines

a The ROSAT and EUVE instrument responses are taken from the IRAF/PROS and I R A F / E U V software packages, respectively. The ROSAT/WFC data are from the WFC Master Calibration File of the ASTERIX software. They are regularly updated according to the most recent releases.

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

pendent or time-independent, c o m p u t e d over spatial domains of one or more dimensions. To manage any of such models, the results are stored in a direct-access binary file whose content and structure is described in a companying ASCII file, called "(data) file descriptor"; the data file and the file descriptor have the same root name and two different suffixes, and both represent the standard ASAP data format. A dictionary of standard "(parameter) keywords" is used to code the required information within the file descriptor. This information allows any user to determine which physical quantities have been computed in the model, whether they are time variable or not, the number of dimensions and the size of the spatial domain, and to retrieve a number of scalar parameter values characteristic of the model. In the phase of reading model results within ASAP, the structure of the data files is transparent to the user: by means of two general ASAP routines, it is sufficient to load the file descriptor into an IDL "structure" and, in turn, the data array. A third routine allows to extract the distribution of any physical quantity over the spatial domain, at any selected time, by providing as an input the data array, the file descriptor, the keyword (a string) of the desired parameter, and an optional record number in case of time-variable parameters.

3.2. Plasma models The numerical models are computed by running Fortran codes which solve the specific equations, and write the results in one or more data files (cf. Section 3.1). ASAP routines can directly manage the execution of programs for models whose computation takes a reasonably short time, and therefore can be computed interactively. Models whose computations can take hours, such as the 2D FCT code, are typically run asynchronously outside ASAP (we have taken care that their output conforms to ASAP standard); the only reason to avoid driving the computation within ASAP is the inconvenience of an extremely long interactive session. The file descriptors may be produced either by the Fortran code itself, or by the ASAP driver. Note that only a

109

standardization of the data output format has been required to allow ASAP to manage any of the available models, whose list is reported in Table 1. The model results consist of the spatial distributions of relevant plasma quantities, namely temperature, density ( a n d / o r pressure), and velocity. A sequence of snapshots of these distributions are present in case of time-dependent models. ASAP is also capable to manage the results of two-fluid plasma models, so to match the ongoing development of this class of models in our team. ASAP allows to present easily any plasma quantity vs. any other, or vs. time, or to display its distribution over any one- or two-dimensional cross section of the spatial domain. This is made first extracting data as in Section 3.1, and then through the built-in visualization facilities of IDL, fully exploitable either interactively or by means of a number of specific graphic routines in the ASAP library (cf. appendix).

3.3. Emissivity models The models of emissivity of an optically thin plasma in thermal equilibrium, typically are functions of the plasma temperature, T, the photon energy, E, (or wavelength), and the element abundances, aE, and include contributions from b o u n d - b o u n d , free, free-bound, free-free, and two-photon processes. The plasma emissivity per unit emission measure, P(T, E, aE), when multiplied by the proton density times the electron density, provides the emission spectrum per unit volume of a plasma at temperature T. The plasma emissivity can be computed in real time after choosing the appropriate energy band, energy resolution, temperature resolution, and plasma abundances, by running external emissivity codes with an ASAP driver - as it is the case with the Raymond thermal model [28,29] - or simply by loading a precomputed emissivity table (e.g. the emissivity tables of Landini and Monsignori Fossi [16]). As mentioned in Section 3.1, all emissivity tables, either precomputed or computed in real time, are available on-line in the same format as the plasma model results, i.e. in

110

A. Maggioet al. /Computer PhysicsCommunications81 (1994) 105-119

the standard ASAP data format (cf. Section 3.1), so they can be loaded with the same transparent I / O procedure described above. In Table 1 we list the emissivity models currently available. We intend to keep them updated, if possible, by incorporating revisions which may be released in the future, and to include new emissivity models if they are made available. At the moment ASAP includes emissivity tables for some X-ray lines or bands observed with specific instruments. Finally, we plan to include into ASAP a series of subroutines to synthesize the spectral region around the resonance line of helioid ions of calcium and iron (Ca XIX 3.18 A, and Fe XXV 1.85 A) at a high level of resolution ( E / A E ~ 6000). The diagnostics based on the high resolution spectroscopy of these ions, in fact, can play a fundamental role in ascertaining the energetics and dynamics of high temperature astrophysical plasma whose features range from those of solar flares to those of AGN.

3.4. Spectral synthesis The synthesis of the spectrum emitted by an optically thin, multi-temperature plasma is performed with the following integration:

S( E) = f no(s)np( s)e( T( s), E)

dV,

where n e is the electron density, np the proton density, P(T, E) is the emissivity over the band of interest, and V is the plasma volume; both ne, n o, and the temperature T are typically functions of the spatial coordinates, s. This integral is computed after interpolation of the P(T, E) function at the temperature values for any point in the model. The function P(T, E) may represent either the total power emitted by the plasma per unit energy, or the line emission in a series of b o u n d - b o u n d transitions. Consequently, the ASAP routine which performs the spectral synthesis returns either the full source spectrum over the energy range defined by the P(T, E) function, or a selected set of line emission fluxes at the source. This routine may optionally return the integrand

o'( E, s) = n~( S)np( s ) P ( T( s), E ) ,

which provides the spatial distribution of plasma emission per unit volume, as a function of energy, computed with the spatial resolution of the model. The display of tr(E, s) is useful per se, but also for a better interpretation of those observations with limited or no spatial resolution. The same routine allows to include a weight function in the integration, a feature useful to evaluate quantities averaged over the plasma volume and over the instrument response, such as the effective temperature of the plasma observed in a particular band (see Section 3.6 for more details).

3.5. ISM absorption To take into account the interstellar absorption, the source emission spectrum must be multiplied by the absorption coefficient F ( E , alSM) = e - ° / ~ ( e ' ~tsM) where lz(E, ais M) is the ISM attenuation length, dependent on the photon energy, E, and the ISM element abundances, a~SM, and D is the source distance. The factor F(E, OqsM) can be currently computed by two ASAP routines: the first one, parameterizes the total absorption cross section for an ISM with cosmic abundances according to Brown and Gould [5] or Morrison and McCammon [17]; the second one - appropriate only to the XUV region below the carbon edge at 0.28 keV - computes the absorption coefficient due to neutral hydrogen, and neutral and ionized helium, according to the prescription of Rumph et al. [35].

3. 6. Instrument folding The emission spectrum incident on the telescope, I(E), is computed dividing the absorbed spectrum by 4TrD 2, where D is the source distance. Then, the synthetic observed spectrum is computed by folding the incident spectrum through the instrument effective area, A(E), itself a function of photon energy (or wavelength). This area is usually tabulated with a spectral resolution higher than that achievable by the detector, and is usually made of two multiplica-

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

tive components, one due to the telescope, and the other due to the detector proper. A third component may be present if the focused photons go through a filter before being detected. Within ASAP we have built an extensive library of instrument effective areas, listed in Table 1, all available in the standard ASAP data format, so to ensure ease of table handling. Other effective areas can be easily included, so that ASAP can be applied to future instruments. An ASAP routine is dedicated to the folding of the input spectrum with the selected effective area function, optionally performing an integration over one or more selected energy ranges, to compute the source flux or the count rate (if a photon spectrum is provided) as measured by a given instrument. A direct comparison with the observations, or even with laboratory measurements, is possible after having converted the synthetic spectrum to the detector units. ASAP can take into account specific features of many detectors, such as the energy binning or the scale of the spatial pixels. For some detectors, like the R O S A T Position Sensitive Proportional Counter (PSPC) or the Soft Imaging Spectrograph of ASTRO-D, the socalled energy redistribution (or response) matrix is used; it describes the probability that one incident photon of a given energy is detected in any of the detector energy channels. If we call this matrix R(J, E), where J runs over the detector channels and E is the energy of the incident photons, then the observed spectrum, in units of counts per second, is computed by an ASAP routine as follows:

I(J) = f

A (E)R(j, E) de,

where all terms have been already defined above. In practice, this expression is applicable to detectors with constant gain, and is valid for simulations of point sources observed at the center of the field of view, and where the photon-collecting area is much larger than the F W H M of the detector point spread function. In fact, the variation of the effective area with the off-axis angle in the detector plane is currently not taken into account, and no correction is applied for photons

111

lost outside the source detection cell, but both effects can be easily introduced when needed. An alternative approach is to use the same ASAP routine to convolve the plasma emissivity table with the instrument response, in order to obtain the quantity

f P( r, E) G(T, J) = j -~ A(E)R(J, E)

dE,

which can be supplied to the spectral synthesis routine to get the observed spectrum, already binned in instrument channels. This approach, although mathematically equivalent, has the advantage of computing the synthesized "observed" spectra directly with the resolution required, rather than performing many integrations over an highly resolved spectral domain, as it might occur in the previous case. In the case of detectors which provide measurements in selected emission lines, like the flat crystal spectrometer on board of the solar maximum mission, or in selected spectral bands, e.g. the Yohkoh soft X-ray telescope, the instrument response might be coded into an effective emissivity table, G(T), which represents the plasma emissivity per unit volume, directly in detector units, for each observed line or spectral band. This table can be also supplied to the ASAP spectral synthesis routine (see Section 3.4) to get either volume-integrated plasma emission, or the spatial distribution of the plasma emission, depending on the needs. The latter can be compared with spatially-resolved images of the source, e.g. a solar coronal loop or a supernova remnant, after performing an integration of the spatial distribution of the model emission over a scale length equal to the instrument resolution. Effective quantities (e.g. effective temperature of the plasma) averaged over the instrument band response (and the volume, if needed) can be computed with the instrument-folded emissivity table G(T), through the spectral synthesis routine (Section 3.4) as follows:

(Q) = fQ(s)ne(s)np(s)G(T(s)) dV fne( s)np( s)G( T( s) ) dV ' where Q(s) is the selected quantity to

be averaged. This kind of computation is useful for sev-

112

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

eral tasks, such as to explore the sensitivity of the instrument to variations of the plasma temperature, or to estimate the expected value of velocity which might be reflected in the Doppler shift in a line.

10 7

........

i

........

I

.

.

.

.

.

.

.

.

011 106 ~ ,'~

010 \

4. Sample applications

"--.

10 5

In order to enlighten the large range of applications possible with this tool, both in terms of physical problems of interest and of possible presentations of results, in the following we show a few representative examples.

4.1. Solar type coronal loops We have modelled a coronal loop of magnetically confined hot plasma, of the kind commonly observed on the sun with X-ray telescopes such as S-054 on Skylab [42], N I X T [10] and Y o h k o h / S X T [41]. We have used the model of Serio et al. [37] to compute the hydrostatic configuration of plasma in steady energy balance under the effect of heat conduction, radiative losses and a phenomenological heating term. The loop base pressure is 1 dyne cm -2, and the arch semi-length is 3 × 101° cm. It has been shown [34,37] that, under the reasonable assumption of negligible heat flux in chromosphere, these two parameters univocally determine the loop plasma stratification. Solar abundances are assumed but different abundances could easily be taken into account. The numerical code is used to generate the density and temperature distribution along the arch, shown in Fig. 2. From the computed values and the table of emissivity vs. wavelength and vs. temperature derived from the Raymond code, we synthesize the spectrum emitted at any location in the loop, with the appropriate choice of energy band, spectral resolution, and temperature resolution. The synthesized global loop spectrum, as shown in Fig. 3a with a resolution of 1 eV, is obtained by integrating along all the loop. The same model can be used to compute stellar emission, as done already for instance by Schmitt et

IE ,0

09 \

lO 4

108

.

.

.

.

.

.

.

.

i

.

.

.

.

.

.

.

.

i

.

.

.

.

.

109 1010 S [cm]

.

.

10 1

Fig. 2. Temperature and density distribution along field line coordinates in a symmetric static coronal arch, as computed with the model by Serio et al. [37]. The model computes only half of the arch, because of its symmetry with respect to its apex. The maximum temperature, located at the loop apex, is 4.3 10 6 K.

x

al. [36] with an early version of some of the ASAP routines, by considering the effect of stellar gravity on plasma stratification, the fraction of a star's surface covered by loop feet and the interstellar absorption. Folding the loop spectrum through the instrument response allows to obtain a realistic prediction of an observation. Here we report the result of a hypothetical observation of the loop reported above by means of the R O S A T / P S P C whose effective area is reported in Fig. 3b. The synthesized observed spectrum is presented in Fig. 3c.

4.2. Observations fitting with a model A model of the emitting object can be used within ASAP to perform a fitting of real observations. In Fig. 4 we present a specific example of this kind. We can, for instance, perform the X z fitting of real observations of a corona - in this case the PSPC spectrum of VB34 [27] - with a grid of loop models. Using the loop model emission synthesis discussed above we have synthesized the spectrum, in the focal plane of the PSPC, for a grid of values of base pressure and loop semi-length, taking into account the value of

.4. Maggio et aL / Computer Physics Communications 81 (1994) 105-119 .

I ,o-'°~

.

.

.

.

113

i

(o) il Jt

1°-1~

, m

e

0_14

10-15 I

L

,

,

,

i

1.0

0.1

E [k,V] .

i

.

.

.

.

.

I

(b) 1 O0

E ,_E, u

10

=o

0.1

I.O E [keY]

1000

(c) I 100

L_

/

L__

/

8 o

10

L

L

I

,

0.1

1.0

E [k,v]

L L 1

Fig. 3. (a) Global spectrum of the loop in Fig. 2, with a spectral resolution of 1 eV. The spectrum has been normalized assuming that 10% of the visible stellar hemisphere is covered by loops with the same characteristics, and the source distance from Earth is 3 pc. (b) ROSAT telescope + PSPC effective area vs. photon energy. (c) Final synthesized spectrum in the focal plane of ROSAT + PSPC, computed by folding the spectrum of panel (a), through the effective area in panel (b), and the instrument response matrix. An exposure time of 104 s has been assumed.

114

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

brightness along a representative loop of semilength 3 x 109 cm, in the first 38 seconds of evolution since the beginning of the impulsive phase, in a format less detailed than a standard graph but more like real observations. We have assumed a loop cross section radius of 1 / 1 0 of the loop half-length, a value consistent with most observations, corresponding to an area of 2.8 x 1017 cm 2, and an exposure time of 50 ms, representative of observing modes with Yohkoh SXT. The loop has been divided along field line coordinates into sections of ~ 2.5 arcsec, the angular resolution of SXT and all the emission from each section has been used to compute the instrumental counts, reported in the Figure. The flare is triggered by a heating pulse at the two footpoints of the loop. The main visible effect is a progressive brightening of the loop upwards from the footpoints, due to the rapid filling up of the loop

the star's gravity, its distance and the interstellar absorption. As shown in Fig. 4, in this case the model is able to find a region of minimum X 2 with a minimum reduced X 2 ~ 1.3. Any of the crosses in the graph correspond to a computed loop model spectrum, the quasi-horizontal lines identify the locus of loop semi-length equal to the pressure scale height (lower line) and ten times this height (upper line). It is worth emphasizing that this approach provides an analysis of observations as straightforward as that typically pursued with single or double temperature models of the coronal plasma but with a more realistic model of the coronal structure.

4.3. Evolution of evaporation fronts during a solar flare This example is taken from Reale and Peres [32], a work performed with the P a l e r m o - H a r v a r d code and centered on the evolution of the surface brightness distribution along a coronal arch during a flare. We refer to Peres et al. [22,23], Pallavicini et al. [21] and Peres [24], for previous modeling of flares with this code, and the interpretation of results with a preliminary version of some routines present in ASAP. The present task is to diagnose the chromospheric evaporation front, i.e. to find the physical conditions which characterize the expansion of the overheated upper layers of the chromosphere into the corona, during the flare impulsive phase, and the possibility of direct imaging of this front and its evolution. ASAP has played a crucial role in this work. This front has been imaged for the first time by Hudson et al. [12] by means of the imaging Soft X-ray Telescope (SXT) on board Yohkoh, and had been predicted by Peres and Reale [26]. In order to model the observations from Y o h k o h / S X T [41], we have folded the results of the code through the instrument response in the filters reported in Tsuneta et al. [41]. Among the five X-ray filters available, we pay particular attention to the hardest one, the Be 119 txm filter, appropriate for flare observations. Fig. 5 reports the evolution of the surface

7.6

"'

' ~

,

~ ~ ~ ~ ~

,

[[

,

,

,

+ /

-

-

! i

+ +

7.4 +

\

7.2 x

~ 7.O )

/ +

6.8

6.6 +

6.4 0.00

~

v

~

~

~

+

~

I

I

0.05

0.10 log

pO

+

+

L~

•

~

,

I 0

15

•

,

•

• 020

[dyn/cm2]

Fig. 4. Fitting of an observation of VB34, made with ROSAT + PSPC [27]. Each cross marks the plasma base pressure and the maximum temperature of a loop model whose focal plane spectrum has been computed for the X 2 fitting. The quasihorizontal lines identify the locus of loop semi-length equal to the pressure scale height at the loop apex (lower line) and ten times such a height (upper line). The closed curves are isolevel curves of reduced X 2 values corresponding to lo--, 2tr-, and 3o--values in the space of the two parameters of interest [2]. The minimum, identified by a small square, corresponds to a reduced X 2 ~ 1.3.

115

A. Maggio et aL/ Computer Physics Communications 81 (1994) 105-119

volume by dense plasma coming up from the chromosphere (the "chromospheric evaporation"). 4.4. S h o c k f r o n t across a bubble in the interstellar medium

Here we present an example computed with the Palermo 2D hydrodynamic code [30]. X-ray observations of supernova remnants often show, within the typical approximately circular emitting region, inhomogeneities and more

.

m

f

,oo

intense knots of emission. The latter are interpreted as overdense regions of interstellar medium, presumably shocked by the passage of the supernova front. In order to explore these structures with a high level of physical insight, Reale et al. [31] are studying.the problem of shock fronts running in the interstellar medium and impinging on a local condensation of plasma, by means of the aforementioned 2D hydro code. The task is to study, among other things, the physical conditions conducive to the bubble disruption (or survival) along

.i

,.o

ii

,,o

i

~

26.0

30.0

0 304 607 911 Fig. 5. Grey scale images of the surface brightness along the loop as detected with Yohkoh SXT in the Be 119 ixm filter, during the early phase of the thermal flare. Time in seconds since the beginning of flare energy release labels each image of the loop. Every section of the loop has the same length along the field line coordinates as a resolution element of SXT/Yohkoh (2.5 arcsec). The grey scales, reported at the panel bottom, are in units of SXT's detector digital numbers.

116

A. Maggio et a L / Computer Physics Communications 81 (1994) 105-119

with a detailed study of time-dependent nonequilibrium ionization [3]. Fig. 6 reports the evolution of a bubble in the interstellar medium, subject to the impact of a shock front. The environment has the same characteristics as that of the Vela SNR. In the figure we report the evolving structure of density, temperature and emission in the O VIII (A 18.97 ,~) resonance line on a section across the shock-bub-

Log(Density)

m

0.5 II 0.0 -0.5 -1.0

Temperature

ble system. The shock propagates upwards, and the model has been computed assuming cylindrical symmetry around the vertical axis.

5. Conclusions and future prospects

The approach pursued with ASAP has already given a significant set of results: we have an

1.0

Log(OVlll)

0.5

0.8

-0.5

0.6

-1.0 -1.5 -2.0

t=2560 yrs

"

~

$

0

0

yrs

yrs

yrs Fig. 6. Bubble in the interstellar medium subject to the impact of a shock front. The figure reports the structure of number density (cm-3), temperature (107 K), and O VIII (A 18.97 ,~) line emission per unit volume (10 -23 erg c m - 3 s -1) in a cross section of the shock-bubble system during the evolution, sampled at the times shown. The length of each panel side is 2.25 × 1019 cm. Initially, the bubble is 14 times denser than the surrounding medium. The grey levels are reported at the top of the figure.

A. Maggio et aL / Computer Physics Communications 81 (1994) 105-119

integrated software package which allows to present, analyze and interpret the results of models of astrophysical plasmas. The approach aims also at filling the traditional gap between theory and experiment, by directly connecting the model of an astrophysical object with a system to synthesize its emission at a highly realistic level of detail, and by folding its emission with the response of existing (or planned) instruments. As a matter of fact, with ASAP we already go from the plasma fluid equations, down to the determination of numbers measured in the focal plane of the instrument. Another point worth emphasizing is that, for studies in which the dynamic evolution is a fundamental aspect of the problem, such as in the last two examples shown, the possibility given by IDL of making movies with model results allows a considerable insight. We already plan to enlarge such a system to include interfacing with two-fluid, and 3D hydrocodes of astrophysical plasmas, with codes for the detailed computation of non-equilibrium ionization in complex systems, such as the shock/bubble system briefly discussed above. On the software side we plan to interface the system with other traditional systems such as IRAF/PROS, MIDAS, and XSPEC, in order to enlarge its scope and to allow an easy exchange of results among such software environments. We plan to include the use of IDL widgets (i.e. an easy mouse-driven menu) to enhance the flexibility and ease of use of the system. One of the advantages of ASAP is to be totally open and expandable also simply by "accretion", as long as the subroutines entering the system just respect a few elementary protocol interfaces within IDL. Therefore our scope is that such a system will keep evolving and enriching, providing more and more a common workbench for the theorist and the observer as well.

Appendix In the following we give a brief description of the main routines currently available in ASAP, grouped in categories.

117

A.1. General data I / 0 rddata: Reads a direct-access binary data file

(standard format), according to a file descriptor. Time-independent ("static") quantities are loaded into an array. In case of time-dependent models, a pointer to the section of the data file containing time-dependent ("dynamic") quantities is also returned. rdhead: Reads an ASCII file descriptor. select: Gets the

distribution of the selected plasma parameter from the input array, according to a file descriptor. The "dynamic" file pointer (see rddata) may be also provided instead of the "static" array. In the latter case, the order number of the desired sampling time must be provided.

selpar: Extracts scalar parameter values from the

input file descriptor. write_des: Writes a file descriptor interactively. wr_stand _fi: Writes an array into a standard for-

mat binary file. A.2. Drivers for external Fortran codes stardy_exe: Driver for the Serio et al. [37] static

coronal loop model. syphon _exe: Driver for the Orlando et al. [19]

siphon flow loop model. gostardy: Creates a set of 1D hydrostatic loop

models interactively. ray _ exe: Driver for the Raymond thermal emis-

sion code. Several ranges of energy bins, with variable bin sizes, can be selected. Emissivities are computed at each temperature in the input set, for a fixed plasma density. Optionally, the plasma density can be changed along with the temperature. Plasma abundances can be also modified. By default, separate contributions from free-free, free-bound, and two-photon processes remain avail-

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

118

able in the output file (in standard format), together with the total emissivity values. Alternatively, a set of line emissivities is computed and saved. The total emissivities are optionally rebinned in wavelengths bins by the code, and the results are returned by the routine.

A.3. Emissivity models get_radia: Performs the integration of the input emissivity table P(T, E) in photon energy, to compute total radiative losses as a function of temperature. The resuiting function, P(T), is optionally downloaded in standard format. rd_etab: Reads the required emissivity table saved in standard format. A. 4. Spectral synthesis emiss: Produces synthetic spectra for 1D optically thin plasma models. Returns both the spatial distribution of emission and the integrated emission flux. Optional space-dependent weights can be included. absism: Computes the ISM absorption coefficient, as a function of energy and neutral hydrogen column density along the line of sight, assuming cosmic abundances, according to Morrison and McCammon [17] or Brown and Gould [5]. ism: Computes the ISM absorption coefficient, as a function of wavelength, for given values of the neutral hydrogen column density, and the number density ratios nne/n H and nHei/n H [35]. A.5. Instrument data input

A.6. Instrument folding afold: Folding of source spectrum or emissivity table with instrument effective area, and energy integration in one or more energy ranges. Optionally, the integrand is returned. Optional energy-dependent weights can be included. afoldbin: Folding of source spectrum or emissivity table with instrument effective area and response matrix. A range of detector channels can be selected. Optional energy-dependent weights can provided. A. 7. Instrument specific routines pspc _resp: Generates a table of Raymond-Smith thermal spectra folded through the ROSAT/PSPC effective area and response matrix. sun_emiss: Computes the emission in selected energy bands of the solar experiments Yohkoh/SXT, NIXT, or the integrated intensity in the SMM/FCS emission lines, from 1D loop plasma models. A.8. Graphic routines eloop: Calculates results from the output of the Palermo-Harvard code, and plots a loop emission distribution folded with the response of several solar X-ray observation instrument, in various representations. plot2d: Plots distributions of plasma quantities, such as density or temperature, of Palermo-Harvard hydrodynamic models at selected times.

rd_area: Reads the effective area for the specified detector, and optional filter, saved in standard format.

plot3d: 3D plots of selected plasma quantities computed by the Palermo-Harvard code, using time as x-axis coordinate and field line coordinate as y-axis.

rd_resp: Reads the specified response matrix, saved in standard format.

pltpalet: Plots a color palette and the related scale.

A. Maggio et al. / Computer Physics Communications 81 (1994) 105-119

psimgseq: Produces a picture with a sequence of images in a PostScript file.

tvloop: Produces loop-shaped images of distributions of plasma quantities.

wimgseq: Stores a sequence of images in a file. Acknowledgements The authors acknowledge partial support from Ministero dell'Universitfi e della Ricerca Scientifica e Tecnologica, the Italian Space Agency, and GNA-CNR.

References [1] L.W. Acton et al., Sol. Phys. 65 (1980) 53. [2] Y. Avni, Astrophys. J. 210 (1976) 642. [3] F. Bocchino, S. Orlando, F. Reale and G. Peres, in preparation (1993). [4] S. Bowyer and R.F. Malina, in: Extreme Ultraviolet Astronomy, eds. R.F. Malina and S. Bowyer (Pergamon Press, New York, 1991) p. 397. [5] R.L. Brown and R.J. Gould, Phys. Rev. D 1 (1970) 2252. [6] A. Ciaravella, G. Peres and S. Serio, Sol. Phys. 145 (1993) 45. [7] J.L. Culhane et al., Sol. Phys. 136 (1991) 89. [8] I.J.D. Craig and J.C. Brown, Astron. Astrophys. 49 (1976) 239. [9] V. Domingo, in: The SOHO Mission, ESA SP-1104 (1988) p. 49. [10] L. Golub et al., Nature 344 (1990) 842. [11] L. Golub, in: Advances in Stellar and Solar Coronal Physics, eds. J. Linsky and S. Serio (Kluwer, Dordrecht, 1993) p. 71. [12] H.S. Hudson, K.T. Strong, B.R. Dennis, D. Zarro, M. Inda, T. Kosugi and T. Sakao, Astrophys. J. Lett. 422 (1994) L25. [13] L. Golub and M. Herant, Proc. SPIE 1160 (1989) 629. [14] W, Jeffrey and R. Rosner, Astrophys. J. 310 (1986) 463. [15] F.R. Harnden, D.G. Fabricant, D.E. Harris and J. Schwartz, Smithsonian Obs. Spec. Rept. 393 (1984). [16] M. Landini and B.C. Monsignori-Fossi, Astron. Astrophys. Suppl. 82 (1990) 229.

119

[17] R. Morrison and D. McCammon, Astrophys. J. 270 (1983) 119. [18] Y. Ogawara, T. Takano, T. Kato, T. Kosugi, S. Tsuneta, T. Watanabe, I. Kondo, and Y. Uchida, Sol. Phys. 136 (1991) 1. [19] S. Orlando, G. Peres and S. Serio, in preparation (1993). [20] L.E. Orwig, K.J. Frost and B.R. Dennis, Sol. Phys. 65 (1980) 25. [21] R. Pallavicini, G. Peres, S. Serio, G.S. Vaiana, L. Acton, L. Leibacher and R. Rosner, Astrophys. J. 270 (1983) 270. [22] G. Peres, R. Rosner, S. Serio and G.S. Vaiana, Astrophys. J. 252 (1982) 791. [23] G. Peres, F. Reale, S. Serio and R. Pallavicini, Astrophys. J. 312 (1987) 895. [24] G. Peres, Sol. Phys. 121 (1989) 289. [25] G. Peres, R. Rosner and G. Bodo, Mem. Soc. Astron. Ital. 139 (1989) 60. [26] G. Peres and F. Reale, Astron. Astrophys. Lett. 275 (1993) L13. [27] J. Pye and S. Hodgkin, private communication (1993). [28] J.C. Raymond and B.W. Smith, Astrophys. J. Suppl. 35 (1977) 419. [29] J.C. Raymond, private communication (1989). [30] F. Reale, G. Peres and S. Serio, Nuovo Cimento B 105 (1990) 1235. [31] F. Reale et al., in preparation (1993). [32] F. Reale and G. Peres, in preparation (1993). [33] E.M. Reeves, M.C.E. Huber and J,C. Thimothy, Appl. Optics 16 (1977) 837. [34] R. Rosner, W. Tucker and G.S. Vaiana, Astrophys. J., 220 (1978) 643. [35] T. Rumph, S. Bowyer and S. Vennes, Astrophys. J., submitted (1993). [36] J.H.M.M. Schmitt, F.R. Harnden, G. Peres, R. Rosner and S. Serio, Astrophys. J. 288 (1985) 751. [37] S. Serio, G. Peres, G.S. Vaiana, L. Golub and R. Rosner, Astrophys. J. 243 (1981) 288. [38] L. Stella et al., Proc. SPIE, in press (1993). [39] Y. Tanaka, in Observatories in Earth Orbit and Beyond, ed. Y. Kondo (Kluwer, Dordrecht, 1990) p. 81. [40] J. Trumper et al., Nature 349 (1991) 579. [41] S. Tsuneta, L. Acton, M. Bruner, J. Lemen, W. Brown, R. Caravalho, R. Catura, S. Freeland, B. Jurcevich, M. Morrison, Y. Ogawara, T. Hirayama and J. Owens, Sol. Phys. 136 (1991) 37. [42] G.S. Vaiana, L. VanSpeibroeck, M.V. Zombek, A.S. Krieger, J.K. Silk and A.P. Timothy, Space Sci. Instrum. 3 (1977) 19.