diff --git a/README.md b/README.md index 895c510..8aad93f 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ This is the official release of the *nuspacesim* simulator tool! This package simulates upward-going extensive air showers caused by neutrino interactions with the atmosphere. It calculates the tau neutrino acceptance for the -Optical Cherenkov technique. The simulation is parameterized by an input XML +Optical Cherenkov technique. The simulation is parameterized by an input TOML configuration file, with settings for detector characteristics and global parameters. The package also provides a python3 API for programatic access. @@ -31,18 +31,18 @@ Tau propagation is interpolated using included data tables from [nupyprop](https ![NuSpaceSim Usage](https://raw.githubusercontent.com/NuSpaceSim/nuSpaceSim/main/docs/_static/run.svg) -### Create an XML configuration script +### Create an TOML configuration script -The command line simulator uses an XML file to store configuration parameters. To +The command line simulator uses an TOML file to store configuration parameters. To generate a default configuration file run the following, with your choice of file name. -`nuspacesim create-config my_config_file.xml` +`nuspacesim create-config my_config_file.toml` ### Run simulator Simulate neutrino interactions and save the results to a named fits file. -`nuspacesim run my_config_file.xml -o my_nss_sim.fits` +`nuspacesim run my_config_file.toml -o my_nss_sim.fits` # Documentation diff --git a/docs/conf.py b/docs/conf.py index 5993a69..1f81e1e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -92,7 +92,6 @@ "click", "dask", "h5py", - "lxml", "matplotlib", "mpl_toolkits", "numpy", diff --git a/docs/index.rst b/docs/index.rst index 0808c42..b30071e 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -13,7 +13,7 @@ This is the public release of the nuspacesim simulator tool! This package simulates upward-going electromagnetic air showers caused by neutrino interactions with the atmosphere. It calculates the tau neutrino acceptance for the -Optical Cherenkov technique. The simulation is parameterized by an input XML +Optical Cherenkov technique. The simulation is parameterized by an input TOML configuration file, with settings for detector characteristics and global parameters. The package also provides a python3 API for programatic access. diff --git a/docs/quickstart.rst b/docs/quickstart.rst index 20da08c..9862c47 100644 --- a/docs/quickstart.rst +++ b/docs/quickstart.rst @@ -27,13 +27,13 @@ Read the nuspacesim help docstring python3 -m nuspacesim --help ---------------------------------- -Create the XML configuration file +Create the TOML configuration file ---------------------------------- Create a configuration file with the ``create-config`` command. This is editable by the user for defining different simulation parameters. :: - nuspacesim create-config --numtrajs 1e6 --monospectrum 10.25 my_config_file.xml + nuspacesim create-config --numtrajs 1e6 --monospectrum 10.25 my_config_file.toml ----------------- Run the simulator @@ -42,11 +42,11 @@ Run the simulator Simulate neutrino interactions, and extensive air showers, then save the results to a FITS file. :: - nuspacesim run my_config_file.xml -output my_nss_sim.fits + nuspacesim run my_config_file.toml -output my_nss_sim.fits Optionally, override the configuration file on the command line. :: - nuspacesim run my_config_file.xml 1e5 --powerspectrum 2 6 12 -o my_nss_sim.fits + nuspacesim run my_config_file.toml 1e5 --powerspectrum 2 6 12 -o my_nss_sim.fits .. raw:: html diff --git a/docs/tutorial/command_line.rst b/docs/tutorial/command_line.rst index 4a8c822..62a2150 100644 --- a/docs/tutorial/command_line.rst +++ b/docs/tutorial/command_line.rst @@ -12,12 +12,12 @@ in application or module mode. nuspacesim --help -Create an XML configuration file. +Create an TOML configuration file. :: nuspacesim create-config --help - nuspacesim create-config my_config_file.xml + nuspacesim create-config my_config_file.TOML Run a simulation. diff --git a/docs/tutorial/configuration.rst b/docs/tutorial/configuration.rst index fdf654c..ee7c80b 100644 --- a/docs/tutorial/configuration.rst +++ b/docs/tutorial/configuration.rst @@ -7,14 +7,14 @@ Configuration File .. toctree:: :hidden: -Simulation configuration is governed by an XML file. The File structure is +Simulation configuration is governed by a TOML file. The File structure is separated into 2 primary sections, the DetectorCharacteristics and the -SimulationParameters. Both XML sections map 1-to-1 to the nuspacesim.config +SimulationParameters. Both TOML sections map 1-to-1 to the nuspacesim.config dataclass objects of the same names. -DetectorCharacteristics +Detector *********************** This is a dataclass holding the Detector Characteristics for a given simulation. @@ -34,7 +34,7 @@ The member attributes are as follows: * **det_gain**: Antenna gain in dB: Default = 1.8 -SimulationParameters +Simulation ******************** This is a dataclass holding the Detector Characteristics for a given simulation. diff --git a/docs/tutorial/simulation.rst b/docs/tutorial/simulation.rst index 2b7b7fe..926259a 100644 --- a/docs/tutorial/simulation.rst +++ b/docs/tutorial/simulation.rst @@ -12,7 +12,7 @@ variables determined by the physics modules. Configuration File ****************** -Simulation configuration is governed by an XML file. +Simulation configuration is governed by a TOML file. Results File ************ diff --git a/recipe/meta.yaml b/recipe/meta.yaml index 962c366..37dea72 100644 --- a/recipe/meta.yaml +++ b/recipe/meta.yaml @@ -37,7 +37,6 @@ requirements: - importlib_metadata # [py<38] - importlib_resources # [py<39] - h5py - - lxml - matplotlib - numpy >=1.20 - rich diff --git a/sample_input_file.toml b/sample_input_file.toml new file mode 100644 index 0000000..aee8f6c --- /dev/null +++ b/sample_input_file.toml @@ -0,0 +1,45 @@ +title = "NuSpaceSim" + +[detector] +name = "Default Name" + +[detector.initial_position] +altitude = "525.0 km" +latitude = "0.0 deg" +longitude = "0.0 deg" + +[detector.optical] +telescope_effective_area = "2.5 m2" +quantum_efficiency = 0.2 +photo_electron_threshold = 10 + +[detector.radio] +low_frequency = "30.0 MHz" +high_frequency = "300.0 MHz" +snr_threshold = 5.0 +nantennas = 10 +gain = "1.8 dB" + +[simulation] +mode = "Diffuse" +thrown_events = 10000 +max_cherenkov_angle = "3.0000000000000004 deg" +max_azimuth_angle = "360.0 deg" +angle_from_limb = "7.0 deg" +cherenkov_light_engine = "Default" + +[simulation.ionosphere] +total_electron_content = 10.0 +total_electron_error = 0.1 + +[simulation.tau_shower] +id = "nupyprop" +etau_frac = 0.5 +table_version = "3" + +[simulation.spectrum] +id = "monospectrum" +log_nu_energy = 8.0 + +[simulation.cloud_model] +id = "no_cloud" diff --git a/sample_input_file.xml b/sample_input_file.xml deleted file mode 100644 index a128bdf..0000000 --- a/sample_input_file.xml +++ /dev/null @@ -1,39 +0,0 @@ - - - - 0.2 - 2.5 - - 10.0 - - 525.0 - 0.0 - 0.0 - 30.0 - 300.0 - 5.0 - 10 - 1.8 - - - 0.05235987755982989 - 0.12217304763960307 - - 0.5 - - - - 8.0 - - - - - - 6.283185307179586 - 100 - 0 - 10.0 - 0.1 - 3 - - diff --git a/setup.cfg b/setup.cfg index f58e95d..264aee3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -31,15 +31,16 @@ install_requires = dask dask[distributed] h5py - lxml matplotlib numpy>=1.21 + pydantic rich scipy - cached-property;python_version<"3.8" + tomli-w configparser;python_version<"3.8" importlib-metadata;python_version<"3.8" importlib-resources;python_version<"3.9" + tomli;python_version<"3.11" python_requires = >=3.8 ext_package = nuspacesim package_dir = diff --git a/src/nuspacesim/__init__.py b/src/nuspacesim/__init__.py index 50e156d..33b52c9 100644 --- a/src/nuspacesim/__init__.py +++ b/src/nuspacesim/__init__.py @@ -74,7 +74,7 @@ nuspacesim.NssConfig nuspacesim.DetectorCharacteristics nuspacesim.SimulationParameters - nuspacesim.xml_config + nuspacesim.config ***************** Simulate Function @@ -131,18 +131,17 @@ """ -from . import constants, data, utils, xml_config +from . import constants, data, utils from ._version import version, version_tuple from .compute import compute -from .config import DetectorCharacteristics, NssConfig, SimulationParameters +from .config import NssConfig, config_from_toml from .simulation import eas_optical, geometry, taus __all__ = [ # Core "constants", "NssConfig", - "DetectorCharacteristics", - "SimulationParameters", + "config_from_toml", "compute", # modules "geometry", @@ -151,7 +150,6 @@ # other "data", "utils", - "xml_config", # version "version", "version_tuple", diff --git a/src/nuspacesim/apps/__init__.py b/src/nuspacesim/apps/__init__.py index f4fc4ec..7165e22 100644 --- a/src/nuspacesim/apps/__init__.py +++ b/src/nuspacesim/apps/__init__.py @@ -15,6 +15,6 @@ """ -__all__ = ["cli"] +__all__ = ["cli", "run", "create_config", "show_plot", "utils"] -from . import cli +from . import cli, create_config, run, show_plot, utils diff --git a/src/nuspacesim/apps/__main__.py b/src/nuspacesim/apps/__main__.py new file mode 100644 index 0000000..4cafccb --- /dev/null +++ b/src/nuspacesim/apps/__main__.py @@ -0,0 +1,3 @@ +from .cli import cli + +cli() diff --git a/src/nuspacesim/apps/cli.py b/src/nuspacesim/apps/cli.py index febc6e2..fca786e 100644 --- a/src/nuspacesim/apps/cli.py +++ b/src/nuspacesim/apps/cli.py @@ -46,341 +46,16 @@ """ -import configparser - import click -from astropy.table import Table as AstropyTable - -from .. import NssConfig, SimulationParameters, simulation -from ..compute import compute -from ..config import FileSpectrum, MonoSpectrum, PowerSpectrum -from ..types.cloud_types import MonoCloud, NoCloud, PressureMapCloud -from ..utils import plots -from ..utils.plot_function_registry import registry -from ..xml_config import config_from_xml, create_xml -__all__ = ["create_config", "run", "show_plot"] +from . import create_config, run, show_plot @click.group() -# @click.option("--debug/--no-debug", default=False) def cli(): pass - # def cli(ctx, debug): - # # ctx.ensure_object(dict) - # # ctx.obj["DEBUG"] = debug - - -@cli.command() -@click.option( - "-o", "--output", type=click.Path(exists=False), default=None, help="Output file." -) -@click.option( - "-p", - "--plot", - type=click.Choice(list(registry), case_sensitive=False), - multiple=True, - default=[], - help="Available plotting functions. Select multiple plots with multiple uses of -p", -) -@click.option( - "-P", - "--plotconfig", - type=click.Path( - exists=True, - dir_okay=False, - readable=True, - ), - help="Read selected plotting functions and options from the specified INI file", -) -@click.option("--plotall", is_flag=True, help="Show all result plots.") -@click.option( - "-w", - "--write-stages", - is_flag=True, - help="Write intermediate values after each simulation stage.", -) -@click.option( - "-n", - "--no-result-file", - is_flag=True, - help="Do not save the results to an output file.", -) -@click.option( - "--monospectrum", - type=float, - default=None, - help="Mono Energetic Spectrum Log Energy.", -) -@click.option( - "--powerspectrum", - nargs=3, - type=click.Tuple([float, float, float]), - default=None, - help="Power Spectrum index, lower_bound, upper_bound.", -) -@click.option( - "--nocloud", - is_flag=True, - default=None, - help="No Cloud Model. [Default]", -) -@click.option( - "--monocloud", - type=float, - default=None, - help="Uniform (mono) Height Cloud Model (km).", -) -@click.option( - "--pressuremapcloud", - type=click.DateTime(["%m", "%B", "%b"]), - default=None, - help="Pressure Map Cloud Model (built in and included with NuSpaceSim). This map is" - "an instance of a global cloud top pressure map sampled from a model of all days in" - "the given month over a 10-year time period from 2011 to 2020. Data provided by the" - "MERRA-2 dataset." - "User should provide a month name, abbreviation, or number.", -) -@click.argument( - "config_file", - default=None, - type=click.Path(exists=True, dir_okay=False, readable=True), -) -@click.argument("count", type=float, default=0.0) -def run( - config_file: str, - count: float, - no_result_file: bool, - output: str, - plot: list, - plotconfig: str, - plotall: bool, - write_stages: bool, - monospectrum, - powerspectrum, - nocloud, - monocloud, - pressuremapcloud, -) -> None: - """Perform the full nuspacesim simulation. - - Main Simulator for nuspacesim. Given a XML configuration file, and - optionally a count of simulated nutrinos, runs nutrino simulation. - - \f - - Parameters - ---------- - config_file: str - XML configuration file for particular simulation particular. - count : int, optional - Number of thrown trajectories. Optionally override value in config_file. - spectrum_type : str, optional - Type of - output: str, optional - Name of the output file holding the simulation results. - plot: list, optional - Plot the simulation results. - plotconfig: str, optional - INI file to select plots for each module, as well as to specifiy global plot settings. - plotall: bool, optional - Plot all the available the simulation results plots. - no_result_file: bool, optional - Disables writing results to output files. - - - Examples - -------- - Command line usage of the run command may work with the following invocation. - - `nuspacesim run sample_input_file.xml 1e5 8 -o my_sim_results.fits` - """ - - # User Inputs - config = config_from_xml(config_file) - - config.simulation.N = int(config.simulation.N if count == 0.0 else count) - - # Spectra - if monospectrum is not None and powerspectrum is not None: - raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") - if monospectrum is not None: - config.simulation.spectrum = MonoSpectrum(monospectrum) - if powerspectrum is not None: - config.simulation.spectrum = PowerSpectrum(*powerspectrum) - - # Clouds - is_nc = nocloud is not None - is_mc = monocloud is not None - is_pc = pressuremapcloud is not None - any_cloud_flags_present = is_nc | is_mc | is_pc - exactly_one_cloud_flag = (is_nc ^ is_mc ^ is_pc) and not (is_nc & is_mc & is_pc) - if any_cloud_flags_present and not exactly_one_cloud_flag: - raise RuntimeError( - "Only one of --nocloud, --monocloud or --pressuremapcloud may be used." - ) - if is_nc: - config.simulation.cloud_model = NoCloud() - if is_mc: - config.simulation.cloud_model = MonoCloud(monocloud) - if is_pc: - config.simulation.cloud_model = PressureMapCloud(pressuremapcloud.month) - - # Compute - plot = ( - list(registry) - if plotall - else read_plot_config(plotconfig) - if plotconfig - else plot - ) - simulation = compute( - config, - verbose=True, - to_plot=plot, - output_file=output, - write_stages=write_stages, - ) - - if not no_result_file: - simulation.write(output, format="fits", overwrite=True) - - -@cli.command() -@click.option( - "-n", "--numtrajs", type=float, default=100, help="number of trajectories." -) -@click.option( - "--monospectrum", - type=float, - default=None, - help="Mono Energetic Spectrum Log Energy.", -) -@click.option( - "--powerspectrum", - nargs=3, - type=click.Tuple([float, float, float]), - default=None, - help="Power Spectrum index, lower_bound, upper_bound.", -) -@click.argument("filename") -def create_config(filename: str, numtrajs: float, monospectrum, powerspectrum) -> None: - """Generate a configuration file from the given parameters. - - \f - - Parameters - ---------- - filename: str - Name of output xml configuration file. - numtrajs: float, optional - Number of thrown trajectories. Optionally override value in config_file. - - Examples - -------- - Command line usage of the create_config command may work with the following invocation. - - `nuspacesim create_config -n 1e5 sample_input_file.xml` - """ - if monospectrum is not None and powerspectrum is not None: - raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") - - spec = MonoSpectrum() - - if monospectrum is not None: - spec = MonoSpectrum(monospectrum) - if powerspectrum is not None: - spec = PowerSpectrum(*powerspectrum) - # spec = FileSpectrum() - - simulation = SimulationParameters(N=int(numtrajs), spectrum=spec) - - create_xml(filename, NssConfig(simulation=simulation)) - - -@cli.command() -@click.option( - "-p", - "--plot", - type=click.Choice(list(registry), case_sensitive=False), - multiple=True, - default=[], - help="Available plotting functions. Select multiple plots with multiple uses of -p", -) -@click.option( - "-P", - "--plotconfig", - type=click.Path( - exists=True, - dir_okay=False, - readable=True, - ), - help="Read selected plotting functions and options from the specified INI file", -) -@click.option("--plotall", is_flag=True, help="Show all result plots.") -@click.argument( - "simulation_file", - default=None, - type=click.Path(exists=True, dir_okay=False, readable=True), -) -def show_plot( - simulation_file: str, - plot: list, - plotconfig: str, - plotall: bool, -) -> None: - """Show predefined plots of results in simulation file. - - \f - - Parameters - ---------- - simulation_file: str - input nuspacesim results file AstropyTable fits file. - plot: list, optional - Plot the simulation results. - plotconfig: str, optional - INI file to select plots for each module, as well as to specifiy global plot settings. - plotall: bool, optional - Plot all the available the simulation results plots. - - - Examples - -------- - Command line usage of the run command may work with the following invocation. - - `nuspacesim show_plot my_sim_results.fits -p taus_overview` - """ - - sim = AstropyTable.read(simulation_file) - - plot = ( - list(registry) - if plotall - else read_plot_config(plotconfig) - if plotconfig - else plot - ) - - simulation.geometry.region_geometry.show_plot(sim, plot) - simulation.taus.taus.show_plot(sim, plot) - simulation.eas_optical.eas.show_plot(sim, plot) - plots.show_plot(sim, plot) - - -def read_plot_config(filename): - plot_list = [] - cfg = configparser.ConfigParser() - cfg.read(filename) - for sec in cfg.sections()[1:]: - for key in cfg[sec]: - try: - if cfg[sec].getboolean(key): - plot_list.append(key) - except Exception as e: - print(e, "Config file contains non-valid option") - return plot_list -if __name__ == "__main__": - cli() +cli.add_command(run.run) +cli.add_command(create_config.create_config) +cli.add_command(show_plot.show_plot) diff --git a/src/nuspacesim/apps/create_config.py b/src/nuspacesim/apps/create_config.py new file mode 100644 index 0000000..0a8b3cd --- /dev/null +++ b/src/nuspacesim/apps/create_config.py @@ -0,0 +1,85 @@ +import click + +from ..config import NssConfig, create_toml +from .utils import parse_cloud_options, parse_spectra_options + + +@click.command() +@click.option( + "-n", "--numthrown", type=float, default=100, help="number of thrown events." +) +@click.option( + "--monospectrum", + type=float, + default=None, + help="Mono Energetic Spectrum Log Energy.", +) +@click.option( + "--powerspectrum", + nargs=3, + type=click.Tuple([float, float, float]), + default=None, + help="Power Spectrum index, lower_bound, upper_bound.", +) +@click.option( + "--nocloud", + is_flag=True, + default=None, + help="No Cloud Model. [Default]", +) +@click.option( + "--monocloud", + type=float, + default=None, + help="Uniform (mono) Height Cloud Model (km).", +) +@click.option( + "--pressuremapcloud", + type=click.DateTime(["%m", "%B", "%b"]), + default=None, + help="Pressure Map Cloud Model (built in and included with NuSpaceSim). This map is" + "an instance of a global cloud top pressure map sampled from a model of all days in" + "the given month over a 10-year time period from 2011 to 2020. Data provided by the" + "MERRA-2 dataset." + "User should provide a month name, abbreviation, or number.", +) +@click.argument("filename") +def create_config( + filename: str, + numthrown: float, + monospectrum: float, + powerspectrum: click.Tuple, + nocloud: bool, + monocloud: float, + pressuremapcloud: click.DateTime, +) -> None: + """Generate a configuration file from the given parameters. + + \f + + Parameters + ---------- + filename: str + Name of output toml configuration file. + numthrown: float, optional + Number of thrown trajectories. Optionally override value in config_file. + + Examples + -------- + Command line usage of the create_config command may work with the following invocation. + + `nuspacesim create_config -n 1e5 sample_input_file.toml` + """ + config = NssConfig() + + overwrite_spectrum = parse_spectra_options(monospectrum, powerspectrum) + if overwrite_spectrum: + config.simulation.spectrum = overwrite_spectrum + + overwrite_cloud = parse_cloud_options(nocloud, monocloud, pressuremapcloud) + if overwrite_cloud: + config.simulation.cloud_model = overwrite_cloud + + config.simulation.thrown_events = int(numthrown) + + create_toml(filename, config) diff --git a/src/nuspacesim/apps/run.py b/src/nuspacesim/apps/run.py new file mode 100644 index 0000000..09c5191 --- /dev/null +++ b/src/nuspacesim/apps/run.py @@ -0,0 +1,210 @@ +# The Clear BSD License +# +# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted (subject to the limitations in the disclaimer +# below) provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY +# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND +# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +""" Command line client source code. + +.. _cli: + +***** + CLI +***** + +.. currentmodule:: nuspacesim.apps.cli + +.. click:: nuspacesim.apps.cli:cli + :prog: nuspacesim + :show-nested: + +""" + +import click + +from ..compute import compute +from ..config import config_from_toml +from ..results_table import output_filename +from ..utils.plot_function_registry import registry +from .utils import parse_cloud_options, parse_spectra_options, read_plot_config + + +@click.command() +@click.option( + "-o", "--output", type=click.Path(exists=False), default=None, help="Output file." +) +@click.option( + "-p", + "--plot", + type=click.Choice(list(registry), case_sensitive=False), + multiple=True, + default=[], + help="Available plotting functions. Select multiple plots with multiple uses of -p", +) +@click.option( + "-P", + "--plotconfig", + type=click.Path( + exists=True, + dir_okay=False, + readable=True, + ), + help="Read selected plotting functions and options from the specified INI file", +) +@click.option("--plotall", is_flag=True, help="Show all result plots.") +@click.option( + "-w", + "--write-stages", + is_flag=True, + help="Write intermediate values after each simulation stage.", +) +@click.option( + "-n", + "--no-result-file", + is_flag=True, + help="Do not save the results to an output file.", +) +@click.option( + "--monospectrum", + type=float, + default=None, + help="Mono Energetic Spectrum Log Energy.", +) +@click.option( + "--powerspectrum", + nargs=3, + type=click.Tuple([float, float, float]), + default=None, + help="Power Spectrum index, lower_bound, upper_bound.", +) +@click.option( + "--nocloud", + is_flag=True, + default=None, + help="No Cloud Model. [Default]", +) +@click.option( + "--monocloud", + type=float, + default=None, + help="Uniform (mono) Height Cloud Model (km).", +) +@click.option( + "--pressuremapcloud", + type=click.DateTime(["%m", "%B", "%b"]), + default=None, + help="Pressure Map Cloud Model (built in and included with NuSpaceSim). This map is" + "an instance of a global cloud top pressure map sampled from a model of all days in" + "the given month over a 10-year time period from 2011 to 2020. Data provided by the" + "MERRA-2 dataset." + "User should provide a month name, abbreviation, or number.", +) +@click.argument( + "config_file", + default=None, + type=click.Path(exists=True, dir_okay=False, readable=True), +) +@click.argument("count", type=float, default=0.0) +def run( + config_file: str, + count: float, + no_result_file: bool, + output: str, + plot: list, + plotconfig: str, + plotall: bool, + write_stages: bool, + monospectrum: float, + powerspectrum: click.Tuple, + nocloud: bool, + monocloud: float, + pressuremapcloud: click.DateTime, +) -> None: + """Perform the full nuspacesim simulation. + + Main Simulator for nuspacesim. Given a TOML configuration file, and + optionally a count of simulated nutrinos, runs nutrino simulation. + + \f + + Parameters + ---------- + config_file: str + TOML configuration file for particular simulation particular. + count : int, optional + Number of thrown trajectories. Optionally override value in config_file. + spectrum_type : str, optional + Type of + output: str, optional + Name of the output file holding the simulation results. + plot: list, optional + Plot the simulation results. + plotconfig: str, optional + INI file to select plots for each module, as well as to specifiy global plot settings. + plotall: bool, optional + Plot all the available the simulation results plots. + no_result_file: bool, optional + Disables writing results to output files. + + + Examples + -------- + Command line usage of the run command may work with the following invocation. + + `nuspacesim run sample_input_file.toml 1e5 8 -o my_sim_results.fits` + """ + + # User Inputs + config = config_from_toml(config_file) + + config.simulation.thrown_events = int( + config.simulation.thrown_events if count == 0.0 else count + ) + + overwrite_spectrum = parse_spectra_options(monospectrum, powerspectrum) + if overwrite_spectrum: + config.simulation.spectrum = overwrite_spectrum + + overwrite_cloud = parse_cloud_options(nocloud, monocloud, pressuremapcloud) + if overwrite_cloud: + config.simulation.cloud_model = overwrite_cloud + + plot = read_plot_config(registry, plotall, plotconfig, plot) + + output = output_filename(output) + simulation = compute( + config, + verbose=True, + to_plot=plot, + output_file=output, + write_stages=write_stages, + ) + + if not no_result_file: + simulation.write(output, format="fits", overwrite=True) diff --git a/src/nuspacesim/apps/show_plot.py b/src/nuspacesim/apps/show_plot.py new file mode 100644 index 0000000..8ccfd9b --- /dev/null +++ b/src/nuspacesim/apps/show_plot.py @@ -0,0 +1,73 @@ +import click +from astropy.table import Table as AstropyTable + +from .. import simulation +from ..utils import plots +from ..utils.plot_function_registry import registry +from .utils import read_plot_config + +__all__ = ["show_plot"] + + +@click.command() +@click.option( + "-p", + "--plot", + type=click.Choice(list(registry), case_sensitive=False), + multiple=True, + default=[], + help="Available plotting functions. Select multiple plots with multiple uses of -p", +) +@click.option( + "-P", + "--plotconfig", + type=click.Path( + exists=True, + dir_okay=False, + readable=True, + ), + help="Read selected plotting functions and options from the specified INI file", +) +@click.option("--plotall", is_flag=True, help="Show all result plots.") +@click.argument( + "simulation_file", + default=None, + type=click.Path(exists=True, dir_okay=False, readable=True), +) +def show_plot( + simulation_file: str, + plot: list, + plotconfig: str, + plotall: bool, +) -> None: + """Show predefined plots of results in simulation file. + + \f + + Parameters + ---------- + simulation_file: str + input nuspacesim results file AstropyTable fits file. + plot: list, optional + Plot the simulation results. + plotconfig: str, optional + INI file to select plots for each module, as well as to specifiy global plot settings. + plotall: bool, optional + Plot all the available the simulation results plots. + + + Examples + -------- + Command line usage of the run command may work with the following invocation. + + `nuspacesim show_plot my_sim_results.fits -p taus_overview` + """ + + sim = AstropyTable.read(simulation_file) + + plot = read_plot_config(registry, plotall, plotconfig, plot) + + simulation.geometry.region_geometry.show_plot(sim, plot) + simulation.taus.taus.show_plot(sim, plot) + simulation.eas_optical.eas.show_plot(sim, plot) + plots.show_plot(sim, plot) diff --git a/src/nuspacesim/apps/utils.py b/src/nuspacesim/apps/utils.py new file mode 100644 index 0000000..c0ab353 --- /dev/null +++ b/src/nuspacesim/apps/utils.py @@ -0,0 +1,46 @@ +import configparser + +from ..config import Simulation + +__all__ = ["parse_spectra_options", "parse_cloud_options", "read_plot_config"] + + +def parse_spectra_options(monospectrum, powerspectrum): + if monospectrum and powerspectrum: + raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") + if monospectrum: + return Simulation.MonoSpectrum(monospectrum) + if powerspectrum: + return Simulation.PowerSpectrum(*powerspectrum) + + +def parse_cloud_options(nocloud, monocloud, pressuremapcloud): + if sum([1 if x else 0 for x in (nocloud, monocloud, pressuremapcloud)]) > 1: + raise RuntimeError( + "Only one of --nocloud, --monocloud or --pressuremapcloud may be used." + ) + if nocloud: + return Simulation.NoCloud() + if monocloud: + return Simulation.MonoCloud(monocloud) + if pressuremapcloud: + return Simulation.PressureMapCloud(month=pressuremapcloud.month) + + +def read_plot_config(registry, plotall, plotconfig, plot): + if plotall: + return list(registry) + elif plotconfig: + plot_list = [] + cfg = configparser.ConfigParser() + cfg.read(plotconfig) + for sec in cfg.sections()[1:]: + for key in cfg[sec]: + try: + if cfg[sec].getboolean(key): + plot_list.append(key) + except Exception as e: + print(e, "Config file contains non-valid option") + return plot_list + else: + return plot diff --git a/src/nuspacesim/compute.py b/src/nuspacesim/compute.py index a491ee9..c430ede 100644 --- a/src/nuspacesim/compute.py +++ b/src/nuspacesim/compute.py @@ -70,7 +70,7 @@ def compute( config: NssConfig, verbose: bool = False, - output_file: str = None, + output_file: str | None = None, to_plot: list = [], write_stages=False, ) -> AstropyTable: @@ -124,7 +124,10 @@ def compute( console = Console(width=80, log_path=False) - FreqRange = (config.detector.low_freq, config.detector.high_freq) + FreqRange = ( + config.detector.radio.low_frequency, + config.detector.radio.high_frequency, + ) def logv(*args): """optionally print descriptive messages.""" @@ -144,7 +147,6 @@ def mc_logv(mcint, mcintgeo, numEvPass, mcunc, method): logv(f"\t[blue]Stat uncert of MC Integral [/][magenta][{method}][/]:", mcunc) sim = results_table.init(config) - output_file = results_table.output_filename(output_file, now=sim.meta["simTime"][0]) geom = RegionGeom(config) cloud = CloudTopHeight(config) spec = Spectra(config) @@ -176,9 +178,11 @@ def add_meta(self, name: str, value: Any, comment: str): logv(f"Running NuSpaceSim with Energy Spectrum ({config.simulation.spectrum})") logv("Computing [green] Geometries.[/]") - beta_tr, thetaArr, pathLenArr = geom(config.simulation.N, store=sw, plot=to_plot) + beta_tr, thetaArr, pathLenArr = geom( + config.simulation.thrown_events, store=sw, plot=to_plot + ) logv( - f"\t[blue]Threw {config.simulation.N} neutrinos. {beta_tr.size} were valid.[/]" + f"\t[blue]Threw {config.simulation.thrown_events} neutrinos. {beta_tr.size} were valid.[/]" ) init_lat, init_long = geom.find_lat_long_along_traj(np.zeros_like(beta_tr)) sw( @@ -200,7 +204,8 @@ def add_meta(self, name: str, value: Any, comment: str): logv("Computing [green] Decay Altitudes.[/]") altDec, lenDec = eas.altDec(beta_tr, tauBeta, tauLorentz, store=sw) - if config.detector.method == "Optical" or config.detector.method == "Both": + # if config.detector.method == "Optical" or config.detector.method == "Both": + if config.detector.optical: logv("Computing [green] EAS Optical Cherenkov light.[/]") numPEs, costhetaChEff = eas( @@ -219,7 +224,7 @@ def add_meta(self, name: str, value: Any, comment: str): numPEs, costhetaChEff, tauExitProb, - config.detector.photo_electron_threshold, + config.detector.optical.photo_electron_threshold, mc_spec_norm, spec_weights_sum, ) @@ -231,7 +236,7 @@ def add_meta(self, name: str, value: Any, comment: str): mc_logv(mcint, mcintgeo, passEV, mcunc, "Optical") - if config.detector.method == "Radio" or config.detector.method == "Both": + if config.detector.radio: logv("Computing [green] EAS Radio signal.[/]") EFields = eas_radio( @@ -241,17 +246,17 @@ def add_meta(self, name: str, value: Any, comment: str): snrs = calculate_snr( EFields, FreqRange, - config.detector.altitude, - config.detector.det_Nant, - config.detector.det_gain, + config.detector.initial_position.altitude, + config.detector.radio.nantennas, + config.detector.radio.gain, ) logv("Computing [green] Radio Monte Carlo Integral.[/]") mcint, mcintgeo, passEV, mcunc = geom.mcintegral( snrs, - np.cos(config.simulation.theta_ch_max), + np.cos(config.simulation.max_cherenkov_angle), tauExitProb, - config.detector.det_SNR_thres, + config.detector.radio.snr_threshold, mc_spec_norm, spec_weights_sum, ) diff --git a/src/nuspacesim/config.py b/src/nuspacesim/config.py index 962d3f3..9e3038a 100644 --- a/src/nuspacesim/config.py +++ b/src/nuspacesim/config.py @@ -1,6 +1,6 @@ # The Clear BSD License # -# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team +# Copyright (c) 2023 Alexander Reustle and the NuSpaceSim Team # All rights reserved. # # Redistribution and use in source and binary forms, with or without @@ -36,238 +36,281 @@ from __future__ import annotations -from dataclasses import dataclass, field -from typing import Any, Union +from datetime import datetime +from typing import Literal, Optional, Union + +import astropy.units as u +import numpy as np +from astropy.units import Quantity +from pydantic import ( # ValidationError, + BaseModel, + ConfigDict, + Field, + field_serializer, + field_validator, + model_validator, +) try: - from functools import cached_property -except ImportError: - from cached_property import cached_property + import tomllib +except ModuleNotFoundError: + import tomli as tomllib -from numpy import nan, radians, sin - -from . import constants as const -from .types.cloud_types import MonoCloud, NoCloud, PressureMapCloud +import tomli_w __all__ = [ - "DetectorCharacteristics", - "SimulationParameters", + "config_from_toml", + "create_toml", "NssConfig", - "FileSpectrum", - "MonoSpectrum", - "PowerSpectrum", + "Detector", + "Simulation", ] -@dataclass -class DetectorCharacteristics: +def config_from_toml(filename: str) -> NssConfig: + with open(filename, "rb") as f: + c = tomllib.load(f) + return NssConfig(**c) + + +def create_toml(filename: str, c: NssConfig): + with open(filename, "wb") as f: + tomli_w.dump(c.model_dump(), f) + + +def parse_units(value: Union[Quantity, float, str], unit: u.Unit) -> float: + if isinstance(value, (Quantity, str)): + return Quantity(value).to(unit).value + else: + return Quantity(value, unit).value + + +class Detector(BaseModel): r"""Dataclass holding Detector Characteristics.""" - method: str = "Optical" - """ Type of detector: Default = Optical """ - altitude: float = 525.0 - """ Altitude from sea-level. Default = 525 KM """ - lat_start: float = 0.0 - """ Right Ascencion (Degrees). Default = 0.0 """ - long_start: float = 0.0 - """ Declination (Degrees). Default = 0.0 """ - telescope_effective_area: float = 2.5 - '"" Effective area of the detector telescope. Default = 2.5 sq.meters ""' - quantum_efficiency: float = 0.2 - """ Quantum Efficiency of the detector telescope. Default = 0.2 """ - photo_electron_threshold: float = 10.0 - """Photo Electron Threshold, Number Photo electrons = 10""" - low_freq: float = 30.0 - """ Low end for radio band in MHz: Default = 30 """ - high_freq: float = 300.0 - """ High end of radio band in MHz: Default = 300 """ - det_SNR_thres: float = 5.0 - """ SNR threshold for radio triggering: Default = 5 """ - det_Nant: int = 10 - """ Number of radio antennas: Default = 10 """ - det_gain: float = 1.8 - """ Antenna gain in dB: Default = 1.8 """ - - def __call__(self) -> dict[str, tuple[Any, str]]: - r"""Dictionary representation of DetectorCharacteristics instance. - - Groups the data member values with descriptive comments in a tuple. Adds - keyword names no longer than eight characters. This is useful for setting the - FITS Header Keywords in the simulation ouput file. Descriptive comments are - shorter than 80 characters, so as to conform to FITS standards. - - Returns - ------- - dict - Representation of the data members with comments. - """ - - return { - "method": (self.method, "Detector: Method (Optical, radio or both"), - "detAlt": (self.altitude, "Detector: Altitude"), - "latStart": (self.lat_start, "Detector: Initial Latitude"), - "lonStart": (self.long_start, "Detector: Initial Longitude"), - "telEffAr": ( - self.telescope_effective_area, - "Detector: Telescope Effective Area", - ), - "quantEff": (self.quantum_efficiency, "Detector: Quantum Efficiency"), - "phEthres": ( - self.photo_electron_threshold, - "Detector: Photo-Electron Threshold", - ), - "lowFreq": (self.low_freq, "Detector: Antenna Low Frequency"), - "highFreq": (self.high_freq, "Detector: Antenna High Frequency"), - "SNRthres": (self.det_SNR_thres, "Detector: Radio SNR threshold"), - "detNant": (self.det_Nant, "Detector: Number of Antennas"), - "detGain": (self.det_gain, "Detector: Antenna Gain"), - } - - -@dataclass -class MonoSpectrum: - log_nu_tau_energy: float = 8.0 - """Log Energy of the tau neutrinos in GeV.""" - - def __call__(self) -> dict: - return { - "specType": ("Mono", "Simulation: nutau energy spectrum"), - "specPara": (self.log_nu_tau_energy, "Simulation: Log Energy (GeV)"), - } - - -@dataclass -class PowerSpectrum: - index: float = 2.0 - """Power Law Log Energy of the tau neutrinos in GeV.""" - - lower_bound: float = 6.0 - """Lower Bound Log nu_tau Energy GeV.""" - - upper_bound: float = 12.0 - """Upper Bound Log nu_tau Energy GeV.""" - - def __call__(self) -> dict: - return { - "specType": ("Power", "Simulation: nutau Power Law energy spectrum"), - "specPara": (self.index, "Simulation: Power Law Index"), - "specLow": (self.lower_bound, "Simulation: Energy Lower Bound"), - "specHigh": (self.upper_bound, "Simulation: Energy Upper Bound"), - } - - -@dataclass -class FileSpectrum: - path: str = "" - """File path to user defined spectrum file""" - - def __call__(self) -> dict: - return { - "specType": ("File", "Simulation: Nutau User Defined energy spectrum"), - "specFile": (self.path, "Simulation: FilePath"), - } - - -@dataclass -class SimulationParameters: - """Dataclass holding Simulation Parameters.""" - - N: int = 10000 - """Number of thrown trajectories. Default = 1000""" - theta_ch_max: float = radians(3.0) - """Maximum Cherenkov Angle in radians. Default = π/60 radians (3 degrees).""" - spectrum: Union[MonoSpectrum, PowerSpectrum, FileSpectrum] = field( - default_factory=MonoSpectrum + model_config = ConfigDict(arbitrary_types_allowed=True) + + class InitialPos(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + altitude: float = Quantity(525.0, u.km).value + """ Altitude from sea-level (KM). """ + latitude: float = Quantity(0.0, u.rad).value + """ Right Ascencion (Radians). """ + longitude: float = Quantity(0.0, u.rad).value + """ Declination (Radians). """ + + @field_validator("altitude", mode="before") + @classmethod + def valid_distkm(cls, x: Union[Quantity, float, str]) -> float: + return parse_units(x, u.km) + + @field_validator("latitude", "longitude", mode="before") + @classmethod + def valid_anglerad(cls, x: Union[Quantity, float, str]) -> float: + return parse_units(x, u.rad) + + @field_serializer("altitude") + def serialize_km(self, altitude: float) -> str: + return str(Quantity(altitude, u.km)) + + @field_serializer("latitude", "longitude") + def serialize_rad(self, x: float) -> str: + return str(Quantity(x, u.rad).to(u.deg)) + + class Optical(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + telescope_effective_area: float = 2.5 # Quantity(2.5, u.m**2) + """ Effective area of the detector telescope (sq.meters). """ + quantum_efficiency: float = 0.2 + """ Quantum Efficiency of the detector telescope. """ + photo_electron_threshold: float = 10 + """ Photo Electron Threshold, Number Photo electrons. """ + + @field_validator("telescope_effective_area", mode="before") + @classmethod + def valid_aream2(cls, x: Union[Quantity, float, str]) -> float: + return parse_units(x, u.m**2) + + @field_serializer("telescope_effective_area") + def serialize_aream2(self, x: float) -> str: + return str(Quantity(x, u.m**2)) + + class Radio(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + low_frequency: float = Quantity(30.0, u.MHz).value + """ Low end for radio band in MHz: Default = 30 """ + high_frequency: float = Quantity(300.0, u.MHz).value + """ High end of radio band in MHz: Default = 300 """ + snr_threshold: float = 5.0 + """ SNR threshold for radio triggering: Default = 5 """ + nantennas: int = 10 + """ Number of radio antennas: Default = 10 """ + gain: float = Quantity(1.8, u.dB).value + """ Antenna gain in dB: Default = 1.8 """ + + @field_validator("low_frequency", "high_frequency", mode="before") + @classmethod + def valid_freqMHz(cls, x: Union[Quantity, float, str]) -> float: + return parse_units(x, u.MHz) + + @field_validator("gain", mode="before") + @classmethod + def valid_powerdB(cls, x: Union[Quantity, float, str]) -> float: + return parse_units(x, u.dB) + + @field_serializer("low_frequency", "high_frequency") + def serialize_freqMHz(self, x: float) -> str: + return str(Quantity(x, u.MHz)) + + @field_serializer("gain") + def serialize_dB(self, x: float) -> str: + return str(Quantity(x, u.dB)) + + @model_validator(mode="after") + def validate_high_frequency(self): + if self.high_frequency <= self.low_frequency: + raise ValueError("High frequency must be greater than low frequency") + return self + + name: str = "Default Name" + initial_position: InitialPos = InitialPos() + """Initial conditions for detector""" + optical: Optional[Optical] = Optical() + """Characteristics of the optical detector""" + radio: Optional[Radio] = Radio() + """Characteristics of the radio detector""" + + +################################################################### + + +class Simulation(BaseModel): + """Model holding Simulation Parameters.""" + + model_config = ConfigDict(arbitrary_types_allowed=True) + + ################ Radio Ionosphere classes ################ + + class Ionosphere(BaseModel): + total_electron_content: float = 10.0 + """Total Electron Content for ionospheric propagation. """ + total_electron_error: float = 0.1 + """Error for TEC reconstruction. """ + + ################ tau_shower classes ################ + + class NuPyPropShower(BaseModel): + id: Literal["nupyprop"] = "nupyprop" + etau_frac: float = 0.5 + """Fraction of ETau in Shower. Default = 0.5.""" + table_version: str = "3" + """Version of tau conversion tables.""" + + ################ spectrum classes ################ + + class MonoSpectrum(BaseModel): + id: Literal["monospectrum"] = "monospectrum" + log_nu_energy: float = 8.0 + """Log Energy of the tau neutrinos in GeV.""" + + class PowerSpectrum(BaseModel): + id: Literal["powerspectrum"] = "powerspectrum" + index: float = 2.0 + """Power Law Log Energy of the tau neutrinos in GeV.""" + lower_bound: float = 6.0 + """Lower Bound Log nu_tau Energy GeV.""" + upper_bound: float = 12.0 + """Upper Bound Log nu_tau Energy GeV.""" + + ################ Cloud Model classes ################ + + class NoCloud(BaseModel): + id: Literal["no_cloud"] = "no_cloud" + + class MonoCloud(BaseModel): + id: Literal["monocloud"] = "monocloud" + altitude: float = float("-inf") + """Altitude of monoheight cloud.""" + + class PressureMapCloud(BaseModel): + id: Literal["pressure_map"] = "pressure_map" + month: int = 1 + """Cloud Map Month integer 1-12 inclusive.""" + version: Union[str, int] = 0 + """Cloud Map File Version.""" + + @field_validator("month", mode="before") + @classmethod + def valid_month(cls, date: str | int | datetime) -> int: + if isinstance(date, datetime): + return date.month + if isinstance(date, int): + if date < 1 or date > 12: + raise ValueError(f"Provided month {date} is invalid") + return date + if isinstance(date, str): + try: + return (datetime.strptime(date, "%m")).month + except ValueError: + pass + try: + return (datetime.strptime(date, "%B")).month + except ValueError: + pass + try: + return (datetime.strptime(date, "%b")).month + except ValueError: + pass + raise ValueError( + f"date {date} does not match valid month patterns (%m, %B, %b)" + ) + + ################################################################################ + + mode: Literal["Diffuse", "ToO"] = "Diffuse" + """ The Simulation Mode """ + thrown_events: int = 1000 + """ Number of thrown event trajectories. """ + max_cherenkov_angle: float = np.radians(3) + """ Maximum Cherenkov Angle (Radians). """ + max_azimuth_angle: float = np.radians(360) + """ Maximum Azimuthal Angle (Radians). """ + angle_from_limb: float = np.radians(7) + """ Angle From Limb. Default (Radians). """ + cherenkov_light_engine: Literal["Default"] = "Default" # "CHASM", "EASCherSim" + ionosphere: Optional[Ionosphere] = Ionosphere() + tau_shower: NuPyPropShower = NuPyPropShower() + """ Tau Shower Generator. """ + spectrum: Union[MonoSpectrum, PowerSpectrum] = Field( + default=MonoSpectrum(), discriminator="id" ) - """Distribution from which to draw nu_tau energies.""" - cloud_model: Union[NoCloud, MonoCloud, PressureMapCloud] = field( - default_factory=NoCloud + """ Distribution from which to draw nu_tau energies. """ + cloud_model: Union[NoCloud, MonoCloud, PressureMapCloud] = Field( + default=NoCloud(), discriminator="id" ) - e_shower_frac: float = 0.5 - """Fraction of ETau in Shower. Default = 0.5.""" - ang_from_limb: float = radians(7.0) - """Angle From Limb. Default = π/25.714 radians (7 degrees).""" - max_azimuth_angle: float = radians(360.0) - """Maximum Azimuthal Angle. Default = 2π radians (360 degrees).""" - model_ionosphere: int = 0 - """Model ionosphere for radio propagation?. Default = 0 (false).""" - TEC: float = 10.0 - """Total Electron Content for ionospheric propagation. Default = 10.""" - TECerr: float = 0.1 - """Error for TEC reconstruction. Default = 0.1""" - tau_table_version: str = "3" - """Version of tau conversion tables.""" - - @cached_property - def log_nu_tau_energy(self) -> float: - """log10 of nu_tau_energy.""" - if isinstance(self.spectrum, MonoSpectrum): - return self.spectrum.log_nu_tau_energy - else: - return nan - - @cached_property - def nu_tau_energy(self) -> float: - """10 ^ log_nu_tau_energy.""" - if isinstance(self.spectrum, MonoSpectrum): - return 10**self.log_nu_tau_energy - else: - return nan - - @cached_property - def sin_theta_ch_max(self) -> float: - """sin of theta_ch_max.""" - return sin(self.theta_ch_max) - - def __call__(self) -> dict[str, tuple[Union[int, float, str], str]]: - r"""Dictionary representation of SimulationParameters instance. - - Groups the data member values with descriptive comments in a tuple. Adds - keyword names no longer than eight characters. This is useful for setting the - FITS Header Keywords in the simulation ouput file. Descriptive comments are - shorter than 80 characters, so as to conform to FITS standards. - - Returns - ------- - dict - Representation of the data members with comments. - """ - d = { - "N": (self.N, "Simulation: thrown neutrinos"), - "thChMax": (self.theta_ch_max, "Simulation: Maximum Cherenkov Angle"), - "specUnit": ("log(GeV)", "Simulation: Energy spectrum units"), - "eShwFrac": (self.e_shower_frac, "Simulation: Fraction of Etau in Shower"), - "angLimb": (self.ang_from_limb, "Simulation: Angle From Limb"), - "maxAzAng": (self.max_azimuth_angle, "Simulation: Maximum Azimuthal Angle"), - "ionosph": ( - self.model_ionosphere, - "Simulation: Radio ionosphere model flg", - ), - "TEC": (self.TEC, "Simulation: Actual slant TEC value"), - "TECerr": (self.TECerr, "Simulation: Uniform distr. err: TEC est."), - } - - d.update(self.spectrum()) - d.update(self.cloud_model()) - - return d - - -@dataclass -class NssConfig: + + @field_validator( + "max_cherenkov_angle", "max_azimuth_angle", "angle_from_limb", mode="before" + ) + @classmethod + def valid_anglerad(cls, x: Union[Quantity, float, str]) -> float: + return parse_units(x, u.rad) + + @field_serializer("max_cherenkov_angle", "max_azimuth_angle", "angle_from_limb") + def serialize_rad(self, x: float) -> str: + return str(Quantity(x, u.rad).to(u.deg)) + + +class NssConfig(BaseModel): r"""Necessary Configuration Data for NuSpaceSim. An :class:`NssConfig` is a container object holding all of the other nuSpaceSim configuration objects for a simplified access API. Instances of :class:`NssConfig` - objects can be serialized to XML. - - .. :ref:`XML`. - + objects can be serialized to TOML. """ - detector: DetectorCharacteristics = field(default_factory=DetectorCharacteristics) + title: str = "NuSpaceSim" + detector: Detector = Detector() """The Detector Characteristics.""" - simulation: SimulationParameters = field(default_factory=SimulationParameters) + simulation: Simulation = Simulation() """The Simulation Parameters.""" - constants: const.Fund_Constants = field(default_factory=const.Fund_Constants) - """The Fudimental physical constants.""" diff --git a/src/nuspacesim/results_table.py b/src/nuspacesim/results_table.py index 846e3ab..526a94d 100644 --- a/src/nuspacesim/results_table.py +++ b/src/nuspacesim/results_table.py @@ -46,6 +46,7 @@ from astropy.table import Table as AstropyTable from .config import NssConfig +from .utils.misc import flatten_dict __all__ = ["init", "output_filename"] @@ -59,10 +60,8 @@ def init(config: NssConfig | None = None): now = f"{datetime.datetime.now():%Y%m%d%H%M%S}" return AstropyTable( meta={ - **config.detector(), - **config.simulation(), - **config.constants(), "simTime": (now, "Start time of Simulation"), + **flatten_dict(config.model_dump(), "HIERARCH Config", sep=" "), } ) elif isinstance(config, AstropyTable): @@ -75,5 +74,5 @@ def init(config: NssConfig | None = None): def output_filename(filename: str | None, now: str | None = None) -> str: if filename is not None: return filename - now = now if not now else f"{datetime.datetime.now():%Y%m%d%H%M%S}" + now = now if now else f"{datetime.datetime.now():%Y%m%d%H%M%S}" return f"nuspacesim_run_{now}.fits" diff --git a/src/nuspacesim/simulation/atmosphere/clouds.py b/src/nuspacesim/simulation/atmosphere/clouds.py index 5d238ba..c3de54a 100644 --- a/src/nuspacesim/simulation/atmosphere/clouds.py +++ b/src/nuspacesim/simulation/atmosphere/clouds.py @@ -45,8 +45,9 @@ import numpy.typing as npt from astropy.io import fits -from ...config import NssConfig -from ...types import cloud_types +from ...config import NssConfig, Simulation + +# from ...types import cloud_types from ..eas_optical import atmospheric_models as atm @@ -57,26 +58,26 @@ def __init__(self, config: NssConfig): if self.cloud_model is None: self.altitude = mono_altitude() - elif isinstance(self.cloud_model, cloud_types.NoCloud): + elif isinstance(self.cloud_model, Simulation.NoCloud): self.altitude = mono_altitude() - elif isinstance(self.cloud_model, cloud_types.MonoCloud): + elif isinstance(self.cloud_model, Simulation.MonoCloud): self.altitude = mono_altitude(self.cloud_model.altitude) - elif isinstance(self.cloud_model, cloud_types.PressureMapCloud): + elif isinstance(self.cloud_model, Simulation.PressureMapCloud): self.map = extract_fits_cloud_pressure_map_v0(self.cloud_model) self.altitude = altitude_from_pressure_map_v0(self.map) else: RuntimeError(f"Unrecognized Cloud Model Type: {type(self.cloud_model)}!") - def __call__(self, *args, **kwargs) -> npt.single: + def __call__(self, *args, **kwargs) -> np.single: return self.altitude(*args, **kwargs) def mono_altitude(altitude=0.0): def f(*args, **kwargs) -> np.single: - return altitude + return np.single(altitude) return f @@ -94,7 +95,7 @@ def f(lat: float, long: float, *args, **kwargs) -> np.single: return f -def extract_fits_cloud_pressure_map_v0(cloud_model: cloud_types.PressureMapCloud): +def extract_fits_cloud_pressure_map_v0(cloud_model: Simulation.PressureMapCloud): month = cloud_model.month version = cloud_model.version with as_file( diff --git a/src/nuspacesim/simulation/eas_optical/eas.py b/src/nuspacesim/simulation/eas_optical/eas.py index c44ed4b..4640fd5 100644 --- a/src/nuspacesim/simulation/eas_optical/eas.py +++ b/src/nuspacesim/simulation/eas_optical/eas.py @@ -32,9 +32,12 @@ # POSSIBILITY OF SUCH DAMAGE. import numpy as np +from astropy import units +from astropy.constants import R_earth, c from ...config import NssConfig from ...utils import decorators +from ..taus.taus import mean_Tau_life from .cphotang import CphotAng from .local_plots import eas_optical_density, eas_optical_histogram @@ -50,7 +53,7 @@ class EAS: def __init__(self, config: NssConfig): self.config = config - self.CphotAng = CphotAng(self.config.detector.altitude) + self.CphotAng = CphotAng(self.config.detector.initial_position.altitude) @decorators.nss_result_store("altDec", "lenDec") def altDec(self, beta, tauBeta, tauLorentz, u=None, *args, **kwargs): @@ -60,17 +63,17 @@ def altDec(self, beta, tauBeta, tauLorentz, u=None, *args, **kwargs): u = np.random.uniform(0, 1, len(beta)) if u is None else u - tDec = (-1.0 * tauLorentz / self.config.constants.inv_mean_Tau_life) * np.log(u) + tDec = -tauLorentz * mean_Tau_life * np.log(u) - lenDec = tDec * tauBeta * self.config.constants.c + lenDec = tDec * tauBeta * c.value altDec = np.sqrt( - self.config.constants.earth_radius**2 + R_earth.to(units.km).value ** 2 + lenDec**2 - + 2.0 * self.config.constants.earth_radius * lenDec * np.sin(beta) + + 2.0 * R_earth.to(units.km).value * lenDec * np.sin(beta) ) - altDec -= self.config.constants.earth_radius + altDec -= R_earth.to(units.km).value return altDec, lenDec @@ -112,12 +115,12 @@ def __call__( numPEs = ( dphots - ###### * showerEnergy # Scaling no longer necessary. Being done in cphotang. - * self.config.detector.telescope_effective_area - * self.config.detector.quantum_efficiency + # * showerEnergy # Scaling no longer necessary. Being done in cphotang. + * self.config.detector.optical.telescope_effective_area + * self.config.detector.optical.quantum_efficiency ) - enhanceFactor = numPEs / self.config.detector.photo_electron_threshold + enhanceFactor = numPEs / self.config.detector.optical.photo_electron_threshold logenhanceFactor = np.empty_like(enhanceFactor) efMask = enhanceFactor > 2.0 logenhanceFactor[efMask] = np.log(enhanceFactor[efMask]) diff --git a/src/nuspacesim/simulation/eas_optical/local_plots.py b/src/nuspacesim/simulation/eas_optical/local_plots.py index b0905d3..8aa69ca 100644 --- a/src/nuspacesim/simulation/eas_optical/local_plots.py +++ b/src/nuspacesim/simulation/eas_optical/local_plots.py @@ -40,7 +40,7 @@ def eas_optical_density(inputs, results, *args, **kwargs): r"""Plot some density plots""" - _, betas, altDec, showerEnergy = inputs + _, betas, altDec, showerEnergy, *_ = inputs numPEs, costhetaChEff = results fig, ax = plt.subplots(2, 3, figsize=(15, 8), constrained_layout=True) diff --git a/src/nuspacesim/simulation/eas_radio/radio.py b/src/nuspacesim/simulation/eas_radio/radio.py index 65ebf0d..4c15153 100644 --- a/src/nuspacesim/simulation/eas_radio/radio.py +++ b/src/nuspacesim/simulation/eas_radio/radio.py @@ -1,5 +1,7 @@ import astropy.io.misc.hdf5 as hf import numpy as np +from astropy.constants import R_earth +from astropy.units import km from ...config import NssConfig from ...utils import decorators @@ -41,7 +43,10 @@ def __call__( """ EAS radio output from ZHAires lookup tables """ - FreqRange = (self.config.detector.low_freq, self.config.detector.high_freq) + FreqRange = ( + self.config.detector.radio.low_frequency, + self.config.detector.radio.high_frequency, + ) radioParams = RadioEFieldParams(FreqRange) mask = (altDec < 0.0) | ( altDec > 10.0 @@ -55,11 +60,11 @@ def __call__( nssDist = distance_to_detector( beta[mask], altDec[mask], - self.config.detector.altitude, - self.config.constants.earth_radius, + self.config.detector.initial_position.altitude, + R_earth.to(km).value, ) zhairesDist = distance_to_detector( - beta[mask], altDec[mask], 525.0, self.config.constants.earth_radius + beta[mask], altDec[mask], 525.0, R_earth.to(km).value ) EFields = np.zeros_like(beta) @@ -74,17 +79,17 @@ def __call__( distScale = np.abs(zhairesDist / nssDist) EFields[mask] = (EFields[mask].T * distScale).T # no ionosphere if the detector is below 90km - if self.config.detector.altitude > 90.0: - if self.config.simulation.model_ionosphere: - if self.config.simulation.TEC < 0: + if self.config.detector.initial_position.altitude > 90.0: + if self.config.simulation.ionosphere: + if self.config.simulation.ionosphere.total_electron_content < 0: print( "TEC should be positive!! continuing without ionospheric dispersion" ) else: ionosphere = IonosphereParams( FreqRange, - self.config.simulation.TECerr, - self.config.simulation.TEC, + self.config.simulation.ionosphere.total_electron_error, + self.config.simulation.ionosphere.total_electron_content, ) ionosphereScaling = ionosphere(EFields[mask]) EFields[mask] *= ionosphereScaling @@ -95,7 +100,7 @@ def __call__( # geomagnetic is only full strength when perpendicular to Earth B-field # here i apply a scaling for vxB for an orbit close to equatorial - Re = self.config.constants.earth_radius + Re = R_earth.to(km).value B_angle = np.ones(altDec[mask].shape) B_angle *= np.pi / 2.0 - np.arccos( (lenDec[mask] ** 2.0 + (altDec[mask] + Re) ** 2.0 - Re**2.0) diff --git a/src/nuspacesim/simulation/geometry/region_geometry.py b/src/nuspacesim/simulation/geometry/region_geometry.py index 61943d7..3b5cae9 100644 --- a/src/nuspacesim/simulation/geometry/region_geometry.py +++ b/src/nuspacesim/simulation/geometry/region_geometry.py @@ -32,6 +32,8 @@ # POSSIBILITY OF SUCH DAMAGE. import numpy as np +from astropy import units as u +from astropy.constants import R_earth from ...utils import decorators from .local_plots import geom_beta_tr_hist @@ -44,20 +46,18 @@ class RegionGeom: def __init__(self, config): self.config = config - - self.earth_rad_2 = self.config.constants.earth_radius**2 + self.earth_radius: np.float64 = R_earth.to(u.km).value + self.earth_rad_2: np.float64 = self.earth_radius**2 self.core_alt = ( - self.config.constants.earth_radius + self.config.detector.altitude + self.earth_radius + self.config.detector.initial_position.altitude ) - self.detLat = np.radians(config.detector.lat_start) - self.detLong = np.radians(config.detector.long_start) + self.detLat = config.detector.initial_position.latitude + self.detLong = config.detector.initial_position.longitude - alphaHorizon = np.pi / 2 - np.arccos( - self.config.constants.earth_radius / self.core_alt - ) - alphaMin = alphaHorizon - config.simulation.ang_from_limb + alphaHorizon = np.pi / 2 - np.arccos(self.earth_radius / self.core_alt) + alphaMin = alphaHorizon - config.simulation.angle_from_limb minChordLen = 2 * np.sqrt( self.earth_rad_2 - (self.core_alt * np.sin(alphaMin)) ** 2 ) @@ -66,10 +66,10 @@ def __init__(self, config): self.minLOSpathLen = self.core_alt * np.cos(alphaMin) - minChordLen / 2 self.maxLOSpathLen = np.sqrt(self.core_alt**2 - self.earth_rad_2) - self.sinOfMaxThetaTrSubV = np.sin(config.simulation.theta_ch_max) + self.sinOfMaxThetaTrSubV = np.sin(config.simulation.max_cherenkov_angle) - self.maxPhiS = config.simulation.max_azimuth_angle / 2 - self.minPhiS = -config.simulation.max_azimuth_angle / 2 + self.maxPhiS = config.simulation.max_azimuth_angle * 0.5 + self.minPhiS = config.simulation.max_azimuth_angle * -0.5 normThetaTrSubV = 2 / self.sinOfMaxThetaTrSubV**2 normPhiTrSubV = np.reciprocal(2 * np.pi) @@ -154,12 +154,12 @@ def throw(self, u=None): rvsqrd = self.losPathLen * self.losPathLen costhetaS = (self.core_alt**2 + self.earth_rad_2 - rvsqrd) / ( - 2 * self.config.constants.earth_radius * self.core_alt + 2 * self.earth_radius * self.core_alt ) self.thetaS = np.arccos(costhetaS) self.costhetaNSubV = (self.core_alt**2 - self.earth_rad_2 - rvsqrd) / ( - 2 * self.config.constants.earth_radius * self.losPathLen + 2 * self.earth_radius * self.losPathLen ) thetaNSubV = np.arccos(self.costhetaNSubV) @@ -289,11 +289,11 @@ def find_lat_long_along_traj(self, dist_along_traj): yPath_v = dist_along_traj * np.sin(self.thetas()) * np.sin( self.phis() - ) + self.config.constants.earth_radius * np.cos(self.valid_elevAngVSubN()) + ) + self.earth_radius * np.cos(self.valid_elevAngVSubN()) - zPath_v = dist_along_traj * np.cos( - self.thetas() - ) + self.config.constants.earth_radius * np.sin(self.valid_elevAngVSubN()) + zPath_v = dist_along_traj * np.cos(self.thetas()) + self.earth_radius * np.sin( + self.valid_elevAngVSubN() + ) # Compute xyz-coordinates in the spot's ENU frame (accomplished by transforming from # the los ENU frame to the spot's ENU frame) diff --git a/src/nuspacesim/simulation/spectra/spectra.py b/src/nuspacesim/simulation/spectra/spectra.py index 38aa61a..91aebd2 100644 --- a/src/nuspacesim/simulation/spectra/spectra.py +++ b/src/nuspacesim/simulation/spectra/spectra.py @@ -36,7 +36,7 @@ import numpy as np from numpy.typing import NDArray -from ...config import FileSpectrum, MonoSpectrum, NssConfig, PowerSpectrum +from ...config import NssConfig, Simulation from ...utils import decorators from .local_plots import spectra_histogram @@ -45,16 +45,16 @@ @decorators.nss_result_store("log_e_nu") def energy_spectra( N: int, - spectra: Union[MonoSpectrum, PowerSpectrum, FileSpectrum, Callable], + spectra: Union[Simulation.MonoSpectrum, Simulation.PowerSpectrum, Callable], *args, **kwargs, ) -> NDArray[Any]: """Energy Spectra of thrown Neutrinos""" - if isinstance(spectra, MonoSpectrum): - return np.full(shape=(N), fill_value=spectra.log_nu_tau_energy) + if isinstance(spectra, Simulation.MonoSpectrum): + return np.full(shape=(N), fill_value=spectra.log_nu_energy) - if isinstance(spectra, PowerSpectrum): + if isinstance(spectra, Simulation.PowerSpectrum): p = spectra.index a = 10**spectra.lower_bound b = 10**spectra.upper_bound @@ -71,14 +71,14 @@ def energy_spectra( def spec_norm( - spectra: Union[MonoSpectrum, PowerSpectrum, FileSpectrum, Callable], + spectra: Union[Simulation.MonoSpectrum, Simulation.PowerSpectrum, Callable], *args, **kwargs, ) -> float: - if isinstance(spectra, MonoSpectrum): + if isinstance(spectra, Simulation.MonoSpectrum): return 1.0 - if isinstance(spectra, PowerSpectrum): + if isinstance(spectra, Simulation.PowerSpectrum): p = spectra.index a = 10**spectra.lower_bound b = 10**spectra.upper_bound @@ -89,14 +89,14 @@ def spec_norm( def sum_spec_weights( - spectra: Union[MonoSpectrum, PowerSpectrum, FileSpectrum, Callable], + spectra: Union[Simulation.MonoSpectrum, Simulation.PowerSpectrum, Callable], *args, **kwargs, ) -> float: - if isinstance(spectra, MonoSpectrum): + if isinstance(spectra, Simulation.MonoSpectrum): return 1.0 - if isinstance(spectra, PowerSpectrum): + if isinstance(spectra, Simulation.PowerSpectrum): p = spectra.index a = 10**spectra.lower_bound b = 10**spectra.upper_bound diff --git a/src/nuspacesim/simulation/taus/taus.py b/src/nuspacesim/simulation/taus/taus.py index 5a1ba1b..2a555f2 100644 --- a/src/nuspacesim/simulation/taus/taus.py +++ b/src/nuspacesim/simulation/taus/taus.py @@ -48,7 +48,10 @@ from importlib_resources import as_file, files -__all__ = ["Taus", "show_plot"] +__all__ = ["Taus", "show_plot", "massTau", "mean_Tau_life"] + +massTau = 1.77686 # GeV/c^2 +mean_Tau_life = 2.903e-13 # seconds class Taus(object): @@ -82,14 +85,14 @@ def __init__(self, config: NssConfig): # grid of pexit table with as_file( files("nuspacesim.data.nupyprop_tables") - / f"nu2tau_pexit.{config.simulation.tau_table_version}.h5" + / f"nu2tau_pexit.{config.simulation.tau_shower.table_version}.h5" ) as file: self.pexit_grid = NssGrid.read(file, path="/", format="hdf5") # grid of tau_cdf tables with as_file( files("nuspacesim.data.nupyprop_tables") - / f"nu2tau_cdf.{config.simulation.tau_table_version}.h5" + / f"nu2tau_cdf.{config.simulation.tau_shower.table_version}.h5" ) as file: self.tau_cdf_grid = NssGrid.read(file, format="hdf5") @@ -172,9 +175,9 @@ def __call__(self, betas, log_e_nu, *args, **kwargs): tauEnergy = self.tau_energy(betas, log_e_nu) # in units of 100 PeV - showerEnergy = self.config.simulation.e_shower_frac * tauEnergy / 1e8 + showerEnergy = self.config.simulation.tau_shower.etau_frac * tauEnergy / 1e8 - tauLorentz = tauEnergy / self.config.constants.massTau + tauLorentz = tauEnergy / massTau tauBeta = np.sqrt(1.0 - np.reciprocal(tauLorentz**2)) return tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb diff --git a/src/nuspacesim/state.py b/src/nuspacesim/state.py new file mode 100644 index 0000000..61b6628 --- /dev/null +++ b/src/nuspacesim/state.py @@ -0,0 +1,40 @@ +from dataclasses import dataclass +from datetime import datetime +from typing import Any, Dict, Sequence + +from astropy.io import fits +from astropy.table import Table as AstropyTable + +from .config import NssConfig + + +@dataclass +class SimulationState: + config: NssConfig = NssConfig() + event_table: AstropyTable = AstropyTable() + results: Dict[str, float] = {} + plots: Sequence[Any] = [] + sim_time: str = f"{datetime.now():%Y%m%d%H%M%S}" + output_file = f"nuspacesim_run_{sim_time}.fits" + + +class StagedWriter: + """Optionally write intermediate values to file""" + + def __init__(self, sim: SimulationState): + self.sim: SimulationState = sim + + def __call__( + self, + col_names: Sequence[str], + columns, # : Sequence[ArrayLike], + *args, + **kwargs, + ): + self.sim.event_table.add_columns(columns, names=col_names, *args, **kwargs) + fits.write(self.sim.event_table, self.sim.output_file, overwrite=True) + + # def add_meta(self, name: str, value: Any, comment: str): + # sim.meta[name] = (value, comment) + # if write_stages: + # sim.write(output_file, format="fits", overwrite=True) diff --git a/src/nuspacesim/utils/misc.py b/src/nuspacesim/utils/misc.py index f9121c5..71c1db1 100644 --- a/src/nuspacesim/utils/misc.py +++ b/src/nuspacesim/utils/misc.py @@ -1,3 +1,5 @@ +from collections.abc import MutableMapping + import numpy as np @@ -7,3 +9,20 @@ def cartesian_product(*arrays): for i, a in enumerate(np.ix_(*arrays)): prod[..., i] = a return prod.reshape(-1, sz) + + +def _flat(d, parent_key, sep): + for k, v in d.items(): + new_key = parent_key + sep + k if parent_key else k + if isinstance(v, MutableMapping): + yield from flatten_dict(v, new_key, sep=sep).items() + else: + yield new_key, v + + +def flatten_dict(d: MutableMapping, parent_key: str = "", sep: str = "."): + return dict(_flat(d, parent_key, sep)) + + +def output_filename(state, filename) -> str: + return filename if filename else f"nuspacesim_run_{state.sim_time}.fits" diff --git a/src/nuspacesim/xml_config/__init__.py b/src/nuspacesim/xml_config/__init__.py deleted file mode 100644 index 9926526..0000000 --- a/src/nuspacesim/xml_config/__init__.py +++ /dev/null @@ -1,79 +0,0 @@ -# The Clear BSD License -# -# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team -# All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted (subject to the limitations in the disclaimer -# below) provided that the following conditions are met: -# -# * Redistributions of source code must retain the above copyright notice, -# this list of conditions and the following disclaimer. -# -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# -# * Neither the name of the copyright holder nor the names of its -# contributors may be used to endorse or promote products derived from this -# software without specific prior written permission. -# -# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY -# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND -# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR -# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR -# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -# POSSIBILITY OF SUCH DAMAGE. - -r"""XML file utilities for configuration objects - --------------------------------------- -NuSpaceSim Configuration XML Interface --------------------------------------- - -Configuration object can be serialized to an XML file, in whole or in part. These -configuration files can be read back into nuspacesim with -:func:`config_from_xml`. All XML config files -are validated with an XSD Schema to by -is_valid_xml ensure correctness. -:func:`is_valid_xml` ensure correctness. - -.. autosummary:: - :toctree: - :nosignatures: - - create_xml - config_from_xml - is_valid_xml - parseXML - parse_detector_chars - parse_simulation_params - config_xml_schema - -""" - -__all__ = [ - "is_valid_xml", - "parse_config", - "parse_detector_chars", - "parse_simulation_params", - "parseXML", - "config_from_xml", - "create_xml", -] - -from . import parse_config -from .parse_config import ( - config_from_xml, - create_xml, - is_valid_xml, - parse_detector_chars, - parse_simulation_params, - parseXML, -) diff --git a/src/nuspacesim/xml_config/config_xml_schema.py b/src/nuspacesim/xml_config/config_xml_schema.py deleted file mode 100644 index e612cf6..0000000 --- a/src/nuspacesim/xml_config/config_xml_schema.py +++ /dev/null @@ -1,213 +0,0 @@ -# The Clear BSD License -# -# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team -# All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted (subject to the limitations in the disclaimer -# below) provided that the following conditions are met: -# -# * Redistributions of source code must retain the above copyright notice, -# this list of conditions and the following disclaimer. -# -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# -# * Neither the name of the copyright holder nor the names of its -# contributors may be used to endorse or promote products derived from this -# software without specific prior written permission. -# -# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY -# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND -# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR -# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR -# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -# POSSIBILITY OF SUCH DAMAGE. -"""XMLSchema string file. For validating XML configuration files.""" - -from io import StringIO - -xsd = StringIO( - """\ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - """ -) -"""XMLSchema string file. For validating XML configuration files.""" diff --git a/src/nuspacesim/xml_config/parse_config.py b/src/nuspacesim/xml_config/parse_config.py deleted file mode 100644 index 2c14b38..0000000 --- a/src/nuspacesim/xml_config/parse_config.py +++ /dev/null @@ -1,423 +0,0 @@ -# The Clear BSD License -# -# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team -# All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted (subject to the limitations in the disclaimer -# below) provided that the following conditions are met: -# -# * Redistributions of source code must retain the above copyright notice, -# this list of conditions and the following disclaimer. -# -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# -# * Neither the name of the copyright holder nor the names of its -# contributors may be used to endorse or promote products derived from this -# software without specific prior written permission. -# -# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY -# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND -# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR -# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR -# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -# POSSIBILITY OF SUCH DAMAGE. -""" -Module contains functions for parsing and interacting with XML configuration files. -""" - -import lxml.etree as ET -import numpy as np - -from .. import constants as const -from ..config import ( - DetectorCharacteristics, - FileSpectrum, - MonoSpectrum, - NssConfig, - PowerSpectrum, - SimulationParameters, -) - -# from ..types import cloud_types -from ..types.cloud_types import MonoCloud, NoCloud, PressureMapCloud, parse_month -from . import config_xml_schema - -__all__ = [ - "is_valid_xml", - "parse_detector_chars", - "parse_simulation_params", - "parseXML", - "config_from_xml", - "create_xml", -] - - -def is_valid_xml(xmlfile: str) -> bool: - r"""Check that the given xml file is valid. - - Confirm that the given xml document is valid by validating it with the XMLSchema. - - Parameters - ---------- - xmlfile: str - The input configuration xml file. - - Returns - ------- - bool - Whether the file is valid. - """ - - xmlschema_doc = ET.parse(config_xml_schema.xsd) - xmlschema = ET.XMLSchema(xmlschema_doc) - return xmlschema.validate(ET.parse(xmlfile)) - - -def parse_detector_chars(xmlfile: str) -> DetectorCharacteristics: - r"""Parse the XML file into a DetectorCharacteristics object. - - Parameters - ---------- - xmlfile: str - The input configuration xml file. - - Returns - ------- - DetectorCharacteristics - The detector characteristics object. - """ - - detchar: dict[str, str] = {} - tree = ET.parse(xmlfile) - root = tree.getroot() - eldetchar = root.find("DetectorCharacteristics") - detchar["Method"] = eldetchar.attrib["Method"] - for node in tree.find("./DetectorCharacteristics"): - if node.tag == "PhotoElectronThreshold": - detchar[node.tag] = str(node.attrib["Preset"]) - if node.attrib["Preset"] == "true": - detchar["NPE"] = str(node.find("NPE").text) - else: - detchar[node.tag] = str(node.text) - - # Convert Degrees to Radians - if "Unit" in node.attrib: - if node.tag in [ - "InitialDetectorLatitude", - "InitialDetectorLongitude", - ]: - x = float(node.text) - detchar[node.tag] = ( - x if node.attrib["Unit"] == "Degrees" else np.radians(x) - ) - - np.degrees(float(node.text)) - elif node.attrib["Unit"] == "Degrees": - detchar[node.tag] = np.radians(float(node.text)) - - return DetectorCharacteristics( - method=detchar["Method"], - altitude=float(detchar["DetectorAltitude"]), - lat_start=float(detchar["InitialDetectorLatitude"]), - long_start=float(detchar["InitialDetectorLongitude"]), - telescope_effective_area=float(detchar["TelescopeEffectiveArea"]), - quantum_efficiency=float(detchar["QuantumEfficiency"]), - photo_electron_threshold=float(detchar["NPE"]), - low_freq=float(detchar["LowFrequency"]), - high_freq=float(detchar["HighFrequency"]), - det_SNR_thres=float(detchar["SNRThreshold"]), - det_Nant=int(detchar["NAntennas"]), - det_gain=float(detchar["AntennaGain"]), - ) - - -def parse_simulation_params(xmlfile: str) -> SimulationParameters: - r"""Parse the XML file into a SimulationParameters object. - - Parameters - ---------- - xmlfile: str - The input configuration xml file. - - Returns - ------- - SimulationParameters - The Simulation Parameters object. - """ - - simparams = {} - tree = ET.parse(xmlfile) - root = tree.getroot() - elsimparams = root.find("SimulationParameters") - simparams[elsimparams.tag] = elsimparams.attrib["DetectionMode"] - - for node in tree.find("./SimulationParameters"): - if node.tag == "TauShowerType": - simparams[node.tag] = node.attrib["Preset"] - if node.attrib["Preset"] == "true": - simparams["FracETauInShower"] = str(node.find("FracETauInShower").text) - - elif node.tag == "NuTauEnergySpecType": - for spectrum_type in node: - if "MonoSpectrum" == spectrum_type.tag: - simparams["Spectrum"] = MonoSpectrum( - log_nu_tau_energy=float(spectrum_type.find("LogNuEnergy").text), - ) - if "PowerSpectrum" == spectrum_type.tag: - simparams["Spectrum"] = PowerSpectrum( - index=float(spectrum_type.find("PowerLawIndex").text), - lower_bound=float(spectrum_type.find("LowerBound").text), - upper_bound=float(spectrum_type.find("UpperBound").text), - ) - if "FileSpectrum" == spectrum_type.tag: - simparams["Spectrum"] = FileSpectrum( - path=str(node.spectrum_type("FilePath").text) - ) - elif node.tag == "CloudModelType": - for cloud_model in node: - if "NoCloudModel" == cloud_model.tag: - simparams["CloudModel"] = NoCloud() - if "MonoCloudModel" == cloud_model.tag: - simparams["CloudModel"] = MonoCloud( - altitude=float(cloud_model.find("CloudTopHeight").text), - ) - if "PressureMapCloudModel" == cloud_model.tag: - simparams["CloudModel"] = PressureMapCloud( - month=parse_month(cloud_model.find("Month").text), - version=str(cloud_model.find("Version").text), - ) - else: - simparams[node.tag] = str(node.text) - - # Convert Degrees to Radians - if "Unit" in node.attrib: - if node.attrib["Unit"] == "Degrees": - simparams[node.tag] = np.radians(float(node.text)) - - return SimulationParameters( - N=int(simparams["NumTrajs"]), - theta_ch_max=float(simparams["MaximumCherenkovAngle"]), - spectrum=simparams["Spectrum"], - cloud_model=simparams["CloudModel"], - e_shower_frac=float(simparams["FracETauInShower"]), - ang_from_limb=float(simparams["AngleFromLimb"]), - max_azimuth_angle=float(simparams["AzimuthalAngle"]), - model_ionosphere=bool(int(simparams["ModelIonosphere"])), - TEC=float(simparams["TEC"]), - TECerr=np.abs(float(simparams["TECerr"])), - tau_table_version=simparams["TauTableVer"], - ) - - -def parseXML(xmlfile: str) -> tuple: - r"""Parse the XML file into a pair of configuration objects. - - If the xml file is valid, parse the file into a DetectorCharacteristics and - SimulationParameters tuple. - - Parameters - ---------- - xmlfile: str - The input configuration xml file. - - Returns - ------- - tuple: - Tuple of [DetectorCharacteristics, SimulationParameters] objects. - - Raises - ------ - RuntimeError - If the xml file is invalid, an exception is raised. - """ - - if is_valid_xml(xmlfile): - return parse_detector_chars(xmlfile), parse_simulation_params(xmlfile) - else: - raise RuntimeError("Invalid XML File!") - - -def config_from_xml( - xmlfile: str, fundcon: const.Fund_Constants = const.Fund_Constants() -) -> NssConfig: - r"""Parse the XML file into an NssConfig object. - - If the xml file is valid, parse the file into a DetectorCharacteristics and - SimulationParameters tuple. - - Parameters - ---------- - xmlfile: str - The input configuration xml file. - fundcon: Fund_Constants - A fundimental constants object. Defaults to default constructed object. - - Returns - ------- - NssConfig: - Parsed NssConfig object. - """ - - d, s = parseXML(xmlfile) - return NssConfig(d, s, fundcon) - - -def create_xml(filename: str, config: NssConfig = NssConfig()) -> None: - r"""Create an XML configuration file. - - Parameters - ---------- - filename: str - The name of the output xml file. - config: NssConfig - A NssConfig object from which to build the XML file. - """ - - nuspacesimparams = ET.Element("NuSpaceSimParams") - - detchar = ET.SubElement(nuspacesimparams, "DetectorCharacteristics") - detchar.set("Type", "Satellite") - detchar.set("Method", "Both") - - qeff = ET.SubElement(detchar, "QuantumEfficiency") - qeff.text = str(config.detector.quantum_efficiency) - - telaeff = ET.SubElement(detchar, "TelescopeEffectiveArea") - telaeff.set("Unit", "Sq.Meters") - telaeff.text = str(config.detector.telescope_effective_area) - - pethres = ET.SubElement(detchar, "PhotoElectronThreshold") - pethres.set("Preset", "true") - - detalt = ET.SubElement(detchar, "DetectorAltitude") - detalt.set("Unit", "km") - detalt.text = str(config.detector.altitude) - - detlat = ET.SubElement(detchar, "InitialDetectorLatitude") - detlat.set("Unit", "Degrees") - detlat.text = str(config.detector.lat_start) - - detlong = ET.SubElement(detchar, "InitialDetectorLongitude") - detlong.set("Unit", "Degrees") - detlong.text = str(config.detector.long_start) - - npe = ET.SubElement(pethres, "NPE") - npe.text = str(config.detector.photo_electron_threshold) - - detlow_freq = ET.SubElement(detchar, "LowFrequency") - detlow_freq.set("Unit", "MHz") - detlow_freq.text = str(config.detector.low_freq) - - dethigh_freq = ET.SubElement(detchar, "HighFrequency") - dethigh_freq.set("Unit", "MHz") - dethigh_freq.text = str(config.detector.high_freq) - - detSNRthres = ET.SubElement(detchar, "SNRThreshold") - detSNRthres.text = str(config.detector.det_SNR_thres) - - detNant = ET.SubElement(detchar, "NAntennas") - detNant.text = str(config.detector.det_Nant) - - detGain = ET.SubElement(detchar, "AntennaGain") - detGain.text = str(config.detector.det_gain) - - simparams = ET.SubElement(nuspacesimparams, "SimulationParameters") - simparams.set("DetectionMode", "Diffuse") - - cherangmax = ET.SubElement(simparams, "MaximumCherenkovAngle") - cherangmax.set("Unit", "Radians") - cherangmax.text = str(config.simulation.theta_ch_max) - - limbang = ET.SubElement(simparams, "AngleFromLimb") - limbang.set("Unit", "Radians") - limbang.text = str(config.simulation.ang_from_limb) - - eshowtype = ET.SubElement(simparams, "TauShowerType") - eshowtype.set("Preset", "true") - - fraceshow = ET.SubElement(eshowtype, "FracETauInShower") - fraceshow.text = str(config.simulation.e_shower_frac) - - nutauspectype = ET.SubElement(simparams, "NuTauEnergySpecType") - - if isinstance(config.simulation.spectrum, MonoSpectrum): - mono = ET.SubElement(nutauspectype, "MonoSpectrum") - nutauen = ET.SubElement(mono, "LogNuEnergy") - nutauen.text = str(config.simulation.spectrum.log_nu_tau_energy) - - elif isinstance(config.simulation.spectrum, PowerSpectrum): - power = ET.SubElement(nutauspectype, "PowerSpectrum") - sp1 = ET.SubElement(power, "PowerLawIndex") - sp2 = ET.SubElement(power, "LowerBound") - sp3 = ET.SubElement(power, "UpperBound") - sp1.text = str(config.simulation.spectrum.index) - sp2.text = str(config.simulation.spectrum.lower_bound) - sp3.text = str(config.simulation.spectrum.upper_bound) - - elif isinstance(config.simulation.spectrum, FileSpectrum): - filespec = ET.SubElement(nutauspectype, "FileSpectrum") - sp1 = ET.SubElement(filespec, "FilePath") - sp1.text = str(config.simulation.spectrum.path) - - cloudmodeltype = ET.SubElement(simparams, "CloudModelType") - - if isinstance(config.simulation.cloud_model, NoCloud): - ET.SubElement(cloudmodeltype, "NoCloudModel") - elif isinstance(config.simulation.cloud_model, MonoCloud): - mono = ET.SubElement(cloudmodeltype, "MonoCloudModel") - cth = ET.SubElement(mono, "CloudTopHeight") - cth.text = str(config.simulation.cloud_model.altitude) - elif isinstance(config.simulation.cloud_model, PressureMapCloud): - pmcm = ET.SubElement(cloudmodeltype, "PressureMapCloudModel") - m = ET.SubElement(pmcm, "Month") - m.text = str(config.simulation.cloud_model.month) - v = ET.SubElement(pmcm, "Version") - v.text = str(config.simulation.cloud_model.version) - - azimuthang = ET.SubElement(simparams, "AzimuthalAngle") - azimuthang.set("Unit", "Radians") - azimuthang.text = str(config.simulation.max_azimuth_angle) - - numtrajs = ET.SubElement(simparams, "NumTrajs") - numtrajs.text = str(config.simulation.N) - - ionosphere = ET.SubElement(simparams, "ModelIonosphere") - ionosphere.text = str(config.simulation.model_ionosphere) - - tec = ET.SubElement(simparams, "TEC") - tec.text = str(config.simulation.TEC) - - tecerr = ET.SubElement(simparams, "TECerr") - tecerr.text = str(config.simulation.TECerr) - - TauTableVer = ET.SubElement(simparams, "TauTableVer") - TauTableVer.text = str(config.simulation.tau_table_version) - - def indent(elem, level=0): - i = "\n" + level * " " - if len(elem): - if not elem.text or not elem.text.strip(): - elem.text = i + " " - if not elem.tail or not elem.tail.strip(): - elem.tail = i - for elem in elem: - indent(elem, level + 1) - if not elem.tail or not elem.tail.strip(): - elem.tail = i - else: - if level and (not elem.tail or not elem.tail.strip()): - elem.tail = i - - indent(nuspacesimparams) - - tree = ET.ElementTree(nuspacesimparams) - tree.write(filename, encoding="utf-8", xml_declaration=True, method="xml") diff --git a/test/core/test_config.py b/test/core/test_config.py index a50bff9..e28b7ba 100644 --- a/test/core/test_config.py +++ b/test/core/test_config.py @@ -31,31 +31,384 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +import tempfile +from datetime import datetime + +import astropy.units as u import numpy as np +import pytest +from pydantic import ValidationError + +from nuspacesim.config import ( + Detector, + NssConfig, + Simulation, + config_from_toml, + create_toml, +) + + +def test_detector_initial_position_default(): + ip = Detector.InitialPos() + assert ip.altitude == 525.0 + assert ip.latitude == 0 + assert ip.longitude == 0 + + assert ip.model_dump() == { + "altitude": "525.0 km", + "latitude": "0.0 deg", + "longitude": "0.0 deg", + } + + +def test_detector_initial_position_units(): + ip = Detector.InitialPos( + altitude=50000 * u.m, latitude=10.0 * u.degree, longitude=np.radians(10.0) + ) + assert ip.altitude == 50 + assert ip.latitude == np.radians(10) + assert ip.longitude == np.radians(10) + + assert ip.model_dump() == { + "altitude": "50.0 km", + "latitude": "10.0 deg", + "longitude": "10.0 deg", + } + + +def test_detector_initial_position_bad_units(): + with pytest.raises(ValidationError): + Detector.InitialPos( + altitude=50000 * u.s, latitude=10.0 * u.s, longitude=10.0 * u.s + ) + + +def test_detector_optical_default(): + a = Detector.Optical() + assert a.telescope_effective_area == 2.5 + assert a.quantum_efficiency == 0.2 + assert a.photo_electron_threshold == 10 + + assert a.model_dump() == { + "telescope_effective_area": "2.5 m2", + "quantum_efficiency": 0.2, + "photo_electron_threshold": 10.0, + } + + +def test_detector_optical_units(): + a = Detector.Optical( + telescope_effective_area=16, + quantum_efficiency=0.5, + photo_electron_threshold=100, + ) + assert a.telescope_effective_area == 16 + assert a.quantum_efficiency == 0.5 + assert a.photo_electron_threshold == 100 + + assert a.model_dump() == { + "telescope_effective_area": "16.0 m2", + "quantum_efficiency": 0.5, + "photo_electron_threshold": 100.0, + } + + b = Detector.Optical(telescope_effective_area=2.5e4 * u.cm * u.cm) + assert b.telescope_effective_area == 2.5 + assert b.quantum_efficiency == 0.2 + assert b.photo_electron_threshold == 10 + + assert b.model_dump() == { + "telescope_effective_area": "2.5 m2", + "quantum_efficiency": 0.2, + "photo_electron_threshold": 10.0, + } + + +def test_detector_optical_bad_units(): + with pytest.raises(ValidationError): + Detector.Optical(telescope_effective_area=2.5 * u.cm) + + +def test_detector_radio_default(): + radio = Detector.Radio() + assert radio.low_frequency == 30.0 + assert radio.high_frequency == 300.0 + assert radio.snr_threshold == 5.0 + assert radio.nantennas == 10 + assert radio.gain == 1.8 + + assert radio.model_dump() == { + "low_frequency": "30.0 MHz", + "high_frequency": "300.0 MHz", + "snr_threshold": 5.0, + "nantennas": 10, + "gain": "1.8 dB", + } + + +def test_detector_radio_units(): + radio = Detector.Radio( + low_frequency=40.0, + high_frequency=200.0, + snr_threshold=7.0, + nantennas=15, + gain=2.5, + ) + assert radio.low_frequency == 40.0 + assert radio.high_frequency == 200.0 + assert radio.snr_threshold == 7.0 + assert radio.nantennas == 15 + assert radio.gain == 2.5 + + assert radio.model_dump() == { + "low_frequency": "40.0 MHz", + "high_frequency": "200.0 MHz", + "snr_threshold": 7.0, + "nantennas": 15, + "gain": "2.5 dB", + } + + radio_with_units = Detector.Radio( + low_frequency=1 * u.kHz, high_frequency=1 * u.GHz, gain=3.0 * u.dB + ) + assert radio_with_units.low_frequency == 0.001 + assert radio_with_units.high_frequency == 1000 + assert radio_with_units.gain == 3.0 + + assert radio_with_units.model_dump() == { + "low_frequency": "0.001 MHz", + "high_frequency": "1000.0 MHz", + "snr_threshold": 5.0, + "nantennas": 10, + "gain": "3.0 dB", + } + + +def test_detector_radio_bad_units(): + with pytest.raises(ValidationError): + Detector.Radio(low_frequency=2.5 * u.cm) + + +def test_detector_radio_invalid(): + with pytest.raises(ValidationError): + Detector.Radio(low_frequency=2, high_frequency=1) + + +def test_detector_default(): + detector = Detector() + assert detector.name == "Default Name" + assert detector.initial_position.altitude == 525.0 + assert detector.initial_position.latitude == 0.0 + assert detector.initial_position.longitude == 0.0 + assert detector.optical.telescope_effective_area == 2.5 + assert detector.optical.quantum_efficiency == 0.2 + assert detector.optical.photo_electron_threshold == 10 + assert detector.radio.low_frequency == 30.0 + assert detector.radio.high_frequency == 300.0 + assert detector.radio.snr_threshold == 5.0 + assert detector.radio.nantennas == 10 + assert detector.radio.gain == 1.8 + + +def test_detector_custom_values(): + a = Detector( + name="Custom Detector", + initial_position=Detector.InitialPos( + altitude=600.0, latitude=1.0, longitude=2.0 + ), + optical=Detector.Optical( + telescope_effective_area=3.0, + quantum_efficiency=0.5, + photo_electron_threshold=20, + ), + radio=Detector.Radio( + low_frequency=40.0, + high_frequency=200.0, + snr_threshold=7.0, + nantennas=15, + gain=2.5, + ), + ) + assert a.name == "Custom Detector" + assert a.initial_position.altitude == 600.0 + assert a.initial_position.latitude == 1.0 + assert a.initial_position.longitude == 2.0 + assert a.optical.telescope_effective_area == 3.0 + assert a.optical.quantum_efficiency == 0.5 + assert a.optical.photo_electron_threshold == 20 + assert a.radio.low_frequency == 40.0 + assert a.radio.high_frequency == 200.0 + assert a.radio.snr_threshold == 7.0 + assert a.radio.nantennas == 15 + assert a.radio.gain == 2.5 + + +def test_detector_serialization(): + assert Detector().model_dump() == { + "name": "Default Name", + "initial_position": { + "altitude": "525.0 km", + "latitude": "0.0 deg", + "longitude": "0.0 deg", + }, + "optical": { + "telescope_effective_area": "2.5 m2", + "quantum_efficiency": 0.2, + "photo_electron_threshold": 10.0, + }, + "radio": { + "low_frequency": "30.0 MHz", + "high_frequency": "300.0 MHz", + "snr_threshold": 5.0, + "nantennas": 10, + "gain": "1.8 dB", + }, + } + + +################################################################################ + + +def test_default_simulation(): + a = Simulation() + assert a.mode == "Diffuse" + assert a.thrown_events == 1000 + assert a.max_cherenkov_angle == np.radians(3.0) + assert a.max_azimuth_angle == np.radians(360.0) + assert a.angle_from_limb == np.radians(7.0) + assert a.cherenkov_light_engine == "Default" + assert a.ionosphere is not None + assert a.ionosphere.total_electron_content == 10.0 + assert a.ionosphere.total_electron_error == 0.1 + assert a.tau_shower.id == "nupyprop" + assert a.tau_shower.etau_frac == 0.5 + assert a.tau_shower.table_version == "3" + assert isinstance(a.spectrum, Simulation.MonoSpectrum) + assert isinstance(a.cloud_model, Simulation.NoCloud) + + assert a.model_dump() == { + "mode": "Diffuse", + "thrown_events": 1000, + "max_cherenkov_angle": "3.0000000000000004 deg", + "max_azimuth_angle": "360.0 deg", + "angle_from_limb": "7.0 deg", + "cherenkov_light_engine": "Default", + "ionosphere": {"total_electron_content": 10.0, "total_electron_error": 0.1}, + "tau_shower": {"id": "nupyprop", "etau_frac": 0.5, "table_version": "3"}, + "spectrum": {"id": "monospectrum", "log_nu_energy": 8.0}, + "cloud_model": {"id": "no_cloud"}, + } + + +def test_custom_simulation(): + a = Simulation( + mode="ToO", + thrown_events=500, + max_cherenkov_angle=3.5 * u.deg, + max_azimuth_angle=270.0 * u.deg, + angle_from_limb=10.0 * u.deg, + ionosphere=Simulation.Ionosphere(total_electron_content=15.0), + tau_shower=Simulation.NuPyPropShower(etau_frac=0.6), + spectrum=Simulation.PowerSpectrum(index=2.5, lower_bound=7.0, upper_bound=11.0), + cloud_model=Simulation.MonoCloud(altitude=20.0), + ) + + assert a.mode == "ToO" + assert a.thrown_events == 500 + assert a.max_cherenkov_angle == np.radians(3.5) + assert a.max_azimuth_angle == np.radians(270.0) + assert a.angle_from_limb == np.radians(10.0) + assert a.ionosphere is not None + assert a.ionosphere.total_electron_content == 15.0 + assert a.tau_shower.etau_frac == 0.6 + assert isinstance(a.spectrum, Simulation.PowerSpectrum) + assert isinstance(a.cloud_model, Simulation.MonoCloud) + assert a.spectrum.index == 2.5 + assert a.spectrum.lower_bound == 7.0 + assert a.spectrum.upper_bound == 11.0 + assert a.cloud_model.altitude == 20.0 + + +def test_invalid_cherenkov_angle(): + with pytest.raises(TypeError): + Simulation(max_cherenkov_angle="invalid_angle") + + +def test_invalid_ionosphere_content(): + with pytest.raises(ValidationError): + Simulation( + ionosphere=Simulation.Ionosphere(total_electron_content="invalid_content") + ) + + +def test_invalid_spectrum_type(): + with pytest.raises(ValidationError): + Simulation(spectrum="invalid_spectrum_type") + + +def test_invalid_spectrum_values(): + with pytest.raises(ValidationError): + Simulation(spectrum=Simulation.MonoSpectrum(log_nu_energy="invalid_value")) + + +def test_invalid_cloud_model_type(): + with pytest.raises(ValidationError): + Simulation(cloud_model="invalid_cloud_model_type") + + +def test_invalid_cloud_model_values(): + with pytest.raises(ValidationError): + Simulation(cloud_model=Simulation.MonoCloud(altitude="invalid_value")) + + +def test_pressure_map_cloud_default(): + cloud = Simulation.PressureMapCloud() + assert cloud.id == "pressure_map" + assert cloud.month == 1 + assert cloud.version == 0 + + +def test_pressure_map_cloud_custom_values(): + cloud = Simulation.PressureMapCloud(month=6, version="2023") + assert cloud.id == "pressure_map" + assert cloud.month == 6 + assert cloud.version == "2023" + + +def test_pressure_map_cloud_serialization(): + cloud = Simulation.PressureMapCloud(month=8, version=2) + assert cloud.model_dump() == {"id": "pressure_map", "month": 8, "version": 2} -import nuspacesim as nss +def test_pressure_map_cloud_valid_month(): + assert Simulation.PressureMapCloud.valid_month(6) == 6 + assert Simulation.PressureMapCloud.valid_month("06") == 6 + assert Simulation.PressureMapCloud.valid_month("June") == 6 + assert Simulation.PressureMapCloud.valid_month("Jun") == 6 + assert ( + Simulation.PressureMapCloud.valid_month( + datetime.strptime("2023-06-01", "%Y-%m-%d") + ) + == 6 + ) -def test_detector_characteristics(): - dc = nss.DetectorCharacteristics() - assert dc.altitude > 0.0 - dc = nss.DetectorCharacteristics(altitude=200.1) - assert dc.altitude == 200.1 +def test_pressure_map_cloud_invalid_month(): + with pytest.raises(ValueError): + Simulation.PressureMapCloud.valid_month(13) + with pytest.raises(ValueError): + Simulation.PressureMapCloud.valid_month("Invalid") -def test_simulation_params(): - sp = nss.SimulationParameters() - assert sp.N > 0 - assert np.log10(sp.nu_tau_energy) == sp.log_nu_tau_energy - assert np.sin(sp.theta_ch_max) == sp.sin_theta_ch_max +####################################################################################### +def test_config_serialization(): + with tempfile.NamedTemporaryFile( + mode="w+", suffix=".toml", delete=False + ) as tmpfile: + tmpfile_name = tmpfile.name + a = NssConfig() + create_toml(tmpfile_name, a) + loaded_config = config_from_toml(tmpfile_name) -def test_nss_config(): - nc1 = nss.NssConfig() - dc = nss.DetectorCharacteristics() - assert nc1.detector == dc - sp = nss.SimulationParameters() - assert nc1.simulation == sp - nc2 = nss.NssConfig(detector=dc, simulation=sp) - assert nc1 == nc2 + assert loaded_config.model_dump() == a.model_dump()