#!/usr/bin/python3
"""Definition of site model."""
import logging
import numpy as np
from astropy import units as u
from simtools.model.model_parameter import ModelParameter
from simtools.simulation.configuration import get_model_writer
from simtools.utils import names
[docs]
class SiteModel(ModelParameter):
"""
Representation of an observatory site model.
The site model includes (among others):
- Reference point coordinates
- Array elements
- Geomagnetic field parameters
- Atmospheric parameters
- NSB parameters
Parameters
----------
site: str
Site name (e.g., South or North).
model_version: str or list
Model version or list of model versions (in which case only the first one is used).
label: str, optional
Instance label.
overwrite_model_parameter_dict: dict, optional
Dictionary to overwrite model parameters from the model repository with provided values.
ignore_software_version: bool, optional
If True, ignore software version checks for deprecated parameters.
model_directory: pathlib.Path or str, optional
Directory for generated model assets and sim_telarray configuration files.
parameter_names: iterable of str, optional
If supplied, load only these site parameters.
load_simulation_software_parameters: bool
If False, do not load CORSIKA and sim_telarray configuration parameters.
"""
def __init__(
self,
site,
model_version,
label=None,
overwrite_model_parameter_dict=None,
ignore_software_version=False,
model_reader=None,
context=None,
model_directory=None,
parameter_names=None,
load_simulation_software_parameters=True,
):
"""Initialize SiteModel."""
self._logger = logging.getLogger(__name__)
self._logger.debug("Init SiteModel for site %s", site)
super().__init__(
site=site,
model_version=model_version,
label=label,
collection="sites",
overwrite_model_parameter_dict=overwrite_model_parameter_dict,
ignore_software_version=ignore_software_version,
model_reader=model_reader,
context=context,
model_directory=model_directory,
parameter_names=parameter_names,
load_simulation_software_parameters=load_simulation_software_parameters,
)
[docs]
def get_reference_point(self):
"""
Get reference point coordinates.
Ground coordinates are calculated relative to this point.
Returns
-------
dict
Reference point coordinates.
"""
return {
"center_altitude": self.get_parameter_value_with_unit("reference_point_altitude"),
"center_northing": self.get_parameter_value_with_unit("reference_point_utm_north"),
"center_easting": self.get_parameter_value_with_unit("reference_point_utm_east"),
"epsg_code": self.get_parameter_value("epsg_code"),
}
[docs]
def get_corsika_site_parameters(self, config_file_style=False, model_directory=None):
"""
Get site-related CORSIKA parameters.
Parameters are returned with units wherever possible ('config_file_style=False')
or in CORSIKA-expected coordinates.
Parameters
----------
config_file_style: bool
Return using CORSIKA config file syntax
model_directory: Path, optional
Model directory to use for file paths
Returns
-------
dict
Site-related CORSIKA parameters.
"""
if config_file_style:
# pylint: disable-next=import-outside-toplevel
from simtools.corsika.corsika_config import get_site_config_parameters
return get_site_config_parameters(self, model_directory)
return {
"corsika_observation_level": self.get_parameter_value_with_unit(
"corsika_observation_level"
),
"geomag_horizontal": self.get_parameter_value_with_unit("geomag_horizontal"),
"geomag_vertical": self.get_parameter_value_with_unit("geomag_vertical"),
"geomag_rotation": self.get_parameter_value_with_unit("geomag_rotation"),
}
[docs]
def get_array_elements_for_layout(self, layout_name):
"""
Return list of array elements for a given array layout.
If ``layout_name`` is not found in the model repository but is a valid,
concrete telescope name
(e.g., ``MSTN-05``) belonging to this site, a single-telescope layout is returned
automatically.
Parameters
----------
layout_name: str
Name of the array layout or a single telescope name (e.g., ``MSTN-05``).
Returns
-------
list
List of array elements.
"""
layouts = self.get_parameter_value("array_layouts")
for layout in layouts:
if layout["name"].lower() == layout_name.lower():
return layout["elements"]
validated_name = self._validate_as_single_telescope(layout_name)
if validated_name is not None:
self._logger.debug(
f"Array layout '{layout_name}' not found in the model repository; "
"treating as single-telescope layout."
)
return [validated_name]
raise ValueError(f"Array layout '{layout_name}' not found in '{self.site}' site model.")
def _validate_as_single_telescope(self, layout_name):
"""
Validate that layout_name is a concrete telescope belonging to this site.
Parameters
----------
layout_name: str
Candidate telescope name.
Returns
-------
str or None
Validated and normalized telescope name, or None if not a valid single telescope.
"""
try:
validated_name = names.validate_array_element_name(layout_name)
if names.get_collection_name_from_array_element_name(validated_name) != "telescopes":
return None
if names.is_design_type(validated_name):
return None
if names.get_site_from_array_element_name(validated_name) != self.site:
return None
if validated_name not in self.get_array_elements_of_type(
names.get_array_element_type_from_name(validated_name)
):
return None
except ValueError:
return None
return validated_name
[docs]
def get_array_elements_of_type(self, array_element_type):
"""
Get all array elements of a given type.
Parameters
----------
array_element_type : str
Type of the array element (e.g. LSTN, MSTS)
Returns
-------
list
List with array elements.
"""
try:
collection = names.get_collection_name_from_array_element_name(array_element_type)
except ValueError:
return []
return self.model_reader.get_array_elements_of_type(
array_element_type=array_element_type,
model_version=self.model_version,
collection=collection,
)
[docs]
def get_list_of_array_layouts(self):
"""
Get list of available array layouts.
Returns
-------
list
List of available array layouts
"""
return [layout["name"] for layout in self.get_parameter_value("array_layouts")]
[docs]
def export_atmospheric_transmission_file(self, model_directory):
"""
Export the atmospheric profile source and sim_telarray table files.
Parameters
----------
model_directory: Path
Model directory to export the file to.
"""
return get_model_writer(self, context=self.context).export_atmospheric_transmission_file(
model_directory=model_directory
)
[docs]
def get_nsb_integrated_flux(self, wavelength_min=300 * u.nm, wavelength_max=650 * u.nm):
"""
Get the integrated flux for the NSB (Night Sky Background) model.
Returns
-------
float
Integrated flux value.
"""
table = self.get_parameter_table("nsb_spectrum")
table.sort("wavelength")
wavelength_column = table["wavelength"]
rate_column = table["differential_photon_rate"]
wl = getattr(wavelength_column, "quantity", wavelength_column).to(u.nm)
rate = getattr(rate_column, "quantity", rate_column).to(1 / (u.nm * u.cm**2 * u.ns * u.sr))
mask = (wl >= wavelength_min) & (wl <= wavelength_max)
integral_cm2 = np.trapezoid(rate[mask], wl[mask])
self._logger.debug(
f"NSB integral between {wavelength_min} and {wavelength_max}: {integral_cm2}"
)
return float(getattr(integral_cm2, "value", integral_cm2))
