# -*- coding: utf-8 -*-
""" Wrapper for SME C library """
import logging
import os, sys
from ctypes import cdll
from importlib import import_module
from os.path import normpath
from .smelib import libtools
import numpy as np
from .smelib.libtools import get_full_datadir, get_full_libfile, _parse_needed_arch_from_error, download_smelib
logger = logging.getLogger(__name__)
_smelib = None
CURRENT_LIB = None
[docs]
def ensure_smelib_ready(libfile=None):
"""Load and initialize the native SME library on first use."""
global _smelib
global CURRENT_LIB
libfile = libfile or get_full_libfile()
if _smelib is not None and CURRENT_LIB == libfile:
return _smelib
if not os.path.exists(libfile):
download_smelib()
try:
cdll.LoadLibrary(libfile)
except OSError as e:
msg = str(e)
if sys.platform == "darwin" and (
"incompatible architecture" in msg or "mach-o file" in msg
):
need = _parse_needed_arch_from_error(msg)
print("Detected arch mismatch; error says need:", need or "(unknown)")
if need is not None:
download_smelib(force_arch=need)
libtools.compile_interface()
cdll.LoadLibrary(libfile)
else:
raise
else:
raise
try:
_smelib = import_module(".smelib._smelib", package="pysme")
except Exception:
libtools.compile_interface()
_smelib = import_module(".smelib._smelib", package="pysme")
CURRENT_LIB = libfile
return _smelib
[docs]
def reload_lib(libfile):
global _smelib
global CURRENT_LIB
target = libfile or get_full_libfile()
if target != CURRENT_LIB:
_smelib = None
ensure_smelib_ready(target)
[docs]
class SME_DLL:
"""Object Oriented interface for the SME C library"""
def __init__(self, libfile=None, datadir=None):
self.libfile = libfile
reload_lib(libfile)
if hasattr(_smelib, "SetHlinopWarningMode"):
_smelib.SetHlinopWarningMode(1)
self.SetLibraryPath(datadir)
self.check_data_files_exist()
@property
def datadir(self):
"""str: Expected directory of the data files"""
return self.GetLibraryPath()
@datadir.setter
def datadir(self, value):
self.SetLibraryPath(value)
@property
def file(self):
"""str: (Expected) Location of the library file"""
return self.libfile or CURRENT_LIB or get_full_libfile()
[docs]
def check_data_files_exist(self):
"""
Checks if required data files for the SME library exist.
If they dont exist, SME will just segfault, without any hint.
Raises
------
FileNotFoundError
If any of the files don't exist
"""
names = self.GetDataFiles()
directory = self.GetLibraryPath()
for name in names:
n = os.path.join(directory, name)
if not os.path.exists(n):
raise FileNotFoundError(
"Could not find required data file {name} in library directory {directory}".format(
name=name, directory=directory
)
)
[docs]
def SMELibraryVersion(self):
"""Return SME library version"""
return _smelib.LibraryVersion()
[docs]
def GetLibraryPath(self):
"""Get the data file directory"""
return _smelib.GetLibraryPath()
[docs]
def GetDataFiles(self):
"""Get the required data files"""
files = _smelib.GetDataFiles()
return files.split(";")
[docs]
def SetLibraryPath(self, libpath=None):
"""Set the path to the library"""
if libpath is None:
libpath = get_full_datadir()
libpath = normpath(libpath) + os.sep
_smelib.SetLibraryPath(libpath)
[docs]
def SetVWscale(self, gamma6):
"""
Set van der Waals scaling factor
Parameters
----------
gamma6 : float
van der Waals scaling factor
"""
_smelib.SetVWscale(gamma6)
[docs]
def SetH2broad(self, h2_flag=True):
"""Set flag for H2 molecule"""
if h2_flag:
_smelib.SetH2broad()
else:
_smelib.ClearH2broad()
[docs]
def ClearH2broad(self):
"""Clear flag for H2 molecule"""
self.SetH2broad(False)
[docs]
def SetHlinopWarningMode(self, mode):
"""Set HLINPROF->HLINOP warning mode (0=stderr, 1=record-only, 2=off)."""
if hasattr(_smelib, "SetHlinopWarningMode"):
_smelib.SetHlinopWarningMode(int(mode))
[docs]
def GetHlinopWarnings(self):
"""Return and clear the last HLINPROF->HLINOP warning summary, if any."""
if hasattr(_smelib, "GetHlinopWarnings"):
return _smelib.GetHlinopWarnings()
return ""
def _log_hlinop_warnings(self):
msg = self.GetHlinopWarnings()
if msg:
logger.warning("%s", msg)
[docs]
def SetLineInfoMode(self, mode):
"""Set handling mode for precomputed line info (0=internal, 1=use_if_valid, 2=trust)."""
_smelib.SetLineInfoMode(int(mode))
[docs]
def OutputLineList(self):
"""
Return line list
Returns
-------
atomic : array of size (nlines, 6)
relevant data of the linelist
wlcent, excit, gflog, gamrad, gamqst, gamvw
"""
return _smelib.OutputLineList()
[docs]
def UpdateLineList(self, atomic, species, index):
"""
Change line list parameters
Parameters
---------
atomic : array of size (nlines, 8)
atomic linelist data for each line
fields are: atom_number, ionization, wlcent, excit, gflog, gamrad, gamqst, gamvw
species : array(string) of size (nlines,)
names of the elements (with Ionization level)
index : array(int) of size (nlines,)
indices of the lines to update relative to the overall linelist
"""
index = np.asarray(index, dtype=np.int16)
species = species[index].astype("S8")
atomic = atomic[index].T
_smelib.UpdateLineList(
species,
atomic,
index,
)
[docs]
def Opacity(self):
"""Calculate opacities"""
_smelib.Opacity()
[docs]
def GetOpacity(self, switch, species=None, key=None):
"""
Returns specific cont. opacity, different output depending on the input
Parameters
----------
switch : str
one of [COPSTD, COPRED, COPBLU, AHYD, AH2P, AHMIN, SIGH, AHE1, AHE2,
AHEMIN, SIGHE, ACOOL, ALUKE, AHOT, SIGEL, SIGH2]
key : str, optional
for ACOOL, one of [new, old, fraction]
species : str, optional
for ACOOL and ALUKE it specifies the element
ACOOL: C1, Mg1, Al1, Si1, Fe1, CH, NH, OH
ALUKE: N1, O1, Mg2, Si2, Ca2
"""
kwargs = {}
if key is not None:
kwargs["key"] = key
if species is not None:
kwargs["species"] = species
return _smelib.GetOpacity(switch, **kwargs)
[docs]
def Ionization(self, ion=0):
"""
Calculate ionization balance for current atmosphere and abundances.
Ionization state is stored in the external library.
Ion is a bit flag with values (add them together to use multiple):
:1: adopt particle number densities from EOS
:2: adopt electron number densities from EOS
:4: adopt gas densities (g/cm^3) from EOS
instead of using values from model atmosphere. Different abundance patterns
in the model atmosphere (usually scaled solar) and SME (may be non-solar)
can affect line shape, e.g. shape of hydrogen lines.
Parameters
----------
ion : int
flag that determines the behaviour of the C function
"""
_smelib.Ionization(ion)
[docs]
def GetDensity(self):
"""
Retrieve density in each layer
Returns
-------
density : array of size (ndepth,)
Density of the atmosphere in each layer
"""
return _smelib.GetDensity()
[docs]
def GetNatom(self):
"""
Get XNA
Returns
-------
XNA : array of size (ndepth,)
XNA in each layer
"""
return _smelib.GetNatom()
[docs]
def GetNelec(self):
"""
Get XNE (Electron number density) for each layer in the atmosphere
Returns
-------
XNE : array of size (ndepth,)
XNE in each layer
"""
return _smelib.GetNelec()
[docs]
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)
[docs]
def Transf(
self,
mu,
wave=None,
nwmax=400000,
accrt=1e-3,
accwi=3e-3,
keep_lineop=False,
long_continuum=True,
):
"""
Radiative Transfer Calculation
Perform the radiative transfer calculation thorugh the atmosphere
Requires that all parameters have been set beforehand
Parameters
---------
mu : array of shape (nmu,)
mu angles (1 - cos(phi)) of different limb points along the stellar surface
accrt : float
accuracy of the radiative transfer integration
accwi : float
accuracy of the interpolation on the wavelength grid
keep_lineop : bool, optional
if True do not recompute the line opacities (default: False)
long_continuum : bool, optional
if True the continuum is calculated at every wavelength (default: True)
nwmax : int, optional
maximum number of wavelength points if wavelength grid is not set with wave (default: 400000)
wave : array, optional
wavelength grid to use for the calculation,
if not set will use an adaptive wavelength grid with no constant step size (default: None)
Returns
-------
nw : int
number of actual wavelength points, i.e. size of wint_seg
wint_seg : array of shape (nwave,)
wavelength grid, the number of wavelengthpoints is equal to the number of lines * 2 - 1
One point in the center of each line + plus one between the next line
sint_seg : array of shape (nmu, nwave)
spectrum for each mu point
cint_seg : array of shape (nmu, nwave)
continuum for each mu point
"""
keep_lineop = 1 if keep_lineop else 0
long_continuum = 1 if long_continuum else 0
# keywords = {"mu", "wave", "nwmax", "accrt", "accwi", "keep_lineop", "long_continuum"}
nw, wave, sint, cint = _smelib.Transf(
mu, wave, nwmax, accrt, accwi, keep_lineop, long_continuum
)
self._log_hlinop_warnings()
# Resize the arrays
wave = wave[:nw]
sint = sint[:nw, :].T
cint = cint[:nw, :].T
sint = np.nan_to_num(sint, copy=False)
cint = np.nan_to_num(cint, copy=False)
return nw, wave, sint, cint
[docs]
def CentralDepth(self, mu, accrt):
"""
This subroutine explicitly solves the transfer equation
for a set of nodes on the star disk in the centers of spectral
lines. The results are specific intensities.
Parameters
----------
mu : array of size (nmu,)
mu values along the stellar disk to calculate
accrt : float
precision of the radiative transfer calculation
Returns
-------
table : array of size (nlines,)
Centeral depth (i.e. specific intensity) of each line
"""
table = _smelib.CentralDepth(mu, accrt)
self._log_hlinop_warnings()
return table
[docs]
def ALMAXRange(self, accrt=1e-4):
"""Compute first-stage ALMAX and line ranges from SMElib preselection logic.
Parameters
----------
accrt : float, optional
Opacity-ratio threshold (same role as accrt in Transf preselection).
Returns
-------
almax : array of size (nlines,)
Line-center opacity-ratio proxy from the first-stage screening.
linerange : array of size (nlines, 2)
Preselection range for each line.
"""
result = _smelib.ALMAXRange(accrt=accrt)
self._log_hlinop_warnings()
return result
[docs]
def GetLineOpacity(self, wave):
"""
Retrieve line opacity data from the C library
Parameters
----------
wave : float
Wavelength of the line opacity to retrieve
Returns
---------
lop : array
line opacity
cop : array
continuum opacity including scatter
scr : array
Scatter
tsf : array
Total source function
csf : array
Continuum source function
"""
lop, cop, scr, tsf, csf = _smelib.GetLineOpacity(wave)
return lop, cop, scr, tsf, csf
[docs]
def GetLineRange(self):
"""Get the effective wavelength range for each line
i.e. the wavelengths for which the line has significant impact
Returns
-------
linerange : array of size (nlines, 2)
lower and upper wavelength for each spectral line
"""
return _smelib.GetLineRange()
[docs]
def GetDepartureCoefficients(self, line):
"""Get the NLTE departure coefficients as stored in the C library
Parameters
----------
line : int
requested line index, i.e. between 0 and number of lines
Returns
-------
bmat : array of size (2, nrhox)
departure coefficients for the given line index
"""
return _smelib.GetDepartureCoefficients(line)
[docs]
def ResetDepartureCoefficients(self):
"""Reset departure coefficients from any previous call, to ensure LTE as default"""
_smelib.ResetDepartureCoefficients()
[docs]
def GetNLTEflags(self):
"""Get an array that tells us which lines have been used with NLTE correction
Returns
-------
nlte_flags : array(bool) of size (nlines,)
True if line was used with NLTE, False if line is only LTE
"""
nlte_flags = _smelib.GetNLTEflags()
return nlte_flags.astype(bool)
[docs]
def GetContributionfunction(self, mu, wave, nwmax=400000, long_continuum=True):
"""Get the contribution function for a wavelength grid.
Returns
-------
"""
long_continuum_val = 1 if long_continuum else 0
res = _smelib.ContributionFunctions(mu, wave, nwmax, long_continuum_val)
table, ctable = res
if long_continuum:
return table, ctable
else:
return table
