THÉMIS/MTR Basic Reduction Software Package

New stuff about a "standard observing procedure & reduction" based on spatio-temporal modulation (e.g., Trujillo Bueno et al. 2001, in "Advanced Solar Polarimetry", ASP Conf. Serie 236, 141) for 2nd solar spectrum data should be put on-line sometime at the begining of 2003. Tarbes MTR, Avril 2003 However, I release two IDL functions for the computation of the phase chromatism of the 400-700 nm and 600-900 nm achromatic Fichou quarter-wave plates used at THÉMIS: - - Please report any problem!
This page aims to be useful for anyone starting to work with THÉMIS spectropolarimetric data obtained with the MTR observing mode. It describes basic and "standard" procedures which have to be applied step by step for data reduction. I also considered in priority users who have none or little experience with the IDL® language. Before getting started, you need first to get FITS files reading subroutines such as MRDFITS that I personnally recommend. So far, all the routines in which it is needed to read raw data use by defaults MRDFITS. It is easy to substitute calls to MRDFITS with others (READFITS, your own...). I also assume that you have been recording your data in the FITS 3 format which consists in "header + NAXIS3 × images", i.e. a "cube of images". The latest version of the MRDFITS package can be found @ as well as at BASS 2000 (OMP/Tarbes site). I would be grateful if you could send me an e-mail message with subject "THEMIS/MTR/BRSP" in case you would pick & use any of the material available from this page. Please also report any troubles and... improvements you would suggest. Thanks!
1. THÉMIS/MTR file format Basically you will have to deal with 3 type of files with respective keywords value (for OBS_MODE): - flat-field sequence of images [RFLAT] - sequence of observation [SCAN] - dark current calibration file [RDARK] It is quite easy to recognize them among the large set of files that you got after your observing run. Indeed, the THÉMIS file format is such that: OBS_MODE value termination ------------------------------------- RFLAT *y3.fts SCAN *b3.fts RDARK *x3.fts For instance, t002_b0300_sp_19991022_08263184_b3.fts is a file which corresponds to a SCAN taken @ 08:26 UT on October 22, 1999. The _sp_ characters indicate that the MTR mode was used while producing the file. The two digits after _b is the reference number of the CCD camera that was used to take the images. The later reference is also directly related to (a) the spectral domain and (b) the output channel (+ ou - Stokes). This information has to be given to you by the local staff. For more details upon the files/data format please consult BASS 2000.
2. Detection of bad images My current estimate of bad images is of 0.5% based on reduction made out of several 104 images taken during Campaign '99. Let us assume that you have downloaded your data on directory /data/your_name/obs_ddmmyy/ and you want to check whether or not you may need to correct or suppress any corrupted image in file t002_b0300_sp_19991022_08263184_b3.fts. The subroutine gives you a graphical output of the "flux" vs. the image # in the datafile and points out to you the rank of any "bad" image. IDL> chemin='/data/your_name/obs_ddmmyy/' IDL> nom_de_fichier='t002_b0300_sp_19991022_08263184_b3.fts' IDL> verification, chemin, nom_de_fichier
Bad Image @ #99!
The figure above shows that image #99 has to be corrected or removed before going to the next step of the reduction process. A correction procedure will be proposed soon: image #p detected as bad by verification it will be changed to the mean of images #(p-1) and #(p+1) [likely to be OK].
3. Summing-up images Once bad images have been taken care of, you will now have to make mean images from the calibration (flat-field and dark current) raw images. The subroutine does the job and the use of it is explained below. IDL> chemin='/data/your_name/obs_ddmmyy/' IDL> nom_de_fichier='t002_b0300_sp_19991022_08263184_b3.fts' IDL> mean_image_f, chemin, nom_de_fichier, a_single_mean_image In case you would have read and modified a large array of size (naxis1,naxis2,naxis3) containing the images to sum-up, use instead IDL> chemin='/data/your_name/obs_ddmmyy/' IDL> nom_de_fichier='t002_b0300_sp_19991022_08263184_b3.fts' IDL> raw=mrdfits(chemin+nom_de_fichier) IDL> ... modifications of array "raw" of dimensions (naxis1,naxis2,naxis3) IDL> ... IDL> mean_image, raw, a_single_mean_image
4. How bended are my spectra? After having performed the previous steps, you should end up with an image of this kind. The example shown here corresponds to the FeI @ 630.1 nm spectral domain. The next step of the reduction will be to correct this flat- field image for its curvature. In this particular case, we may evaluate the curvature parameters - a parabolic branch - from the central telluric spectral line. For this, you need to define a spectral-pixel domain "of use" around the line of interest. This is very easy to do interac- tively, using the profiles procedure [provided by IDL®]. Assume the mean flat-field image is in array "flap", do the following: IDL> tvscl, flap IDL> profiles, flap This is going to open automatically an extra window displaying rows or columns profiles of the image [see your IDL® reference guide]. After inspection of the image [note that you have to do this for each state of polarisation e.g. the (I+V) and (I-V) images], pick the values [x1,x2] delimiting the spectral pixel range containing the line of interest. The procedure will determine for the image passed as an argument, the positions of the spectral line for each row xplus (dimension: naxis2), and the coefficients for a 2nd degree polynomial fit for the positions; [AP y2 + BP y + CP] - with y=findgen(naxis2) - fits [xplus vs. y]. IDL> bending, flap, x1, x2, 0, naxis2-1, AP, BP, CP, xplus % Compiled module: BENDING. % Compiled module: CENT_LINE. % Compiled module: NINT. % Compiled module: AVY_X. % Compiled module: RUN_AV. % Compiled module: POLY_FIT. % Compiled module: POLY_EXTR. % Compiled module: POLY. % Compiled module: SVDFIT. % Compiled module: REVERSE. % Compiled module: SVDFUNCT. BENDING A= -2.51941e-05 BENDING B= -0.00183213 BENDING C= 70.3273 IDL> Note that, for instance, line res=svdfit(yyyy,xcent,3) could be changed to res=poly_fit(yyyy,xcent,2). bending uses the KIS_LIB (courtesy Kiepenheuer-Institut für Sonnenphysik, Freiburg) package routine cent_line. Download here extra routines called by bending: - nint - avy_x - run_av - poly_extr - cent_line Please contact Dr H. Schleicher ( for more information. Another subroutine for determining the curvature A faster and sometimes more reliable way of doing the same tasks as bending is to use another subroutine called curvature. uses the bisector method for the determination of the "line positions", in that case defined by the position of the bisector at a given depth in the line (Rayrole J., 1967, Ann. Astrophys. 30, 257). Hence you need also to download The use of curvature follows as: IDL> curvature, flap, aprox_pos, width_1, width_2, ap, bp, cp, xplus % Compiled module: CURVATURE. % Compiled module: POS_BISECTEUR. CURVATURE A= -2.85022e-05 CURVATURE B= -0.00107062 CURVATURE C= 70.3715 IDL> aprox_pos is the approximate location in pixel of the line, width_1,2 values depending on the line width and depth [typically 5 to 10 pixels; for a 1st use we recommend to pass 2 identical values of 9 pixels], ap, bp, cp, xplus identical than those of bending. Please contact Dr Véronique BOMMIER ( for more information.
5. Correcting for the curvature It is critical to re-align the spectra line by line in order to build a flat-field image by "removing" [see next section] an average spectrum from the mean raw image (after having removed the dark current level of course). Usually, this is done by "collapsing" all the lines of the image down to one single profile. We are going to use the 3 parameters [AP, BP, CP] as input parameters for the correction for curvature procedure IDL> straight, flap, AP, BP, CP % Compiled module: STRAIGHT. % Compiled module: ECLAT. IDL> tvscl, flap
before... after...

Download here:

   - eclat

The procedure straight uses eclat (courtesy G. Molodij, THÉMIS,

6. The flat-field image The mean profile from the flat-field image is found very easily using the total procedure. IDL> mpro=total(flap,2)/naxis2 IDL> plot, mpro, /xst, xtitle='spectral pixel',ytitle='counts'
mean profile from flat
Note that the right part of the plot correspond to the blue(!) part of the spectrum.
before... after...

The procedure while make the rightmost image directly from
the curvature-corrected image. Assuming that array "flap" contains the
later image, the command-line to make the flat image is just:

IDL> flat, flap
IDL> tvscl, flap

The telluric line of reference used for the curvature parameters
determination has completely disappeared in the new image. Remain the
effects of dust on the slit, on the chip window... and fringes.

Always keep in mind that:

                     « fringes in polarization are the
                         enemies of spectropolarimetry » 

                                              Meir Semel

7. Processing the "observation images" Then starts a more automatic treatement of the observation images since (a) we know how to correct the curvature for each spectral line and each channel ([I+V] and [I-V], say) and (b) we have already computed a flat-field correction image (or matrix) for each channel.
Quiet region Small spots
quite quiet small spots

The corrected images shown above correspond to 2 different steps of
a scan across a group of small sunspots taken on August 29, 1999.

Given that scanp is a 2-dimensional array to which the dark
current level has been removed, the basic operations needed to get
corrected images are the following:

IDL> straight, scanp, AP, BP, CP
IDL> scanp=scanp/flap

In general, you may store the SCAN images directly as a 3-dimensional
array when reading the relevant file(s) with MRDFITS. In such a case,
you may need to use loops of this kind:

IDL> for ii=0, naxis3-1 do begin straight, scanp(*,*,ii), AP, BP, CP
IDL> for ii=0, naxis3-1 do begin scanp(*,*,ii)=scanp(*,*,ii)/flap(*,*)

Do not forget that the coefficients [AP, BP, CP] have to be stored
in memory at the stage of the evaluation of line bending from the
analysis of the flat-field image.

8. Training with real data It is time for some training for you to get acquainted with the numerical tools that you found on this page. First download an IDL save file which contains for each channel, raw flat-field and dark current mean images. Start an IDL session and when you get the IDL> prompt, restore the arrays and check their size with help, and substract the dark current levels for one of the channels: IDL> restore, filename='flatdark.sav' IDL> help % At $MAIN$ DRKM FLOAT = Array[140, 286] DRKP FLOAT = Array[140, 286] FLAM FLOAT = Array[140, 286] FLAP FLOAT = Array[140, 286] Compiled Procedures: $MAIN$ Compiled Functions: IDL> flap=flap-drkp Check that the central O2 telluric line is located @ spectral pixel 68. Now, let us go for curvature evaluation followed by the correction of it: IDL> bending,flap, 60, 78, 0, 285, AP, BP, CP, xcp % Compiled module: BENDING. % Compiled module: CENT_LINE. % Compiled module: NINT. % Compiled module: AVY_X. % Compiled module: RUN_AV. % Compiled module: POLY_FIT. % Compiled module: POLY_EXTR. % Compiled module: POLY. % Compiled module: SVDFIT. % Compiled module: REVERSE. % Compiled module: SVDFUNCT. 70.3273 -0.00183213 -2.51941e-05 BENDING A= -2.51941e-05 BENDING B= -0.00183213 BENDING C= 70.3273 IDL> tvscl, flap IDL> straight,flap, AP, BP, CP % Compiled module: STRAIGHT. % Compiled module: ECLAT. IDL> tvscl, flap And finally, remove the mean profile followed by a visualization of the result: IDL> flat,flap % Compiled module: FLAT. IDL> tvscl, flap You should recover the same images as the ones displayed down to §6.
9. Quicklook @ the scanned area Assuming that you have now at hand two (at least) arrays of reduced images of your scan (e.g. scanp of dimensions [naxis1,naxis2,naxis3]), it may be useful to visualize the actual area and structures which have been scanned. This can also be very useful for BASS 2000 in order to document, with the help of small GIF images for instance, the content of your sequence of observations. Using the nearby continuum intensity from the spectra, gives back a [naxis3, naxis2] array containing crude information upon the scanned area. If you also wish to build a GIF image out of it, follow those lines: IDL> map_ic, scanp, px_value_for_continuum, icont IDL> tvwin, icont IDL> write_gif,'whatever_filename.gif',tvrd()
Quicklook Somehow improved...
bunch of sunspots same plus cosmetics

This group of small sunspots was observed on August 29, 1999. On the
leftmost image, the abscissa correspond there to 200 positions by
successive steps of 0.8" of the 0.5" wide spectrograph entrance slit;
the ordinate dimension is about 2' on the surface of the Sun.

The rightmost image experienced some "cosmetics" in order to
(visually) recover the 2D spatial aspect/extension of the scan.
However, you need to feed subroutine with further
technical information such as (a) the "spatial pixel" [arcsec/px]
and (b) the scan step [arcsec].

IDL> cute_map_ic, scanp, pixcont, spatpix, scanstep, icont, cute_icont

Above is given the subroutine usage, quite similar to the one of
map_ic plus the improved image as an output cute_icont.

10. Quicklook for B// and v// (sort of...) Longitudinal components of B and v can be derived from [I+V] and [I-V] images of circular polarization, by considering the positions of the "spectral line" (more precisely it is profiles of I±V). Very crudely, if pos(+) and pos(-) are the respective positions of the [I±V] profiles, B// is proportional to [pos(+) - pos(-)] and, v// to [pos(+) + pos(-)].
Longitudinal Magnetic Field Radial Velocities
B long. map V rad. map makes maps of the (a) difference and (b) half-sum of the
respective positions pos(+) and pos(-) from a set of two arrays of size
(naxis1,naxis2,naxis3) corrected for flat-field and dark current, as @
end of §7.

IDL> restore,filename='/data/paletou/FeI_6301/scans.sav'
IDL> b_v_long, scanp, scanm, 77, 117, blong, vlong
% Compiled module: B_V_LONG.
% Compiled module: POS_LINE.
% Compiled module: CENT_LINE.
% Compiled module: NINT.
% Compiled module: AVY_X.
% Compiled module: RUN_AV.
% Compiled module: POLY_FIT.
% Compiled module: POLY_EXTR.
% Compiled module: POLY.
IDL> tvscl, blong
IDL> tvscl, vlong uses

Training with real data can be done by downloading a file scans.sav
[64+ Megs!] as follows:

          ftp -i
          (login: anonymous; passwd: your e-mail address)
          cd /pub/paletou/opensoft
          get scans.sav

And then, restore the arrays with an IDL session:

IDL> restore,filename='/data/paletou/FeI_6301/scans.sav'
IDL> help
SCANM           FLOAT     = Array[140, 286, 200]
SCANP           FLOAT     = Array[140, 286, 200]

Next chapters (in preparation): - Extraction of the Stokes parameters
- Vector magnetograms... I would then suggest that you consult the Community Inversion Codes page hosted by HAO/NCAR

Author: Frédéric Paletou
Last updated: April 10, 2003