tp-en.md 14 KB
Newer Older
eric buchlin's avatar
eric buchlin committed
1 2 3
# Instrumentation, diagnostics, signal processing
# Analysis of the Sun by UV spectroscopy

eric buchlin's avatar
eric buchlin committed
4
[Éric Buchlin](mailto:eric.buchlin@ias.u-psud.fr), 8 October 2018
eric buchlin's avatar
eric buchlin committed
5 6 7 8 9

The aim of this practical work is to determine different parameters of the solar corona plasma, using spectroscopy of lines emitted in the UV. We will use observations by
* **slit spectrometers**: a grating disperses light at each position along a slit, and we get an intensity as a function of wavelength. When the slit scans the field-of-view, we can get spectral images.
* **narrow-band imagers**: an image is obtained by a telescope by selecting some wavelength band; this band is generally chosen so that it is dominated, as much as possible, by only one spectral line.

eric buchlin's avatar
eric buchlin committed
10
The questions are noted **[Qn]**, the other action items are noted [An]. Only questions **[Qn]** have to be answered in the report, which can be sent in PDF format by e-mail or by a file hosting service such as [http://dl.free.fr](http://dl.free.fr) until 14 October 2018 at midnight.
eric buchlin's avatar
eric buchlin committed
11 12 13 14 15 16

## First contact with the data

You have an account on the `edu-calcul1` and `edu-calcul1` servers, on which IDL and Python are available (as well as Matlab on the former).
For Python, Python3 is preferred (`python3` command, or `ipython3` for an interactive shell).

eric buchlin's avatar
eric buchlin committed
17
The below-mentioned paths are relative to the `/home/ebuchlin/2018tp-plasmas` directory. Please access these data directly there instead of copying them to your account.
eric buchlin's avatar
eric buchlin committed
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105

The `data/` directory contains two sub-directories:
* `AIA/` for the data from the [AIA](http://aia.lmsal.com/) instrument on board the SDO satellite (launched in 2010). AIA takes images of the Sun in narrow wavelength bands (identified by a number corresponding to a wavelength in ångströms) around different UV lines emitted in the solar atmosphere. These data a "level 1.5" data, meaning that basic calibration routines (dark current, flat field...) have already been applied, but the intensities are still in «data numbers» (DN) instead of being in physical units.
* `EIS/` for the data from the [EIS](http://msslxr.mssl.ucl.ac.uk:8080/SolarB/Solar-B.jsp) instrument on-board the Hinode satellite (launched in 2006). EIS is a slit spectrometer, in the  extreme UV (EUV). The `L1/` sub-directory contains "level 1" data (after dark current, flat field... corrections), for 5 successive scans of some field-of-view.

The data files are in [FITS](https://fits.gsfc.nasa.gov/) format, a format standardised by NASA and frequently used in astronomy; they include a header with metadata (in text format, on 80 columns, then one or several data tables (with their own header, also in text format). FITS images can be displayed with the [`ds9`](http://ds9.si.edu/site/Home.html) software.


**[Q1]** According to the EIS `EIS/L1/eis_l1_20110102_125317.fits` FITS file headers, what is the time interval during which the observation has been done (during that time, the spectrograph slit has moved over the field-of-view), the pixel size, the field-of-view, some of the available wavelength ranges?

> Indications: use `less` in a terminal, and look for these keywords: `DATE_OBS`, `DATE_END`, `XCEN`, `FOVX`, `CDELT1`, `CTYPE1`, `TTYPE1`....

**[Q2]** Determine the AIA FITS file (in the 193Å band) taken at the closest to the middle of this time interval (give its name).

[A3] Display AIA FITS files, at several times and in different wavelength bands (193Å, 304Å...), in the `ds9` software. These data can also be explored, more conveniently but in lower quality, using [HelioViewer](http://helioviewer.ias.u-psud.fr/).

**[Q4]** Display the EIS field-of-view [Q1] on the AIA image selected at [Q2].

> Indications: compute the position in pixels on the AIA image corresponding to the centre of the EIS field-of-view determined at [Q1], knowing that pixel (*i*,*j*) in an image has coordinates (`CRVAL1`+(*i*-`CRPIX1`)×`CDELT1`, `CRVAL2`+(*j*-`CRPIX2`)×`CDELT2`); open the AIA image with `ds9` and locate the position of centre of the EIS field-of-view.


## EIS intensity map

The specific libraries needed to process the original Hinode/EIS FITS files are not available on the servers where you have accounts, so we will use data saved by IDL (proprietary "IDL save" format, `.sav` extension) which are the result of a preprocessing by these specific libraries.

These files can be read from Python using the `readsav` function from the `scipy.io.idl` module.

The Hinode/EIS level-1 data are then available in `data/EIS/L1/`, in files whose name are in the form `eis_l1_20110102_125317_09.sav`: after `l1`, the file name is the date and time, as well as the index of the wavelength band, according to the following table

| N. | Main line        | log T<sub>max</sub> |
|----|------------------|---------------------|
| 0  | Fe X 184.54Å     | 6.0                 |
| 1  | Fe VIII 185.21Å  | 5.6                 |
| 2  | Fe XII 186.88Å   | 6.1                 |
| 3  | Fe XI 188.23Å    | 6.1                 |
| 4  | Fe XXIV 191.93Å  | 7.3                 |
| 5  | Ca XVII 192.82Å  | 6.7                 |
| 6  | Fe XII 195.12Å   | 6.1                 |
| 7  | Fe XIII 196.54Å  | 6.2                 |
| 8  | Fe XIII 202.04Å  | 6.2                 |
| 9  | Fe XIII 203.83Å  | 6.2                 |
| 10 | Fe XXIV 255.10Å  | 7.3                 |
| 11 | He II 256.32Å    | 4.7                 |
| 12 | Fe XVI 262.98Å   | 6.4                 |
| 13 | Fe XXIII 263.79Å | 6.6                 |
| 14 | Fe XIV 264.68Å   | 6.3                 |
| 15 | Fe XIV 274.20Å   | 6.3                 |
| 16 | Si VII 275.35Å   | 5.8                 |
| 17 | Fe XV 284.16Å    | 6.3                 |


[A5] Read one of these files in Python or IDL, for window number 6 (Fe XII 195.12Å). One gets the `d` variable. Explore the structure of this variable.

> In Python (in a iPython shell):
> ```python
> from scipy.io.idl import readsav
> d = readsav('..._06.sav')
> d.wininfo      # informations on all windows
> d.d.wvl[0]     # wavelength axis [Å]
> d.d.time[0]    # time axis [s] since the start of the scan
> d.d.solar_x[0] # x axis (heliocentric coordinates) [arcsec]
> d.d.solar_y[0] # y axis (heliocentric coordinates [arcsec]
> d.d.int[0]     # intensity as a function of position and wavelength
> ...
> ```

> In IDL: `help, d, /struct`, then `help, d.int, /struct`, etc.

> Please note: at each time, the slit is at some position *x*, the time axis and the *x* axis are in then in correspondance.


**[Q6]** Compute a map of the integrated intensity in the Fe XII 19.512nm line. What solar coronal structures (as seen with AIA, see [Q4]) can be seen in the EIS field-of-view?

> The idea is to compute at each position the integrated intensity in the line.
> A simple way of doing this is to sum along the wavelength axis.

> Python code skeleton:
> ```python
> import numpy as np
> import matplotlib.pyplot as plt
> from scipy.io.idl import readsav
>
> d = readsav('..._06.sav')
> i = d.d.int[0]    # spectral intensity cube
> x = d.d.solar_x[0]
> y = d.d.solar_y[0]
> ii = np.sum(..., axis=...)  # sum over the λ axis
>
106
> plt.imshow(ii, origin='lower',
eric buchlin's avatar
eric buchlin committed
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
>            vmin=..., vmax=...,
>            extent=[x[0], x[-1], y[0], y[-1]])
> plt.colorbar()
> plt.show()
> ```


## EIS density map

For a given ion, the intensity emitted in each line is proportional to its "emissivity", which is a fuction of the electronic density *n<sub>e</sub>*.
The emissivities in 8377 Fe XII lines, computed using the [CHIANTI](http://www.chiantidatabase.org/) atomic physics database, are stored in the `fe12e` variable that can be read in the `data/EIS/fe12e.sav` file.

In Python, the emissivities for the line of index `n` are then available using:
```python
e = readsav('fe12e.sav')
dens = e.fe12e.density[0]         # density table
em = e.fe12e.emiss[0][n].em[:,0]  # emissivity table (as a function of density)
```

The lines we are interested in now are:
127
* Fe XII 195.119Å, index *n*=873 in the line list
eric buchlin's avatar
eric buchlin committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
* Fe XII 186.887Å, index *n*=750 in the line list

**[Q7]** Plot the ration of emissivities for these two lines as a function of electronic density.

> Indications: Plot the ratio in log-log (Python: `plt.loglog(...)`); put labels on axes (Python: `plt.xlabel(...)`, `plt.ylabel(...)`).

**[Q8]** Compute a density map in the EIS field-of-view.

> Indications: read the EIS data corresponding to both lines; create an integrated intensity map (see [Q6]) for each of them; use `np.interp()` in Python or `interpol()` in IDL to interpolate the relationship between the intensity ratio (equal to the emissivity ratio) and density.

> Beware: the second argument of `np.interp()` must be strictly increasing.


## EIS Doppler shifts

**[Q9]** Compute Doppler velocity maps in the Fe XII 195.119Å and Ca XVII 192.82Å lines, from the first EIS observation.

145
> Indications: to make this simple, we can compute at each position the gravity centre of the line in the wavelength window, subtract the wavelength of the line at rest, and compute the corresponding Doppler velocity (*v*/*c*=Δλ/λ<sub>0</sub>).
eric buchlin's avatar
eric buchlin committed
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227

> In Python, the first step (gravity centre) can be compute very easily on all positions simultaneously using the `np.dot()` and `np.sum()` functions.


A difficulty in such a Doppler velocity determination is that the wavelength axis calibration is actually sensitive to different effects, depending in particular on the position of the Hinode satellite on its orbit (Low Earth Orbit, with about 90min period), including:
* the Doppler effect due to the motion of the satellite around the Earth,
* the variations of the light (in particular in infrared) coming from the Earth (producing thermal effects on the mechanical structure of the spectroscope).

This gives rise to temporal variations of the position of each line, that come in addition to the variations that have a solar origin.

**[Q10]** Show that there is such a temporal variation, by plotting the Doppler velocity in the Fe XII 195.119Å line, average along the spectrometer slit, as a function of time.

> Indications: we can open successively the 5 data files for this line, compute the Doppler shift map for each of them (as in the previous question), average along the slit (Python: `np.mean`), and plot this average as a function of time.

> The difficulty is to have a common time axis for all 5 observations; this can be done by starting from the list of file names and then computing the duration in seconds between 12:00:00 and the start of the observation:

> ```python
> files = ['125317', '131848', '134422', '140956', '143530']
> for f in files:
>    d = readsav(...)
>    t0 = int(f[4:6]) + int(f[2:4]) * 60 + (int(f[0:2]) - 12) * 3600
>    ...
>    dopp = ...
>    plt.plot(t0 + d.d.time[0], dopp)
> ...
> ```

**[Q11]** Estimate the variation of the Doppler velocity originating from the motion of the satellite around the Earth. Is this motion sufficient to explain the observed Doppler velocity variations?

**[Q12]** A wavelength correction (corresponding to known instrumental effects) as a function of time is proposed in the Hinode/EIS data files (`d.wave_corr_t` variable in the `.sav` files): can it effectively correct the time variation measured in [Q10]? Use it to produce corrected Doppler velocity maps. Is the average shift now a blue shift or a red shift?

> Indication: `d.wave_corr_t` has to be subtracted from the measured Doppler shifts.

> More corrections are available:
> * `d.wave_corr_tilt`, depending on the position along the slit
> * `d.wave_corr_t`, for each pixel; this is the sum of `d.wave_corr_t` on the time (or *x*) axis and of `d.wave_corr_tilt` along the slit (*y* axis).


## Determination of coronal temperature

The table below gives typical temperatures of the plasma emitting the main lines contained in the different AIA bands (the name of channels correspond to a wavelength in ångströms).

| Channel | Main line        | log T<sub>max</sub> |
|---------|------------------|---------------------|
| 304     | He II            | 4.7                 |
| 171     | Fe IX            | 5.8                 |
| 193     | Fe XII, XXIV     | 6.1, 7.3            |
| 211     | Fe XIV           | 6.3                 |
| 335     | Fe XVI           | 6.4                 |
|  94     | Fe XVIII         | 6.8                 |
| 131     | Fe VI, XX, XXIII | 5.6, 7.0, 7.2       |


The measured intensity in each band actually depends on the *wavelength response* of the instrument and on the intensities received in all spectral lines, and these intensities depend (in the coronal approximation) on the plasma quantity at some density and temperature.
These informations can be combined to give the instrument *temperature response*, that represents for each temperature the intensity measured by the instrument in one of its band, for a given plasma quantity at this temperature, in the volume constituted by a pixel and the line of sight.

The `data/AIA/aia_response.sav` file contains the AIA wavelength response (`wvlresp` variable) and temperature response (`tempresp` variable), for each AIA channel, as a function of the wavelength (in ångströms) and temperature (base 10 logarithm, kelvins) respectively.

**[Q13]** Plot the temperature response curves for the different SDO/AIA bands.

> To get the useful data in Python, for example for the 211Å channel:
> ```python
> r = readsav('aia_response.sav')
> logt = r.tempresp.logte[0]
> r211 = r.tempresp.a211[0].tresp[0]
> ```

**[Q14]** Assuming that the plasma is isothermal, and using the ratio between the intensities in two well-chosen AIA channels, get an estimate of the plasma temperature.

> Indications: for an isothermal plasma (at temperature *T*), the intensity measured in each channel is proportional to the temperature response at temperature *T* in this channel.

**[Q15]** Using two SDO/AIA intensity ratios, show that the plasma is *not* isothermal everywhere.


A more precise determination of the temperature (and more generally of the distribution of the quantity of plasma as a function of temperature) could be obtained by using a large number of line intensities from a spectrometer like EIS.



*[AIA]: Atmospheric Imaging Assembly
*[SDO]: Solar Dynamics Observatory
*[EIS]: Extreme-UV Imaging Spectrometer
*[DEM]: Differential Emission Measure