search

rhdtownsend / msg

Multidimensional Spectral Grids
GNU General Public License v3.0
11 stars 3 forks source link

Possible unit error in PyMSG #18

Closed MartianColonist closed 1 year ago

MartianColonist commented 1 year ago

First, thank you for creating MSG! I'm excited to explore applications of your package.

When comparing the output of PyMSG (following the tutorial in your documentation) with Pysynphot I've noticed what could be a unit misspecification in the PyMSG output. Below is a comparison for a star with Teff = 4803, [Fe/H] = 0.0, and log_g = 4.57 (both spectra are surface flux, not the observed flux at Earth).

Converting the PyMSG output on the right (erg / cm^2 / s / A) to W / m^2 / m requires multiplying by 1e-7 1e4 1e10 = 1e7, which produces a spectrum 8 orders of magnitude greater than the Pysynphot output on the left. The peak flux on the left matches a black body with the same Teff, so the Pysynphot output looks correct.

I'm wondering if the output units for PyMSG are actually erg / cm^2 / s / cm ? Then the conversion factor would be 1e-7 1e4 1e2 = 1e-1 and PyMSG would agree with Pysynphot after converting to W / m^2 / m.

Thanks in advance!

Pysynphot_PyMSG_comparison

rhdtownsend commented 1 year ago

Hi Ryan --

Thanks for the issue. To help track down what's going on, can you tell me what steps you're using in Pysynphot to calculate the spectrum? Also, I presume you're using the demo grid with MSG, right?

Best wishes,

Rich

MartianColonist commented 1 year ago

Hi Rich,

With Pysynphot I'm using the following code:

import pysynphot as psyn

T_eff = 4803.0
Met = 0.0
log_g = 4.57

sp = psyn.Icat('ck04models', T_eff, Met, log_g)

sp.convert("um")               # Convert wavelengths to micron
sp.convert('flam')             # Convert to flux units (erg/s/cm^2/A)

wl_s = sp.wave                      # Stellar wavelength array (um)
F_s = (sp.flux*1e-7*1e4*1e10)       # Convert to W/m^2/m

I'm using the Goettingen-HiRes grid with MSG.

Thanks,

Ryan

rhdtownsend commented 1 year ago

Hi Ryan --

I can confirm your problem -- looks like I had assumed the wrong units when converting the Goettingen grids to MSG format. I'm just finished uploading new versions of the grids to our webserver -- please go ahead and download the revised versions. You can check that you've got a new version by e.g. running the command

h5dump -a label sg-Goettingen-HiRes.h5

In the output, you should see the line

(0): "Goetingen HiRes Grid / v2 / 2023-02-19 / Rich Townsend"

cheers,

Rich

MartianColonist commented 1 year ago

Thanks, it looks like that did the trick for the units.

On an unrelated issue, I'm hitting many "unavailable data" exceptions in regions of parameter space that fall well within the listed photospheric parameter range.

Examples:

Both throw the attached error, while the supported range for the Göttingen grid is:

I wonder if this could be related to the grid dimension being different as a function of the alpha element abundance. Husser+2013, Table 1, notes that "Alpha element abundances [alpha/Fe] =/= 0 are only available for 3500 K ≤ Teff ≤ 8000 K and −3 ≤ [Fe/H] ≤ 0." The two examples above should satisfy this requirement (since [alpha/Fe] = 0), but maybe that data isn't included in the 'sg-Goettingen-HiRes.h5' file?

pymsg_unavailable_data

rhdtownsend commented 1 year ago

Hi Ryan --

Apologies for the delay in getting back to you about this. What seems to be happening is that the extent of the Goettingen grids (in Teff, logg, etc) is significantly larger at solar [alpha/Fe] than for alpha-enhanced/depleted cases. You are trying to interpolate in the extended part of the grid, but at solar [alpha/Fe] -- so one would think your calculations should still work.

However, MSG is only able to interpolate in N-dimensional cells where all 2**N corners of the cell are defined. So, even though you are requesting an interpolation at solar [alpha/Fe], there are no data at the adjacent higher/lower [alpha/Fe], and so MSG cannot form interpolation cells with all corners defined.

The solution to this is to split the Goettingen grids into two: an N dimensional grid that doesn't include the extended data at solar [alpha/Fe], and an N-1 dimensional grid that includes only solar [alpha/Fe] data.

You can grab these grids from the following URLs:

http://user.astro.wisc.edu/~townsend/resource/download/msg/grids/Goettingen/HiRes/v3/sg-Goettingen-HiRes.h5 http://user.astro.wisc.edu/~townsend/resource/download/msg/grids/Goettingen/HiRes-alpha/v3/sg-Goettingen-HiRes-alpha.h5

http://user.astro.wisc.edu/~townsend/resource/download/msg/grids/Goettingen/MedRes-R/v3/sg-Goettingen-MedRes-R.h5 http://user.astro.wisc.edu/~townsend/resource/download/msg/grids/Goettingen/MedRes-R-alpha/v3/sg-Goettingen-MedRes-R-alpha.h5

http://user.astro.wisc.edu/~townsend/resource/download/msg/grids/Goettingen/MedRes-A/v3/sg-Goettingen-MedRes-A.h5 http://user.astro.wisc.edu/~townsend/resource/download/msg/grids/Goettingen/MedRes-A-alpha/v3/sg-Goettingen-MedRes-A-alpha.h5

(the ones with alpha in their name include [alpha/Fe] as an interpolation axis, the ones without have fixed, solar [alpha/Fe]).

For your purposes, try switching to the sg-Goettingen-Hires.h5 grid.

cheers,

Rich

MartianColonist commented 1 year ago

Hi Rich,

Sorry for the delay getting back to you - I just returned from vacation.

I can confirm that everything is now working great. Feel free to close this issue.

Thank you so much for tracking down these issues (and remaking all these HDF5 files!). I'm excited to use PyMSG in my research.

Thanks,

Ryan

rhdtownsend commented 1 year ago

You're welcome!