I know that this has come up in conversation at some point, but I think it would be really useful to use astropy units throughout this package, particularly at the user interface layer.

I fully realize that introduction of units is not a trivial addition and can also involve some performance penalties as well. However, the units on some inputs/outputs are sufficiently ambiguous to warrant this. For example, I was modifying the z component of the velocity on the FAL atmosphere and I was not sure what the components on the velocity were (I believe they are SI, m/s).

Have a kwarg on LwContext.formal_sol_gamma_matrices that tells it not to recalculate collisional rates. (Maybe make sure they've been computed once...?)

This doesn't work if we simply want to do something of the form for a in radSet. Need to handle more like AtomicStateTable.

This probably also applies to SpeciesStateTable too.

Seems that all 3 directions need to specified or a buffer error happens when the Context is being set up.

i.e.
atmos.rays(muz=1.0, wmu=0.0)
doesn't work, but
atmos.rays(muz=[1.0], muy=[0.0], mux=[0.0], wmu=[1.0]) does.

If an int is passed into atmos.rays instead of a float, e.g. due to the case of ctx.compute_rays(wave, 1), it is not caught by the non-iterable check, as this checks isinstance(muz, float). This cascades onwards causing a number of issues, resulting in failure to create a new Context.
Further if one instead puts a list with a singular int, different but similar issues arise due to types getting carried over somewhere, and being incompatible with the backend memory views (e.g. ctx.compute_rays(wave, [1])).

The best fix here is likely to be using something like np.asarray(muz, dtype=np.float64) instead of the current checks (and then ensuring that it's 1- and not 0-dimensional).

For now this issue can be worked around by only passing floats to the mu argument: e.g. ctx.compute_rays(wave, 1.0).

This would be a nice QoL improvement when experimenting wityh different timesteppers like the Theta method in MsLightweaver

We can damp the Newton-Raphson correction if it's too large -- it's still self-consistent, obviously the process isn't converged till we're done with damping; like the Armijo line search.

Due to rewrite of Bezier3 Scalar FS, the Full Stokes FS currently doesn't work due to some re-ordering of terms. This isn't a big deal with Jaime's derivation and should only a take a few hours to fix next week.

Need to check if there are pre-existing height/cmass/tau arrays in the atmos object, and reuse if there is, to avoid breaking the link to the Context.

Probably do this through a function pointer? They have the same interface anyway.

The Bezier3 solver, even rederived, seems to have issues on very fine grids (RADYN gives us < 5cm in places) and a standard linear seems better behaved here.

Another alternative is to check the size of the step in Bezier and just use linear info when distance between points less than say, 50-100 m?

For install purposes, it would be useful to have a list of which packages are required. I'm setting up a new virtual environment and keep missing which ones are necessary. . .

This is a hassle for reprocessing mu=1.0 snapshots otherwise, as wmu is ill-defined.

The lightweaver directory should contain a __init__.py file so that modules can be imported using the syntax from lightweaver.<module_name> import <function_or_class_name>.

The underlying implementation is relatively slow (not compared to the RH15d implementation), and recreates the entire thread pool which is very expensive (and shouldn't be necessary now due to the change to the use of a View on the atom object which should automatically update).

Maybe see if we can use the thread pool (if already present) to accelerate the calculation of these terms.

Currently install_requires has numpy>=1.19, however, I get the following RuntimeError on import

This error was resolved by upgrading to Numpy v1.24.4.

It would be nice to have the option to do this.

Currently numpy 1.24 isn't supported by numba, and triggers a cryptic error on import:

SystemError: initialization of _internal failed without raising an exception

Can be fixed by installing a numpy from the 1.23 series (pip install --upgrade 'numpy<1.24'), but the version should be pinned for safety for now. Track the issue below for updates.

numba/numba#8464

Test12.py fails because the file `pf_Kurucz.input' is missing from the main distro. It's on one of the branches, though.

If ctx.Nthreads > 1 on fork, then only main thread is copied, and none of the workers are cloned (as per the definition of fork(2)). If Nthreads is set to one and then increased in each child process then all is well. This happens with a Context in global state on creating a ProcessPoolExecutor.

Interestingly if ctx.Nthreads is left > 1, everything still runs (albeit only 1 thread per process). This is due to the way the task library works. It spawns N-1 threads, with the Nth being the main thread. When we call scheduler_join in the main thread, it jumps in and does the work (as if the other workers were busy).

This is all surprisingly close to "as one would/should expect" IMO, however it should be documented somewhere. For the most part if a job can use inter-process concurrency, it is unlikely to need the sub-process multi-threaded formal solutions that ctx.Nthreads exists for, so in most places ctx.Nthreads should be set to 1 before spawning. If multiple threads are needed in the children then there are two simple options off the top of my head: set ctx.Nthreads > 1 and the scheduler will happily spin up threads, or don't rely on the global state fork behaviour, and pass the context through the initialiser options available in ProcessPoolExecutor.

This works transparently for arrays, there's no reason not to have the same for scalars.

Change this to be something like if background and any([...]):, or maybe we should just go straight to an interface that allows different background implementations?

Need to make sure this is run.

This variable isn't used usefully anymore and we waste memory and cpu allocating it and clearing it.

Numba is incompatible with numpy > 1.20.3, but by default the pyproject.toml sources 1.22. This should be fixed with numba 0.55 in a few weeks and all the ABIs should match again. Until then we limit Lightweaver to numpy==1.20.3

I'm trying to adapt the simple example in the gallery to calculate the Stokes parameters for a constant magnetic field (as a function of z). Besides passing stokes=True to compute_rays, all of the code from the example is the same. I've then modified the atmosphere as follows,

atmos_const_b = Falc82()
atmos_const_b.B = 0.05 * np.ones(atmos_const_b.z.shape)
atmos_const_b.gammaB = 0.78539819 * np.ones(atmos_const_b.B.shape)
atmos_const_b.chiB = np.zeros(atmos_const_b.B.shape)

Running synth8542 with this atmosphere, I then get the following exception,

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-147-63eb88a56916> in <module>
----> 1 ctx, IQUV_B_const = synth_8542(
      2     atmos_const_b,
      3     conserve=False,
      4     useNe=False,
      5     wave=wave_caII,

<ipython-input-91-f7bfabc22d02> in synth_8542(atmos, conserve, useNe, wave, stokes)
     73     # compute the final intensity for mu=1 on the provided wavelength grid.
     74     eqPops.update_lte_atoms_Hmin_pops(atmos)
---> 75     Iwave = ctx.compute_rays(wave, [atmos.muz[-1]], stokes=stokes)
     76     return ctx, Iwave
     77 

Source/LwMiddleLayer.pyx in lightweaver.LwCompiled.LwContext.compute_rays()

ValueError: could not broadcast input array from shape (3,1000,1,1) into shape (3,1000)

I'm not really sure where this dimension mismatch is coming from. It looks like some extra dimensions need to be added or squeezed somewhere to account for the fact that this is a 1D atmosphere?

It probably makes sense to allow PRD sub-iterations on these lines. If the atom with fixed populations was computed out of PRD then the LineTypes should probably be converted to LineType.Crd

The imports in the Test2.py script in Samples/do not reflect the current structure of the package. Doing the following gets me most of the way there:

• Created lightweaver/__init__.py (See #2)
• Prepended imports with `lightweaver.<module_name>

However, there's one module that I cannot seem to find anywhere else in the repo:

from PyProto import ComputationalAtom, background, gamma_matrices, stat_equil

I noticed that on some of the other branches, there is a PyProto.py file in the root of the repo. However, I can't seem to find the imported classes/functions (i.e. ComputationalAtom, background, gamma_matrices, stat_equil) on the master branch. Has the workflow for doing these calculations changed?

I'm happy to submit a PR to fix the examples.

Currently the installation process is a little involved, and I find myself needing to set environment variables for most systems. This can hopefully be improved.
Additionally, see if we can use a cloud system to build a nice set of wheels so most users won't need to compile.

more memes plz

This is perhaps less of an issue and more of a feature/docs request. My ignorance here is likely due to my unfamiliarity with Lightweaver, RH, and radiative transfer so I will apologize in advance 😅

I'm working through this toy example which uses RH (from the command line) + a combination of scripts to demonstrate the effect of the velocity on the shape of the absorption line profile around 656.25 nm. For simplicity, I'll focus just on the first part which leaves the velocity unmodified

Starting from the simple gallery example in the docs, my understanding is that I need to modify the H_6 atom that is passed to RadiativeSet in the following way:

h6 = H_6_atom()
h6.lines[4].quadrature.qCore = 25.0
h6.lines[4].quadrature.Nlambda = 200

(the 4 comes from the fact that I want to modify the first line in the Balmer series, j=2, i=1). However, my resulting line profile is

compared to the example produced in RH,

I've included the complete code to reproduce my figure below. I'm sure I'm likely doing something wrong or not completely understanding what I'm doing here!

from lightweaver.fal import Falc82
from lightweaver.rh_atoms import (H_6_atom, CaII_atom, O_atom, C_atom, Si_atom, Al_atom, Fe_atom,
                                  He_atom, Mg_atom, N_atom, Na_atom, S_atom)
import lightweaver as lw
import matplotlib.pyplot as plt
import numpy as np
import time

def synth_8542(atmos, conserve, useNe, wave):
    '''
    Synthesise a spectral line for given atmosphere with different
    conditions.

    Parameters
    ----------
    atmos : lw.Atmosphere
        The atmospheric model in which to synthesise the line.
    conserve : bool
        Whether to start from LTE electron density and conserve charge, or
        simply use from the electron density present in the atomic model.
    useNe : bool
        Whether to use the electron density present in the model as the
        starting solution, or compute the LTE electron density.
    wave : np.ndarray
        Array of wavelengths over which to resynthesise the final line
        profile for muz=1.

    Returns
    -------
    ctx : lw.Context
        The Context object that was used to compute the equilibrium
        populations.
    Iwave : np.ndarray
        The intensity at muz=1 for each wavelength in `wave`.
    '''
    # Configure the atmospheric angular quadrature
    atmos.quadrature(5)
    # Configure the set of atomic models to use.
    # Set some things specific to our example.
    h6 = H_6_atom()
    h6.lines[4].quadrature.qCore = 25.0
    h6.lines[4].quadrature.Nlambda = 200
    aSet = lw.RadiativeSet([
        h6,
        C_atom(),
        O_atom(),
        Si_atom(),
        Al_atom(),
        CaII_atom(),
        Fe_atom(),
        He_atom(),
        Mg_atom(),
        N_atom(),
        Na_atom(),
        S_atom(),
    ])
    # Set H and Ca to "active" i.e. NLTE, everything else participates as an
    # LTE background.
    aSet.set_active('H', 'Ca')
    # Compute the necessary wavelength dependent information (SpectrumConfiguration).
    spect = aSet.compute_wavelength_grid()

    # Either compute the equilibrium populations at the fixed electron density
    # provided in the model, or iterate an LTE electron density and compute the
    # corresponding equilibrium populations (SpeciesStateTable).
    if useNe:
        eqPops = aSet.compute_eq_pops(atmos)
    else:
        eqPops = aSet.iterate_lte_ne_eq_pops(atmos)

    # Configure the Context which holds the state of the simulation for the
    # backend, and provides the python interface to the backend.
    # Feel free to increase Nthreads to increase the number of threads the
    # program will use.
    ctx = lw.Context(atmos, spect, eqPops, conserveCharge=conserve, Nthreads=1)
    start = time.time()
    # Iterate the Context to convergence
    iterate_ctx(ctx)
    end = time.time()
    print('%.2f s' % (end - start))
    # Update the background populations based on the converged solution and
    # compute the final intensity for mu=1 on the provided wavelength grid.
    eqPops.update_lte_atoms_Hmin_pops(atmos)
    Iwave = ctx.compute_rays(wave, [atmos.muz[-1]], stokes=False)
    return ctx, Iwave


def iterate_ctx(ctx, Nscatter=3, NmaxIter=500):
    '''
    Iterate a Context to convergence.
    '''
    for i in range(NmaxIter):
        # Compute the formal solution
        dJ = ctx.formal_sol_gamma_matrices()
        # Just update J for Nscatter iterations
        if i < Nscatter:
            continue
        # Update the active populations under statistical equilibrium,
        # conserving charge if this option was set on the Context.
        delta = ctx.stat_equil()

        # If we are converged in both relative change of J and populations,
        # then print a message and return
        # N.B. as this is just a simple case, there is no checking for failure
        # to converge within the NmaxIter. This could be achieved simpy with an
        # else block after this for.
        if dJ < 3e-3 and delta < 1e-3:
            print('%d iterations' % i)
            print('-'*80)
            return


atmos_rest = Falc82()
wave = np.linspace(656.1, 656.5, 1000)
ctx, Iwave_rest = synth_8542(atmos_rest, conserve=False, useNe=False, wave=wave)
plt.plot(wave, Iwave_rest)

Hi Chris,
I noticed that in your example you just show that Lightweaver can give specific intensity along with their corresponding quadrature and integration weights. And I know some other RT codes like Turbospectrum and MULTI can output flux so I am curious if LW can give that as well? If not, can user calculate flux based on current LW output?

Thanks!
Yangyang

This is important for MsLightweaver style stuff, and is 80+% done.