diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 4e0466b6..dad9febe 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -14,7 +14,7 @@ jobs: strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-15, macos-13] # macos-15/14=arm64，macos-13=x86_64 + os: [ubuntu-latest, macos-15] # macos-15/14=arm64，macos-13=x86_64 py: ['3.9','3.10','3.11','3.12','3.13'] name: py${{ matrix.py }} on ${{ matrix.os }} runs-on: ${{ matrix.os }} diff --git a/.readthedocs.yaml b/.readthedocs.yaml index e1abd260..f7e3f303 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -27,3 +27,6 @@ sphinx: python: install: - requirements: docs/requirements.txt + # Install the project itself so version metadata resolves from tags. + - method: pip + path: . diff --git a/changelog.md b/changelog.md new file mode 100644 index 00000000..d997f01c --- /dev/null +++ b/changelog.md @@ -0,0 +1,16 @@ +# Changelog + +## 2026-02-10 + +- Added `derived_param` as the preferred name for derived-parameter callbacks in `solve()`. +- Kept backward compatibility with `dynamic_param`: + - If `dynamic_param` is used, a `DeprecationWarning` is emitted. + - If both `derived_param` and `dynamic_param` are provided (and differ), `ValueError` is raised. +- Updated internal solver logic and user-facing messages to use the "derived parameter" terminology. +- Added `smelib_lineinfo_mode` passthrough in `solve()` call paths (`_residuals` and `_jacobian`) so fitting runs can use SMElib precomputed line-info modes. +- Updated `linelist_mode` naming: + - Preferred values are now `"all"` and `"dynamic"`. + - `"auto"` is kept as a deprecated compatibility alias and maps to `"dynamic"` with a `DeprecationWarning`. +- Added segment-aware optional input `sme.wint` for synthesis transfer grids. + - Priority is now `sme.wint[segment]` first, then internal cached grids (when enabled), then SMElib adaptive grid generation. +- Updated user docs accordingly (`sme_struct`, `quickstart`, and `how-to`). diff --git a/docs/_static/custom.css b/docs/_static/custom.css new file mode 100644 index 00000000..b8f0d867 --- /dev/null +++ b/docs/_static/custom.css @@ -0,0 +1,52 @@ +.home-cards { + display: grid; + grid-template-columns: repeat(3, minmax(0, 1fr)); + gap: 1rem; + margin: 1.25rem 0 1.75rem; +} + +.home-card { + display: flex; + flex-direction: column; + justify-content: space-between; + min-height: 13.5rem; + padding: 1.1rem 1rem; + border: 1px solid var(--pst-color-border, #d0d7de); + border-radius: 0.75rem; + background: var(--pst-color-surface, #f6f8fa); + text-decoration: none; + color: var(--pst-color-text-base, #24292f); + box-shadow: 0 1px 2px rgba(0, 0, 0, 0.05); + transition: transform 120ms ease, box-shadow 120ms ease, border-color 120ms ease; +} + +.home-card:hover { + transform: translateY(-2px); + box-shadow: 0 6px 18px rgba(0, 0, 0, 0.12); + border-color: var(--pst-color-primary, #2b6cb0); + text-decoration: none; +} + +.home-card h3 { + margin: 0 0 0.5rem; + font-size: 1.35rem; + line-height: 1.2; +} + +.home-card p { + margin: 0; + color: var(--pst-color-text-muted, #57606a); +} + +.home-card-cta { + margin-top: 1rem; + align-self: flex-start; + font-weight: 600; + color: var(--pst-color-primary, #2b6cb0); +} + +@media (max-width: 1050px) { + .home-cards { + grid-template-columns: 1fr; + } +} diff --git a/docs/_static/figures/first_spectrum.png b/docs/_static/figures/first_spectrum.png new file mode 100644 index 00000000..73d1f071 Binary files /dev/null and b/docs/_static/figures/first_spectrum.png differ diff --git a/docs/_static/nlte_table.png b/docs/_static/nlte_table.png new file mode 100644 index 00000000..3488feb8 Binary files /dev/null and b/docs/_static/nlte_table.png differ diff --git a/docs/advance/derived_param.md b/docs/advance/derived_param.md new file mode 100644 index 00000000..d8e829b8 --- /dev/null +++ b/docs/advance/derived_param.md @@ -0,0 +1,73 @@ +# Derived Parameters + +`derived_param` lets you define SME parameters as functions of other parameters +during fitting, without adding them to the free-parameter list. + +## What It Does + +In `solve(...)`, PySME: + +1. sets the current free parameters, +2. evaluates each function in `derived_param`, +3. writes the returned value back to `sme`, +4. runs synthesis and computes residuals. + +So derived parameters are updated at every iteration. + +## API + +```py +from pysme.solve import solve + +sme = solve( + sme, + param_names=[...], # free parameters + derived_param={...}, # derived parameters +) +``` + +- `derived_param` is a `dict[str, callable]`. +- Key: SME parameter name (for example `"vmic"` or `"abund Mg"`). +- Value: function `f(sme) -> float`. + +## Rules and Notes + +- A parameter cannot be both free and derived in the same run. +- `dynamic_param` is still accepted as a legacy alias, but it is deprecated. +- For abundance keys (for example `"abund Mg"`), return the **final abundance** + you want in the usual abundance scale. PySME applies the internal + `monh` conversion for you when writing into `sme.abund`. + +## Example 1: Tie `vmic` to `teff` and `logg` + +```py +import numpy as np +from pysme.solve import solve + +derived = { + "vmic": lambda s: np.clip(1.1 + 1e-4 * (s.teff - 5500.0) - 0.3 * (s.logg - 4.0), 0.2, 5.0) +} + +fit = ["teff", "logg", "monh", "vsini"] +sme = solve(sme, fit, derived_param=derived) +``` + +Here `vmic` is never fitted directly; it is recomputed from the current model +state at each iteration. + +## Example 2: Enforce fixed `[Mg/Fe]` + +```py +from pysme.abund import Abund +from pysme.solve import solve + +solar_mg = Abund.solar()["Mg"] + +# Target relation: [Mg/Fe] = +0.20 -> A(Mg) = A_sun(Mg) + [M/H] + 0.20 +derived = { + "abund Mg": lambda s: solar_mg + s.monh + 0.20 +} + +sme = solve(sme, ["teff", "logg", "monh"], derived_param=derived) +``` + diff --git a/docs/usage/faq.md b/docs/advance/faq.md similarity index 97% rename from docs/usage/faq.md rename to docs/advance/faq.md index 58213afb..61c32de3 100644 --- a/docs/usage/faq.md +++ b/docs/advance/faq.md @@ -6,7 +6,7 @@ Call `util.start_logging(filename)`. ```py -from SME import util +from pysme import util util.start_logging("your_log_file.log") ``` diff --git a/docs/usage/fordev.md b/docs/advance/fordev.md similarity index 90% rename from docs/usage/fordev.md rename to docs/advance/fordev.md index a5698ecd..e78ce2c7 100644 --- a/docs/usage/fordev.md +++ b/docs/advance/fordev.md @@ -70,7 +70,8 @@ The follwing table shows the version matching between SMElib and PySME. |PySME version|SMElib release version|SMElib version| |:--:|:--:|:--:| -|v0.5.0-|latest (v6.13.x)|6.13 (June 2025)| +|v0.7.0-|latest (v6.13.x)|6.13 (June 2025)| +|v0.6.23|v6.13.12 (freezed)|6.13 (June 2025)| |v0.4.199|v6.0.6 (freezed)|6.03 (July 2019)| |v0.4.167-v0.4.198|v6.0.6 (not freezed)|6.03 (July 2019)| @@ -85,6 +86,21 @@ The versions are controlled by `git tag`. For PySME, setting a new tag will update the version to the latest tag. For SMElib, setting a new release will based on a tag. +### Docs versioning policy (RTD) + +PySME docs use a tag-driven versioning strategy. + +1. Create release tags in `vX.Y.Z` format (for example `v0.6.23`). +2. `docs/conf.py` resolves `release` from `src/pysme/_version.py` (versioneer), so docs title/version follow git tags. +3. `.readthedocs.yaml` installs both docs requirements and the project itself, so tag metadata is available during RTD build. + +RTD project settings must also be configured once (in the RTD web UI): + +- Keep both `latest` (branch docs) and `stable` (latest release tag) enabled. +- Set default docs version to `stable`. +- Add an automation rule to activate SemVer tags (`v*`) as RTD versions. +- Hide or deactivate outdated branch versions if needed. + ## NLTE departure coefficients What happens when the NLTE grid is added into PySME? @@ -150,4 +166,4 @@ Some extra test is needed to clear the situation. ### Installation -TBD. \ No newline at end of file +TBD. diff --git a/docs/usage/how-to.md b/docs/advance/how-to.md similarity index 97% rename from docs/usage/how-to.md rename to docs/advance/how-to.md index f5214779..f0a371d2 100644 --- a/docs/usage/how-to.md +++ b/docs/advance/how-to.md @@ -245,7 +245,7 @@ The synthesis of long spectra, i.e., a spectra covering a wide wavelength range The deafult synthesis mode of PySME takes all the lines in the `linelist` into account to synthesize the spectra, thus the synthesis takes a long time. To make it even worse, the synthesis is repeated for each segment, thus if we have 10 segment covering from 3000-9000AA, all the lines (including those in ~9000AA) will be used for the first segment (3000-3600AA), and the total time cost is 10 times longer than just using one segment to cover the whole range. -This situation can be improved by using the central depth and line range parameters to select the relevant lines for each segment, by simply adding `linelist_mode='auto'` in `synthesize_spectrum`: +This situation can be improved by using the central depth and line range parameters to select the relevant lines for each segment, by simply adding `linelist_mode='dynamic'` in `synthesize_spectrum`: ```py from pysme.synthesize import synthesize_spectrum @@ -256,10 +256,13 @@ sme.wave = np.arange(3000, 6000, 0.02).reshape(-1, 100) sme.linelist = line_list sme.nlte.set_nlte('Na') sme.cdr_depth_thres = 0.1 -sme = synthesize_spectrum(sme, linelist_mode='auto') +sme = synthesize_spectrum(sme, linelist_mode='dynamic') ``` -The detailed explanaiton of what `linelist_mode='auto'` triggers is as follow: +`linelist_mode='auto'` is still accepted as a deprecated compatibility alias +and maps to `linelist_mode='dynamic'`. + +The detailed explanaiton of what `linelist_mode='dynamic'` triggers is as follow: 1. If the `sme.linelist.cdr_paras` is None (means `update_cdr` hasn't run for it), or either the Teff, logg, monh and vmic (optional) in `sme.linelist.cdr_paras` differ from the input parameters by 250K, 0.5, 0.5 or 1 (you can change these parameter in `sme.linelist.cdr_paras_thres` as a dictionary), then run `update_cdr` under current input stellar parameters. 2. Find the beginning and ending wavelength of each segment, wbeg and wend. @@ -289,7 +292,7 @@ sme.spec = obs_flux sme.uncs = obs_flux_err sme.linelist = line_list sme.cdr_depth_thres = 0.1 -sme = solve(sme, linelist_mode='auto') +sme = solve(sme, linelist_mode='dynamic') ``` -It only triggers the `linelist_mode` variable in `synthesize_spectrum`. \ No newline at end of file +It only triggers the `linelist_mode` variable in `synthesize_spectrum`. diff --git a/docs/advance/index.rst b/docs/advance/index.rst new file mode 100644 index 00000000..9eb2fea2 --- /dev/null +++ b/docs/advance/index.rst @@ -0,0 +1,14 @@ +Advanced usage +======== + +This section collects advanced PySME functions, including derived-parameter +fitting, line-filtering strategies for performance, and developer-oriented notes. + +.. toctree:: + :maxdepth: 1 + + derived_param.md + line_filtering.md + how-to.md + fordev.md + faq.md \ No newline at end of file diff --git a/docs/advance/line_filtering.md b/docs/advance/line_filtering.md new file mode 100644 index 00000000..8c249af1 --- /dev/null +++ b/docs/advance/line_filtering.md @@ -0,0 +1,67 @@ +# Line Filtering + +For wide wavelength coverage (or many segments), using the full line list in +every segment is expensive. PySME provides dynamic line filtering to keep only +relevant lines per segment (see Jian et al. in prep). + +## Core Option + +Use `linelist_mode` in synthesis or solve: + +- `"all"`: use all lines (default). +- `"dynamic"`: filter lines by precomputed line properties (recommended for long spectra). +- `"auto"`: legacy alias of `"dynamic"` (deprecated). + +## How Dynamic Filtering Works + +When `linelist_mode="dynamic"`: + +1. PySME ensures line metadata exists (`central_depth`, `line_range_s`, `line_range_e`). + If missing or stale, it updates them via `update_cdr(...)`. +2. For each segment, PySME keeps lines that overlap the segment range + (with broadening margin) and pass a strength threshold. +3. Only this reduced line subset is sent to SMElib for that segment. + +This can significantly reduce runtime for long or segmented spectra. + +## Main Controls + +- `sme.cdr_depth_thres`: minimum line-strength threshold used in filtering. +- `sme.cdr_N_line_chunk`: chunk size used in `update_cdr`. +- `sme.cdr_parallel`: enable/disable parallel `update_cdr`. +- `sme.cdr_n_jobs`: number of parallel jobs. +- `cdr_database` / `cdr_create` (function args): reuse or build a CDR grid on disk. + +## Example 1: Dynamic Filtering in Synthesis + +```py +from pysme.synthesize import Synthesizer, synthesize_spectrum + +synth = Synthesizer() +sme = synth.update_cdr(sme) # populate central_depth / line_range_* once +sme.cdr_depth_thres = 0.02 # keep only stronger lines + +sme = synthesize_spectrum(sme, linelist_mode="dynamic") +``` + +## Example 2: Dynamic Filtering in Solve + +```py +from pysme.solve import solve + +fit = ["teff", "logg", "monh", "vmic"] +sme = solve( + sme, + fit, + linelist_mode="dynamic", + cdr_database="path/to/cdr_grid", # optional on-disk CDR cache/grid + cdr_create=False, # set True to force regeneration +) +``` + +## Practical Guidance + +- Start with `sme.cdr_depth_thres = 0.0` and increase gradually if needed. +- Use `"all"` for short, narrow windows where filtering overhead may not help. +- Use `"dynamic"` for wide ranges or many segments. + diff --git a/docs/usage/abundance.md b/docs/concepts/abundance.md similarity index 96% rename from docs/usage/abundance.md rename to docs/concepts/abundance.md index 4bd86f97..77198a75 100644 --- a/docs/usage/abundance.md +++ b/docs/concepts/abundance.md @@ -58,7 +58,7 @@ values of H, He, and Li are approximately 3.16e4, 2.69e3, and ## Solar metallicity -PySME contains three pre defined sets of solar abundances, +PySME contains multiple pre defined sets of solar abundances for you to choose from. They are: From solar photosphere: @@ -66,7 +66,7 @@ From solar photosphere: - `grevesse1996`: [Grevesse, Noels & Sauval (1996, ASPC)](https://ui.adsabs.harvard.edu/abs/1996ASPC...99..117G) - `grevesse1998`: [Grevesse, & Sauval (1998, SSRv)](https://ui.adsabs.harvard.edu/abs/1998SSRv...85..161G) - `asplund2005`: [Asplund, Grevesse & Sauval (2005, ASPC)](https://ui.adsabs.harvard.edu/abs/2005ASPC..336...25A) -- `grevesse2007`: [Grevesse, Asplund & Sauval (2007, SSRv)](https://ui.adsabs.harvard.edu/abs/2007SSRv..130..105G) +- `grevesse2007` (alias: `solar`): [Grevesse, Asplund & Sauval (2007, SSRv)](https://ui.adsabs.harvard.edu/abs/2007SSRv..130..105G) - `asplund2009`: [Asplund, Grevesse & Sauval (2009, ARA&A)](https://ui.adsabs.harvard.edu/abs/2009ARA&A..47..481A) - `asplund2021`: [Asplund, Amarsi & Grevesse (2021, A&A)](https://ui.adsabs.harvard.edu/abs/2021A&A...653A.141A) @@ -116,4 +116,4 @@ thus - $K_\mathrm{H} = \left(1 + \sum_\mathrm{X}10^{H_\mathrm{X}-12} \right)^{-1}$ - $K_\mathrm{He} = \left(10^{12-H_\mathrm{He}} + 1 + 10^{12-H_\mathrm{He}}\times\sum_\mathrm{X}10^{H_\mathrm{X}-12} \right)^{-1}$ -- $K_\mathrm{X} = H_\mathrm{X}-12-\log{(1+10^{H_\mathrm{He}-12})}$ \ No newline at end of file +- $K_\mathrm{X} = H_\mathrm{X}-12-\log{(1+10^{H_\mathrm{He}-12})}$ diff --git a/docs/usage/atmosphere.md b/docs/concepts/atmosphere.md similarity index 87% rename from docs/usage/atmosphere.md rename to docs/concepts/atmosphere.md index 669889fc..c3b56ee6 100644 --- a/docs/usage/atmosphere.md +++ b/docs/concepts/atmosphere.md @@ -47,7 +47,27 @@ The atmopshere object has the following fields: |`wlstd`|Wavelength for continuum optical depth scale.|Å, Default value: 5000Å| |`opflag`|Flags that indicate whether to enable various opacity packages during the radiative transfer calculation|| -## Grid atmospheres +## Atmosphere grids: + +- recommended: + - marcs2012.sav [(Gustafsson et al. 2008)](https://ui.adsabs.harvard.edu/abs/2008A%26A...486..951G) + - marcs2012p_t0.0.sav + - marcs2012p_t1.0.sav + - marcs2012p_t2.0.sav + - marcs2012s_t1.0.sav + - marcs2012s_t2.0.sav + - marcs2012s_t5.0.sav + - marcs2012t00cooldwarfs.sav + - marcs2012t01cooldwarfs.sav + - marcs2012t02cooldwarfs.sav + +- deprecated: + - atlas12.sav + - atlas9_vmic0.0.sav + - atlas9_vmic2.0.sav + - ll_vmic2.0.sav + +### Grid plots ![](../img/atmosphere/marcs2012_grid.png) ![](../img/atmosphere/marcs2012p_t0.0_grid.png) diff --git a/docs/concepts/data_files.md b/docs/concepts/data_files.md new file mode 100644 index 00000000..c96d3d4b --- /dev/null +++ b/docs/concepts/data_files.md @@ -0,0 +1,20 @@ +Data Files You Need +=================== + +PySME uses local cache data under `~/.sme/`. + +Typical files: + +- Atmosphere grids +- NLTE grids +- Optional custom resources + +```{admonition} Accessing data files +The atmosphere and nlte data files should be downloaded from the server automatically when used, so network connection is required when using PySME (not only during installation). +These files can also be downloaded manually by: +- Download data files as part of IDL SME from [here](http://www.stsci.edu/~valenti/sme.html). +- copy them into their respective storage locations in ~/.sme/atmospheres and ~/.sme/nlte_grids + - atmospheres: everything from SME/atmospheres +- Download the nlte_grids in [zenodo](https://doi.org/10.5281/zenodo.3888393). +- The files (mainly atmosphere models and NLTE departure coefficnent grids) required by PySME will be saved inside `~/.sme/`. These files can be large thus if your home directory is small, we recommend to create a softlink for `~/.sme`. +``` diff --git a/docs/usage/fitparameters.md b/docs/concepts/fitparameters.md similarity index 94% rename from docs/usage/fitparameters.md rename to docs/concepts/fitparameters.md index 85af5122..09e47987 100644 --- a/docs/usage/fitparameters.md +++ b/docs/concepts/fitparameters.md @@ -45,3 +45,4 @@ It is therefore recommended to use `'cscale_flag'` and `'vrad_flag'` to specify the desired fitting method. If `'vrad'` is passed as a fitparameter it is equivalent to `'vrad_flag' = 'each'`, and if `'cont'` is passed it is the same as `'cscale_flag' = 'linear'`. +See the [radial velocity](../fundamentals/rv.md) and [continuum](../fundamentals/continuum.md) page for more information. \ No newline at end of file diff --git a/docs/usage/fitresults.md b/docs/concepts/fitresults.md similarity index 99% rename from docs/usage/fitresults.md rename to docs/concepts/fitresults.md index 90670116..803d0a02 100644 --- a/docs/usage/fitresults.md +++ b/docs/concepts/fitresults.md @@ -35,7 +35,8 @@ We estimate the cumulative distribution function of the **generalized** normal d #normalize y /= y[-1] ``` -![](https://upload.wikimedia.org/wikipedia/commons/thumb/c/ca/Normal_Distribution_CDF.svg/500px-Normal_Distribution_CDF.svg.png) +![](https://upload.wikimedia.org/wikipedia/commons/thumb/c/ca/Normal_Distribution_CDF.svg/500px-Normal_Distribution_CDF.svg.png). + (Example of the cumulative distribution function for various values of sigma) diff --git a/docs/concepts/index.rst b/docs/concepts/index.rst new file mode 100644 index 00000000..93d419d2 --- /dev/null +++ b/docs/concepts/index.rst @@ -0,0 +1,17 @@ +Concepts +======== + +Core model and data concepts behind PySME. + +.. toctree:: + :maxdepth: 1 + + sme_struct.md + abundance.md + atmosphere.md + nlte.md + fitparameters.md + fitresults.md + lsf.md + data_files.md + system_info.md diff --git a/docs/concepts/lfs.md b/docs/concepts/lfs.md new file mode 100644 index 00000000..9171cb2f --- /dev/null +++ b/docs/concepts/lfs.md @@ -0,0 +1,6 @@ +# Large File Server + +PySME does not come with all of the large atmosphere or NLTE grids +as part of the distribution. Instead Uppsala University provides +a server that serves the files when needed. Simply specify one of +the available filenames and PySME will fetch the file on your next run. \ No newline at end of file diff --git a/docs/concepts/nlte.md b/docs/concepts/nlte.md new file mode 100644 index 00000000..59e0ee68 --- /dev/null +++ b/docs/concepts/nlte.md @@ -0,0 +1,128 @@ +# NLTE grids + +Non Local Thermal Equilibrium (NLTE) calculations are important to accuarately fit certain lines. +PySME supports them using pre-computed grids of NLTE departure coefficients, which need to be created for every element. +For common elements PySME provides grids (see below) via the LFS (see [](lfs.md)). +If any of these grids are used, please kindly take care to cite the papers describing the NLTE models and departure coefficient calculations. + +NLTE calculations need to be specified for each element they are supposed to be used for individually using `sme.nlte.set_nlte(el, grid)` (the `grid` can be omitted if there is a grid in lfs). +Similarly they can be disabled for each element using `sme.nlte.remove_nlte(el)`, where sme is your SME structure. +If no element is set to NLTE in the structure PySME will perform +LTE calculations only. + +## Fields of NLTE object + +- `elements`: The elements for which NLTE has been activated +- `grids`: The grid file that is used for each active element +- `subgrid_size`: A small segment of the NLTE grid will be cached in memory + to speed up calculations. This sets the size of that cache + by defining the number of points in each + axis (rabund, teff, logg, monh). +- `flags`: After the synthesis all lines are flaged if they used NLTE + +## Grid interpolation + +The grid has 6 dimensions. + +- teff: Effective Temperature +- logg: Surface Gravity +- monh: Overall Metallicity +- rabund: relative abundance of that element +- depth: optical depth in the atmosphere +- departure coefficients: + The NLTE departure coefficients describing how much + it varies from the LTE calculation + +We then perform linear interpolation to the stellar parameters +we want to model. And we then perform a cubic spline fit to the depth scale +of the model atmosphere we specified (See [](atmosphere)). + +We then use the linelist to find only the relevant transitions in the grid, +and pass the departure coefficients for each line to the C library. + +## NLTE flags in line list + +PySME provides information on whehter a line is synthesized in NLTE through the `nlte_flag` column in the line list. + +- `1`: this line was synthesized with in NLTE. +- `0`: this line was synthesized in LTE. +- `-1`: this line was not included in the current synthesis pass (see [](../advance/line_filtering.md)). + +## Recommended and default grids + +| Element | Grid name | Zenodo version | Citation | +|:---:|:---:|:---:|:---:| +| H | `nlte_H_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Li | `nlte_Li_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| C | `nlte_C_pysme.grd` | 2 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| N | `nlte_N_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| O | `nlte_O_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Na | `nlte_Na_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Mg | `nlte_Mg_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Al | `nlte_Al_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Si | `nlte_Si_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| S | `nlte_Si_pysme.grd` | 7 | [Amarsi et al. 2025](https://ui.adsabs.harvard.edu/abs/2025A%26A...703A..35A/abstract) | +| K | `nlte_K_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Ca | `nlte_Ca_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Ti | `nlte_Ti_pysme.grd` | 5 | [Mallinson et al. 2024](https://ui.adsabs.harvard.edu/abs/2024A%26A...687A...5M) | +| Mn | `nlte_Mn_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Fe | `nlte_Fe_pysme.grd` | 4 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | +| Cu | `nlte_Cu_pysme.grd` | 6 | [Caliskan et al. 2025](https://ui.adsabs.harvard.edu/abs/2025A%26A...696A.210C/abstract) | +| Ba | `nlte_Ba_pysme.grd` | 3 | [Amarsi et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) | + +All recommended/default grids are available in [Zenodo](https://zenodo.org/records/17064337). + +## Deprecated grids + + - H + - marcs2012_H2018.grd + - Li + - marcs2012_Li.grd [(Lind et al. 2009)](https://ui.adsabs.harvard.edu/abs/2009A%26A...503..541L) + - marcs2012_Li2009.grd [(Lind et al. 2009)](https://ui.adsabs.harvard.edu/abs/2009A%26A...503..541L) + - nlte_Li_multi.grd + - C + - marcs2012_C.grd + - N + - marcs2012_N.grd + - O + - marcs2012p_t1.0_O.grd [(Sitnova et al. 2013)](https://ui.adsabs.harvard.edu/abs/2013AstL...39..126S) + - marcs2012_O2015.grd [(Amarsi et al. 2016b)](https://ui.adsabs.harvard.edu/abs/2016MNRAS.455.3735A) + - marcs2012_O.grd + - Na + - marcs2012p_t1.0_Na.grd [(Mashonkina et al. 2001)](https://ui.adsabs.harvard.edu/abs/2000ARep...44..790M) + - marcs2012_Na.grd [(Lind et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011A%26A...528A.103L) + - marcs2012_Na2011.grd [(Lind et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011A%26A...528A.103L) + - nlte_Na_multi_full.grd + - nlte_Na_multi_sun.grd + - Al + - marcs2012_Al.grd + - marcs2012_Al2017.grd + - Mg + - marcs2012_Mg2016.grd [(Osorio et al. 2016)](https://ui.adsabs.harvard.edu/abs/2016A%26A...586A.120O) + - marcs2012_Mg.grd + - Si + - marcs2012_Si2016.grd [(Amarsi & Asplund 2017)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.464..264A) + - marcs2012_Si.grd + - K + - marcs2012_K.grd + - Ca + - marcs2012s_t2.0_Ca.grd [(Mashonkina et al. 2007)](https://ui.adsabs.harvard.edu/abs/2007A%26A...461..261M) + - marcs2012p_t1.0_Ca.grd [(Mashonkina et al. 2017)](https://ui.adsabs.harvard.edu/abs/2007A%26A...461..261M) + - marcs2012_Ca.grd + - Ti + - marcs2012s_t2.0_Ti.grd [(Sitnova et al. 2016)](https://ui.adsabs.harvard.edu/abs/2016MNRAS.461.1000S) + - Mn + - marcs2012_Mn.grd + - Fe + - marcs2012_Fe2016.grd [(Amarsi et al. 2016a)](https://ui.adsabs.harvard.edu/abs/2016MNRAS.463.1518A) + - marcs2012s_t2.0_Fe.grd [(Mashonkina et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011A%26A...528A..87M) + - nlte_Fe_multi_full.grd + - marcs2012_Fe.grd + - Ba + - marcs2012p_t1.0_Ba.grd [(Mashonkina et al. 1999)](https://ui.adsabs.harvard.edu/abs/1999A%26A...343..519M) + - Eu + - nlte_Eu.grd + + + + diff --git a/docs/usage/sme_struct.md b/docs/concepts/sme_struct.md similarity index 81% rename from docs/usage/sme_struct.md rename to docs/concepts/sme_struct.md index 5ada7f1a..8e567cfc 100644 --- a/docs/usage/sme_struct.md +++ b/docs/concepts/sme_struct.md @@ -19,7 +19,6 @@ Stellar parameters describe the star in general and are usually what we want to - `vmac`: The macro-turbulence velocity in km s⁻¹. Describes turbulence on scales **larger** than the photon mean free path, also contributing to broadening. - `mu`: Limb-angle values ($\mu = \cos\theta$) at which the radiative-transfer calculation is performed. $\mu = 1$ corresponds to disk center, $0$ to the limb. - ## Radial velocity and Continuum The radial velocity and continuum shift the continuum @@ -43,17 +42,10 @@ PySME has many options to determine them. - each: Fit each wavelength segment individally - whole: Fit the whole spectrum at once -- `normalize_by_continuum`: - A flag that determines, whether the synthetic flux should be normalized - by the continous intensities or not. As long as you have a normalized - observation this should be True, but if you have a flux calibrated - spectrum this should be set to False. Note that even if this is False, - you can still fit a continuum normally using cscale_flag. -Spectra -------- +## Spectra -Spectra are given as a list of arrays[^#], where each array represents +Spectra are given as a list of arrays[^iliffe], where each array represents one wavelength segment of the spectrum. If there is only one segment, the list will only have one element. For legacy reasons there is also an interface to the 'old' system and names (e.g. smod instead of synth) @@ -74,10 +66,18 @@ from IDL SME. It is recommend however to use the new variables. this if you dont have an observation and dont want to specify the exact wavelength grid of the synthetic observation. Note that this is not an Illiffe vector. +:wint: + Optional wavelength grid passed to SMElib `Transf` as adaptive grid for each segment. + Accepts the same segment-aware input styles as `wave`: + `numpy.ndarray` (single segment), `list` of arrays (multi-segment), + or `Iliffe_vector`. + If provided, `sme.wint[segment]` is used directly for synthesis. + If not provided, PySME can reuse an internal cached adaptive grid + when `reuse_wavelength_grid=True`; otherwise SMElib computes a new + adaptive grid. -Abundance ---------- +## Abundance The individal abundances are stored in a seperate Abundance object, which shares the same metallicity as the overall structure. @@ -85,13 +85,12 @@ For more detailed information see [](abundance). :abund: The abundance object -Linelist --------- +## Linelist The sme structure does contain the whole linelist in the linelist property. For legacy reasons, it also provides direct access to the 'species' and 'atomic' arrays. They refer directly to the linelist however. -For more detailed information see :ref:`linelist`. +For more detailed information see [](../fundamentals/linelist.md). :linelist: The linelist object :species: Names of the species of each spectral line @@ -99,26 +98,23 @@ For more detailed information see :ref:`linelist`. Atomic linelist data with columns "atom_number", "ionization", "wlcent", "excit", "gflog", "gamrad", "gamqst", "gamvw" -Atmosphere ----------- +## Atmosphere Unlike the linelist the atmosphere is stored in an external file, that is only referenced by name in the structure. -For more detailed information see :ref:`atmosphere`. +For more detailed information see [](atmosphere). :atmo: The atmosphere object -NLTE ----- +## NLTE Unlike the linelist, but similar to the atmosphere, the NLTE parameters are stored in external tables, which are only referenced -by name. For more detailed information see :ref:`nlte`. +by name. For more detailed information see [](nlte). :nlte: The NLTE object -Instrument Parameters ---------------------- +## Instrument Parameters PySME can also model instrumental broadening as part of the spectral synthesis. For this you need to specify the resolution @@ -135,32 +131,30 @@ and the broadening method to use. The y points of the instrument profile. Only relevant if iptype is 'table'. -Fitresults ----------- +## Fitresults :fitparameters: The fitparameters used for the fitting. - See :ref:`fitparameters`. -:fitresults: The fitresults object. See :ref:`fitresults`. + See [fitparameters](../concepts/fitparameters.md). +:fitresults: The fitresults object. See [fitresults](../concepts/fitresults.md). -System Information ------------------- +## System Information The sme structure does contain information about the host system. E.g. which operating system was used. This is mostly for legacy reasons, and potential debugging information. -For more information see :ref:`system_info`. +For more information see [system_info](../concepts/system_info.md). :system_info: The system information object. It replaces the idlver object. -Other Parameters ----------------- +## Other Parameters :gam6: van der Waals scaling factor (usually 1) :h2broad: flag determing whether to use H2 broadening or not (usually True) :accrt: Minimum accuracy for synthethized spectrum at wavelength grid - points in sme.wave. Values below 1e-4 are not meaningful + points in `sme.wint` (or SMElib adaptive grid if `sme.wint` is not set). + Values below 1e-4 are not meaningful. :accwi: Minimum accuracy for linear spectrum interpolation vs. wavelength. Values below 1e-4 are not meaningful. @@ -169,4 +163,4 @@ Other Parameters The date and time when this structure or the last synthetic spectrum was created -[^#] They are called Illiffe vectors in the code, and they were that in IDL. But they are technically not Illiffe vectors anymore, but just lists of individal numpy arrays. \ No newline at end of file +[^iliffe]: They are called Illiffe vectors in the code, and they were that in IDL. But they are technically not Illiffe vectors anymore, but just lists of individal numpy arrays. diff --git a/docs/usage/system_info.rst b/docs/concepts/system_info.md similarity index 68% rename from docs/usage/system_info.rst rename to docs/concepts/system_info.md index 961687a1..c07eb179 100644 --- a/docs/usage/system_info.rst +++ b/docs/concepts/system_info.md @@ -1,12 +1,10 @@ -.. _system_info: +# System Info -System Info -=========== The system info object is mostly for debugging, so you probably -dont need to worry about it. +don't need to worry about it. -Still here are the fields +Still here are the fields: :arch: System architecture :os: Operating system @@ -15,5 +13,4 @@ Still here are the fields :release: Python version :build_date: build date of the python version :memory_bits: Platform architecture bit size (32bit or 64bit) -:file_offset_bits: ? -:host: Name of the machine +:host: Name of the machine \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py index 15992621..43dfdb10 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -14,8 +14,10 @@ # import os import sys +import importlib.util -src_dir = os.path.abspath("..") +docs_dir = os.path.abspath(os.path.dirname(__file__)) +src_dir = os.path.abspath(os.path.join(docs_dir, "..")) sys.path.insert(0, src_dir) sys.path.insert(0, src_dir + "/src") @@ -26,10 +28,35 @@ copyright = "2025, Jeff Valenti, Nikolai Piskunov, Mingjie Jian, Ansgar Wehrhahn" author = "Jeff Valenti, Nikolai Piskunov, Mingjie Jian, Ansgar Wehrhahn" -# The short X.Y version -version = "" -# The full version, including alpha/beta/rc tags -release = "0.6.20" +def _resolve_release(): + """Resolve docs version from versioneer without importing pysme package. + + Importing ``pysme`` would trigger heavy runtime side effects in ``__init__``. + We therefore load ``src/pysme/_version.py`` directly. + """ + version_file = os.path.join(src_dir, "src", "pysme", "_version.py") + try: + spec = importlib.util.spec_from_file_location("pysme_version", version_file) + if spec is None or spec.loader is None: + raise RuntimeError(f"Could not load version spec from: {version_file}") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module.get_versions().get("version", "unknown") + except Exception: + # Fallbacks keep docs buildable even without full VCS metadata. + return ( + os.environ.get("READTHEDOCS_GIT_IDENTIFIER") + or os.environ.get("READTHEDOCS_VERSION") + or "unknown" + ) + + +# The full version, including local build metadata (e.g. +g.dirty) +release = _resolve_release() +# The short X.Y version shown by Sphinx. +short_release = release.split("+")[0] +version_parts = short_release.split(".") +version = ".".join(version_parts[:2]) if len(version_parts) >= 2 else short_release # -- General configuration --------------------------------------------------- @@ -82,7 +109,7 @@ # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. -language = None +language = "en" # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. @@ -118,6 +145,7 @@ "image_dark": "_static/pysme_logo-dark.png", } } +html_css_files = ["custom.css"] # Custom sidebar templates, must be a dictionary that maps document names # to template names. diff --git a/docs/usage/changelog.md b/docs/dev/changelog.md similarity index 81% rename from docs/usage/changelog.md rename to docs/dev/changelog.md index 018a27f1..b830fe70 100644 --- a/docs/usage/changelog.md +++ b/docs/dev/changelog.md @@ -5,6 +5,9 @@ This page stores the change log for pysme since May 2024. ## In-development - (model) modify the atmosphere geometry to '' for auto desicion instead of 'PP'. +- (sme structure) add `sme.wint` as a segment-aware optional transfer grid input for synthesis. +- (synthesize) when `sme.wint` is provided, synthesis now prefers user-provided `wint` over internal cached wavelength grids. +- (docs) update long-spectrum examples to prefer `linelist_mode='dynamic'` (keep `auto` as deprecated compatibility alias). ## In github-repo @@ -42,4 +45,4 @@ This page stores the change log for pysme since May 2024. - Add line depth result to sme_structure. - Add line range result to sme_structure. - Support python 3.12. -- Support new VALD3 line format. \ No newline at end of file +- Support new VALD3 line format. diff --git a/docs/usage/changes.rst b/docs/dev/changes.rst similarity index 97% rename from docs/usage/changes.rst rename to docs/dev/changes.rst index 0ea7a288..e9e5e84a 100644 --- a/docs/usage/changes.rst +++ b/docs/dev/changes.rst @@ -36,8 +36,7 @@ SME structure: Abundance: * Python Class object - * set of three default solar abundances available - ("asplund2009", "grevesse2007", "lodders2003") + * set of multiple solar abundances available * can output abundances in different formats ('sme', 'n/nTot', 'n/nH', and 'H=12') * internal format is 'H=12', i.e. log(X/H) + 12 diff --git a/docs/dev/index.rst b/docs/dev/index.rst new file mode 100644 index 00000000..97361388 --- /dev/null +++ b/docs/dev/index.rst @@ -0,0 +1,12 @@ +Development +=========== + +Developer-facing guides and maintenance notes. + +.. toctree:: + :maxdepth: 1 + + changelog.md + changes + old_docs.md + \ No newline at end of file diff --git a/docs/dev/old_docs.md b/docs/dev/old_docs.md new file mode 100644 index 00000000..931701cc --- /dev/null +++ b/docs/dev/old_docs.md @@ -0,0 +1,4 @@ +# Legacy SME Documentations + +- [IDLSME structure document (v5)](../_static/IDLSME_592_Struct-v5.pdf) +- [SME abundances (2012)](../_static/SME_abundances_2012.pdf) diff --git a/docs/fundamentals/abund.md b/docs/fundamentals/abund.md new file mode 100644 index 00000000..50b91189 --- /dev/null +++ b/docs/fundamentals/abund.md @@ -0,0 +1,69 @@ +# Abundance + +`sme.abund` stores elemental abundances, while `sme.monh` stores the global metallicity offset. +In PySME, these two are combined at synthesis time. + +## Core rule in PySME + +For elements heavier than He, the abundance used by SMElib is: + +``` +effective abundance = pattern abundance + sme.monh +``` + +Hydrogen and helium are not shifted by `monh`. + +## What this means when you modify abundances + +When you update one element in `sme.abund`, you are writing the **pattern** +value (before the global metallicity shift). +So if you want a target final abundance `A_target` (for example in `H=12`), use: + +``` +pattern_value = A_target - sme.monh +``` + +then assign that pattern value. + +```{warning} +Note that when printing out `abund`, `monh` is applied to the abundance values. +``` + +## Practical examples + +```py +from pysme.abund import Abund + +sme.abund = Abund.solar() +sme.monh = -0.20 +``` + +### 1. Set a target absolute abundance (H=12) + +If you want final `A(Mg) = 7.40`: + +```py +target = 7.40 +sme.abund["Mg"] = target - sme.monh +``` + +Because PySME adds `monh` later, this produces the intended final Mg abundance. + +### 2. Keep scaled-solar composition, only change metallicity + +```py +sme.monh = -0.50 +``` + +This shifts all metals (Z > 2) together by `-0.50 dex`, while keeping abundance ratios fixed. + +### 3. Set an element enhancement at fixed metallicity + +If you want `[Mg/Fe] = +0.30` at current `monh`, add `+0.30` to Mg in the +pattern relative to your base pattern. + + diff --git a/docs/fundamentals/continuum.md b/docs/fundamentals/continuum.md new file mode 100644 index 00000000..a90ef9c1 --- /dev/null +++ b/docs/fundamentals/continuum.md @@ -0,0 +1,12 @@ +# Continuum + +In PySME, the continuum level of the synthetic spectrum can be adjusted to match the observed spectrum using `cscale` and `cscale_flag`. + +- `cscale`: Polynomial coefficients (per segment) applied to the synthetic spectrum. + The polynomial is evaluated on the observed wavelength grid after shifting the first point to zero, i.e., `f(wave - wave[0])`. +- `cscale_flag`: Controls whether and how continuum correction is applied. + - `none`: No continuum correction. + - `fix`: Use the current `cscale` values without updating them. + - `constant`: Scale by a constant factor. + - `linear`: First-order polynomial (straight line). + - `quadratic`: Second-order polynomial. diff --git a/docs/fundamentals/flux_inten.md b/docs/fundamentals/flux_inten.md new file mode 100644 index 00000000..708e06bb --- /dev/null +++ b/docs/fundamentals/flux_inten.md @@ -0,0 +1,37 @@ +# Flux and intensity + +This page summarizes two switches that control how PySME outputs synthetic spectra. + +## `normalize_by_continuum` + +`normalize_by_continuum` controls whether the synthetic spectrum is divided by +the synthetic continuum: + +- `True` (default): output a continuum-normalized spectrum (usually for normalized observations) +- `False`: keep flux-level output (usually for flux-calibrated observations) + +Even when this is `False`, continuum fitting can still be enabled through +`cscale_flag`. + +## `specific_intensities_only` + +`specific_intensities_only` controls whether PySME returns angle-dependent +specific intensities or disk-integrated flux: + +- `False` (default): integrate specific intensities over the stellar disk and return flux +- `True`: return specific intensities directly, without flux-level post-processing + - Note that in this case, `synthesize_spectrum` will only return the specific intensities, not the `sme` object. + +This is useful when you want the radiative-transfer output itself (as a function +of angle), rather than only the final integrated spectrum. + +## Typical usage + +- Normalized stellar spectrum fitting: + - `normalize_by_continuum = True` + - `specific_intensities_only = False` +- Flux-calibrated analysis: + - `normalize_by_continuum = False` + - `specific_intensities_only = False` +- Intensity-level diagnostics / custom integration: + - `specific_intensities_only = True` diff --git a/docs/fundamentals/index.rst b/docs/fundamentals/index.rst new file mode 100644 index 00000000..4a265ade --- /dev/null +++ b/docs/fundamentals/index.rst @@ -0,0 +1,17 @@ +Fundamentals and usage +================ + +This section introduces the fundamentals you need before running synthesis or fitting in PySME. +It covers the key data structures and controls that define spectra, line selection, and model behavior. + +.. toctree:: + :maxdepth: 1 + + mask.md + rv.md + continuum.md + linelist.md + abund.md + segment.md + flux_inten.md + io.md diff --git a/docs/fundamentals/io.md b/docs/fundamentals/io.md new file mode 100644 index 00000000..2f38a99a --- /dev/null +++ b/docs/fundamentals/io.md @@ -0,0 +1,10 @@ +# Save and Load + +The SME structure can be loaded with +```py +sme = SME_Structure.load("in.sme") +``` +or saved with +```py +sme.save("out.sme") +``` \ No newline at end of file diff --git a/docs/fundamentals/linelist.md b/docs/fundamentals/linelist.md new file mode 100644 index 00000000..8ca4972f --- /dev/null +++ b/docs/fundamentals/linelist.md @@ -0,0 +1,142 @@ +# Linelist + +`linelist` is one of the core inputs for spectral synthesis. It stores the +physical parameters of each transition (element/species, wavelength, `log gf`, +damping parameters, etc.). + +## Getting a Linelist from VALD + +PySME is designed to work well with all types of line lists in [VALD](https://vald.astro.uu.se/) database. +A common workflow is: + +1. Download a line list from VALD (typically from `extract stellar`). +2. Load it in PySME with `ValdFile`. + +```py +from pysme.linelist.vald import ValdFile + +vald = ValdFile("my_linelist.lin") +sme.linelist = vald +``` + +## Linelist is essentially a DataFrame + +Internally, PySME stores `LineList` data in a `pandas.DataFrame` +(`sme.linelist._lines`). +So you can use DataFrame-style operations directly: filtering, sorting, +editing columns, adding flags, and so on. + +```py + +# Filter by wavelength range +sub = sme.linelist[(sme.linelist["wlcent"] > 6436.0) & (sme.linelist["wlcent"] < 6444.0)] + +# Select one species (for example, Fe I) +fe1 = sme.linelist[sme.linelist["species"] == "Fe 1"] +``` + +## Short vs Long Format + +VALD provides two line list formats: + +- short: fewer fields, suitable for LTE. +- long: includes additional upper/lower level and quantum-number information. + +If you plan to run NLTE, you should use **VALD long format**. +In PySME, NLTE relies on the extra level information available in long format; +short format does not provide enough information for that workflow. + +## Line parameters + +The short format fields are + +:`species`: + A string identifier including the element + and ionization state or the molecule +:`wlcent`: The central wavelength of the line in Angstrom +:`gflog`: + The log of the product of the statistical weight of + the lower level and the oscillator strength for the transition. +:`excit`: The excitation energy in the lower level +:`ionization`: The ionization state of the species, where 1 is neutral +:`gamrad`: The radiation broadening parameter +:`gamqst`: Stark broadening parameter +:`gamvw`: van der Waals broadening parameter +:`lande`: The lande factor +:`depth`: An arbitrary depth estimation of the line +:`reference`: A citation where this data came from +:`atom_number`: + Identifies the species by the atomic number + (i.e. the number of protons) + +In addition the long format has the following fields + +:`lande_lower`: The lower Lande factor +:`lande_upper`: The upper Lande factor +:`j_lo`: The spin of the lower level +:`j_up`: The spin of the upper level +:`e_upp`: The energy of the upper level +:`term_lower`: The electron configuration of the lower level +:`term_upper`: The electron configuration of the upper level +:`error`: An uncertainty estimate for this linedata + +### Runtime columns + +:`nlte_flag`: Per-line NLTE usage flag written after synthesis. See [NLTE flags in line list](../concepts/nlte.md#nlte-flags-in-line-list). +:`central_depth`: Estimated central line depth used for line selection/filtering. See [Line Filtering](../advance/line_filtering.md). +:`line_range_s`: Start wavelength of the estimated line-contribution range. See [Line Filtering](../advance/line_filtering.md). +:`line_range_e`: End wavelength of the estimated line-contribution range. See [Line Filtering](../advance/line_filtering.md). + + +### Important Note + +As far as the radiative transfer code is concerned the ionization is defined as part of the species term. +I.e. a line with species = "Fe 2" will be calculated as ionization = 2. +If no number is set within the species field, SME will use an ionization of 1. + +Also atom_number is ignored in the radiative transfer calculations and therefore does not need to be set. + +## Example Linelist + +Below is a typical `sme.linelist` excerpt (middle rows omitted): + +```text + species wlcent gflog excit j_lo e_upp j_up lande_lower \ +0 Fe 1 6436.4055 -2.460 4.186364 2.0 6.112128 1.0 0.68 +1 Eu 2 6437.6106 -1.242 1.319612 5.0 3.245008 5.0 1.73 +2 Eu 2 6437.6121 -1.280 1.319636 5.0 3.245027 5.0 1.73 +3 Eu 2 6437.6135 -2.473 1.319612 5.0 3.245008 5.0 1.73 +4 Eu 2 6437.6194 -2.511 1.319636 5.0 3.245027 5.0 1.73 +.. ... ... ... ... ... ... ... ... +62 Mn 1 6443.4689 -2.818 3.772307 1.5 5.695960 2.5 1.21 +63 Mn 1 6443.4726 -3.538 3.772307 1.5 5.695960 2.5 1.21 +64 Mn 1 6443.4821 -2.275 3.772307 1.5 5.695960 2.5 1.21 +65 Mn 1 6443.4887 -2.964 3.772307 1.5 5.695960 2.5 1.21 +66 Mn 1 6443.4938 -3.918 3.772307 1.5 5.695960 2.5 1.21 + + lande_upper lande ... term_lower couple_upper \ +0 0.56 0.73 ... 3d8 c3F LS +1 1.79 1.76 ... 4f7.(8S).5d a9D* JJ +2 1.79 1.76 ... 4f7.(8S).5d a9D* JJ +3 1.79 1.76 ... 4f7.(8S).5d a9D* JJ +4 1.79 1.76 ... 4f7.(8S).5d a9D* JJ +.. ... ... ... ... ... +62 1.37 1.49 ... 3d5.4s2 b4D LS +63 1.37 1.49 ... 3d5.4s2 b4D LS +64 1.37 1.49 ... 3d5.4s2 b4D LS +65 1.37 1.49 ... 3d5.4s2 b4D LS +66 1.37 1.49 ... 3d5.4s2 b4D LS + + term_upper error nlte_flag atom_number ionization \ +0 3d6.(3F2).4s.4p.(3P*) v3D* 0.50 0 1.0 1.0 +1 4f7.(8S<7/2>).6p<3/2> (7/2,3/2) 1.00 0 1.0 2.0 +2 4f7.(8S<7/2>).6p<3/2> (7/2,3/2) 1.00 0 1.0 2.0 +3 4f7.(8S<7/2>).6p<3/2> (7/2,3/2) 1.00 0 1.0 2.0 +4 4f7.(8S<7/2>).6p<3/2> (7/2,3/2) 1.00 0 1.0 2.0 +.. ... ... ... ... ... +62 3d6.(5D).4p z4D* 0.25 0 1.0 1.0 +63 3d6.(5D).4p z4D* 0.25 0 1.0 1.0 +64 3d6.(5D).4p z4D* 0.25 0 1.0 1.0 +65 3d6.(5D).4p z4D* 0.25 0 1.0 1.0 +66 3d6.(5D).4p z4D* 0.25 0 1.0 1.0 +``` diff --git a/docs/fundamentals/mask.md b/docs/fundamentals/mask.md new file mode 100644 index 00000000..a0452891 --- /dev/null +++ b/docs/fundamentals/mask.md @@ -0,0 +1,18 @@ +# Mask + +Masking the observed spectrum is supported in `solve`. +The code will then allocate the corresponding pixel to different mask, and perform different manupulation to them. + +The mask, `sme.mask` is an array with `dtype=int` with the same length with `sme.wave`. +The types of mask are: + +|mask value|mask name|Description| +|:--:|:--:|:--:| +|0|bad pixel|Pixels excluded| +|1|line pixel|Pixels included in the parameter fitting| +|2|continuum pixel|Pixels included in continuum fitting| +|4|vrad pixel|Pixels included in radial velocity fitting| + +The masks are additive, i.e., you can set mask value to 5 for line and vrad pixel. Only the good pixels will be included in the fit, but the synthetic spectrum will still be calculated. + +In default, mask is all 1. \ No newline at end of file diff --git a/docs/fundamentals/rv.md b/docs/fundamentals/rv.md new file mode 100644 index 00000000..71a6f230 --- /dev/null +++ b/docs/fundamentals/rv.md @@ -0,0 +1,15 @@ +# Radial velocity (vrad) + +`vrad` in PySME is segment-aware and has two functions (with `vrad_flag`): + +- fit RV from observation vs. synthetic spectrum +- apply RV shift to synthetic spectrum before comparison + +## Core parameters + +- `vrad`: radial velocity in km/s for each segment +- `vrad_flag`: controls how RV is determined + - `none`: no RV fitting + - `fix`: use the current `sme.vrad` to shift the synthetic spectra and do not change it + - `each`: fit one RV per segment + - `whole`: fit one shared RV from all selected segments \ No newline at end of file diff --git a/docs/fundamentals/segment.md b/docs/fundamentals/segment.md new file mode 100644 index 00000000..0edb977b --- /dev/null +++ b/docs/fundamentals/segment.md @@ -0,0 +1,71 @@ +# Spectral segment + +A segment is one wavelength chunk of the spectrum. +PySME is segment-aware by design: most spectral arrays are stored as a list-like +object with one entry per segment. + +## Why segments exist + +Segments are useful when: + +- your observation is naturally split into orders/chunks +- different chunks have different wavelength sampling/resolution/radial velocity. +- you want to process only part of the spectrum in synthesis/solve + +## How segments are defined + +PySME can get segment information from either: + +- `sme.wave`: explicit wavelength arrays (recommended when available) +- `sme.wran`: only segment boundaries `[[w0_start, w0_end], [w1_start, w1_end], ...]` + +If `sme.wave` exists, it effectively defines segment boundaries. +If only `sme.wran` is given, PySME synthesizes each segment within those ranges. + +## Segment-aware inputs + +Typical segment-aware fields are: + +- `wave` +- `spec` +- `uncs` +- `mask` +- `synth` +- `cont` +- `wint` (optional transfer grid) +- `vrad` +- `cscale` + +You can pass these as: + +- a single 1D array for one segment +- a list of arrays for multiple segments + +## Parameters that are not segmented + +Some model parameters are global (single value for all segments), e.g.: + +- `vmic` +- `vmac` +- `vsini` + +Some are scalar-or-per-segment, e.g.: + +- `ipres`: can be one value or one value per segment +- `vrad`: handled per segment, depending on `vrad_flag` + +## Choosing which segments to run + +Both `synthesize_spectrum(...)` and `solve(...)` accept a `segments` argument. + +- `segments="all"`: run all segments +- `segments=[0, 2, ...]`: run selected segments only + +Invalid indices raise an error. +Segments that are fully masked as bad pixels are skipped automatically. + +## Practical notes + +- Keep segment ordering consistent across all segment-aware inputs. +- If you provide per-segment arrays, lengths must match within each segment. +- For simple single-chunk workflows, using one 1D `wave` array is enough. diff --git a/docs/getting_started/first_fit.md b/docs/getting_started/first_fit.md new file mode 100644 index 00000000..aa2ac15b --- /dev/null +++ b/docs/getting_started/first_fit.md @@ -0,0 +1,27 @@ +# Fit an observed spectrum + +Assuming that we have an observed spectrum, with its wavelength, normalized flux, and uncertainties, as an array of `wave`, `flux` and `uncertainties`. + +They can be inserted into the SME structure with: +```py +sme = SME_Structure() +sme.teff, sme.logg, sme.monh, sme.vmic, sme.vmac, sme.vsini = 5777, 4.4, 0, 1.09, 4.19, 1.60 +sme.abund = Abund.solar() +sme.linelist = vald +sme.iptype = 'gauss' +sme.ipres = 42000 +sme.wave = wave +sme.spec = flux +sme.uncs = uncertainties +``` + +The new inputs are the [instrument resolution](../concepts/sme_struct.md#instrument-parameters) and [observed spectra](../concepts/sme_struct.md#spectra). + +Then the `solve` function can be used to find the best fit solution: +```py +from pysme.solve import solve +fitparameters = ["teff", "logg", "monh", "abund Mg"] +sme = solve(sme, fitparameters) +``` + +The [fitresults](../concepts/sme_struct.md#fitresults) are stored in `sme.fitresults`. \ No newline at end of file diff --git a/docs/getting_started/first_spectrum.md b/docs/getting_started/first_spectrum.md new file mode 100644 index 00000000..33012188 --- /dev/null +++ b/docs/getting_started/first_spectrum.md @@ -0,0 +1,58 @@ +# Synthesize your first spectrum + +Three parts of information is needed to synthesize a spectrum: +- Stellar atmosphere model; +- Stellar parameters; +- Line list. + +The first two can be handeled in PySME: + +```py +from pysme.sme import SME_Structure +from pysme.abund import Abund + +sme = SME_Structure() +sme.teff, sme.logg, sme.monh = 5700, 4.4, 0 +sme.abund = Abund.solar() +``` + +Here we defined effective temperature `teff`, surface gravity `logg`, metallicity `monh`, and solar [chemical abundance](../concepts/abundance.md) `abund`. +For a full list of stellar parameters, see [](../concepts/sme_struct.md). + +LineList can be download from [VALD database](https://vald.astro.uu.se/). +Here we provide an example line list: [sun.lin](https://raw.githubusercontent.com/MingjieJian/PySME/master/examples/sun.lin). +Load the line list using: +```py +from pysme.linelist.vald import ValdFile +vald = ValdFile("linelist.lin") +sme.linelist = vald +``` +Define wavelength grid or wavelength range of the synthetic spectrum: +```py +import numpy as np +sme.wave = np.arange(6436, 6440, 0.1) +# Or +sme.wran = [6436, 6440] # pysme will choose sampling automatically +``` + +Then use the `synthesize_spectrum` function: +```py +from pysme.synthesize import synthesize_spectrum +sme = synthesize_spectrum(sme) +``` + +The synthesized spectra are stored in `sme.wave` and `sme.synth`: +```py +Iliffe_vector([array([6436. , 6436.1, 6436.2, ..., 6439.8, 6439.9])]) +Iliffe_vector([array([0.99986606, 0.99984309, 0.99980888, ..., 0.99754268, 0.99807395])]) +``` + +and can be plot by (for example): + +```py +import matplotlib.pyplot as plt + +plt.plot(sme.wave[0], sme.synth[0]) +``` + +`[0]` is required because PySME stores spectra as segments; here the result is recognized as segment 0. \ No newline at end of file diff --git a/docs/getting_started/index.rst b/docs/getting_started/index.rst new file mode 100644 index 00000000..36780003 --- /dev/null +++ b/docs/getting_started/index.rst @@ -0,0 +1,14 @@ +Getting Started +=============== + +Start here if you are new to PySME. The goal is to run a full minimal workflow quickly. + +.. toctree:: + :maxdepth: 1 + + installation.md + first_spectrum.md + first_fit.md + nlte.md + .. data_files + .. troubleshooting diff --git a/docs/getting_started/installation.md b/docs/getting_started/installation.md new file mode 100644 index 00000000..9f806b0d --- /dev/null +++ b/docs/getting_started/installation.md @@ -0,0 +1,88 @@ +# Installation + +PySME can be installed through PyPI (recommended; stable release) or from github (the latest version in `develop` branch) directly. + +```{admonition} Supported environments: +- Platforms: Linux, macOS (arm platform) + - The x86/Intel platform macOS is not supported anymore, but you can still install v0.6.23, the last supported PySME version, or download [SMElib](https://github.com/MingjieJian/SMElib) and compile it manually. +- Windows: supported via WSL2 (install/run PySME inside the Linux subsystem) +- Python versions: + - 3.9–3.13. +``` + +## Set up virtual environment + +This step is optional but recommended. + +### conda + +```bash +conda create -n pysme python=3.12 +conda activate pysme +``` + +### venv (alternative) + +```bash +python -m venv .venv +source .venv/bin/activate +``` + +## Install PySME + +### Stable release (recommended) + +```bash +pip install pysme-astro +``` + +### From Github + +#### Clone the repository + +```bash +git clone https://github.com/MingjieJian/SME.git +cd SME +``` + +#### Install PySME from source + +```bash +pip install -U pip +pip install . +``` + +## Verify installation + +```bash +python -c "import pysme; print('PySME version:', pysme.__version__)" +``` + +You should see an output of `PySME version: [version]`. + +## Uninstall + +You can uninstall PySME by: +```sh +pip uninstall pysme-astro +``` + +```{warning} +Note that several files (data file, SMElib file etc) will remain after the uninstall. +They are all list in the output of the pip command, and it is recommended to remove them manually. +``` + +The content below to be removed. + +## Running SME +- An simple minimum example is provided in the [examples directory](https://github.com/MingjieJian/SME/tree/master/examples). Make sure to also download the provided input structure. +- You can then run it with: `python minimum.py` + +```{warning} +PySME requires a **pre-compiled C++/Fortran SME library**. +Wheels are currently provided for Linux and macOS. +For Windows, we recommend using **WSL2** (install PySME inside the Linux environment). + +- PySME requires the pre-compled C++/Fortran SME library to run. Currently we deliver SME library with Linux and Mac version; for Windows users, we recommend to use WSL and hence the Linux version. + +``` \ No newline at end of file diff --git a/docs/getting_started/nlte.md b/docs/getting_started/nlte.md new file mode 100644 index 00000000..b4177500 --- /dev/null +++ b/docs/getting_started/nlte.md @@ -0,0 +1,30 @@ +## NLTE correction + +The non-Local Thermal Equilibrium (NLTE) correction can be applied to the spectrum by setting the `nlte` attribute of the SME structure, before `synthesize_spectrum` or `solve`: +```py +sme.nlte.set_nlte(el) # such as "Ca" +# Or provide your own grid: +sme.nlte.set_nlte(el, "path_to_your_grid") +``` + +The current elements with NLTE correction includes: + +![](../_static/nlte_table.png) + +```{warning} +- Long format VALD linelist is required for NLTE calculations. +If only a short format has been given, then the calculations will only be in LTE. (See [](../fundamentals/linelist.md)) +- The NLTE grids are only compatible with the MARCS model atmosphere. +``` + +```{warning} + +``` + +```{note} +Upon the first use of the NLTE correction, the NLTE grid will be downloaded from the server and this may takes a while, depending on the size of the NLTE grid. +``` + +Similarly they can be disabled for each element using `sme.nlte.remove_nlte(el)`. + +More details on NLTE can be found in the [concept section](../concepts/nlte.md). \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index 4f24a6dd..e2d3cc68 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,7 +1,7 @@ PySME documentation =================== -More than two decades ago `Valenti & Piskunov (1996) `_ developed SME – Spectroscopy Made Easy, a high-precision stellar-spectra synthesis/analysis engine that has powered hundreds of studies. +More than two decades ago `Valenti & Piskunov (1996) `_ developed SME - Spectroscopy Made Easy, a high-precision stellar-spectra synthesis/analysis engine that has powered hundreds of studies. PySME is its modern Python front-end: a wrapper around the original C++/Fortran core that lets you (1) compute accurate, high-resolution synthetic spectra from a linelist + model atmosphere, (2) invert observed spectra to derive stellar parameters, and (3) explore NLTE corrections — all from an interactive notebook or scripted pipeline. The same capabilities make PySME invaluable for exoplanet work, where characterising the host star is essential for understanding its planets. 🎉 PySME (version 0.6.14+) is now fully compatible with the arm64 MacOS! Feel free to try it out. 🎉 @@ -11,29 +11,42 @@ PySME is its modern Python front-end: a wrapper around the original C++/Fortran * Plane-parallel and spherical radiative-transfer engine * LTE & 1-D NLTE line formation with pre-computed grids * Automatic :math:`\chi^2` fitting for :math:`T_\mathrm{eff}`, :math:`\log{g}`, :math:`v_\mathrm{mic}`, [X/Fe] … - * Seamless use of ATLAS, MARCS, Phoenix and PINN model atmospheres and VALD line lists - -.. note:: - - * If you are new to PySME: follow the :doc:`usage/installation` and :doc:`usage/quickstart` to get started. - * If you want to get familiar with PySME: read the :doc:`usage/sme_struct` for detail information on using the code. - * If you are familiar with PySME: check out the :doc:`usage/how-to` and :doc:`usage/changelog` to see the new functions. + * Seamless use of ATLAS and MARCS model atmospheres and VALD line lists + +.. raw:: html + +
+ +

