Plot spectra example

This example demonstrates how to extract a spectra directly from a grid and plots all the available spectra.

NOTE: this only works on 2D grids at the moment

----------
SUMMARY OF SED
Number of wavelength points: 9244
Wavelength range: [0.00 Å,             299293000000.00 Å]
log10(Peak luminosity/erg/(Hz*s)):             21.49
log10(Bolometric luminosity/erg/s):36.84609875911521----------
/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/unyt/array.py:1824: RuntimeWarning: divide by zero encountered in log10
  out_arr = func(np.asarray(inp), out=out_func, **kwargs)
----------
SUMMARY OF SED
Number of wavelength points: 9244
Wavelength range: [0.00 Å,             299293000000.00 Å]
log10(Peak luminosity/erg/(Hz*s)):             24.41
log10(Bolometric luminosity/erg/s):36.166199101822286----------
----------
SUMMARY OF SED
Number of wavelength points: 9244
Wavelength range: [0.00 Å,             299293000000.00 Å]
log10(Peak luminosity/erg/(Hz*s)):             24.44
log10(Bolometric luminosity/erg/s):36.59752277341082----------
----------
SUMMARY OF SED
Number of wavelength points: 9244
Wavelength range: [0.00 Å,             299293000000.00 Å]
log10(Peak luminosity/erg/(Hz*s)):             21.39
log10(Bolometric luminosity/erg/s):36.487843672084175----------
----------
SUMMARY OF SED
Number of wavelength points: 9244
Wavelength range: [0.00 Å,             299293000000.00 Å]
log10(Peak luminosity/erg/(Hz*s)):             24.44
log10(Bolometric luminosity/erg/s):36.84716642614972----------
----------
SUMMARY OF SED
Number of wavelength points: 9244
Wavelength range: [0.00 Å,             299293000000.00 Å]
log10(Peak luminosity/erg/(Hz*s)):             23.26
log10(Bolometric luminosity/erg/s):36.39658434503608----------

import argparse

import matplotlib.pyplot as plt
import numpy as np
from synthesizer.grid import Grid

if __name__ == "__main__":
    # Get the location of this script, __file__ is the absolute path of this
    # script, however we just want to directory
    # script_path = os.path.abspath(os.path.dirname(__file__))

    # define the test grid dir
    # test_grid_dir = script_path + "/../../tests/test_grid/"
    test_grid_dir = "../../tests/test_grid/"

    # initialise argument parser
    parser = argparse.ArgumentParser(
        description=(
            "Create a plot of all spectra types for a given metallicity and \
            age"
        )
    )

    # The name of the grid. Defaults to the test grid.
    parser.add_argument(
        "-grid_name",
        "--grid_name",
        type=str,
        required=False,
        default="test_grid",
    )

    # The path to the grid directory. Defaults to the test grid directory.
    parser.add_argument(
        "-grid_dir",
        "--grid_dir",
        type=str,
        required=False,
        default=test_grid_dir,
    )

    # The target metallicity. The code function will find the closest
    # metallicity and report it back. The rationale behind this is
    # that this code can easily be adapted to explore other grids.
    parser.add_argument(
        "-metallicity", type=float, required=False, default=0.01
    )

    # The target log10(age/yr). The code function will find the closest
    # metallicity and report it back. The rationale behind this is that
    # this code can easily be adapted to explore other grids.
    parser.add_argument("-log10age", type=float, required=False, default=6.0)

    # Get dictionary of arguments
    args = parser.parse_args()

    # initialise grid
    grid = Grid(args.grid_name, grid_dir=args.grid_dir)

    # get the grid point for this log10age and metallicity
    grid_point = grid.get_grid_point((args.log10age, args.metallicity))

    # loop over available spectra and plot
    for spec_name in grid.available_spectra:
        # get Sed object
        sed = grid.get_spectra(grid_point, spectra_id=spec_name)
        # print summary of SED object
        print(sed)
        plt.plot(
            np.log10(sed.lam),
            np.log10(sed.lnu),
            lw=1,
            alpha=0.8,
            label=spec_name,
        )

    plt.xlim([2.0, 4.0])
    plt.ylim([18.0, 23])
    plt.legend(fontsize=8, labelspacing=0.0)
    plt.xlabel(r"$\rm log_{10}(\lambda/\AA)$")
    plt.ylabel(r"$\rm log_{10}(L_{\nu}/erg\ s^{-1}\ Hz^{-1} M_{\odot}^{-1})$")
    plt.show()