Error while importing .stab IndexedSEDFamily

[francesco-chiti-tegli]

Good morning,
I hope this message finds you well.

I am working with SKIRT to monochromatically (so in oligochromatic regime) simulate the dust scattering of a Halpha gas emission from a simulated galaxy.
With Phyton I previously computed the luminosity of each gas cell (which are 8388608).

In order to do so, I thought of importing the gas cells as a CellSource Class from a .txt file (containing and index "ID" that identified each cell), and then importing the Halpha luminosity as a .stab file through a IndexedSEDFamily Class.

The two files are respectively called Halfa_Ridotto.txt and SEDfamily_Ridotto.stab.

I tried to write a function that could create a .stab file from a (8388608, 2) numpy array whose first column is just the wavelenght of Halpha (0.65628 micron) and the second are the luminosities of the cells.
Reading said file, after saving it, returns me this text:

{'axisNames': ['lambda', 'ID'],
'axisUnits': ['m', '1'],
'axisScales': ['lin', 'lin'],
'quantityNames': ['Llambda'],
'quantityUnits': ['W/m'],
'quantityScales': ['lin'],
'lambda': <Quantity [0.65628] m>,
'ID': <Quantity [0.000000e+00, 1.000000e+00, 2.000000e+00, ..., 8.388605e+06,
8.388606e+06, 8.388607e+06]>,
'Llambda': <Quantity [[2.27619808e+40, 2.27619808e+40, 2.00511234e+40, ...,
5.37732446e+39, 5.20180433e+39, 5.20180433e+39]] W / m>}

So I believe the file has been written correctly.
On the other hand, the .txt file containing the gas cells is structured with the following header (with # before) and data:

Halfa_Ridotto.txt: import file for cell source -- Halfa
Column 1: xmin (pc)
Column 2: ymin (pc)
Column 3: zmin (pc)
Column 4: xmax (pc)
Column 5: ymax (pc)
Column 6: zmax (pc)
Column 7: mass density (g/cm3)
Column 8: ID (1)

-512 -512 -256 -508 -508 -252 7.32388e-28 0
-512 -512 -252 -508 -508 -248 7.32388e-28 1
-512 -512 -248 -508 -508 -244 6.95805e-28 2
-512 -512 -244 -508 -508 -240 6.95805e-28 3
...

When trying to run SKIRT with an sources-only oligochromatic (0.65628 micron) with the previous input files, I get the following error:

Error: Axis 1 does not have name index(1) in stored table SEDfamily_Ridotto.stab

This error realy puzzles me, since the .stab file can be correctly read on Python through PTS readStoredTable function.

So I was wondering if you had any advice on how to proceed, as I am completely stuck on this error.

I thank you very much in advance for your availability,
Best regards,

F

[petercamps]

The SKIRT code is quite picky about the axis names and units given in the stab file. In this case, SKIRT expects the second axis to be named "index", while in your stab file the second axis is named "ID". The "(1)" in the error message indicates the expected units - dimensionless - for the index axis. This seems to be correct in your stab file.

As further guidance, below is a short python script that I used some time ago to generate a trivial stab file for testing. I hope this helps.

# script that generates an artifial emission line spectrum file
# the file contains two significantly different spectra:
#  - a single narrow line spectrum
#  - a flat spectrum over a significantly larger wavelength range

import numpy as np
import pts.storedtable as stab

lam = np.array([0.09, 0.1, 0.14, 0.15, 0.16, 0.17, 10.0, 10.01]) * 1.0e-6
index = np.array([0.0, 1.0])
Llambda = np.zeros((lam.shape[0], index.shape[0]))
# first spectrum
Llambda[3:5, 0] = 1.0e30
# second spectrum
Llambda[1:7, 1] = 1.0e28

# write to file
stab.writeStoredTable(
    "spectrum.stab",
    ["lambda", "index"],
    ["m", "1"],
    ["lin", "lin"],
    [lam, index],
    ["Llambda"],
    ["W/m"],
    ["lin"],
    [Llambda],
)

[francesco-chiti-tegli]

Good morning,
thank you very much for your prompt and very precise answer: after said correction, now I can successfully import the Halpha luminosities of the ISM cells.
I just had a second, very brief, question concerning the SEDInstrument output: is the flux unit output in Jy?
Thank you very much again,
Best regards,

F

[petercamps]

The physical quantity and corresponding units for each column are listed in the header of the text output files. So you can check for yourself.

Also, the units used will depend on your selections in the ski file, as explained at the start of the first SKIRT tutorial: https://skirt.ugent.be/root/_tutorial_basics_mono.html.