Code, data, and plots for the publication: Characterising the surface magnetic fields of T Tauri stars with high-resolution near-infrared spectroscopy. ADS, arXiv, DOI.
The content of this repository is licensed under MIT license (see the included LICENSE).
This code comes with no warranty whatsoever and you may find it poorly written. You are very welcome to email astro@lavail.net if you have questions, but keep in mind that I will not guarantee assistance if you experience problems installing or running this code.
Here are listed the python packages needed to run the python scripts, and their version as output by pipreqs.
- numpy==1.14.3
- corner==2.0.1
- scipy==1.1.0
- h5py==2.8.0
- matplotlib==2.2.2
- pandas==0.23.0
- requests==2.18.4
- emcee==2.2.1
Additionally, couple of large files in the repository (listed in .gitattributes) are handled with git Large File Storage
analysis includes three codes that fit Zeeman broadening of spectral lines. magn0x corresponds to model 1 in the paper, magn0xx corresponds to model 2 in the paper, and magnotron is the MCMC-based model 4 in the paper (of which model 3 is a subset).
Each code uses grids of synthetic spectra, which have been computed with IDL codes (not included here, not my own), and they are stored in the HDF5 files synth-spectra-2-h5.
I include in this repository the 3 reduced observations of TW Hya, so that the analysis can be run fully for these three spectra. In each directory, there is a bash script run.sh showing the 3 commands to run the analysis for the 3 spectra.
The input observation to each code are the same dataset of reduced, telluric-removed, normalized CRIRES high-resolution near-infrared spectra, which are within the data subfolders. Each .ascii spectrum file is a 2-columns file with:
- wavelength in Å
- normalized intensity
The data for the observation of TW Hya are included in this repository.
magn0x fits the input spectrum with a 2-component magnetic model (one part of the stellar surface is magnetic, the rest is non-magnetic). The array of filling factors used by magn0x.py is generated by running ffgen0x.py.
magn0x.py is run with one argument (integer) which specifies which spectrum shall be input. bash run.sh runs the code on all input spectra.
Results are spit out in the ./analysis/magn0x/results/ directory. One 4-column ascii file is produced for each input spectrum with:
- wavelength in Å
- observed intensity
- best fit synthetic spectrum
- best fit non-magnetic synthetic spectrum
One line is also appended in ./analysis/magn0x/results/magn0x.log, for instance
306726 TWA7 2.34 (4.07)
with the format [OBS ID] [STAR ID] [B in kG] [(chi^2)]
Exact same thing than magn0x but with 2 magnetic components
magnotron.py (models 3--4) performs a MCMC fitting of the observed spectrum. The code starts with the simplest magnetic model (0,2 kG) and fits the spectrum and computes the Bayesian Information Criterion (bIC). Then, iteratively, more components are added (4,6,8,10,12, and ultimately 14 kG) and the process is repeated. The model favoured by the data is the model yielding the smallest BIC. If you want to run magnotron on your computer, make sure that you change the ncpu variable to the number of CPUs available on your machine.
magnotron.py is run with one argument (integer) which specifies which spectrum shall be input. bash runall.sh runs the code on all input spectra.
Results are spilled out in the ./analysis/magn0x/results/ folder, and quality-control plots in ./analysis/magn0x/plots/.
Files for each iteration of the model are named:
- *iter0.ascii for models with 0,2 kG
- *iter1.ascii for models with 0,2,4 kG
- *iter2.ascii for models with 0,2,4,6 kG
- *iter3.ascii for models with 0,2,4,6,8 kG
- *iter4.ascii for models with 0,2,4,6,8,10 kG
- *iter5.ascii for models with 0,2,4,6,8,10,12 kG
- *iter6.ascii for models with 0,2,4,6,8,10,12,14 kG
- Best-fit synthetic spectra for each iteration are output to *-spec*.ascii files, which are 3-columns ascii files with:
- wavelength in Å
- observed intensity
- best fit synthetic spectrum
- The Markov chains (the
samplesvariable) from the MCMC runs are saved into pickle files names samples*.pickle - A log summarizing the results of all iteration for a given observed spectrum is output to log_*, which is an ascii file file with the following columns:
- iter: the model iteration number
- bmin: minimum B value (kG) of the 3-sigma confidence interval
- bbest: median B value (kG)
- bmax: maximum B value (kG) of the 3-sigma confidence interval
- chi2: reduced chi squared value from the fit
- aic: Akaike Information Criterion value
- bic: Bayesian Information Criterion value
- acceptance fraction: the acceptance fraction from the AIE sampler
- ESSmin: the number of independent samples
These plots show the first values of the markov chains for the spectrum [OBSID], the model iteration number N, and variable X (with variables 0 to X-1 being filling factors, and variable X being the scaling factor). The burn-in length is indicated by the vertical line, the median value of the chain after burn-in is indicated by the horizontal line.
These plots show the best fit to the observed spectrum [OBSID] obtained with the model iteration N. x-axis is just the array index, y-axis is normalized intensity, blue points with error bars are observed spectrum, and yellow line is best-fit synthetic spectrum. The [OBSID], name of the star, and median B value are indicated on top of the plot.
These show an histogram of the posterior probability for the magnetic field B, the observed spectrum [OBSID] obtained with the model iteration N.
Luminosity determination: a notebook that I used to determine consistent stellar luminosities of the sample. The requirements to run this notebook as listed by pipreqs are:
- numpy==1.14.3
- astropy==3.0.2
- matplotlib==2.2.2
- astroquery==0.3.9
- scipy==1.1.0