Installation

+

Set up PySME and verify your code before synthesizing.

+ Go to installation + + +

Getting started

+

Run the first spectrum and first fit with a minimal, beginner-friendly path.

+ Open getting started + + +

Advanced usage

+

Learn core usage and structure details for deeper and more controlled use.

+ To advanced usage + +
.. toctree:: :maxdepth: 1 - :caption: Contents: - - Installation - usage/quickstart - usage/sme_struct - Large File Server - PySME how to - usage/faq - usage/system_info - usage/changes - Changelog - For dev - _sources/modules + + getting_started/index + fundamentals/index + advance/index + concepts/index + dev/index + +Citation +~~~~~~~~ + +- Jian et al. (2026; in prep.) +- `Wehrhahn et al. (2023) `_ Indices and tables ~~~~~~~~~~~~~~~~~~ @@ -44,4 +57,4 @@ Indices and tables .. rubric:: Quick links :GitHub repository: https://github.com/MingjieJian/SME -:Issue tracker: https://github.com/MingjieJian/SME/issues \ No newline at end of file +:Issue tracker: https://github.com/MingjieJian/SME/issues diff --git a/docs/reference/api.rst b/docs/reference/api.rst new file mode 100644 index 00000000..140c92a0 --- /dev/null +++ b/docs/reference/api.rst @@ -0,0 +1,7 @@ +API Reference +============= + +.. toctree:: + :maxdepth: 2 + + /_sources/modules diff --git a/docs/reference/cli_or_scripts.rst b/docs/reference/cli_or_scripts.rst new file mode 100644 index 00000000..69a5bd35 --- /dev/null +++ b/docs/reference/cli_or_scripts.rst @@ -0,0 +1,8 @@ +CLI Or Scripts +============== + +PySME is primarily used as a Python library. + +For developer-side scripts and internal workflows, see: + +- :doc:`/usage/fordev` diff --git a/docs/reference/index.rst b/docs/reference/index.rst new file mode 100644 index 00000000..273c55bb --- /dev/null +++ b/docs/reference/index.rst @@ -0,0 +1,10 @@ +Reference +========= + +API and interface reference materials. + +.. toctree:: + :maxdepth: 1 + + api + cli_or_scripts diff --git a/docs/usage/atmosphere.rst b/docs/usage/atmosphere.rst deleted file mode 100644 index 7763650a..00000000 --- a/docs/usage/atmosphere.rst +++ /dev/null @@ -1,53 +0,0 @@ -.. _atmosphere: - -Atmosphere -========== - -For the spectral synthesis PySME needs a model atmosphere -to perform the radiative transfer in. PySME does not come -with a set of atmospheres in each distribution but instead -uses the LFS (See :ref:`lfs`) to fetch only the required -model atmosphere when run. - -If you want to provide your own model atmosphere file, it should be present in `~/.sme/atmospheres/`. - -Each atmosphere model file describes a grid of models, on -which we then linearly interpolate to the desired stellar parameters. -Sometimes we dare extrapolate from this grid as well, but in that case, -we always show a warnning. - -Note that the atmosphere also contains a seperate set of stellar -parameters, which is usually the same as that of the sme structure, -but can be different, if for example the atmosphere is embedded, i.e. -fixed, or has not been calculated yet. - -The atmopshere object has the following fields - -:teff: Effective Temperature in Kelvin -:logg: Surface Gravity in log(cgs) -:monh: Metallicity relative to the individual abundances -:abund: The individual abundances (see :ref:`abund`) -:vsini: Projected Rotational velocity in km/s -:vmic: Microturbulence velocity in km/s -:vmac: Macroturbulence veclocity in km/s -:vturb: Turbulent velocity in km/s -:lonh: ? Metallicity -:source: Filename of the atmosphere grid -:depth: - The depth scale to use for calculations. - Either RHOX or TAU -:interp: - The depth scale to use for interpolation. - Either RHOX or TAU -:geom: - The geometry of the atmopshere. Either Plane - Parallel 'PP' or Spherical 'SPH'. -:method: - The method to use for interpolation. Either - 'grid' for a model grid or 'embedded' if - only a single atmosphere is given. -:rhox: 'Column density' depth scale -:tau: 'Optical depth' depth scale -:temp: Temperature profile -:xna: Number density of atoms, ions, and molecules in each depth -:xne: Number density of electrons in each depth diff --git a/docs/usage/installation.md b/docs/usage/installation.md deleted file mode 100644 index b6cb6275..00000000 --- a/docs/usage/installation.md +++ /dev/null @@ -1,46 +0,0 @@ -# Installation - -PySME can be installed through PyPI (recommended) or from github directly. - -Currently PySME is tested with Python verion 3.9-3.13. - -## (optional) set up virtual environment -- `conda env create --name pysme` - - You can specify python verion using `python=3.12` -- `source activate pysme` - -## Install -- For the "stable" version (recommended): - - `pip install pysme-astro` - -```{warning} -- PySME requires the pre-compled C++/Fortran SME library to run. Currently we deliver SME library with Linux and Mac version; for Windows users, we recommend to use WSL and hence the Linux version. -- The files (mainly atmosphere models and NLTE departure coefficnent grids) required by PySME will be saved inside `~/.sme/`. These files can be large thus if your home directory is small, we recommend to create a softlink for `~/.sme`. -``` - -```{warning} -PySME requires the pre-compled C++/Fortran SME library to run. Currently we deliver SME library with Linux and Mac version; for Windows users, we recommend to use WSL and hence the Linux version. -``` - -## Running SME -- An simple minimum example is provided in the [examples directory](https://github.com/MingjieJian/SME/tree/master/examples). Make sure to also download the provided input structure. -- You can then run it with: `python minimum.py` - -```{admonition} Accessing data files -The atmosphere and nlte data files should be downloaded from the server automatically when used, so network connection is required when using PySME (not only during installation). -These files can also be downloaded manually by: -- Download data files as part of IDL SME from [here](http://www.stsci.edu/~valenti/sme.html). -- copy them into their respective storage locations in ~/.sme/atmospheres and ~/.sme/nlte_grids - - atmospheres: everything from SME/atmospheres -- Download the nlte_grids in [zenodo](https://doi.org/10.5281/zenodo.3888393). -``` - -## Uninstall - -You can uninstall PySME by: -```sh -pip uninstall pysme-astro -``` - -Note that several files (data file, SMElib file etc) will remain after the uninstall. -They are all list in the output of the pip command, and it is recommended to remove them manually. \ No newline at end of file diff --git a/docs/usage/lfs.md b/docs/usage/lfs.md deleted file mode 100644 index ffb69473..00000000 --- a/docs/usage/lfs.md +++ /dev/null @@ -1,112 +0,0 @@ -# Large File Server - -PySME does not come with all of the large atmosphere or NLTE grids -as part of the distribution. Instead Uppsala University provides -a server that serves the files when needed. Simply specify one of -the available filenames and PySME will fetch the file on your next run. - -## Available Files - -### Atmosphere grids: - -- recommended: - - marcs2012.sav [(Gustafsson et al. 2008)](https://ui.adsabs.harvard.edu/abs/2008A%26A...486..951G) - - marcs2012p_t0.0.sav - - marcs2012p_t1.0.sav - - marcs2012p_t2.0.sav - - marcs2012s_t1.0.sav - - marcs2012s_t2.0.sav - - marcs2012s_t5.0.sav - - marcs2012t00cooldwarfs.sav - - marcs2012t01cooldwarfs.sav - - marcs2012t02cooldwarfs.sav - -- deprecated: - - atlas12.sav - - atlas9_vmic0.0.sav - - atlas9_vmic2.0.sav - - ll_vmic2.0.sav - -### NLTE grids: - -- recommended and default grids: - - H - - nlte_H_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Li - - nlte_Li_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - C - - nlte_C_pysme.grd [(Amarsi et al. 2020; version 2)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - N - - nlte_N_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - O - - nlte_O_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Na - - nlte_Na_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Mg - - nlte_Mg_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Al - - nlte_Al_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Si - - nlte_Si_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - K - - nlte_K_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Ca - - nlte_Ca_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Ti - - nlte_Ti_pysme.grd [(Mallinson et al. 2024)](https://ui.adsabs.harvard.edu/abs/2024A%26A...687A...5M) - - Mn - - nlte_Mn_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Fe - - nlte_Fe_pysme.grd [(Amarsi et al. 2020; version 4)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - - Ba - - nlte_Ba_pysme.grd [(Amarsi et al. 2020; version 3)](https://ui.adsabs.harvard.edu/abs/2020A%26A...642A..62A) - -- deprecated: - - H - - marcs2012_H2018.grd - - Li - - marcs2012_Li.grd [(Lind et al. 2009)](https://ui.adsabs.harvard.edu/abs/2009A%26A...503..541L) - - marcs2012_Li2009.grd [(Lind et al. 2009)](https://ui.adsabs.harvard.edu/abs/2009A%26A...503..541L) - - nlte_Li_multi.grd - - C - - marcs2012_C.grd - - N - - marcs2012_N.grd - - O - - marcs2012p_t1.0_O.grd [(Sitnova et al. 2013)](https://ui.adsabs.harvard.edu/abs/2013AstL...39..126S) - - marcs2012_O2015.grd [(Amarsi et al. 2016b)](https://ui.adsabs.harvard.edu/abs/2016MNRAS.455.3735A) - - marcs2012_O.grd - - Na - - marcs2012p_t1.0_Na.grd [(Mashonkina et al. 2001)](https://ui.adsabs.harvard.edu/abs/2000ARep...44..790M) - - marcs2012_Na.grd [(Lind et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011A%26A...528A.103L) - - marcs2012_Na2011.grd [(Lind et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011A%26A...528A.103L) - - nlte_Na_multi_full.grd - - nlte_Na_multi_sun.grd - - Al - - marcs2012_Al.grd - - marcs2012_Al2017.grd - - Mg - - marcs2012_Mg2016.grd [(Osorio et al. 2016)](https://ui.adsabs.harvard.edu/abs/2016A%26A...586A.120O) - - marcs2012_Mg.grd - - Si - - marcs2012_Si2016.grd [(Amarsi & Asplund 2017)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.464..264A) - - marcs2012_Si.grd - - K - - marcs2012_K.grd - - Ca - - marcs2012s_t2.0_Ca.grd [(Mashonkina et al. 2007)](https://ui.adsabs.harvard.edu/abs/2007A%26A...461..261M) - - marcs2012p_t1.0_Ca.grd [(Mashonkina et al. 2017)](https://ui.adsabs.harvard.edu/abs/2007A%26A...461..261M) - - marcs2012_Ca.grd - - Ti - - marcs2012s_t2.0_Ti.grd [(Sitnova et al. 2016)](https://ui.adsabs.harvard.edu/abs/2016MNRAS.461.1000S) - - Mn - - marcs2012_Mn.grd - - Fe - - marcs2012_Fe2016.grd [(Amarsi et al. 2016a)](https://ui.adsabs.harvard.edu/abs/2016MNRAS.463.1518A) - - marcs2012s_t2.0_Fe.grd [(Mashonkina et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011A%26A...528A..87M) - - nlte_Fe_multi_full.grd - - marcs2012_Fe.grd - - Ba - - marcs2012p_t1.0_Ba.grd [(Mashonkina et al. 1999)](https://ui.adsabs.harvard.edu/abs/1999A%26A...343..519M) - - Eu - - nlte_Eu.grd \ No newline at end of file diff --git a/docs/usage/linelist.rst b/docs/usage/linelist.rst deleted file mode 100644 index c62d65a3..00000000 --- a/docs/usage/linelist.rst +++ /dev/null @@ -1,67 +0,0 @@ -.. _linelist: - -Linelist -======== - -The linelist is an important part of every synthetic spectrum, -since it includes all information that is line specific. - -Line format ------------ - -PySME knows to types of linelists. A short format, and a long format. -The difference between the two is the amount of information contained therein. -The short format contains only enough parameters for LTE calculations, -while the long format is required for NLTE calculations. - -Line parameters ---------------- - -The short format fields are - -:species: - A string identifier including the element - and ionization state or the molecule -:atom number: - Identifies the species by the atomic number - (i.e. the number of protons) -:ionization: The ionization state of the species, where 1 is neutral (?) -:wlcent: The central wavelength of the line in Angstrom -:excit: The excitation energy in ? -:gflog: - The log of the product of the statistical weight of - the lower level and the oscillator strength for the transition. -:gamrad: The radiation broadening parameter -:gamqst: A broadening parameter -:gamvw: van der Waals broadening parameter -:lande: The lande factor -:depth: An arbitrary depth estimation of the line -:reference: A citation where this data came from - -In addition the long format has the following fields - -:lande_lower: The lower Lande factor -:lande_upper: The upper Lande factor -:j_lo: The spin of the lower level -:j_up: The spin of the upper level -:e_upp: The energy of the upper level -:term_lower: The electron configuration of the lower level -:term_upper: The electron configuration of the upper level -:error: An uncertainty estimate for this linedata - -Important Note --------------- -As far as the radiative transfer code is concerned the ionization is defined as part of the species term. -I.e. a line with species = "Fe 2" will be calculated as ionization = 2. -If no number is set within the species field, SME will use an ionization of 1. - -Also atom_number is ignored in the radiative transfer calculations and therefore does not need to be set. - -VALD integration ----------------- - -PySME is designed to be used in combination with -VALD3 (https://vald.astro.uu.se/). The easiest way to -get a linelist into PySME is therefore to use VALD -extract stellar, as that can be directly imported -using the ValdFile class. diff --git a/docs/usage/nlte.md b/docs/usage/nlte.md deleted file mode 100644 index 2e7f91d6..00000000 --- a/docs/usage/nlte.md +++ /dev/null @@ -1,66 +0,0 @@ -# NLTE - -Non Local Thermal Equilibrium (NLTE) calculations are important to accuarately fit certain lines. -PySME supports them using pre-computed grids of NLTE departure coefficients, which need to be created for -every element. For common elements PySME provides grids via the LFS -(see [](lfs)). If any of these grids are used, please kindly take care to cite the papers describing the NLTE models -and departure coefficient calculations. - -If you want to provide your own NLTE grid files, they should be present in `~/.sme/nlte_grids`. - -NLTE calculations need to be specified for each element they are -supposed to be used for individually using `sme.nlte.set_nlte(el, grid)` (the `grid` can be omitted if there is a grid in lfs). -Similarly they can be disabled for each element using -`sme.nlte.remove_nlte(el)`, where sme is your SME structure. -If no element is set to NLTE in the structure PySME will perform -LTE calculations only. - -```{warning} -Long format VALD linelist is required for NLTE calculations. -If only a short format has been given, then the calculations will -only be in LTE as well. (See [](linelist)) -``` - -Currently PySME supports the following NLTE grids in default: - -```{raw} html - -``` - -```{warning} -The NLTE grids are only compatible with the MARCS model atmosphere. -``` - -The NLTE object has the following fields - -- `elements`: The elements for which NLTE has been activated -- `grids`: The grid file that is used for each active element -- `subgrid_size`: A small segment of the NLTE grid will be cached in memory - to speed up calculations. This sets the size of that cache - by defining the number of points in each - axis (rabund, teff, logg, monh). -- `flags`: After the synthesis all lines are flaged if they used NLTE - - -## Grid interpolation - - -The grid has 6 dimensions. - -- teff: Effective Temperature -- logg: Surface Gravity -- monh: Overall Metallicity -- rabund: relative abundance of that element -- depth: optical depth in the atmosphere -- departure coefficients: - The NLTE departure coefficients describing how much - it varies from the LTE calculation - -We then perform linear interpolation to the stellar parameters -we want to model. And we then perform a cubic spline fit to the depth scale -of the model atmosphere we specified (See [](atmosphere)). - -We then use the linelist to find only the relevant transitions in the grid, -and pass the departure coefficients for each line to the C library. \ No newline at end of file diff --git a/docs/usage/quickstart.md b/docs/usage/quickstart.md deleted file mode 100644 index 094b4d58..00000000 --- a/docs/usage/quickstart.md +++ /dev/null @@ -1,115 +0,0 @@ -# Quickstart - -## SME Structure - -The first step in each SME project is to create an SME structure - -```py -from pysme.sme import SME_Structure -``` - -This can be done in done in a few different ways: -- assign values manually: `sme = SME_Structure()` -- load an existing SME save file (from Python or IDL): `sme = SME_Structure.load("sme.inp")` -- load an .ech file spectrum: `sme = SME_Structure.load("obs.ech")` - - -## Synthesize a spectrum - -Once the SME structure is created, it can be used to synthesize a spectrum. -First we need to define necessary parameters: -* Stellar parameters (effective temperature `teff`, surface gravity `logg`, metallicity `monh`, chemical abundance `abund`) -```py -from pysme.abund import Abund -sme.teff, sme.logg, sme.monh = 5700, 4.4, -0.1 -sme.abund = Abund.solar() -``` -* LineList (linelist), e.g. from [VALD database](https://vald.astro.uu.se/) -```py -from pysme.linelist.vald import ValdFile -vald = ValdFile("linelist.lin") -sme.linelist = vald -``` -* Wavelength grid or wavelength range -```py -sme.wave = [np.arange(6436, 6440, 0.1), np.arange(6442, 6443, 0.1)] -# Or -sme.wran = [[6436, 6440], [6442, 6443]] -``` - -Then use the `synthesize_spectrum` function: -```py -from pysme.synthesize import synthesize_spectrum -sme = synthesize_spectrum(sme) -``` - -The synthesized spectra are stored in `sme.wave` and `sme.synth`: -```py -Iliffe_vector([array([6436. , 6436.1, 6436.2, ..., 6439.8, 6439.9]), - array([6442. , 6442.1, 6442.2, ... 6442.8, 6442.9])]) -Iliffe_vector([array([0.99986606, 0.99984309, 0.99980888, ..., 0.99754268, 0.99807395]), - array([0.99983848, 0.99985176, 0.99986281, ..., 0.99979064, 0.99471954])]) -``` - -## Fit an observed spectrum - -Assuming that we have an observed spectrum, with its wavelength, normalized flux, and uncertainties, as an array of `wave`, `flux` and `uncertainties`. - -They can be inserted into the SME structure with: -```py -sme.wave = wave -sme.spec = flux -sme.uncs = uncertainties -sme.mask = np.ones(len(Spectrum), dtype=int) -``` - -- The wavelength is always given in Angstrom. -- Note that the observation may be split into segments (orders etc). - - Then Wavelength is a list of arrays [segment1, segment2, ...], and the same applies to spec, uncs, and mask. -- The mask values are: 0: bad pixel, 1: line pixel, 2: continuum pixel, 4: vrad pixel - - The masks are additive, i.e., you can set mask value to 5 for line and vrad pixel. - - Note that the dtype of sme.mask must be int. - -Then the `solve` function can be used to find the best fit solution: -```py -from sme.solve import solve -# for more details on the fitparameter option, see fitparameters -fitparameters = ["teff", "logg", "monh", "abund Mg"] -sme = solve(sme, fitparameters) -``` - -## NLTE correction - -The non-Local Thermal Equilibrium (NLTE) correction can be applied to the spectrum by setting the `nlte` attribute of the SME structure: -```py -# SME also comes with a few NLTE grids, see NLTE section -# The current NLTE grid is only applicable to the MARCS model! -sme.nlte.set_nlte("Ca") -``` -```{note} -Upon the first use of the NLTE correction, the NLTE grid will be downloaded from the server and this may takes a while, depending on the size of the NLTE grid. -``` - -The results in the output sme structure can, for example, be plotted using the gui module (under development). -```py -from gui import plot_plotly -fig = plot_plotly.FinalPlot(sme) -fig.save(filename="sme.html") -``` - -```{raw} html - -``` - -## Save and Load - -The SME structure can be loaded with -```py -sme = SME_Structure.load("in.sme") -``` -or saved with -```py -sme.save("out.sme") -``` diff --git a/docs/usage/quickstart_old.md b/docs/usage/quickstart_old.md deleted file mode 100644 index 5f42de0c..00000000 --- a/docs/usage/quickstart_old.md +++ /dev/null @@ -1,95 +0,0 @@ -# Quickstart - -## SME Structure - -The first step in each SME project is to create an SME structure - -```py -from pysme.sme import SME_Structure -``` - -This can be done in done in a few different ways: -- load an existing SME save file (from Python or IDL): `sme = SME_Structure.load("sme.inp")` -- load an .ech file spectrum: `sme = SME_Structure.load("obs.ech")` -- assign values manually: `sme = SME_Structure()` - -Either way one has to make sure that a few essential properties are set in the object, those are: - * Stellar parameters (teff, logg, monh, abund) - >>> from pysme.abund import Abund - >>> sme.teff, sme.logg, sme.monh = 5700, 4.4, -0.1 - >>> sme.abund = Abund.solar() - * LineList (linelist), e.g. from VALD - >>> from pysme.linelist.vald import ValdFile - >>> vald = ValdFile("linelist.lin") - >>> sme.linelist = vald - * Atmosphere (atmo), the file has to be in PySME/src/sme/atmospheres - >>> # SME comes with a few model atmospheres see Atmosphere section - >>> sme.atmo.source = "marcs2012p_t1.0.sav" - >>> sme.atmo.method = "grid" - >>> sme.atmo.geom = "PP" - -If no wavelength grid (sme.wave) is set, one has to set the wavelength range: - * Wavelength range(s) in Ångstrom - >>> sme.wran = [[4500, 4600], [5200, 5400]] - -Furthermore for fitting to an observation an observation is required: - * Wavelength wave - >>> # The observation may be split into segments (orders etc) - >>> # Then Wavelength is a list of arrays [segment1, segment2, ...] - >>> # The same applies to spec, uncs, and mask - >>> # The wavelength is always given in Angstrom - >>> sme.wave = Wavelength - * Spectrum spec - >>> sme.spec = Spectrum - * Uncertainties uncs - >>> sme.uncs = Uncertainties - * Mask mask - >>> # The mask values are: 0: bad pixel, 1: line pixel, 2: continuum pixel, 4: vrad pixel - >>> # The masks are additive, i.e., you can set mask value to 5 for line and vrad pixel. - >>> # Note that the dtype of sme.mask must be int. - >>> sme.mask = np.ones(len(Spectrum), dtype=int) - * radial velocity and continuum flags - >>> # possible values are: "each", "whole", "fix", "none" - >>> # "each": Each segment is fitted individually - >>> # "whole": All segments are fit with the same radial velocity - >>> # "fix": use the set radial velocity - >>> # "none": Apply no radial velocity offset - >>> sme.vrad_flag = "whole" - >>> # possible values are: "constant", "linear", "fix", "none" - >>> # "constant": Scale the synthetic spectrum by a constant - >>> # "linear": Scale the synthetic spectrum by a linear polynomial - >>> # "fix": Use the set continnuum scale - >>> # "none": apply no continuum correction - >>> sme.cscale_flag = "linear" - >>> # possible values are: "whole", "mask" - >>> # "whole": use MCMC to match the synthetic to the observed spectrum - >>> # "mask": use the continuum mask points to fit the continuum - >>> sme.cscale_type = "mask" - -Optionally the following can be set: - * NLTE nlte for non local thermal equilibrium calculations - >>> # SME also comes with a few NLTE grids, see NLTE section - >>> # The NLTE grid is atmosphere model dependant! - >>> sme.nlte.set_nlte("Ca", "marcs2012p_t1.0_Ca.grd") - -Once the SME structure is prepared, SME can be run in one of its two modes: - 1. Synthesize a spectrum - >>> from sme.solve import synthesize_spectrum - >>> sme = synthesize_spectrum(sme) - 2. Finding the best fit (least squares) solution - >>> from sme.solve import solve - >>> # for more details on the fitparameter option, see fitparameters - >>> fitparameters = ["teff", "logg", "monh", "abund Mg"] - >>> sme = solve(sme, fitparameters) - -The results will be contained in the output sme structure. These can for example be plotted using the gui module. - >>> from gui import plot_plotly - >>> fig = plot_plotly.FinalPlot(sme) - >>> fig.save(filename="sme.html") - -.. raw:: html - :file: ../_static/sun.html - -or saved with - >>> sme.save("out.npy") - diff --git a/src/pysme/abund.py b/src/pysme/abund.py index fa297234..c6f937d7 100644 --- a/src/pysme/abund.py +++ b/src/pysme/abund.py @@ -678,8 +678,10 @@ def set_pattern_by_name(self, pattern_name): Parameters ---------- pattern_name : str - Name of the predefined option to use. One of 'asplund2009', 'grevesse2007', - 'lodders2003', 'empty' + Name of the predefined option to use. One of: + 'grevesse2007', 'asplund2009', 'asplund2021', 'asplund2005', + 'lodders2003', 'anders1989', 'grevesse1996', 'grevesse1998', ,'lodders2010', + 'empty' Raises ------ @@ -719,8 +721,10 @@ def set_pattern_by_name(self, pattern_name): self._pattern = self.empty_pattern() else: raise ValueError( - f"Got abundance pattern name {pattern_name} should be one of" - "'asplund2009', 'asplund2021', 'grevesse2007', 'lodders2003', or 'empty'." + f"Got abundance pattern name {pattern_name}. It should be one of: " + "'grevesse2007' (or 'solar'), 'asplund2009', 'asplund2021', " + "'asplund2005', 'lodders2003', 'lodders2010', 'anders1989', " + "'grevesse1996', 'grevesse1998', or 'empty'." ) def set_pattern_by_value(self, pattern, type): @@ -867,6 +871,6 @@ def _load_v1(file, names, folder=""): @staticmethod def solar(): - """Return solar abundances of asplund 2009""" + """Return solar abundances of Grevesse 2007""" solar = Abund(pattern="solar", monh=0) return solar diff --git a/src/pysme/continuum_and_radial_velocity.py b/src/pysme/continuum_and_radial_velocity.py index 958af6df..eebfe68d 100644 --- a/src/pysme/continuum_and_radial_velocity.py +++ b/src/pysme/continuum_and_radial_velocity.py @@ -322,7 +322,7 @@ def __call__( if sme.cscale is not None: cscale = sme.cscale[segments] else: - cscale = np.zeros(nseg, ndeg + 1) + cscale = np.zeros((nseg, ndeg + 1)) for i, seg in enumerate(segments): cscale[i, -1] = np.nanpercentile(y_obs[seg], 95) @@ -1199,18 +1199,21 @@ def match_rv_continuum(sme, segments, x_syn, y_syn): else: if sme.cscale_type == "mcmc": for s in segments: - ( - vrad[s], - vrad_unc[s], - cscale[s], - cscale_unc[s], - ) = continuum_normalization( + rv_i, rv_unc_i, cscale_i, cscale_unc_i = continuum_normalization( sme, x_syn[s], y_syn[s], s, rvel=vrad[s], ) + # MCMC returns arrays even for one segment; normalize them to row/scalar. + cscale_row = np.asarray(cscale_i, dtype=float).reshape(1, -1)[0] + vrad[s] = np.asarray(rv_i, dtype=float).reshape(-1)[0] + vrad_unc[s] = np.asarray(rv_unc_i, dtype=float).reshape(1, 2)[0] + cscale[s] = cscale_row + cscale_unc[s] = np.asarray(cscale_unc_i, dtype=float).reshape( + 1, cscale_row.size, 2 + )[0] else: for s in segments: cscale[s] = continuum_normalization( diff --git a/src/pysme/linelist/linelist.py b/src/pysme/linelist/linelist.py index 8bad7f4c..7e3e5d79 100644 --- a/src/pysme/linelist/linelist.py +++ b/src/pysme/linelist/linelist.py @@ -220,7 +220,14 @@ def __init__(self, linedata=None, lineformat="short", medium=None, **kwargs): #:pandas.DataFrame: DataFrame that contains all the data self._lines = linedata # should have all the fields (20) self.cdr_paras = None - self.cdr_paras_thres = {'teff':250, 'logg':0.5, 'monh':0.5, 'vmic':1} + self.cdr_paras_thres = { + 'teff': 250, + 'logg': 0.5, + 'monh': 0.5, + 'vmic': 1, + 'strong_depth': 0.001, + 'strong_bin_width': 0.2, + } if medium in ["air", "vac", None]: self._medium = medium else: diff --git a/src/pysme/linelist/vald.py b/src/pysme/linelist/vald.py index 4b9776ee..603a250c 100644 --- a/src/pysme/linelist/vald.py +++ b/src/pysme/linelist/vald.py @@ -142,6 +142,28 @@ def load(filename): """ return ValdFile(filename) + + def __getitem__(self, index): + if isinstance(index, str) and hasattr(self, index): + return getattr(self, index) + if isinstance(index, (list, str)): + if len(index) == 0: + return_list = deepcopy(self) + return_list._lines = self._lines.iloc[[]] + return_list.nlines = len(return_list._lines) + return return_list + values = self._lines[index].values + if index in self.string_columns: + values = values.astype(str) + return values + + if isinstance(index, int): + index = slice(index, index + 1) + return_list = deepcopy(self) + return_list._lines = self._lines.iloc[index] + return_list.nlines = len(return_list._lines) + return return_list + def identify_valdtype(self, lines): """Determines whether the file was created with extract_all, extract_stellar, or extract_element and whether it is in long or short format diff --git a/src/pysme/sme.py b/src/pysme/sme.py index 79e502a5..61af7118 100644 --- a/src/pysme/sme.py +++ b/src/pysme/sme.py @@ -252,6 +252,8 @@ class SME_Structure(Parameters): """), ("wran", None, this, this, "array of size (nseg, 2): beginning and end wavelength points of each segment"), + ("wint", None, vector, this, + "Iliffe_vector of shape (nseg, ...): optional wavelength grid passed to SMElib Transf"), ("wave", None, vector, this, "Iliffe_vector of shape (nseg, ...): wavelength"), ("spec", None, vector, this, @@ -299,7 +301,8 @@ def __init__(self, **kwargs): self.cdr_parallel = True self.cdr_n_jobs = 10 self.cdr_pysme_out = False - self.cdr_depth_thres = 0.0 + self.strong_depth_thres = 0.001 + self.strong_bin_width = 0.2 self.tdnlte_H = False # self.tdnlte_H_new = False diff --git a/src/pysme/sme_synth.py b/src/pysme/sme_synth.py index 83afa59e..2e572e2f 100644 --- a/src/pysme/sme_synth.py +++ b/src/pysme/sme_synth.py @@ -148,6 +148,21 @@ def ClearH2broad(self): """Clear flag for H2 molecule""" self.SetH2broad(False) + def SetLineInfoMode(self, mode): + """Set handling mode for precomputed line info (0=internal, 1=use_if_valid, 2=trust).""" + _smelib.SetLineInfoMode(int(mode)) + + def InputLinePrecomputedInfo(self, line_range_s, line_range_e, strong_mask, central_depth=None): + """Input precomputed line ranges and strong mask to SMElib.""" + range_s = np.ascontiguousarray(line_range_s, dtype=np.float64) + range_e = np.ascontiguousarray(line_range_e, dtype=np.float64) + strong = np.ascontiguousarray(strong_mask, dtype=np.uint8) + if central_depth is None: + _smelib.InputLinePrecomputedInfo(range_s, range_e, strong) + else: + depth = np.ascontiguousarray(central_depth, dtype=np.float64) + _smelib.InputLinePrecomputedInfo(range_s, range_e, strong, depth) + def InputLineList(self, linelist): """ Read in line list @@ -416,7 +431,7 @@ def GetOpacity(self, switch, species=None, key=None): if key is not None: kwargs["key"] = key if species is not None: - kwargs["species "] = species + kwargs["species"] = species return _smelib.GetOpacity(switch, **kwargs) def Ionization(self, ion=0): @@ -473,6 +488,21 @@ def GetNelec(self): """ return _smelib.GetNelec() + def GetFraction(self, species, mode=0): + """ + Get species fraction/partition function/number density. + + Parameters + ---------- + species : str + Species name in SPLIST (e.g., 'Fe', 'Fe+', 'CO') + mode : int, optional + 0: number density (FRACT * PF) + 1: partition function + other: FRACT + """ + return _smelib.GetFraction(species, mode=mode) + def Transf( self, mu, @@ -652,4 +682,4 @@ def GetContributionfunction(self, mu, wave, nwmax=400000, long_continuum=True): if long_continuum: return table, ctable else: - return table \ No newline at end of file + return table diff --git a/src/pysme/smelib/_smelib.cpp b/src/pysme/smelib/_smelib.cpp index 26302bee..c34a3134 100644 --- a/src/pysme/smelib/_smelib.cpp +++ b/src/pysme/smelib/_smelib.cpp @@ -163,6 +163,105 @@ static PyObject *smelib_ClearH2broad(PyObject *self, PyObject *args) Py_RETURN_NONE; } +static char smelib_SetLineInfoMode_docstring[] = "Set handling mode for precomputed line info"; +static PyObject *smelib_SetLineInfoMode(PyObject *self, PyObject *args) +{ + const int n = 1; + const char *result = NULL; + void *args_c[n]; + int mode = 0; + + if (!PyArg_ParseTuple(args, "i", &mode)) + return NULL; + + args_c[0] = &mode; + result = SetLineInfoMode(n, args_c); + if (result != NULL && result[0] != OK_response) + { + PyErr_SetString(PyExc_RuntimeError, result); + return NULL; + } + + Py_RETURN_NONE; +} + +static char smelib_InputLinePrecomputedInfo_docstring[] = "Input precomputed line ranges and strong mask"; +static PyObject *smelib_InputLinePrecomputedInfo(PyObject *self, PyObject *args) +{ + const char *result = NULL; + PyObject *range_s_obj = NULL, *range_e_obj = NULL, *strong_obj = NULL, *depth_obj = NULL; + PyArrayObject *range_s_arr = NULL, *range_e_arr = NULL, *strong_arr = NULL, *depth_arr = NULL; + void *args_c[5]; + int n = 4; + int nlines = 0; + + if (!PyArg_ParseTuple(args, "OOO|O", &range_s_obj, &range_e_obj, &strong_obj, &depth_obj)) + return NULL; + + range_s_arr = (PyArrayObject *)PyArray_FROM_OTF(range_s_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); + range_e_arr = (PyArrayObject *)PyArray_FROM_OTF(range_e_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); + strong_arr = (PyArrayObject *)PyArray_FROM_OTF(strong_obj, NPY_UBYTE, NPY_ARRAY_IN_ARRAY); + if (range_s_arr == NULL || range_e_arr == NULL || strong_arr == NULL) + goto fail; + + if (PyArray_NDIM(range_s_arr) != 1 || PyArray_NDIM(range_e_arr) != 1 || PyArray_NDIM(strong_arr) != 1) + { + PyErr_SetString(PyExc_ValueError, "Expected 1D arrays for range_s, range_e, and strong_mask"); + goto fail; + } + + nlines = (int)PyArray_DIM(range_s_arr, 0); + if ((int)PyArray_DIM(range_e_arr, 0) != nlines || (int)PyArray_DIM(strong_arr, 0) != nlines) + { + PyErr_SetString(PyExc_ValueError, "Expected range_s, range_e, and strong_mask to have same length"); + goto fail; + } + + args_c[0] = &nlines; + args_c[1] = PyArray_DATA(range_s_arr); + args_c[2] = PyArray_DATA(range_e_arr); + args_c[3] = PyArray_DATA(strong_arr); + + if (depth_obj != NULL && depth_obj != Py_None) + { + depth_arr = (PyArrayObject *)PyArray_FROM_OTF(depth_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); + if (depth_arr == NULL) + goto fail; + if (PyArray_NDIM(depth_arr) != 1) + { + PyErr_SetString(PyExc_ValueError, "Expected 1D array for central_depth"); + goto fail; + } + if ((int)PyArray_DIM(depth_arr, 0) != nlines) + { + PyErr_SetString(PyExc_ValueError, "Expected central_depth to have same length as line ranges"); + goto fail; + } + args_c[4] = PyArray_DATA(depth_arr); + n = 5; + } + + result = InputLinePrecomputedInfo(n, args_c); + if (result != NULL && result[0] != OK_response) + { + PyErr_SetString(PyExc_RuntimeError, result); + goto fail; + } + + Py_XDECREF(range_s_arr); + Py_XDECREF(range_e_arr); + Py_XDECREF(strong_arr); + Py_XDECREF(depth_arr); + Py_RETURN_NONE; + +fail: + Py_XDECREF(range_s_arr); + Py_XDECREF(range_e_arr); + Py_XDECREF(strong_arr); + Py_XDECREF(depth_arr); + return NULL; +} + static char smelib_InputLineList_docstring[] = "Read in line list"; static PyObject *smelib_InputLineList(PyObject *self, PyObject *args) { @@ -791,8 +890,11 @@ static PyObject *smelib_GetOpacity(PyObject *self, PyObject *args, PyObject *kwd PyArrayObject *arr; static const char *keywords[] = {"flag", "species", "key", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "s|ss", const_cast(keywords), &choice)) + if (!PyArg_ParseTupleAndKeywords(args, kwds, "s|ss", const_cast(keywords), + &choice, &species, &key)) + { return NULL; + } if (strcmp(choice, "COPSTD") == 0) number = -3; @@ -890,6 +992,48 @@ static PyObject *smelib_GetOpacity(PyObject *self, PyObject *args, PyObject *kwd return (PyObject *)arr; } +static char smelib_GetFraction_docstring[] = "Returns fraction/partition function for a species"; +static PyObject *smelib_GetFraction(PyObject *self, PyObject *args, PyObject *kwds) +{ + int n = 4; + void *args_c[4]; + const char *result = NULL; + char *species = NULL; + short mode = 2; + short length; + IDL_STRING species_idl; + + PyArrayObject *arr; + + static const char *keywords[] = {"species", "mode", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "s|h", const_cast(keywords), + &species, &mode)) + return NULL; + + species_idl.slen = strlen(species); + species_idl.s = species; + species_idl.stype = 0; + + length = GetNRHOX(); + npy_intp dims[] = {length}; + arr = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE); + + args_c[0] = &species_idl; + args_c[1] = &mode; + args_c[2] = &length; + args_c[3] = PyArray_DATA(arr); + + result = GetFraction(n, args_c); + + if (result != NULL && result[0] != OK_response) + { + Py_DECREF(arr); + PyErr_SetString(PyExc_RuntimeError, result); + return NULL; + } + return (PyObject *)arr; +} + static char smelib_Ionization_docstring[] = "Perform EOS calculations"; static PyObject *smelib_Ionization(PyObject *self, PyObject *args) { @@ -1368,6 +1512,8 @@ static PyMethodDef module_methods[] = { {"SetVWscale", smelib_SetVWscale, METH_VARARGS, smelib_SetVWscale_docstring}, {"SetH2broad", smelib_SetH2broad, METH_NOARGS, smelib_SetH2broad_docstring}, {"ClearH2broad", smelib_ClearH2broad, METH_NOARGS, smelib_ClearH2broad_docstring}, + {"SetLineInfoMode", smelib_SetLineInfoMode, METH_VARARGS, smelib_SetLineInfoMode_docstring}, + {"InputLinePrecomputedInfo", smelib_InputLinePrecomputedInfo, METH_VARARGS, smelib_InputLinePrecomputedInfo_docstring}, {"InputLineList", smelib_InputLineList, METH_VARARGS, smelib_InputLineList_docstring}, {"OutputLineList", smelib_OutputLineList, METH_NOARGS, smelib_OutputLineList_docstring}, {"UpdateLineList", smelib_UpdateLineList, METH_VARARGS, smelib_UpdateLineList_docstring}, @@ -1378,6 +1524,7 @@ static PyMethodDef module_methods[] = { {"InputAbund", smelib_InputAbund, METH_VARARGS, smelib_InputAbund_docstring}, {"Opacity", smelib_Opacity, METH_NOARGS, smelib_Opacity_docstring}, {"GetOpacity", (PyCFunction)(void (*)(void))smelib_GetOpacity, METH_VARARGS | METH_KEYWORDS, smelib_GetOpacity_docstring}, + {"GetFraction", (PyCFunction)(void (*)(void))smelib_GetFraction, METH_VARARGS | METH_KEYWORDS, smelib_GetFraction_docstring}, {"Ionization", smelib_Ionization, METH_VARARGS, smelib_Ionization_docstring}, {"GetDensity", smelib_GetDensity, METH_NOARGS, smelib_GetDensity_docstring}, {"GetNatom", smelib_GetNatom, METH_NOARGS, smelib_GetNatom_docstring}, diff --git a/src/pysme/smelib/libtools.py b/src/pysme/smelib/libtools.py index 90696895..9b5b7dfc 100644 --- a/src/pysme/smelib/libtools.py +++ b/src/pysme/smelib/libtools.py @@ -25,6 +25,7 @@ smelib_releases = { "default": "latest/download", + "0.6.23": "download/v6.13.12", "0.4.198": "download/v6.0.6", "0.4.199": "download/v6.0.6", } diff --git a/src/pysme/smelib/sme_synth_faster.h b/src/pysme/smelib/sme_synth_faster.h index 35e421f6..500f2ebd 100644 --- a/src/pysme/smelib/sme_synth_faster.h +++ b/src/pysme/smelib/sme_synth_faster.h @@ -47,6 +47,8 @@ extern "C" const char *SME_DLL InputWaveRange(int n, void *arg[]); /* Read in extern "C" const char *SME_DLL SetVWscale(int n, void *arg[]); /* Set van der Waals scaling factor */ extern "C" const char *SME_DLL SetH2broad(int n, void *arg[]); /* Set flag for H2 molecule */ extern "C" const char *SME_DLL ClearH2broad(int n, void *arg[]); /* Clear flag for H2 molecule */ +extern "C" const char *SME_DLL SetLineInfoMode(int n, void *arg[]); /* Set handling mode for precomputed line info */ +extern "C" const char *SME_DLL InputLinePrecomputedInfo(int n, void *arg[]); /* Input precomputed line ranges/strong mask */ extern "C" const char *SME_DLL InputLineList(int n, void *arg[]); /* Read in line list */ extern "C" const char *SME_DLL OutputLineList(int n, void *arg[]); /* Return line list */ extern "C" const char *SME_DLL UpdateLineList(int n, void *arg[]); /* Change line list parameters */ @@ -58,6 +60,7 @@ extern "C" const char *SME_DLL ResetDepartureCoefficients(int n, void *arg[]); / extern "C" const char *SME_DLL InputAbund(int n, void *arg[]); /* Read in abundances */ extern "C" const char *SME_DLL Opacity(int n, void *arg[]); /* Calculate opacities */ extern "C" const char *SME_DLL GetOpacity(int n, void *arg[]); /* Returns specific cont. opacity */ +extern "C" const char *SME_DLL GetFraction(int n, void *arg[]); /* Returns fraction or partition function */ extern "C" const char *SME_DLL Ionization(int n, void *arg[]); /* Perfrom EOS calculations */ extern "C" const char *SME_DLL GetDensity(int n, void *arg[]); /* Returns density in g/cm^3 */ extern "C" const char *SME_DLL GetNatom(int n, void *arg[]); /* Returns atomic number density */ @@ -66,4 +69,4 @@ extern "C" const char *SME_DLL Transf(int n, void *arg[]); / extern "C" const char *SME_DLL CentralDepth(int n, void *arg[]); /* Computes line central depths */ extern "C" const char *SME_DLL GetLineOpacity(int n, void *arg[]); /* Returns specific line opacity */ extern "C" const char *SME_DLL GetLineRange(int n, void *arg[]); /* Get validity range for every line */ -extern "C" const char *SME_DLL Contribution_functions(int n, void *arg[]); /* Get contribution functions */ \ No newline at end of file +extern "C" const char *SME_DLL Contribution_functions(int n, void *arg[]); /* Get contribution functions */ diff --git a/src/pysme/solve.py b/src/pysme/solve.py index fa791a7a..327181bb 100644 --- a/src/pysme/solve.py +++ b/src/pysme/solve.py @@ -92,6 +92,7 @@ def __init__(self, filename=None, restore=False): self.iteration = 0 self.parameter_names = [] self.update_linelist = False + self.derived_param = None self._latest_residual = None self._latest_jacobian = None self.restore = restore @@ -135,7 +136,20 @@ def backup(self, sme): pass def _residuals( - self, param, sme, spec, uncs, mask, segments="all", isJacobian=False, linelist_mode='all', **_ + self, + param, + sme, + spec, + uncs, + mask, + segments="all", + isJacobian=False, + linelist_mode="all", + smelib_lineinfo_mode=0, + cdr_database=None, + cdr_create=False, + vbroad_expend_ratio=2, + **_, ): """ Calculates the synthetic spectrum with sme_func and @@ -178,20 +192,20 @@ def _residuals( # change parameters for name, value in zip(self.parameter_names, param): - if self.dynamic_param is not None: - if name in self.dynamic_param.keys(): - raise ValueError(f'fitting parameter {name} cannot be also in dynamic parameter.') + if self.derived_param is not None: + if name in self.derived_param.keys(): + raise ValueError(f"fitting parameter {name} cannot be also in derived parameter.") sme[name] = value - # change dynamic parameters - if self.dynamic_param is not None: - for name in self.dynamic_param.keys(): - if 'abund' in name: + # change derived parameters + if self.derived_param is not None: + for name in self.derived_param.keys(): + if "abund" in name: abund_name = name.split()[1] - sme.abund[abund_name] = self.dynamic_param[name](sme) - sme.monh - print(f'Changing dynamic parameter {name} to {sme.abund[abund_name]:.2f}.') + sme.abund[abund_name] = self.derived_param[name](sme) - sme.monh + print(f"Changing derived parameter {name} to {sme.abund[abund_name]:.2f}.") else: - sme[name] = self.dynamic_param[name](sme) - print(f'Changing dynamic parameter {name} to {sme[name]:.2f}.') + sme[name] = self.derived_param[name](sme) + print(f"Changing derived parameter {name} to {sme[name]:.2f}.") # run spectral synthesis try: result = self.synthesizer.synthesize_spectrum( @@ -202,7 +216,11 @@ def _residuals( passLineList=True, updateLineList=self.update_linelist, radial_velocity_mode=radial_velocity_mode, - linelist_mode=linelist_mode + linelist_mode=linelist_mode, + smelib_lineinfo_mode=smelib_lineinfo_mode, + cdr_database=cdr_database, + cdr_create=cdr_create, + vbroad_expend_ratio=vbroad_expend_ratio, ) except AtmosphereError as ae: # Something went wrong (left the grid? Don't go there) @@ -257,6 +275,10 @@ def _jacobian( step_sizes=None, method="2-point", linelist_mode='all', + smelib_lineinfo_mode=0, + cdr_database=None, + cdr_create=False, + vbroad_expend_ratio=2, **_, ): """ @@ -283,7 +305,15 @@ def _jacobian( abs_step=step_sizes, bounds=bounds, args=args, - kwargs={"isJacobian": True, "segments": segments, "linelist_mode": linelist_mode}, + kwargs={ + "isJacobian": True, + "segments": segments, + "linelist_mode": linelist_mode, + "smelib_lineinfo_mode": smelib_lineinfo_mode, + "cdr_database": cdr_database, + "cdr_create": cdr_create, + "vbroad_expend_ratio": vbroad_expend_ratio, + }, ) if not np.all(np.isfinite(g)): @@ -613,7 +643,21 @@ def sanitize_parameter_names(self, sme, param_names): ) return param_names - def solve(self, sme, param_names=None, segments="all", bounds=None, step_sizes=None, dynamic_param=None, linelist_mode='all'): + def solve( + self, + sme, + param_names=None, + segments="all", + bounds=None, + step_sizes=None, + derived_param=None, + dynamic_param=None, + linelist_mode="all", + smelib_lineinfo_mode=0, + cdr_database=None, + cdr_create=False, + vbroad_expend_ratio=2, + ): """ Find the least squares fit parameters to an observed spectrum @@ -637,7 +681,15 @@ def solve(self, sme, param_names=None, segments="all", bounds=None, step_sizes=N assert "wave" in sme, "SME Structure has no wavelength" assert "spec" in sme, "SME Structure has no observation" - self.dynamic_param = dynamic_param + if derived_param is not None and dynamic_param is not None and derived_param is not dynamic_param: + raise ValueError("Specify only one of derived_param or dynamic_param, not both.") + if dynamic_param is not None: + warnings.warn( + "'dynamic_param' is deprecated and will be removed in a future release; use 'derived_param' instead.", + DeprecationWarning, + stacklevel=2, + ) + self.derived_param = derived_param if derived_param is not None else dynamic_param if self.restore and self.filename is not None: fname = self.filename.rsplit(".", 1)[0] @@ -767,7 +819,11 @@ def solve(self, sme, param_names=None, segments="all", bounds=None, step_sizes=N "segments": segments, "step_sizes": step_sizes, "method": sme.leastsquares_jac, - "linelist_mode": linelist_mode + "linelist_mode": linelist_mode, + "smelib_lineinfo_mode": smelib_lineinfo_mode, + "cdr_database": cdr_database, + "cdr_create": cdr_create, + "vbroad_expend_ratio": vbroad_expend_ratio, }, ) # The jacobian is altered by the loss function @@ -793,7 +849,15 @@ def solve(self, sme, param_names=None, segments="all", bounds=None, step_sizes=N # We could try to reuse the already calculated synthetic spectrum (if it already exists) # However it is much lower resolution then the newly synthethized one (usually) # Therefore the radial velocity wont be as good as when redoing the whole thing - sme = self.synthesizer.synthesize_spectrum(sme, segments, linelist_mode=linelist_mode) + sme = self.synthesizer.synthesize_spectrum( + sme, + segments, + linelist_mode=linelist_mode, + smelib_lineinfo_mode=smelib_lineinfo_mode, + cdr_database=cdr_database, + cdr_create=cdr_create, + vbroad_expend_ratio=vbroad_expend_ratio, + ) else: raise ValueError("No fit parameters given") @@ -804,7 +868,21 @@ def solve(self, sme, param_names=None, segments="all", bounds=None, step_sizes=N def solve( - sme, param_names=None, segments="all", filename=None, restore=False, **kwargs + sme, + param_names=None, + segments="all", + filename=None, + restore=False, + derived_param=None, + dynamic_param=None, + **kwargs, ): solver = SME_Solver(filename=filename, restore=restore) - return solver.solve(sme, param_names, segments, **kwargs) + return solver.solve( + sme, + param_names, + segments, + derived_param=derived_param, + dynamic_param=dynamic_param, + **kwargs, + ) diff --git a/src/pysme/synthesize.py b/src/pysme/synthesize.py index 23a4c6bb..fee4747b 100644 --- a/src/pysme/synthesize.py +++ b/src/pysme/synthesize.py @@ -4,6 +4,7 @@ """ import logging import uuid +import warnings import numpy as np from scipy.constants import speed_of_light @@ -606,7 +607,8 @@ def synthesize_spectrum( cdr_create=False, keep_line_opacity=False, vbroad_expend_ratio=2, - contribution_function=False + contribution_function=False, + smelib_lineinfo_mode=0, ): """ Calculate the synthetic spectrum based on the parameters passed in the SME structure @@ -691,28 +693,126 @@ def synthesize_spectrum( wave = [w for w in sme.wave] dll = self.get_dll(dll_id) + dll.SetLineInfoMode(int(smelib_lineinfo_mode)) + + # "dynamic" is the preferred name in publications; keep "auto" as compatibility alias. + if linelist_mode == "auto": + warnings.warn( + "'linelist_mode=\"auto\"' is deprecated; use 'linelist_mode=\"dynamic\"' instead.", + DeprecationWarning, + stacklevel=2, + ) + linelist_mode = "dynamic" + elif linelist_mode not in ("all", "dynamic"): + raise ValueError("linelist_mode must be one of: 'all', 'dynamic'") - # Calculate the line central depth and line range if necessary - if linelist_mode == 'auto': + # Calculate line properties and strong-line flags if necessary. + if linelist_mode == 'dynamic': # logger.info(f'linelist mode: {linelist_mode}') - if sme.linelist.cdr_paras is None or not {'central_depth', 'line_range_s', 'line_range_e'}.issubset(sme.linelist._lines.columns) or np.abs(sme.linelist.cdr_paras[0]-sme.teff) >= sme.linelist.cdr_paras_thres['teff'] or (np.abs(sme.linelist.cdr_paras[1]-sme.logg) >= sme.linelist.cdr_paras_thres['logg']) or (np.abs(sme.linelist.cdr_paras[2]-sme.monh) >= sme.linelist.cdr_paras_thres['monh']) or cdr_create: - logger.info(f'Updating linelist central depth and line range.') - sme = self.update_cdr(sme, cdr_database=cdr_database, cdr_create=cdr_create, show_progress_bars=show_progress_bars) + need_update_cdr = ( + sme.linelist.cdr_paras is None + or not {'central_depth', 'line_range_s', 'line_range_e'}.issubset(sme.linelist._lines.columns) + or np.abs(sme.linelist.cdr_paras[0] - sme.teff) >= sme.linelist.cdr_paras_thres['teff'] + or np.abs(sme.linelist.cdr_paras[1] - sme.logg) >= sme.linelist.cdr_paras_thres['logg'] + or np.abs(sme.linelist.cdr_paras[2] - sme.monh) >= sme.linelist.cdr_paras_thres['monh'] + or cdr_create + ) + if need_update_cdr: + logger.info('Updating linelist central depth and line range.') + sme = self.update_cdr( + sme, + cdr_database=cdr_database, + cdr_create=cdr_create, + show_progress_bars=show_progress_bars, + ) + + strong_depth_prev = sme.linelist.cdr_paras_thres.get('strong_depth') + strong_bin_width_prev = sme.linelist.cdr_paras_thres.get('strong_bin_width') + strong_depth_matches = ( + strong_depth_prev is not None + and np.isclose(float(strong_depth_prev), float(sme.strong_depth_thres)) + ) + strong_bin_width_matches = ( + strong_bin_width_prev is not None + and np.isclose(float(strong_bin_width_prev), float(sme.strong_bin_width)) + ) + need_update_strong = ( + need_update_cdr + or 'strong' not in sme.linelist._lines.columns + or not strong_depth_matches + or not strong_bin_width_matches + ) + + if need_update_strong: + strong_mask = self.flag_strong_lines_by_bins( + sme.linelist['wlcent'], + sme.linelist['central_depth'], + bin_width=sme.strong_bin_width, + threshold=sme.strong_depth_thres, + ) + sme.linelist._lines['strong'] = np.asarray(strong_mask, dtype=bool) + sme.linelist.cdr_paras_thres['strong_depth'] = float(sme.strong_depth_thres) + sme.linelist.cdr_paras_thres['strong_bin_width'] = float(sme.strong_bin_width) # Input Model data to C library dll.SetLibraryPath() if passLineList: - if linelist_mode == 'auto': + linelist_for_smelib = sme.linelist + if linelist_mode == 'dynamic': line_indices = sme.linelist['wlcent'] < 0 v_broad = np.sqrt(sme.vmac**2 + sme.vsini**2 + (clight/sme.ipres)**2) for i in range(sme.nseg): line_indices |= (sme.linelist['line_range_e'] > sme.wran[i][0] * (1 - vbroad_expend_ratio*v_broad/clight)) & (sme.linelist['line_range_s'] < sme.wran[i][1] * (1 + vbroad_expend_ratio*v_broad/clight)) - line_indices &= self.flag_strong_lines_by_bins(sme.linelist['wlcent'], sme.linelist['central_depth']) + line_indices &= np.asarray(sme.linelist['strong'], dtype=bool) sme.linelist._lines['use_indices'] = line_indices - line_ion_mask = dll.InputLineList(sme.linelist[line_indices]) - else: - line_ion_mask = dll.InputLineList(sme.linelist) + linelist_for_smelib = sme.linelist[line_indices] + line_ion_mask = dll.InputLineList(linelist_for_smelib) sme.line_ion_mask = line_ion_mask + + if smelib_lineinfo_mode in (1, 2): + cols = set(linelist_for_smelib._lines.columns) + required_cols = {"line_range_s", "line_range_e"} + missing_cols = sorted(required_cols - cols) + + # Build strong mask from existing column or central depth. + if "strong" in cols: + strong_all = np.asarray(linelist_for_smelib["strong"], dtype=np.uint8) + elif "central_depth" in cols: + strong_all = np.asarray( + self.flag_strong_lines_by_bins( + linelist_for_smelib["wlcent"], + linelist_for_smelib["central_depth"], + bin_width=sme.strong_bin_width, + threshold=sme.strong_depth_thres, + ), + dtype=np.uint8, + ) + else: + strong_all = None + + if strong_all is None: + missing_cols.append("strong/central_depth") + + if missing_cols: + msg = ( + f"SMElib lineinfo mode={smelib_lineinfo_mode} requested but missing " + f"precomputed columns: {', '.join(missing_cols)}" + ) + if smelib_lineinfo_mode == 2: + raise ValueError(msg) + logger.warning("%s. Falling back to mode 0.", msg) + dll.SetLineInfoMode(0) + else: + keep_mask = ~np.asarray(line_ion_mask, dtype=bool) + range_s = np.asarray(linelist_for_smelib["line_range_s"], dtype=np.float64)[keep_mask] + range_e = np.asarray(linelist_for_smelib["line_range_e"], dtype=np.float64)[keep_mask] + strong_mask = np.asarray(strong_all, dtype=np.uint8)[keep_mask] + + if "central_depth" in cols: + depth = np.asarray(linelist_for_smelib["central_depth"], dtype=np.float64)[keep_mask] + dll.InputLinePrecomputedInfo(range_s, range_e, strong_mask, depth) + else: + dll.InputLinePrecomputedInfo(range_s, range_e, strong_mask) if hasattr(updateLineList, "__len__") and len(updateLineList) > 0: # TODO Currently Updates the whole linelist, could be improved to only change affected lines dll.UpdateLineList(sme.atomic, sme.species, updateLineList) @@ -831,7 +931,7 @@ def synthesize_spectrum( sme.vrad = np.asarray(vrad) sme.vrad_unc = np.asarray(vrad_unc) nlte_flags = dll.GetNLTEflags() - if linelist_mode == 'auto': + if linelist_mode == 'dynamic': sme.linelist._lines.loc[sme.linelist._lines['use_indices'], 'nlte_flag'] = nlte_flags.astype(int) else: sme.nlte.flags = nlte_flags @@ -907,7 +1007,7 @@ def synthesize_segment( # ipres_segment = sme.ipres if np.size(sme.ipres) == 1 else sme.ipres[segment] # if ipres_segment != 0: # del_wav += sme.linelist['wlcent'] / ipres_segment - # indices = (~((sme.linelist['line_range_e'] < wbeg - del_wav - line_margin) | (sme.linelist['line_range_s'] > wend + del_wav + line_margin))) & (sme.linelist['central_depth'] > sme.cdr_depth_thres) + # indices = (~((sme.linelist['line_range_e'] < wbeg - del_wav - line_margin) | (sme.linelist['line_range_s'] > wend + del_wav + line_margin))) & (sme.linelist['central_depth'] > sme.strong_depth_thres) # # logger.info(f"There are {len(sme.linelist[indices][np.char.find(sme.linelist[indices]['species'], 'Fe') >= 0])} Fe lines in sub linelist.") # _ = dll.InputLineList(sme.linelist[indices]) # sme.linelist._lines['use_indices'] = indices @@ -928,8 +1028,22 @@ def synthesize_segment( dll.InputWaveRange(wbeg-2, wend+2) dll.Opacity() - # Reuse adaptive wavelength grid in the jacobians - if reuse_wavelength_grid and segment in self.wint.keys(): + # Priority for wavelength grid passed to SMElib: + # 1) user-provided sme.wint for this segment + # 2) internal cache when reuse_wavelength_grid=True + # 3) None (let SMElib compute it) + user_wint_seg = None + if getattr(sme, "wint", None) is not None: + try: + user_wint_seg = sme.wint[segment] + if user_wint_seg is not None and len(user_wint_seg) == 0: + user_wint_seg = None + except (IndexError, KeyError, TypeError): + user_wint_seg = None + + if user_wint_seg is not None: + wint_seg = np.asarray(user_wint_seg, dtype=np.float64) + elif reuse_wavelength_grid and segment in self.wint: wint_seg = self.wint[segment] else: wint_seg = None @@ -957,7 +1071,7 @@ def synthesize_segment( # Store the adaptive wavelength grid for the future # if it was newly created - if wint_seg is None: + if user_wint_seg is None and wint_seg is None: self.wint[segment] = wint if not sme.specific_intensities_only: @@ -984,7 +1098,7 @@ def synthesize_segment( # instrument broadening if "iptype" in sme: logger.debug("Apply detector broadening") - ipres = sme.ipres if np.size(sme.ipres) == 1 else sme.ipres[segment] + ipres = sme.ipres.item() if np.size(sme.ipres) == 1 else sme.ipres[segment] sint = broadening.apply_broadening(ipres, wint, sint, type=sme.iptype, sme=sme) # Apply the correction on Ha, Hb and Hgamma line here. @@ -1026,81 +1140,81 @@ def update_cdr(self, sme, cdr_database=None, cdr_create=False, cdr_grid_overwrit N_line_chunk, parallel, n_jobs, pysme_out = sme.cdr_N_line_chunk, sme.cdr_parallel, sme.cdr_n_jobs, sme.cdr_pysme_out self.update_cdr_switch = True + try: + # Decide how many chunks to be divided + N_chunk = int(np.ceil(len(sme.linelist) / N_line_chunk)) - # Decide how many chunks to be divided - N_chunk = int(np.ceil(len(sme.linelist) / N_line_chunk)) + # Divide the line list to sub line lists + sub_linelist = [sme.linelist[N_line_chunk*i:N_line_chunk*(i+1)] for i in range(N_chunk)] - # Divide the line list to sub line lists - sub_linelist = [sme.linelist[N_line_chunk*i:N_line_chunk*(i+1)] for i in range(N_chunk)] - - if sum(len(item) for item in sub_linelist) != len(sme.linelist): - raise ValueError - - if cdr_database is not None: - # cdr_database is provided, use it to update the central depth and line range - self._interpolate_or_compute_and_update_linelist(sme, cdr_database, cdr_create=cdr_create, cdr_grid_overwrite=cdr_grid_overwrite, mode=mode, dims=dims, show_progress_bars=show_progress_bars) - return sme - - logger.info('[cdr] Using calculation to update central depth and line range.') - sub_sme_init = SME_Structure() - exclude_keys = ['_wave', '_synth', '_spec', '_uncs', '_mask', '_SME_Structure__wran', '_normalize_by_continuum', '_specific_intensities_only', '_telluric', '__cont', '_linelist', '_fitparameters', '_fitresults'] - for key, value in sme.__dict__.items(): - if key not in exclude_keys and 'cscale' not in key and 'vrad' not in key: - setattr(sub_sme_init, key, deepcopy(value)) - sub_sme_init.wave = np.arange(5000, 5010, 1) - - if not parallel: - for i in tqdm(range(N_chunk), disable=not show_progress_bars): - sub_sme_init.linelist = sub_linelist[i] - sub_sme_init = self.synthesize_spectrum(sub_sme_init) - if i == 0: - stack_linelist = deepcopy(sub_sme_init.linelist) - else: - stack_linelist._lines = pd.concat([stack_linelist._lines, sub_sme_init.linelist._lines]) - else: - sub_sme = [] - sub_sme_init.linelist = sme.linelist[:1] - sub_sme_init = self.synthesize_spectrum(sub_sme_init) - for i in range(N_chunk): - sub_sme.append(deepcopy(sub_sme_init)) - sub_sme[i].linelist = sub_linelist[i] - - if pysme_out: - sub_sme = pqdm(sub_sme, self.synthesize_spectrum, n_jobs=n_jobs, disable=not show_progress_bars) - else: - with redirect_stdout(open(f"/dev/null", 'w')): - sub_sme = pqdm(sub_sme, self.synthesize_spectrum, n_jobs=n_jobs, disable=not show_progress_bars) + if sum(len(item) for item in sub_linelist) != len(sme.linelist): + raise ValueError - for i in range(N_chunk): - sub_linelist[i] = sub_sme[i].linelist - stack_linelist = deepcopy(sub_linelist[0]) - stack_linelist._lines = pd.concat([ele._lines for ele in sub_linelist]) - # logger.info(f'{sub_linelist}') - - # Remove - if len(stack_linelist) != len(sme.linelist): - raise ValueError - for column in ['central_depth', 'line_range_s', 'line_range_e']: - if column in sme.linelist.columns: - sme.linelist._lines = sme.linelist._lines.drop(column, axis=1) - for column in ['central_depth', 'line_range_s', 'line_range_e']: - sme.linelist._lines[column] = stack_linelist._lines[column] - # pickle.dump([sme.linelist._lines, stack_linelist._lines], open('linelist.pkl', 'wb')) - - # Manually change the depth of all H 1 lines to 1, to include them back. - sme.linelist._lines.loc[sme.linelist['species'] == 'H 1', 'central_depth'] = 1 - - # Manually change the 2000 line_range to 0.03. - indices = np.isclose(sme.linelist['line_range_e'] - sme.linelist['line_range_s'], 2000, rtol=1e-4, atol=5, equal_nan=False) - sme.linelist._lines.loc[indices, 'line_range_s'] = sme.linelist._lines.loc[indices, 'wlcent']-0.3 - sme.linelist._lines.loc[indices, 'line_range_e'] = sme.linelist._lines.loc[indices, 'wlcent']+0.3 - - # Write the stellar parameters used here to the line list - sme.linelist.cdr_paras = np.array([sme.teff, sme.logg, sme.monh, sme.vmic]) + if cdr_database is not None: + # cdr_database is provided, use it to update the central depth and line range + self._interpolate_or_compute_and_update_linelist(sme, cdr_database, cdr_create=cdr_create, cdr_grid_overwrite=cdr_grid_overwrite, mode=mode, dims=dims, show_progress_bars=show_progress_bars) + return sme + + logger.info('[cdr] Using calculation to update central depth and line range.') + sub_sme_init = SME_Structure() + exclude_keys = ['_wave', '_synth', '_spec', '_uncs', '_mask', '_SME_Structure__wran', '_normalize_by_continuum', '_specific_intensities_only', '_telluric', '__cont', '_linelist', '_fitparameters', '_fitresults'] + for key, value in sme.__dict__.items(): + if key not in exclude_keys and 'cscale' not in key and 'vrad' not in key: + setattr(sub_sme_init, key, deepcopy(value)) + sub_sme_init.wave = np.arange(5000, 5010, 1) + + if not parallel: + for i in tqdm(range(N_chunk), disable=not show_progress_bars): + sub_sme_init.linelist = sub_linelist[i] + sub_sme_init = self.synthesize_spectrum(sub_sme_init) + if i == 0: + stack_linelist = deepcopy(sub_sme_init.linelist) + else: + stack_linelist._lines = pd.concat([stack_linelist._lines, sub_sme_init.linelist._lines]) + else: + sub_sme = [] + sub_sme_init.linelist = sme.linelist[:1] + sub_sme_init = self.synthesize_spectrum(sub_sme_init) + for i in range(N_chunk): + sub_sme.append(deepcopy(sub_sme_init)) + sub_sme[i].linelist = sub_linelist[i] - self.update_cdr_switch = False + if pysme_out: + sub_sme = pqdm(sub_sme, self.synthesize_spectrum, n_jobs=n_jobs, disable=not show_progress_bars) + else: + with redirect_stdout(open(f"/dev/null", 'w')): + sub_sme = pqdm(sub_sme, self.synthesize_spectrum, n_jobs=n_jobs, disable=not show_progress_bars) + + for i in range(N_chunk): + sub_linelist[i] = sub_sme[i].linelist + stack_linelist = deepcopy(sub_linelist[0]) + stack_linelist._lines = pd.concat([ele._lines for ele in sub_linelist]) + # logger.info(f'{sub_linelist}') + + # Remove + if len(stack_linelist) != len(sme.linelist): + raise ValueError + for column in ['central_depth', 'line_range_s', 'line_range_e']: + if column in sme.linelist.columns: + sme.linelist._lines = sme.linelist._lines.drop(column, axis=1) + for column in ['central_depth', 'line_range_s', 'line_range_e']: + sme.linelist._lines[column] = stack_linelist._lines[column] + # pickle.dump([sme.linelist._lines, stack_linelist._lines], open('linelist.pkl', 'wb')) + + # Manually change the depth of all H 1 lines to 1, to include them back. + sme.linelist._lines.loc[sme.linelist['species'] == 'H 1', 'central_depth'] = 1 + + # Manually change the 2000 line_range to 0.03. + indices = np.isclose(sme.linelist['line_range_e'] - sme.linelist['line_range_s'], 2000, rtol=1e-4, atol=5, equal_nan=False) + sme.linelist._lines.loc[indices, 'line_range_s'] = sme.linelist._lines.loc[indices, 'wlcent']-0.3 + sme.linelist._lines.loc[indices, 'line_range_e'] = sme.linelist._lines.loc[indices, 'wlcent']+0.3 + + # Write the stellar parameters used here to the line list + sme.linelist.cdr_paras = np.array([sme.teff, sme.logg, sme.monh, sme.vmic]) - return sme + return sme + finally: + self.update_cdr_switch = False def _interpolate_or_compute_and_update_linelist( self, sme, cdr_database, cdepth_decimals=4, cdepth_thres=0, diff --git a/test/test_atmospheres.py b/test/test_atmospheres.py index f4fed997..241907e2 100644 --- a/test/test_atmospheres.py +++ b/test/test_atmospheres.py @@ -23,7 +23,6 @@ def atmosphere(atmosphere_name): @pytest.fixture -@pytest.mark.usefixtures("lfs_atmo") def interpolator(atmosphere, lfs_atmo): interp = AtmosphereInterpolator( depth=atmosphere.depth, @@ -35,7 +34,6 @@ def interpolator(atmosphere, lfs_atmo): @pytest.fixture -@pytest.mark.usefixtures("lfs_atmo") def atmosphere_grid(atmosphere_name, lfs_atmo): name = lfs_atmo.get(atmosphere_name) atmo = SavFile(name, source=name) diff --git a/test/test_largefilestorage.py b/test/test_largefilestorage.py index f82f75cc..3078048a 100644 --- a/test/test_largefilestorage.py +++ b/test/test_largefilestorage.py @@ -12,7 +12,7 @@ def lfs_available(): return r.status_code == 200 -skipif_lfs = pytest.mark.skipif(lfs_available(), reason="LFS not available") +skipif_lfs = pytest.mark.skipif(not lfs_available(), reason="LFS not available") @pytest.fixture diff --git a/test/test_sme_structure.py b/test/test_sme_structure.py index 5b896c7f..536dac6b 100644 --- a/test/test_sme_structure.py +++ b/test/test_sme_structure.py @@ -5,6 +5,7 @@ import numpy as np import pytest +from pysme.iliffe_vector import Iliffe_vector from pysme.sme import SME_Structure as SME_Struct @@ -140,3 +141,27 @@ def test_fitresults(): sme.fitresults.chisq = 100 sme.fitresults.clear() assert sme.fitresults.chisq is None + + +def test_wint_single_segment_array_input(): + sme = SME_Struct() + wint = np.linspace(5000.0, 5001.0, 11) + + sme.wint = wint + + assert isinstance(sme.wint, Iliffe_vector) + assert sme.wint.nseg == 1 + assert np.allclose(sme.wint[0], wint) + + +def test_wint_multi_segment_list_input(): + sme = SME_Struct() + wint0 = np.linspace(5000.0, 5001.0, 5) + wint1 = np.linspace(6000.0, 6002.0, 7) + + sme.wint = [wint0, wint1] + + assert isinstance(sme.wint, Iliffe_vector) + assert sme.wint.nseg == 2 + assert np.allclose(sme.wint[0], wint0) + assert np.allclose(sme.wint[1], wint1) diff --git a/test/test_synthesis.py b/test/test_synthesis.py index 66695ac9..468ddb5a 100644 --- a/test/test_synthesis.py +++ b/test/test_synthesis.py @@ -4,7 +4,8 @@ import pytest from pysme.iliffe_vector import Iliffe_vector -from pysme.synthesize import synthesize_spectrum +from pysme.sme import SME_Structure as SME_Struct +from pysme.synthesize import Synthesizer, synthesize_spectrum def test_synthesis_simple(sme_2segments): @@ -50,3 +51,82 @@ def test_synthesis_segment(sme_2segments): assert sme2.wave.shape[1][1] != 0 assert np.all(sme2.synth[0] == orig) + + +class _DummyDLL: + def __init__(self): + self.last_wave = "unset" + + def InputWaveRange(self, *_): + return None + + def Opacity(self): + return None + + def Transf(self, mu, accrt, accwi, keep_lineop, wave=None): + self.last_wave = wave + if wave is None: + wint = np.linspace(5000.0, 5001.0, 5) + else: + wint = np.asarray(wave, dtype=float) + sint = np.ones((len(mu), len(wint)), dtype=float) + cint = np.ones((len(mu), len(wint)), dtype=float) + return len(wint), wint, sint, cint + + def CentralDepth(self, mu, accrt): + return np.zeros(0, dtype=float) + + def GetLineRange(self): + return np.zeros((0, 2), dtype=float) + + +def _minimal_sme(): + sme = SME_Struct() + sme.wran = [[5000.0, 5001.0]] + sme.vsini = 0.0 + sme.vmac = 0.0 + sme.vrad_flag = "none" + sme.specific_intensities_only = True + sme.normalize_by_continuum = True + return sme + + +def test_synthesize_segment_prefers_user_wint_over_cache(): + dll = _DummyDLL() + synth = Synthesizer(dll=dll) + sme = _minimal_sme() + + user_wint = np.linspace(5000.0, 5001.0, 7) + cached_wint = np.linspace(5000.0, 5001.0, 9) + sme.wint = user_wint + synth.wint[0] = cached_wint.copy() + + synth.synthesize_segment(sme, 0, reuse_wavelength_grid=True) + + assert np.allclose(dll.last_wave, user_wint) + assert np.allclose(synth.wint[0], cached_wint) + + +def test_synthesize_segment_uses_cache_when_user_wint_missing(): + dll = _DummyDLL() + synth = Synthesizer(dll=dll) + sme = _minimal_sme() + + cached_wint = np.linspace(5000.0, 5001.0, 9) + synth.wint[0] = cached_wint + + synth.synthesize_segment(sme, 0, reuse_wavelength_grid=True) + + assert np.allclose(dll.last_wave, cached_wint) + + +def test_synthesize_segment_populates_cache_when_no_wint_available(): + dll = _DummyDLL() + synth = Synthesizer(dll=dll) + sme = _minimal_sme() + + synth.synthesize_segment(sme, 0, reuse_wavelength_grid=True) + + assert dll.last_wave is None + assert 0 in synth.wint + assert np.allclose(synth.wint[0], np.linspace(5000.0, 5001.0, 5))