diff --git a/docs/changes/2190.feature.md b/docs/changes/2190.feature.md
new file mode 100644
index 0000000000..0bbcd13eb7
--- /dev/null
+++ b/docs/changes/2190.feature.md
@@ -0,0 +1,2 @@
+Major improvements to simulation production configuration preparation with clear separation introduced between observation grid generation and backend submission code.
+Allow NSHOW to follow a pre-defined power law for simulation production grid definition.
diff --git a/docs/source/api-reference/production_configuration.md b/docs/source/api-reference/production_configuration.md
index 7f1744570f..b8cf346fb2 100644
--- a/docs/source/api-reference/production_configuration.md
+++ b/docs/source/api-reference/production_configuration.md
@@ -43,15 +43,16 @@ the calculation of the number of events to be simulated given a pre-determined m
    :members:
 ```
 
-(generate-production-grid)=
+(corsika-limits-lookup)=
 
-## generate_production_grid
+## corsika_limits_lookup
 
 ```{eval-rst}
-.. automodule:: production_configuration.generate_production_grid
+.. automodule:: production_configuration.corsika_limits_lookup
    :members:
 ```
 
+
 (interpolation-handler)=
 
 ## interpolation_handler
@@ -61,6 +62,35 @@ the calculation of the number of events to be simulated given a pre-determined m
    :members:
 ```
 
+(job-grid-io)=
+
+## job_grid_io
+
+```{eval-rst}
+.. automodule:: production_configuration.job_grid_io
+   :members:
+```
+
+(observation-grid)=
+
+## observation_grid
+
+```{eval-rst}
+.. automodule:: production_configuration.observation_grid
+   :members:
+```
+
+(simulation-jobs)=
+
+## simulation_jobs
+
+```{eval-rst}
+.. automodule:: production_configuration.simulation_jobs
+   :members:
+```
+
+(merge-corsika-limits)=
+
 ## merge_corsika_limits
 
 ```{eval-rst}
diff --git a/docs/source/api-reference/visualization.md b/docs/source/api-reference/visualization.md
index a4cb14f40b..6d66c8139b 100644
--- a/docs/source/api-reference/visualization.md
+++ b/docs/source/api-reference/visualization.md
@@ -131,8 +131,6 @@ the visualization module.
    :members:
 ```
 
-(merge-corsika-limits)=
-
 
 (plot-corsika-limits)=
 
diff --git a/src/simtools/applications/plot_production_grid.py b/src/simtools/applications/plot_production_grid.py
index f1500dcae6..73165b4f65 100644
--- a/src/simtools/applications/plot_production_grid.py
+++ b/src/simtools/applications/plot_production_grid.py
@@ -17,7 +17,7 @@
     Model version used to read the site reference coordinates.
 observation_time (str, optional)
     Observation time in UTC ISO format used for Alt/Az <-> RA/Dec transformations.
-    If omitted, the application uses ``metadata.observing_time_utc`` from the
+    If omitted, the application uses ``metadata.time_of_observation_utc`` from the
     grid file when available.
 plot_ra_dec_tracks (flag, optional)
     If provided, plot RA/Dec guide tracks on top of the sky projection. When native
diff --git a/src/simtools/applications/production_generate_grid.py b/src/simtools/applications/production_generate_grid.py
index 7331d3f36f..65684e01cd 100644
--- a/src/simtools/applications/production_generate_grid.py
+++ b/src/simtools/applications/production_generate_grid.py
@@ -1,46 +1,30 @@
 #!/usr/bin/python3
 
 r"""
-Generate a grid of simulation points using flexible axes definitions.
-
-This application generates a grid of simulation points based on the provided axes
-definitions. The axes definitions (range, binning) are specified in a file.
-The viewcone, radius and energy thresholds are provided as a lookup table and
-are interpolated based on the generated grid points. The generated grid points are
-filtered based on the specified telescope IDs and the limits from the lookup table.
-The generated grid points are saved to a file.
-It can also convert the generated points to RA/Dec coordinates if the selected
-coordinate system is 'ra_dec'.
-
-For ``coordinate_system='ra_dec'``, the underlying grid generation supports
-declination-line sampling with hour-angle spacing and applies zenith-angle
-filtering based on the configured zenith range in that mode.
-When explicit ``ra`` / ``dec`` axes are provided, all YAML-defined grid points are
-preserved in the serialized output.
+Generate simulation job grid for production configurations.
+
+This application expands simulation job rows and writes them to disk as ECSV files.
+It supports both:
+
+- explicit cartesian job-grid configuration (primary, zenith, energy range, etc.), and
+- axes-based production-grid configuration with optional ``ra_dec`` coordinate handling
+  and lookup-table interpolation.
 
 Command line arguments
 ----------------------
-axes (str, required)
-    Path to a YAML or JSON file defining the axes of the grid.
-coordinate_system (str, optional, default='zenith_azimuth')
-    The coordinate system for the grid generation ('zenith_azimuth' or 'ra_dec').
-    In ``ra_dec`` mode, observing location/time are used to build sky directions and
-    derive corresponding zenith/azimuth values for interpolation (ICRS/J2000 frame).
-observing_time (str, optional)
+axis (repeatable)
+    Compact axis definition in the form
+    ``--axis <name> <min> <unit> <max> <unit> <binning> [scaling]``.
+    Example: ``--axis azimuth 310 deg 20 deg 3 linear``.
+time_of_observation (str, optional)
     Time of the observation in UTC (format: 'YYYY-MM-DD HH:MM:SS').
-    Used only in ``ra_dec`` mode (for coordinate transforms and sidereal-time
-    sampling). Ignored in ``zenith_azimuth`` mode.
-lookup_table (str, required)
+    Used only if RA/Dec axes are provided (for coordinate transforms and sidereal-time
+    sampling). Ignored otherwise.
+corsika_limits (str, optional)
     Path to the lookup table for simulation limits. The table should contain
-    varying azimuth and/or zenith angles.
-telescope_ids (list of str, optional)
-    List of telescope names used to filter the lookup table rows
-    (e.g. ``MSTN-15``).
-simtel_file (str, optional)
-    Path to a sim_telarray file used only when lookup-table telescope selections
-    are stored as numeric telescope IDs.
-output_file (str, optional, default='grid_output.ecsv')
-    Output file for the generated grid points (default: 'grid_output.ecsv').
+    varying azimuth and/or zenith angles for the selected array layout.
+output_file (str, optional, default='job_grid.ecsv')
+    Output file for the generated executable job grid.
 
 
 Example
@@ -50,11 +34,12 @@
 .. code-block:: console
 
         simtools-production-generate-grid --site North --model_version 6.0.2 \
-            --axes tests/resources/production_grid_generation_axes_definition.yml \
-            --coordinate_system zenith_azimuth \
-            --lookup_table tests/resources/corsika_simulation_limits/
-                merged_corsika_limits_for_test.ecsv \
-            --telescope_ids MSTN-15
+            --array_layout_name alpha \
+            --axis azimuth 310 deg 20 deg 3 linear \
+            --axis zenith 30 deg 40 deg 2 linear \
+            --axis nsb 4 MHz 5 MHz 2 linear \
+            --axis offset 0 deg 10 deg 2 linear \
+            --corsika_limits tests/resources/corsika_simulation_limits/merged_corsika_limits.ecsv
 
 To generate an all-sky RA/Dec direction grid and serialize output in RA/Dec,
 execute:
@@ -62,44 +47,37 @@
 .. code-block:: console
 
         simtools-production-generate-grid --site North --model_version 6.0.2 \
-            --axes tests/resources/production_grid_generation_axes_definition_radec.yml \
-            --coordinate_system ra_dec --observing_time "2017-09-16 00:00:00" \
-            --lookup_table tests/resources/corsika_simulation_limits/
-                merged_corsika_limits_for_test.ecsv \
-            --telescope_ids MSTN-15
+            --array_layout_name alpha \
+            --axis ra 0 deg 360 deg 36 linear \
+            --axis dec -90 deg 90 deg 18 linear \
+            --axis nsb 4 MHz 4 MHz 1 linear \
+            --axis offset 0 deg 10 deg 2 linear \
+            --time_of_observation "2017-09-16 00:00:00" \
+            --corsika_limits tests/resources/corsika_simulation_limits/merged_corsika_limits.ecsv
 """
 
-from pathlib import Path
-
-from astropy.coordinates import EarthLocation
-from astropy.time import Time
-
 from simtools.application_control import build_application
-from simtools.io.ascii_handler import collect_data_from_file
-from simtools.model.site_model import SiteModel
-from simtools.production_configuration.generate_production_grid import GridGeneration
+from simtools.production_configuration.job_grid_io import serialize_job_grid
+from simtools.production_configuration.simulation_jobs import (
+    build_job_grid_metadata,
+    build_simulation_jobs,
+)
 
 
 def _add_arguments(parser):
     """Register application-specific command line arguments."""
     parser.add_argument(
-        "--axes",
-        type=str,
-        required=True,
-        help="Path to a file defining the grid axes.",
-    )
-    parser.add_argument(
-        "--coordinate_system",
-        type=str,
-        default="zenith_azimuth",
+        "--axis",
+        action="append",
+        nargs="+",
+        required=False,
         help=(
-            "Coordinate system ('zenith_azimuth' or 'ra_dec'). "
-            "In 'ra_dec' mode, sky directions are generated using observing"
-            " location/time and converted to zenith/azimuth for interpolation."
+            "Compact axis definition: --axis <name> <min> <unit> <max> <unit> <binning> "
+            "[scaling]. May be repeated."
         ),
     )
     parser.add_argument(
-        "--observing_time",
+        "--time_of_observation",
         type=str,
         required=False,
         help=(
@@ -109,103 +87,63 @@ def _add_arguments(parser):
     parser.add_argument(
         "--output_file",
         type=str,
-        default="grid_output.ecsv",
-        help="Output file for the generated grid points (default: 'grid_output.ecsv').",
+        default="job_grid.ecsv",
+        help="Output file for the generated executable job grid.",
     )
     parser.add_argument(
-        "--telescope_ids",
+        "--corsika_limits",
         type=str,
-        nargs="*",
-        default=None,
-        help=(
-            "List of telescope names used to get specific limits from the lookup table "
-            "(e.g. MSTN-15)."
-        ),
-    )
-    parser.add_argument(
-        "--lookup_table",
-        type=str,
-        required=True,
+        required=False,
         help="Path to the lookup table for simulation limits. "
         "Table required with varying azimuth and or zenith angle. ",
     )
     parser.add_argument(
-        "--simtel_file",
-        type=str,
+        "--number_of_runs",
+        help="Number of runs to be simulated.",
+        type=int,
+        required=False,
+        default=1,
+    )
+    parser.add_argument(
+        "--nshow_power_index",
+        help=(
+            "Power-law index used to scale the baseline nshow with the geometric-mean energy "
+            "of each energy_range entry."
+        ),
+        type=float,
         required=False,
+        default=None,
+    )
+    parser.add_argument(
+        "--nshow_reference_energy",
         help=(
-            "Optional path to a sim_telarray file used to map sim_telarray telescope IDs "
-            "to telescope names when lookup-table selections are numeric IDs."
+            "Reference energy for nshow power-law scaling (for example: '100 GeV'). "
+            "Required together with --nshow_power_index."
         ),
+        type=str,
+        required=False,
+        default=None,
     )
 
 
-def load_axes(file_path: str):
-    """
-    Load axes definitions from a YAML or JSON file.
-
-    Parameters
-    ----------
-    file_path : str
-        Path to the axes YAML or JSON file.
-
-    Returns
-    -------
-    list[dict]
-        List of axes definitions with Quantity values.
-    """
-    if not Path(file_path).exists():
-        raise FileNotFoundError(f"Axes file {file_path} not found.")
-
-    return collect_data_from_file(file_path)
-
-
 def main():
     """See CLI description."""
     app_context = build_application(
         initialization_kwargs={
             "db_config": True,
-            "simulation_model": ["version", "site", "model_version"],
+            "preserve_by_version_keys": ["array_layout_name"],
+            "simulation_model": ["site", "layout", "telescope", "model_version"],
+            "simulation_configuration": {"software": None, "corsika_configuration": ["all"]},
         },
     )
 
-    output_filepath = app_context.io_handler.get_output_file(app_context.args["output_file"])
-
-    axes = load_axes(app_context.args["axes"])
-    site_model = SiteModel(
-        model_version=app_context.args["model_version"],
-        site=app_context.args["site"],
+    job_rows = build_simulation_jobs(app_context.args)
+    serialize_job_grid(
+        job_rows=job_rows,
+        output_file=app_context.io_handler.get_output_file(app_context.args["output_file"]),
+        metadata=build_job_grid_metadata(app_context.args),
     )
 
-    ref_lat = site_model.get_parameter_value_with_unit("reference_point_latitude")
-    ref_long = site_model.get_parameter_value_with_unit("reference_point_longitude")
-    altitude = site_model.get_parameter_value_with_unit("reference_point_altitude")
-
-    observing_location = EarthLocation(lat=ref_lat, lon=ref_long, height=altitude)
-
-    coordinate_system = app_context.args["coordinate_system"]
-    observing_time = None
-    if app_context.args.get("observing_time"):
-        observing_time = Time(app_context.args["observing_time"], scale="utc")
-    elif coordinate_system == "ra_dec":
-        observing_time = Time.now()
-
-    grid_gen = GridGeneration(
-        axes=axes,
-        coordinate_system=coordinate_system,
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=app_context.args["lookup_table"],
-        telescope_ids=app_context.args["telescope_ids"],
-        simtel_file=app_context.args.get("simtel_file"),
-    )
-
-    grid_points = grid_gen.generate_grid()
-
-    if coordinate_system == "ra_dec":
-        grid_points = grid_gen.convert_coordinates(grid_points)
-    grid_gen.serialize_grid_points(grid_points, output_file=output_filepath)
-
 
 if __name__ == "__main__":
     main()
diff --git a/src/simtools/applications/simulate_prod_htcondor_generator.py b/src/simtools/applications/simulate_prod_htcondor_generator.py
index e4dc97f76b..1a7ff0372e 100644
--- a/src/simtools/applications/simulate_prod_htcondor_generator.py
+++ b/src/simtools/applications/simulate_prod_htcondor_generator.py
@@ -3,7 +3,7 @@
 r"""
 Generate a run script and submit file for HT Condor job submission of a simulation production.
 
-This tool generates HTCondor submission files for one or more simulation production grids and
+This tool generates HTCondor submission files from a pre-generated executable job grid and
 supports either a single Apptainer image or a label-to-image mapping for multi-image submissions.
 For each image label, it writes a dedicated ``simulate_prod.submit.<label>.condor`` file and a
 matching ``simulate_prod.submit.<label>.params.txt`` file. When only one default image is used,
@@ -37,6 +37,8 @@
 For the single-image default case, use ``condor_submit simulate_prod.submit.condor``.
 
 Simulation data products are written to the directory controlled by ``--simulation_output``.
+The executable job grid must be prepared first, for example using
+``simtools-production-generate-grid``.
 
 Command line arguments
 ----------------------
@@ -49,6 +51,8 @@
     Directory for HTCondor log files. Defaults to ``output_path/htcondor_logs``.
 simulation_output (str, optional)
     Base directory for simulation output packages passed through as ``pack_for_grid_register``.
+job_grid_file (str, required)
+    Path to the pre-generated executable job grid file in ECSV format.
 priority (int, optional)
     Job priority (default: 1).
 
@@ -63,11 +67,10 @@
 def _add_arguments(parser):
     """Register application-specific command line arguments."""
     parser.add_argument(
-        "--number_of_runs",
-        help="Number of runs to be simulated.",
-        type=int,
+        "--job_grid_file",
+        help="Path to a pre-generated executable job grid file.",
+        type=str,
         required=True,
-        default=1,
     )
     parser.add_argument(
         "--priority",
@@ -90,25 +93,11 @@ def _add_arguments(parser):
         required=False,
         default="./simtools-output",
     )
-    parser.add_argument(
-        "--corsika_limits",
-        help="Path to an ECSV file with CORSIKA limits.",
-        type=str,
-        required=False,
-        default=None,
-    )
 
 
 def main():
     """See CLI description."""
-    app_context = build_application(
-        initialization_kwargs={
-            "db_config": False,
-            "preserve_by_version_keys": ["array_layout_name"],
-            "simulation_model": ["site", "layout", "telescope", "model_version"],
-            "simulation_configuration": {"software": None, "corsika_configuration": ["all"]},
-        },
-    )
+    app_context = build_application()
 
     htcondor_script_generator.generate_submission_script(app_context.args)
 
diff --git a/src/simtools/job_execution/htcondor_script_generator.py b/src/simtools/job_execution/htcondor_script_generator.py
index 705ee7868d..60ca90a4e7 100644
--- a/src/simtools/job_execution/htcondor_script_generator.py
+++ b/src/simtools/job_execution/htcondor_script_generator.py
@@ -8,33 +8,16 @@
 
 """
 
-import ast
-import itertools
 import logging
 import re
 from pathlib import Path
 
 import astropy.units as u
 
-import simtools.version as simtools_version
-from simtools.configuration import defaults
+from simtools.production_configuration.job_grid_io import read_job_grid
 
 _logger = logging.getLogger(__name__)
 
-_GRID_AXES = [
-    "primary",
-    "azimuth_angle",
-    "zenith_angle",
-    "model_version",
-    "corsika_le_interaction",
-    "corsika_he_interaction",
-]
-
-_GRID_AXIS_DEFAULTS = {
-    "corsika_le_interaction": defaults.CORSIKA_LE_INTERACTION,
-    "corsika_he_interaction": defaults.CORSIKA_HE_INTERACTION,
-}
-
 _PARAMS_FIELDS = [
     "apptainer_label",
     "primary",
@@ -44,8 +27,13 @@
     "energy_min_unit",
     "energy_max_value",
     "energy_max_unit",
+    "core_scatter_number",
     "core_scatter_max_value",
     "core_scatter_max_unit",
+    "view_cone_min_value",
+    "view_cone_min_unit",
+    "view_cone_max_value",
+    "view_cone_max_unit",
     "nshow",
     "model_version",
     "array_layout_name",
@@ -55,44 +43,7 @@
     "pack_for_grid_register",
 ]
 
-
-def _normalize_to_list(value):
-    """Normalize scalar values to lists of length one."""
-    if isinstance(value, list):
-        return value
-    if isinstance(value, tuple):
-        return list(value)
-    return [value]
-
-
-def _normalize_grid_axes(args_dict):
-    """Return normalized grid axes for cartesian product expansion."""
-    return {
-        axis: (
-            _normalize_to_list(args_dict[axis])
-            if axis in args_dict and args_dict[axis] is not None
-            else [_GRID_AXIS_DEFAULTS[axis]]
-            if axis in _GRID_AXIS_DEFAULTS
-            else [None]
-        )
-        for axis in _GRID_AXES
-    }
-
-
-def _normalize_energy_ranges(energy_range):
-    """Normalize energy range argument to a list of (e_min, e_max) pairs."""
-    if isinstance(energy_range, tuple) and len(energy_range) == 2:
-        return [energy_range]
-
-    if isinstance(energy_range, list):
-        if len(energy_range) == 2 and all(hasattr(item, "to") for item in energy_range):
-            return [(energy_range[0], energy_range[1])]
-        if all(isinstance(item, (list, tuple)) and len(item) == 2 for item in energy_range):
-            return [tuple(item) for item in energy_range]
-
-    raise ValueError(
-        "energy_range must be one pair (e_min, e_max) or a list of (e_min, e_max) pairs."
-    )
+_REQUIRED_JOB_GRID_METADATA = ("site", "simulation_software")
 
 
 def _resolve_apptainer_images(apptainer_image_arg):
@@ -150,15 +101,19 @@ def _format_param_value(value, field_name):
     if field_name in ("apptainer_label", "pack_for_grid_register"):
         return _sanitize_label_for_params(value)
 
-    if field_name in ("energy_min_value", "energy_max_value"):
-        return _format_quantity(value, default_unit=u.GeV)
+    if field_name == "core_scatter_number":
+        return f"{int(value)}"
 
-    if field_name == "core_scatter_max_value":
-        return _format_quantity(
-            value,
-            default_unit=u.m,
-            convert_to=u.m,
-        )
+    quantity_fields = {
+        "energy_min_value": (u.GeV, None),
+        "energy_max_value": (u.GeV, None),
+        "core_scatter_max_value": (u.m, u.m),
+        "view_cone_min_value": (u.deg, u.deg),
+        "view_cone_max_value": (u.deg, u.deg),
+    }
+    if field_name in quantity_fields:
+        default_unit, convert_to = quantity_fields[field_name]
+        return _format_quantity(value, default_unit=default_unit, convert_to=convert_to)
 
     if field_name in ("azimuth_angle", "zenith_angle"):
         if isinstance(value, u.Quantity):
@@ -179,141 +134,11 @@ def _sanitize_label_for_params(label):
     return re.sub(r"\s+", "_", str(label).strip())
 
 
-def _resolve_array_layout_name(array_layout_name, model_version):
-    """Resolve array layout configuration for a specific model version."""
-    if isinstance(array_layout_name, list) and len(array_layout_name) == 1:
-        array_layout_name = array_layout_name[0]
-
-    # Configurator/_fill_config stringifies dict values when rebuilding argparse arguments.
-    if isinstance(array_layout_name, str) and array_layout_name.strip().startswith("{"):
-        try:
-            parsed_layout = ast.literal_eval(array_layout_name)
-            if isinstance(parsed_layout, dict):
-                array_layout_name = parsed_layout
-        except (SyntaxError, ValueError):
-            return array_layout_name
-
-    if not isinstance(array_layout_name, dict) or list(array_layout_name) != ["by_version"]:
-        return array_layout_name
-
-    resolved = simtools_version.resolve_by_version(
-        {"array_layout_name": array_layout_name}, model_version
-    )
-    return resolved["array_layout_name"]
-
-
-def _get_energy_range_for_zenith_angle(zenith_angle, energy_range_pair, corsika_limits):
-    """
-    Return a zenith-dependent energy range pair or None to skip the simulation step.
-
-    Dummy function - to be implemented.
-    """
-    _ = (zenith_angle, corsika_limits)
-    return energy_range_pair
-
-
-def _get_core_scatter_max_for_zenith_angle(zenith_angle, core_scatter, corsika_limits):
-    """Return zenith-dependent max core-scatter value.
-
-    Dummy function - to be implemented.
-    """
-    _ = (zenith_angle, corsika_limits)
-    return core_scatter[1]
-
-
-def _get_nshow_for_energy_range_and_zenith_angle(
-    energy_range_pair, zenith_angle, nshow, corsika_limits
-):
-    """Return nshow that may depend on energy range and zenith angle.
-
-    Dummy function - to be implemented.
-    """
-    _ = (energy_range_pair, zenith_angle, corsika_limits)
-    return nshow
-
-
-def _build_job_specs(args_dict, image_labels):
-    """Build backend-agnostic job specs from comparison and production grids."""
-    grid_axes = _normalize_grid_axes(args_dict)
-    energy_ranges = _normalize_energy_ranges(args_dict["energy_range"])
-    base_pack_dir = args_dict.get("simulation_output") or "simtools-output"
-    corsika_limits = args_dict.get("corsika_limits")
-    core_scatter = args_dict["core_scatter"]
-    nshow = args_dict["nshow"]
-
-    combinations = list(
-        itertools.product(
-            grid_axes["primary"],
-            grid_axes["azimuth_angle"],
-            grid_axes["zenith_angle"],
-            grid_axes["model_version"],
-            grid_axes["corsika_le_interaction"],
-            grid_axes["corsika_he_interaction"],
-            energy_ranges,
-        )
-    )
-
-    number_of_runs = args_dict.get("number_of_runs", 1)
-    run_number = int(args_dict.get("run_number") or 1)
-
-    job_specs = []
-    for label in image_labels:
-        params_label = _sanitize_label_for_params(label)
-        row_index = 0
-        for (
-            primary,
-            azimuth,
-            zenith,
-            model_version,
-            corsika_le,
-            corsika_he,
-            energy_range_pair,
-        ) in combinations:
-            selected_energy_range_pair = energy_range_pair
-            selected_core_scatter_max = core_scatter[1]
-            selected_nshow = nshow
-
-            if corsika_limits is not None:
-                selected_energy_range_pair = _get_energy_range_for_zenith_angle(
-                    zenith, energy_range_pair, corsika_limits
-                )
-                if selected_energy_range_pair is None:
-                    continue
-                selected_core_scatter_max = _get_core_scatter_max_for_zenith_angle(
-                    zenith, core_scatter, corsika_limits
-                )
-                selected_nshow = _get_nshow_for_energy_range_and_zenith_angle(
-                    selected_energy_range_pair, zenith, nshow, corsika_limits
-                )
-
-            for _ in range(number_of_runs):
-                job_specs.append(
-                    {
-                        "apptainer_label": str(label),
-                        "primary": primary,
-                        "azimuth_angle": azimuth,
-                        "zenith_angle": zenith,
-                        "model_version": model_version,
-                        "array_layout_name": args_dict.get("array_layout_name"),
-                        "corsika_le_interaction": corsika_le,
-                        "corsika_he_interaction": corsika_he,
-                        "energy_min": selected_energy_range_pair[0],
-                        "energy_max": selected_energy_range_pair[1],
-                        "core_scatter_max": selected_core_scatter_max,
-                        "nshow": selected_nshow,
-                        "pack_for_grid_register": f"{base_pack_dir}/{params_label}",
-                        "run_number": run_number + row_index,
-                    }
-                )
-                row_index += 1
-    return job_specs
-
-
 def _group_job_specs_by_label(job_specs):
     """Group job specs by apptainer image label."""
     grouped = {}
     for job_spec in job_specs:
-        label = job_spec["apptainer_label"]
+        label = job_spec["image_label"]
         grouped.setdefault(label, []).append(job_spec)
     return grouped
 
@@ -322,22 +147,27 @@ def _write_params_file(params_file_path, label_job_specs):
     """Write parameter file consumed by HTCondor queue-from syntax."""
     with open(params_file_path, "w", encoding="utf-8") as params_file_handle:
         for job_spec in label_job_specs:
-            array_layout_name = _resolve_array_layout_name(
-                job_spec["array_layout_name"], job_spec["model_version"]
-            )
-
             energy_min_value, energy_min_unit = _format_param_value(
                 job_spec["energy_min"], "energy_min_value"
             )
             energy_max_value, energy_max_unit = _format_param_value(
                 job_spec["energy_max"], "energy_max_value"
             )
+            core_scatter_number = _format_param_value(
+                job_spec["core_scatter_number"], "core_scatter_number"
+            )
             core_scatter_max_value, core_scatter_max_unit = _format_param_value(
                 job_spec["core_scatter_max"], "core_scatter_max_value"
             )
+            view_cone_min_value, view_cone_min_unit = _format_param_value(
+                job_spec["view_cone_min"], "view_cone_min_value"
+            )
+            view_cone_max_value, view_cone_max_unit = _format_param_value(
+                job_spec["view_cone_max"], "view_cone_max_value"
+            )
 
             row = [
-                _format_param_value(job_spec["apptainer_label"], "apptainer_label"),
+                _format_param_value(job_spec["image_label"], "apptainer_label"),
                 _format_param_value(job_spec["primary"], "primary"),
                 _format_param_value(job_spec["azimuth_angle"], "azimuth_angle"),
                 _format_param_value(job_spec["zenith_angle"], "zenith_angle"),
@@ -345,11 +175,16 @@ def _write_params_file(params_file_path, label_job_specs):
                 energy_min_unit,
                 energy_max_value,
                 energy_max_unit,
+                core_scatter_number,
                 core_scatter_max_value,
                 core_scatter_max_unit,
+                view_cone_min_value,
+                view_cone_min_unit,
+                view_cone_max_value,
+                view_cone_max_unit,
                 _format_param_value(job_spec["nshow"], "nshow"),
                 _format_param_value(job_spec["model_version"], "model_version"),
-                _format_param_value(array_layout_name, "array_layout_name"),
+                _format_param_value(job_spec["array_layout_name"], "array_layout_name"),
                 _format_param_value(job_spec["corsika_le_interaction"], "corsika_le_interaction"),
                 _format_param_value(job_spec["corsika_he_interaction"], "corsika_he_interaction"),
                 _format_param_value(job_spec["run_number"], "run_number"),
@@ -368,7 +203,7 @@ def generate_submission_script(args_dict):
         Arguments dictionary.
     """
     apptainer_images = _resolve_apptainer_images(args_dict["apptainer_image"])
-    job_specs = _build_job_specs(args_dict, list(apptainer_images.keys()))
+    job_specs, job_grid_metadata = build_job_specs(args_dict, list(apptainer_images.keys()))
     grouped_job_specs = _group_job_specs_by_label(job_specs)
 
     work_dir = Path(args_dict["output_path"])
@@ -377,14 +212,17 @@ def generate_submission_script(args_dict):
         if args_dict.get("htcondor_log_path")
         else work_dir / "htcondor_logs"
     )
-    log_dir = htcondor_log_path / "log"
-    error_dir = htcondor_log_path / "error"
-    output_dir = htcondor_log_path / "output"
+    htcondor_dirs = {
+        "log": htcondor_log_path / "log",
+        "error": htcondor_log_path / "error",
+        "output": htcondor_log_path / "output",
+    }
     work_dir.mkdir(parents=True, exist_ok=True)
-    for subdir in (log_dir, error_dir, output_dir):
+    for subdir in htcondor_dirs.values():
         subdir.mkdir(parents=True, exist_ok=True)
     submit_file_name = "simulate_prod.submit"
     _logger.info(f"Generating HT Condor submission scripts (path: {work_dir})")
+    submit_args = {**job_grid_metadata, **args_dict}
 
     for label, label_job_specs in grouped_job_specs.items():
         suffix = (
@@ -404,21 +242,17 @@ def generate_submission_script(args_dict):
                     apptainer_images[label],
                     args_dict["priority"],
                     params_file_name,
-                    log_dir=log_dir,
-                    error_dir=error_dir,
-                    output_dir=output_dir,
+                    htcondor_dirs=htcondor_dirs,
                 )
             )
 
     with open(work_dir / f"{submit_file_name}.sh", "w", encoding="utf-8") as submit_script_handle:
-        submit_script_handle.write(_get_submit_script(args_dict))
+        submit_script_handle.write(_get_submit_script(submit_args))
 
     Path(work_dir / f"{submit_file_name}.sh").chmod(0o755)
 
 
-def _get_submit_file(
-    executable, apptainer_image, priority, params_file_name, log_dir, error_dir, output_dir
-):
+def _get_submit_file(executable, apptainer_image, priority, params_file_name, htcondor_dirs):
     """
     Return HTCondor submit file.
 
@@ -434,12 +268,9 @@ def _get_submit_file(
         Priority of the job.
     params_file_name: str
         Name of the params file for queue-from submission.
-    log_dir: Path
-        Directory for HTCondor log files.
-    error_dir: Path
-        Directory for HTCondor error files.
-    output_dir: Path
-        Directory for HTCondor output files.
+    htcondor_dirs: dict
+        Directory mapping with HTCondor files locations. Expected keys are
+        ``log``, ``error``, and ``output``.
 
     Returns
     -------
@@ -454,9 +285,9 @@ def _get_submit_file(
 transfer_container = false
 
 executable = {executable}
-error      = {error_dir}/err.$(cluster)_$(process)
-output     = {output_dir}/out.$(cluster)_$(process)
-log        = {log_dir}/log.$(cluster)_$(process)
+error      = {htcondor_dirs["error"]}/err.$(cluster)_$(process)
+output     = {htcondor_dirs["output"]}/out.$(cluster)_$(process)
+log        = {htcondor_dirs["log"]}/log.$(cluster)_$(process)
 
 priority = {priority}
 arguments = "{arguments_string}"
@@ -486,41 +317,25 @@ def _get_submit_script(args_dict):
         idx = 3 + i
         bash_indices[field] = f"${{{idx}}}"
 
-    core_scatter = args_dict["core_scatter"]
-    n_core_scatter = core_scatter[0]
-    view_cone = args_dict["view_cone"]
-    view_cone_string = f'"{view_cone[0].to(u.deg)} {view_cone[1].to(u.deg)}"'
-
     label = args_dict["label"] if args_dict["label"] else "simulate-prod"
-    run_number_offset_arg = args_dict["run_number_offset"]
-    run_number_offset = 0 if run_number_offset_arg is None else run_number_offset_arg
-
-    azimuth_angle_idx = bash_indices["azimuth_angle"]
-    zenith_angle_idx = bash_indices["zenith_angle"]
-    energy_min_value_idx = bash_indices["energy_min_value"]
-    energy_min_unit_idx = bash_indices["energy_min_unit"]
-    energy_max_value_idx = bash_indices["energy_max_value"]
-    energy_max_unit_idx = bash_indices["energy_max_unit"]
-    core_scatter_max_value_idx = bash_indices["core_scatter_max_value"]
-    core_scatter_max_unit_idx = bash_indices["core_scatter_max_unit"]
-    nshow_idx = bash_indices["nshow"]
-    model_version_idx = bash_indices["model_version"]
-    array_layout_name_idx = bash_indices["array_layout_name"]
-    corsika_le_interaction_idx = bash_indices["corsika_le_interaction"]
-    corsika_he_interaction_idx = bash_indices["corsika_he_interaction"]
-    run_number_idx = bash_indices["run_number"]
-    pack_for_grid_register_idx = bash_indices["pack_for_grid_register"]
+    run_number_offset = args_dict.get("run_number_offset", 0)
 
     energy_range_string = (
-        f'"{energy_min_value_idx} {energy_min_unit_idx} '
-        f'{energy_max_value_idx} {energy_max_unit_idx}"'
+        f'"{bash_indices["energy_min_value"]} {bash_indices["energy_min_unit"]} '
+        f'{bash_indices["energy_max_value"]} {bash_indices["energy_max_unit"]}"'
     )
     core_scatter_string = (
-        f'"{n_core_scatter} {core_scatter_max_value_idx} {core_scatter_max_unit_idx}"'
+        f'"{bash_indices["core_scatter_number"]} {bash_indices["core_scatter_max_value"]} '
+        f'{bash_indices["core_scatter_max_unit"]}"'
+    )
+    view_cone_string = (
+        f'"{bash_indices["view_cone_min_value"]} {bash_indices["view_cone_min_unit"]} '
+        f"{bash_indices['view_cone_max_value']} "
+        f'{bash_indices["view_cone_max_unit"]}"'
     )
     energy_range_tag = (
-        f"erange-{energy_min_value_idx}{energy_min_unit_idx}-"
-        f"{energy_max_value_idx}{energy_max_unit_idx}"
+        f"erange-{bash_indices['energy_min_value']}{bash_indices['energy_min_unit']}-"
+        f"{bash_indices['energy_max_value']}{bash_indices['energy_max_unit']}"
     )
 
     return f"""#!/usr/bin/env bash
@@ -531,12 +346,12 @@ def _get_submit_script(args_dict):
 set -a; source "$2"
 apptainer_label="{bash_indices["apptainer_label"]}"
 primary="{bash_indices["primary"]}"
-model_version="{model_version_idx}"
-array_layout_name="{array_layout_name_idx}"
-corsika_le_interaction="{corsika_le_interaction_idx}"
-corsika_he_interaction="{corsika_he_interaction_idx}"
-run_number="{run_number_idx}"
-pack_for_grid_register="{pack_for_grid_register_idx}"
+model_version="{bash_indices["model_version"]}"
+array_layout_name="{bash_indices["array_layout_name"]}"
+corsika_le_interaction="{bash_indices["corsika_le_interaction"]}"
+corsika_he_interaction="{bash_indices["corsika_he_interaction"]}"
+run_number="{bash_indices["run_number"]}"
+pack_for_grid_register="{bash_indices["pack_for_grid_register"]}"
 energy_range_tag="{energy_range_tag}"
 job_label="{label}_${{corsika_he_interaction}}-${{corsika_le_interaction}}_${{energy_range_tag}}"
 
@@ -547,9 +362,9 @@ def _get_submit_script(args_dict):
     --site {args_dict["site"]} \\
     --array_layout_name "$array_layout_name" \\
     --primary "$primary" \\
-    --azimuth_angle "{azimuth_angle_idx}" \\
-    --zenith_angle "{zenith_angle_idx}" \\
-    --nshow "{nshow_idx}" \\
+    --azimuth_angle "{bash_indices["azimuth_angle"]}" \\
+    --zenith_angle "{bash_indices["zenith_angle"]}" \\
+    --nshow "{bash_indices["nshow"]}" \\
     --energy_range {energy_range_string} \\
     --core_scatter {core_scatter_string} \\
     --view_cone {view_cone_string} \\
@@ -562,3 +377,34 @@ def _get_submit_script(args_dict):
     --log_level {args_dict["log_level"]} \\
     --pack_for_grid_register "$pack_for_grid_register"
 """
+
+
+def build_job_specs(args_dict, image_labels):
+    """Build backend-agnostic job specs from comparison and production grids."""
+    base_pack_dir = args_dict.get("simulation_output") or "simtools-output"
+    normalized_rows, job_grid_metadata = read_job_grid(args_dict["job_grid_file"])
+
+    missing_metadata = [
+        key
+        for key in _REQUIRED_JOB_GRID_METADATA
+        if key not in job_grid_metadata or job_grid_metadata.get(key) in (None, "")
+    ]
+    if missing_metadata:
+        missing_keys = ", ".join(missing_metadata)
+        raise ValueError(
+            "Job grid metadata is missing required field(s): "
+            f"{missing_keys}. Regenerate the job grid with "
+            "simtools-production-generate-grid so metadata includes these values."
+        )
+
+    job_specs = []
+    for label in image_labels:
+        for row in normalized_rows:
+            job_specs.append(
+                {
+                    "image_label": str(label),
+                    **row,
+                    "pack_for_grid_register": f"{base_pack_dir}/{label!s}",
+                }
+            )
+    return job_specs, job_grid_metadata
diff --git a/src/simtools/layout/array_layout_utils.py b/src/simtools/layout/array_layout_utils.py
index e09104774b..b934f4d90d 100644
--- a/src/simtools/layout/array_layout_utils.py
+++ b/src/simtools/layout/array_layout_utils.py
@@ -1,5 +1,6 @@
 """Retrieve, merge, and write layout dictionaries."""
 
+import ast
 import logging
 from pathlib import Path
 
@@ -7,6 +8,7 @@
 from astropy.table import QTable, Table
 
 import simtools.utils.general as gen
+import simtools.version as simtools_version
 from simtools.data_model import data_reader
 from simtools.data_model.model_data_writer import ModelDataWriter
 from simtools.io import ascii_handler, io_handler
@@ -564,30 +566,6 @@ def read_layouts(args_dict):
     return [], background_layout
 
 
-def _get_array_name(array_name):
-    """
-    Return telescope size and number of telescopes from regular array name.
-
-    Finetuned to array names like "4MST", "1LST", etc.
-
-    Parameters
-    ----------
-    array_name : str
-        Name of the regular array (e.g. "4MST").
-
-    Returns
-    -------
-    tel_size : str
-        Telescope size (e.g. "MST").
-    n_tel : int
-        Number of telescopes (e.g. 4).
-    """
-    if len(array_name) < 2 or not array_name[0].isdigit():
-        raise ValueError(f"Invalid array_name: '{array_name}'")
-
-    return array_name[1:], int(array_name[0])
-
-
 def _create_star_array(tel_name, pos_x, pos_y, pos_z, n_telescopes, tel_type, site, distance):
     """Create star-shaped array positions along x and y axes."""
     axis_sequence = ["x", "y", "-x", "-y"]
@@ -818,3 +796,35 @@ def write_array_elements_info_yaml(
         }
 
     ascii_handler.write_data_to_file(data, output_file)
+
+
+def resolve_array_layout_name(array_layout_name, model_version):
+    """
+    Resolve array layout configuration for a specific model version.
+
+    Parameters
+    ----------
+    array_layout_name : str or dict or list of str or list of dict
+        Array layout name(s) or configuration(s) to resolve.
+    model_version : str
+        Model version to resolve for.
+
+    """
+    if isinstance(array_layout_name, list) and len(array_layout_name) == 1:
+        array_layout_name = array_layout_name[0]
+
+    if isinstance(array_layout_name, str) and array_layout_name.strip().startswith("{"):
+        try:
+            parsed_layout = ast.literal_eval(array_layout_name)
+            if isinstance(parsed_layout, dict):
+                array_layout_name = parsed_layout
+        except (SyntaxError, ValueError):
+            return array_layout_name
+
+    if not isinstance(array_layout_name, dict) or list(array_layout_name) != ["by_version"]:
+        return array_layout_name
+
+    resolved = simtools_version.resolve_by_version(
+        {"array_layout_name": array_layout_name}, model_version
+    )
+    return resolved["array_layout_name"]
diff --git a/src/simtools/production_configuration/corsika_limits_lookup.py b/src/simtools/production_configuration/corsika_limits_lookup.py
new file mode 100644
index 0000000000..7305a1245a
--- /dev/null
+++ b/src/simtools/production_configuration/corsika_limits_lookup.py
@@ -0,0 +1,228 @@
+"""Lookup-table access and interpolation for CORSIKA production limits."""
+
+import numpy as np
+from astropy import units as u
+from astropy.table import Table
+from scipy.interpolate import LinearNDInterpolator, griddata
+from scipy.spatial import QhullError  # pylint: disable=no-name-in-module
+
+from simtools.utils.value_conversion import get_value_in_unit
+
+_LOOKUP_FIELDS = (
+    "lower_energy_threshold",
+    "upper_scatter_radius",
+    "viewcone_radius",
+)
+
+_LOOKUP_FIELD_UNITS = {
+    "lower_energy_threshold": "TeV",
+    "upper_scatter_radius": "m",
+    "viewcone_radius": "deg",
+}
+
+_POINT_LIMIT_KEYS = {
+    "lower_energy_threshold": "lower_energy_threshold",
+    "upper_scatter_radius": "scatter_radius",
+    "viewcone_radius": "viewcone_radius",
+}
+
+
+def attach_lookup_limits_to_point(point, limits):
+    """Attach interpolated CORSIKA limits to a grid point."""
+    for lookup_key, point_key in _POINT_LIMIT_KEYS.items():
+        point[point_key] = limits[lookup_key] * u.Unit(_LOOKUP_FIELD_UNITS[lookup_key])
+
+
+class CorsikaLimitsLookup:
+    """Read and interpolate CORSIKA limits for production grids."""
+
+    def __init__(self, lookup_table, array_layout_name=None):
+        """
+        Initialize lookup-table access.
+
+        Parameters
+        ----------
+        lookup_table : str or Path
+            Path to the lookup-table ECSV file.
+        array_layout_name : str, optional
+            Array layout name used to select lookup-table rows.
+        """
+        self.lookup_table = lookup_table
+        self.array_layout_name = array_layout_name
+        self.lookup_points_for_interpolation = None
+        self.lookup_values_for_interpolation = None
+        self.lookup_interpolators_for_point = None
+
+    def load_matching_lookup_arrays(self):
+        """
+        Load and filter lookup-table arrays for the selected array layout.
+
+        Returns
+        -------
+        dict
+            Lookup arrays for interpolation.
+        """
+        lookup_table = Table.read(self.lookup_table, format="ascii.ecsv")
+        if self.array_layout_name is None:
+            matching_rows = list(lookup_table)
+        else:
+            matching_rows = [
+                row for row in lookup_table if str(row["array_name"]) == str(self.array_layout_name)
+            ]
+
+        if not matching_rows:
+            raise ValueError(
+                "No matching rows in the lookup table for "
+                f"array_layout_name: {self.array_layout_name}"
+            )
+
+        def extract_array(field, transform=lambda x: x):
+            return np.array([transform(row[field]) for row in matching_rows])
+
+        zeniths = extract_array("zenith")
+        azimuths = extract_array("azimuth", lambda x: x % 360)
+        nsb_values = extract_array("nsb_level", float)
+
+        return {
+            "points": np.column_stack((zeniths, azimuths, nsb_values)),
+            "lower_energy_threshold": extract_array("lower_energy_limit"),
+            "upper_scatter_radius": extract_array("upper_radius_limit"),
+            "viewcone_radius": extract_array("viewcone_radius"),
+        }
+
+    def prepare_point_interpolators(self):
+        """
+        Prepare lookup arrays for per-point interpolation.
+
+        Returns
+        -------
+        dict
+            Interpolators keyed by lookup quantity.
+        """
+        lookup_arrays = self.load_matching_lookup_arrays()
+        self.lookup_points_for_interpolation = self._build_wrapped_interpolation_points(
+            lookup_arrays["points"]
+        )
+        self.lookup_values_for_interpolation = {
+            key: self._repeat_wrapped_lookup_values(lookup_arrays[key]) for key in _LOOKUP_FIELDS
+        }
+        try:
+            self.lookup_interpolators_for_point = {
+                key: LinearNDInterpolator(
+                    self.lookup_points_for_interpolation,
+                    self.lookup_values_for_interpolation[key],
+                    fill_value=np.nan,
+                )
+                for key in _LOOKUP_FIELDS
+            }
+        except QhullError as exc:
+            raise ValueError(
+                "Lookup table does not contain enough unique points for 3D interpolation "
+                "(zenith, azimuth, nsb_level). Provide a denser lookup table "
+                "or run without --corsika_limits if limits are not required for this use case."
+            ) from exc
+        return self.lookup_interpolators_for_point
+
+    def interpolate_grid_limits(self, target_values):
+        """
+        Interpolate lookup values on a regular zenith/azimuth/NSB grid.
+
+        Parameters
+        ----------
+        target_values : dict
+            Generated target-axis values.
+
+        Returns
+        -------
+        dict
+            Interpolated lookup values on the target grid.
+        """
+        lookup_arrays = self.load_matching_lookup_arrays()
+        points = self._build_wrapped_interpolation_points(lookup_arrays["points"])
+
+        target_grid = (
+            np.array(
+                np.meshgrid(
+                    target_values["zenith_angle"].value,
+                    target_values["azimuth"].value,
+                    target_values["nsb_level"].value,
+                    indexing="ij",
+                )
+            )
+            .reshape(3, -1)
+            .T
+        )
+
+        def interpolate(values):
+            return griddata(
+                points,
+                self._repeat_wrapped_lookup_values(values),
+                target_grid,
+                method="linear",
+                fill_value=np.nan,
+            ).reshape(
+                len(target_values["zenith_angle"]),
+                len(target_values["azimuth"]),
+                len(target_values["nsb_level"]),
+            )
+
+        return {key: interpolate(lookup_arrays[key]) for key in _LOOKUP_FIELDS}
+
+    def interpolate_point(self, zenith, azimuth, nsb=1.0):
+        """
+        Interpolate lookup-table limits for a single point.
+
+        Parameters
+        ----------
+        zenith : float or Quantity
+            Zenith angle.
+        azimuth : float or Quantity
+            Azimuth angle.
+        nsb : float or Quantity, optional
+            Night-sky background level.
+
+        Returns
+        -------
+        dict
+            Interpolated lower-energy threshold, scatter radius, and viewcone radius.
+        """
+        if self.lookup_interpolators_for_point is None:
+            self.prepare_point_interpolators()
+
+        target = np.array(
+            [
+                [
+                    get_value_in_unit(zenith, "deg"),
+                    get_value_in_unit(azimuth, "deg") % 360.0,
+                    get_value_in_unit(nsb),
+                ]
+            ],
+            dtype=float,
+        )
+        return {
+            "lower_energy_threshold": float(
+                self.lookup_interpolators_for_point["lower_energy_threshold"](target)[0]
+            ),
+            "upper_scatter_radius": float(
+                self.lookup_interpolators_for_point["upper_scatter_radius"](target)[0]
+            ),
+            "viewcone_radius": float(
+                self.lookup_interpolators_for_point["viewcone_radius"](target)[0]
+            ),
+        }
+
+    @staticmethod
+    def _repeat_wrapped_lookup_values(values):
+        """Repeat lookup values for azimuth wrapping."""
+        return np.tile(values, 3)
+
+    def _build_wrapped_interpolation_points(self, points):
+        """Return lookup points extended across azimuth wrap boundaries."""
+        azimuths_wrapped = np.concatenate([points[:, 1] + shift for shift in (0, 360, -360)])
+        return np.column_stack(
+            (
+                self._repeat_wrapped_lookup_values(points[:, 0]),
+                azimuths_wrapped,
+                self._repeat_wrapped_lookup_values(points[:, 2]),
+            )
+        )
diff --git a/src/simtools/production_configuration/generate_production_grid.py b/src/simtools/production_configuration/generate_production_grid.py
deleted file mode 100644
index 0d747ed1b7..0000000000
--- a/src/simtools/production_configuration/generate_production_grid.py
+++ /dev/null
@@ -1,729 +0,0 @@
-"""
-Module defines the `GridGeneration` class.
-
-Used to generate a grid of simulation points based on flexible axes definitions such
-azimuth, zenith angle, night-sky background, and camera offset.
-The module handles axis binning, scaling and interpolation of energy thresholds, viewcone,
-and radius limits from a lookup table.
-Additionally, it allows for converting between Altitude/Azimuth and Right Ascension
-Declination coordinates. The resulting grid points are saved to a file.
-"""
-
-import json
-import logging
-from pathlib import Path
-
-import numpy as np
-from astropy import units as u
-from astropy.coordinates import AltAz, EarthLocation, SkyCoord
-from astropy.table import Table
-from astropy.units import Quantity
-from scipy.interpolate import LinearNDInterpolator, griddata
-from scipy.spatial import QhullError  # pylint: disable=no-name-in-module
-
-from simtools.simtel.simtel_io_metadata import (
-    get_sim_telarray_telescope_id_to_telescope_name_mapping,
-)
-
-DEFAULT_SERIALIZATION_ROUND_DECIMALS = 6
-
-
-class GridGeneration:
-    """
-    Defines and generates a grid of simulation points based on flexible axes definitions.
-
-    This class generates a grid of points for a simulation based on parameters such as
-    azimuth, zenith angle, night-sky background, and camera offset,
-    taking into account axis definitions, scaling, and units and interpolating values
-    for simulations from a lookup table.
-    """
-
-    def __init__(
-        self,
-        axes,
-        coordinate_system="zenith_azimuth",
-        observing_location=None,
-        observing_time=None,
-        lookup_table=None,
-        telescope_ids=None,
-        simtel_file=None,
-    ):
-        """
-        Initialize the grid with the given axes and coordinate system.
-
-        Parameters
-        ----------
-        axes : dict
-            Dictionary where each key is the axis name and the value is a dictionary
-            defining the axis properties (range, binning, scaling, etc.).
-        coordinate_system : str, optional
-            The coordinate system for the grid generation (default is 'zenith_azimuth').
-        observing_location : EarthLocation, optional
-            The location of the observation (latitude, longitude, height).
-        observing_time : Time, optional
-            The time of the observation. If None, coordinate conversion to RA/Dec not working.
-        lookup_table : str, optional
-            Path to the lookup table file (ECSV format).
-        telescope_ids : list of str, optional
-            List of telescope IDs to get the limits for.
-        simtel_file : str, optional
-            Path to a sim_telarray file used to map sim_telarray telescope IDs to
-            telescope names when matching lookup-table telescope selections.
-        """
-        self._logger = logging.getLogger(__name__)
-        self.axes = axes["axes"] if "axes" in axes else axes
-        self.coordinate_system = coordinate_system
-        self.observing_location = (
-            observing_location
-            if observing_location is not None
-            else EarthLocation(lat=0.0 * u.deg, lon=0.0 * u.deg, height=0 * u.m)
-        )
-        self.observing_time = observing_time
-        self.lookup_table = lookup_table
-        self.telescope_ids = telescope_ids
-        self.simtel_file = simtel_file
-        self.interpolated_limits = {}
-        self.serialization_round_decimals = DEFAULT_SERIALIZATION_ROUND_DECIMALS
-        self._simtel_id_to_name = (
-            get_sim_telarray_telescope_id_to_telescope_name_mapping(simtel_file)
-            if simtel_file
-            else {}
-        )
-
-        # Store target values for each axis
-        self.target_values = self._generate_target_values()
-
-        if self.lookup_table:
-            if self.coordinate_system == "ra_dec":
-                self._prepare_lookup_table_limits_for_point_interpolation()
-            else:
-                self._apply_lookup_table_limits()
-
-    @staticmethod
-    def _coerce_identifier_container(value):
-        """Coerce identifier input into a list."""
-        if value is None:
-            return []
-        if isinstance(value, str):
-            stripped = value.strip()
-            return json.loads(stripped) if stripped.startswith("[") else [stripped]
-        if isinstance(value, (list, tuple, set)):
-            return list(value)
-        return [value]
-
-    def _normalize_lookup_identifier(self, identifier):
-        """Normalize one telescope identifier and report if it is numeric."""
-        if isinstance(identifier, (int, np.integer)):
-            return self._simtel_id_to_name.get(int(identifier), str(int(identifier))), True
-
-        text = str(identifier).strip()
-        if text.lstrip("+-").isdigit():
-            return self._simtel_id_to_name.get(int(text), text), True
-        return text, False
-
-    def _normalized_identifier_set(self, identifiers):
-        """Return normalized telescope identifiers as a set."""
-        return {
-            self._normalize_lookup_identifier(identifier)[0]
-            for identifier in self._coerce_identifier_container(identifiers)
-        }
-
-    def _lookup_contains_numeric_telescope_ids(self, lookup_table):
-        """Return True when any lookup-table telescope identifier is numeric."""
-        return any(
-            any(
-                self._normalize_lookup_identifier(identifier)[1]
-                for identifier in self._coerce_identifier_container(row["telescope_ids"])
-            )
-            for row in lookup_table
-        )
-
-    def _load_matching_lookup_arrays(self):
-        """Load and filter lookup-table arrays for selected telescope IDs."""
-        lookup_table = Table.read(self.lookup_table, format="ascii.ecsv")
-        selected_telescope_ids = self._normalized_identifier_set(self.telescope_ids)
-
-        matching_rows = [
-            row
-            for row in lookup_table
-            if selected_telescope_ids == self._normalized_identifier_set(row["telescope_ids"])
-        ]
-
-        if not matching_rows:
-            has_numeric_lookup_ids = self._lookup_contains_numeric_telescope_ids(lookup_table)
-
-            if has_numeric_lookup_ids and not self.simtel_file:
-                raise ValueError(
-                    "Lookup table telescope selections contain numeric IDs. "
-                    "Provide --simtel_file to map those IDs to telescope names."
-                )
-
-            raise ValueError(
-                f"No matching rows in the lookup table for telescope_ids: {self.telescope_ids}"
-            )
-
-        def extract_array(field, transform=lambda x: x):
-            return np.array([transform(row[field]) for row in matching_rows])
-
-        zeniths = extract_array("zenith")
-        azimuths = extract_array("azimuth", lambda x: x % 360)
-        nsb_values = extract_array("nsb_level", float)
-
-        return {
-            "points": np.column_stack((zeniths, azimuths, nsb_values)),
-            "lower_energy_threshold": extract_array("lower_energy_limit"),
-            "upper_scatter_radius": extract_array("upper_radius_limit"),
-            "viewcone_radius": extract_array("viewcone_radius"),
-        }
-
-    def _require_observing_time(self):
-        """Return observing time if available, else raise a clear error."""
-        if self.observing_time is None:
-            raise ValueError("Observing time is required for ra_dec grid generation.")
-        return self.observing_time
-
-    def _get_max_zenith_for_radec_mode(self):
-        """Read maximum zenith from axes for RA/Dec direction sampling."""
-        zenith_axis = self.axes.get("zenith_angle")
-        if not zenith_axis or "range" not in zenith_axis or len(zenith_axis["range"]) != 2:
-            raise ValueError(
-                "RA/Dec direction sampling requires 'zenith_angle' axis with a valid "
-                "two-element 'range' in the axes definition."
-            )
-        return float(zenith_axis["range"][1])
-
-    def _prepare_lookup_table_limits_for_point_interpolation(self):
-        """Prepare lookup arrays for per-point interpolation in RA/Dec grid mode."""
-        lookup_arrays = self._load_matching_lookup_arrays()
-
-        points = lookup_arrays["points"]
-        azimuths = points[:, 1]
-        azimuths_wrapped = np.concatenate([azimuths + shift for shift in (0, 360, -360)])
-
-        def repeat_3(arr):
-            """Repeat an array three times."""
-            return np.tile(arr, 3)
-
-        self.lookup_points_for_interpolation = np.column_stack(
-            (
-                repeat_3(points[:, 0]),
-                azimuths_wrapped,
-                repeat_3(points[:, 2]),
-            )
-        )
-        self.lookup_values_for_interpolation = {
-            "lower_energy_threshold": repeat_3(lookup_arrays["lower_energy_threshold"]),
-            "upper_scatter_radius": repeat_3(lookup_arrays["upper_scatter_radius"]),
-            "viewcone_radius": repeat_3(lookup_arrays["viewcone_radius"]),
-        }
-        try:
-            self.lookup_interpolators_for_point = {
-                key: LinearNDInterpolator(
-                    self.lookup_points_for_interpolation,
-                    self.lookup_values_for_interpolation[key],
-                    fill_value=np.nan,
-                )
-                for key in ("lower_energy_threshold", "upper_scatter_radius", "viewcone_radius")
-            }
-        except QhullError as exc:
-            raise ValueError(
-                "Lookup table does not contain enough unique points for 3D interpolation "
-                "(zenith, azimuth, nsb_level). Provide a denser lookup table "
-                "or run without --lookup_table if limits are not required for this use case."
-            ) from exc
-
-    def _has_radec_axes(self):
-        """Return True if axes define a native RA/Dec grid."""
-        return "ra" in self.axes and "dec" in self.axes
-
-    def _generate_target_values(self):
-        """
-        Generate target axis values and store them as Quantities.
-
-        Returns
-        -------
-        dict
-            Dictionary of target values for each axis, stored as Quantity objects.
-        """
-        target_values = {}
-        for axis_name, axis in self.axes.items():
-            axis_range = axis["range"]
-            binning = axis["binning"]
-            scaling = axis.get("scaling", "linear")
-            units = axis.get("units", None)
-
-            if axis_name == "azimuth":
-                # Use circular binning for azimuth
-                values = self.create_circular_binning(axis_range, binning)
-            elif scaling == "log":
-                # Log scaling
-                values = np.logspace(np.log10(axis_range[0]), np.log10(axis_range[1]), binning)
-            elif scaling == "1/cos":
-                # 1/cos scaling
-                cos_min = np.cos(np.radians(axis_range[0]))
-                cos_max = np.cos(np.radians(axis_range[1]))
-                inv_cos_values = np.linspace(1 / cos_min, 1 / cos_max, binning)
-                values = np.degrees(np.arccos(1 / inv_cos_values))
-            else:
-                # Linear scaling
-                values = np.linspace(axis_range[0], axis_range[1], binning)
-
-            if units:
-                values = values * u.Unit(units)
-
-            target_values[axis_name] = values
-
-        return target_values
-
-    def _apply_lookup_table_limits(self):
-        """Apply limits from the lookup table and interpolate values."""
-        lookup_arrays = self._load_matching_lookup_arrays()
-        points_base = lookup_arrays["points"]
-        lower_energy_thresholds = lookup_arrays["lower_energy_threshold"]
-        upper_scatter_radii = lookup_arrays["upper_scatter_radius"]
-        viewcone_radii = lookup_arrays["viewcone_radius"]
-
-        zeniths = points_base[:, 0]
-        azimuths = points_base[:, 1]
-        nsb_values = points_base[:, 2]
-
-        # Wrap azimuths and repeat others
-        azimuths_wrapped = np.concatenate([azimuths + shift for shift in (0, 360, -360)])
-
-        def repeat_3(arr):
-            """Repeat an array three times."""
-            return np.tile(arr, 3)
-
-        points = np.column_stack(
-            (
-                repeat_3(zeniths),
-                azimuths_wrapped,
-                repeat_3(nsb_values),
-            )
-        )
-
-        target_grid = (
-            np.array(
-                np.meshgrid(
-                    self.target_values["zenith_angle"].value,
-                    self.target_values["azimuth"].value,
-                    self.target_values["nsb_level"].value,
-                    indexing="ij",
-                )
-            )
-            .reshape(3, -1)
-            .T
-        )
-
-        def interpolate(values):
-            return griddata(
-                points, repeat_3(values), target_grid, method="linear", fill_value=np.nan
-            ).reshape(
-                len(self.target_values["zenith_angle"]),
-                len(self.target_values["azimuth"]),
-                len(self.target_values["nsb_level"]),
-            )
-
-        self.interpolated_limits = {
-            "lower_energy_threshold": interpolate(lower_energy_thresholds),
-            "upper_scatter_radius": interpolate(upper_scatter_radii),
-            "viewcone_radius": interpolate(viewcone_radii),
-        }
-
-    def _generate_radec_grid_direction_points(self):
-        """Generate direction points from declination lines and hour-angle spacing."""
-        observing_time = self._require_observing_time()
-        max_zenith = self._get_max_zenith_for_radec_mode()
-        lst_deg = observing_time.sidereal_time(
-            "apparent", longitude=self.observing_location.lon
-        ).deg
-
-        direction_points = []
-        for declination in np.arange(-90.0, 91.0, 1.0):
-            cos_dec = np.cos(np.deg2rad(declination))
-            step_ha = 1.0 / cos_dec if cos_dec > 1e-6 else 360.0
-            n_ha = max(1, int(np.ceil(360.0 / step_ha)))
-            hour_angles = np.linspace(-180.0, 180.0, n_ha, endpoint=False)
-            ra_values = (lst_deg - hour_angles) % 360.0
-
-            skycoord = SkyCoord(
-                ra=ra_values * u.deg,
-                dec=np.full_like(ra_values, declination) * u.deg,
-                frame="icrs",
-            )
-            altaz = skycoord.transform_to(
-                AltAz(location=self.observing_location, obstime=observing_time)
-            )
-
-            zenith_values = (90.0 * u.deg - altaz.alt).to(u.deg).value
-            mask = (zenith_values >= 0.0) & (zenith_values <= max_zenith)
-
-            for idx in np.nonzero(mask)[0]:
-                direction_points.append(
-                    {
-                        "zenith_angle": zenith_values[idx] * u.deg,
-                        "azimuth": altaz.az.deg[idx] * u.deg,
-                    }
-                )
-        return direction_points
-
-    def _generate_extra_axis_combinations(self, excluded_keys):
-        """Generate combinations for all axes except the excluded ones."""
-        extra_axes = {
-            key: value for key, value in self.target_values.items() if key not in excluded_keys
-        }
-        if not extra_axes:
-            return list(extra_axes.keys()), [], [np.array([])]
-
-        extra_value_arrays = [value.value for value in extra_axes.values()]
-        extra_units = [value.unit for value in extra_axes.values()]
-        extra_grid = np.meshgrid(*extra_value_arrays, indexing="ij")
-        extra_combinations = np.vstack(list(map(np.ravel, extra_grid))).T
-        return list(extra_axes.keys()), extra_units, extra_combinations
-
-    def _add_lookup_limits_to_point(self, point, zenith, azimuth):
-        """Interpolate and attach lookup-table limits to a grid point."""
-        if not self.lookup_table:
-            return
-
-        nsb_value = point.get("nsb_level", 1)
-        if isinstance(nsb_value, Quantity):
-            nsb_value = nsb_value.value
-        limits = self._interpolate_limits_for_point(
-            zenith=zenith,
-            azimuth=azimuth,
-            nsb=float(nsb_value),
-        )
-        point["lower_energy_threshold"] = limits["lower_energy_threshold"] * u.TeV
-        point["scatter_radius"] = limits["upper_scatter_radius"] * u.m
-        point["viewcone_radius"] = limits["viewcone_radius"] * u.deg
-
-    def _generate_grid_from_radec_axes(self):
-        """Generate grid points from explicit RA/Dec axes definitions.
-
-        All explicit RA/Dec combinations defined by the input axes are preserved,
-        even when their transformed Alt/Az coordinates fall below the local horizon.
-        """
-        observing_time = self._require_observing_time()
-
-        axis_keys = [key for key in self.target_values if key not in ("zenith_angle", "azimuth")]
-        value_arrays = [self.target_values[key].value for key in axis_keys]
-        units = [self.target_values[key].unit for key in axis_keys]
-        grid = np.meshgrid(*value_arrays, indexing="ij")
-        combinations = np.vstack(list(map(np.ravel, grid))).T
-
-        grid_points = []
-        for combination in combinations:
-            grid_point = {
-                key: Quantity(combination[i], units[i]) for i, key in enumerate(axis_keys)
-            }
-
-            skycoord = SkyCoord(
-                ra=grid_point["ra"].to(u.deg),
-                dec=grid_point["dec"].to(u.deg),
-                frame="icrs",
-            )
-            altaz = skycoord.transform_to(
-                AltAz(location=self.observing_location, obstime=observing_time)
-            )
-            zenith = (90.0 * u.deg - altaz.alt).to(u.deg).value
-
-            self._add_lookup_limits_to_point(grid_point, zenith=zenith, azimuth=altaz.az.deg)
-
-            grid_points.append(grid_point)
-
-        return grid_points
-
-    def _interpolate_limits_for_point(self, zenith, azimuth, nsb):
-        """Interpolate lookup-table limits for a single point."""
-        target = np.array([[zenith, azimuth % 360.0, nsb]], dtype=float)
-        return {
-            "lower_energy_threshold": float(
-                self.lookup_interpolators_for_point["lower_energy_threshold"](target)[0]
-            ),
-            "upper_scatter_radius": float(
-                self.lookup_interpolators_for_point["upper_scatter_radius"](target)[0]
-            ),
-            "viewcone_radius": float(
-                self.lookup_interpolators_for_point["viewcone_radius"](target)[0]
-            ),
-        }
-
-    def _generate_grid_radec_mode(self):
-        """Generate grid points for RA/Dec mode.
-
-        If explicit ``ra``/``dec`` axes are provided, use those points directly.
-        Otherwise, sample directions along declination lines with hour-angle spacing
-        and keep only points that satisfy the configured zenith-angle limits.
-        """
-        if self._has_radec_axes():
-            return self._generate_grid_from_radec_axes()
-
-        direction_points = self._generate_radec_grid_direction_points()
-        extra_keys, extra_units, extra_combinations = self._generate_extra_axis_combinations(
-            excluded_keys=("zenith_angle", "azimuth")
-        )
-
-        grid_points = []
-        for direction_point in direction_points:
-            for extra_combination in extra_combinations:
-                point = dict(direction_point)
-                for i, key in enumerate(extra_keys):
-                    point[key] = Quantity(extra_combination[i], extra_units[i])
-
-                self._add_lookup_limits_to_point(
-                    point,
-                    zenith=point["zenith_angle"].value,
-                    azimuth=point["azimuth"].value,
-                )
-
-                grid_points.append(point)
-
-        return grid_points
-
-    def create_circular_binning(self, azimuth_range, num_bins):
-        """
-        Create bin centers for azimuth angles, handling circular wraparound (0 deg to 360 deg).
-
-        Parameters
-        ----------
-        azimuth_range : tuple
-            (min_azimuth, max_azimuth), can wrap around 0 deg.
-        num_bins : int
-            Number of bins.
-
-        Returns
-        -------
-        np.ndarray
-            Array of bin centers.
-        """
-        azimuth_min, azimuth_max = azimuth_range
-        azimuth_min %= 360  # Normalize to [0, 360)
-        azimuth_max %= 360
-
-        clockwise_distance = (azimuth_max - azimuth_min) % 360
-        counterclockwise_distance = (azimuth_min - azimuth_max) % 360
-
-        if clockwise_distance <= counterclockwise_distance:
-            bin_centers = (
-                np.linspace(azimuth_min, azimuth_min + clockwise_distance, num_bins, endpoint=True)
-                % 360
-            )
-        else:
-            bin_centers = (
-                np.linspace(
-                    azimuth_min, azimuth_min - counterclockwise_distance, num_bins, endpoint=True
-                )
-                % 360
-            )
-
-        return bin_centers
-
-    def generate_grid(self):
-        """
-        Generate the grid based on the required axes and include interpolated limits.
-
-        Takes energy threshold, viewcone, and radius from the interpolated lookup table.
-
-        Returns
-        -------
-        list of dict
-            A list of grid points, each represented as a dictionary with axis names
-            as keys and axis values as values. Axis values may include units where defined.
-        """
-        if self.coordinate_system == "ra_dec":
-            return self._generate_grid_radec_mode()
-
-        value_arrays = [value.value for value in self.target_values.values()]
-        units = [value.unit for value in self.target_values.values()]
-        grid = np.meshgrid(*value_arrays, indexing="ij")
-        combinations = np.vstack(list(map(np.ravel, grid))).T
-        grid_points = []
-        for combination in combinations:
-            grid_point = {
-                key: Quantity(combination[i], units[i])
-                for i, key in enumerate(self.target_values.keys())
-            }
-
-            if "lower_energy_threshold" in self.interpolated_limits:
-                zenith_idx = np.searchsorted(
-                    self.target_values["zenith_angle"].value, grid_point["zenith_angle"].value
-                )
-                azimuth_idx = np.searchsorted(
-                    self.target_values["azimuth"].value, grid_point["azimuth"].value
-                )
-                nsb_idx = np.searchsorted(
-                    self.target_values["nsb_level"].value,
-                    grid_point["nsb_level"].value,
-                )
-                energy_lower = self.interpolated_limits["lower_energy_threshold"][
-                    zenith_idx, azimuth_idx, nsb_idx
-                ]
-                grid_point["lower_energy_threshold"] = energy_lower * u.TeV
-
-            if "upper_scatter_radius" in self.interpolated_limits:
-                radius_value = self.interpolated_limits["upper_scatter_radius"][
-                    zenith_idx, azimuth_idx, nsb_idx
-                ]
-                grid_point["scatter_radius"] = radius_value * u.m
-
-            if "viewcone_radius" in self.interpolated_limits:
-                viewcone_value = self.interpolated_limits["viewcone_radius"][
-                    zenith_idx, azimuth_idx, nsb_idx
-                ]
-                grid_point["viewcone_radius"] = viewcone_value * u.deg
-
-            grid_points.append(grid_point)
-
-        return grid_points
-
-    def convert_altaz_to_radec(self, alt, az):
-        """
-        Convert Altitude/Azimuth (AltAz) coordinates to Right Ascension/Declination (RA/Dec).
-
-        Parameters
-        ----------
-        alt : float
-            Altitude angle in degrees.
-        az : float
-            Azimuth angle in degrees.
-
-        Returns
-        -------
-        SkyCoord
-            SkyCoord object containing the RA/Dec coordinates.
-
-        Raises
-        ------
-        ValueError
-            If observing_time is not set.
-        """
-        if self.observing_time is None:
-            raise ValueError(
-                "Observing time is not set. "
-                "Please provide an observing_time to convert coordinates."
-            )
-
-        alt_rad = alt.to(u.rad)
-        az_rad = az.to(u.rad)
-        aa = AltAz(
-            alt=alt_rad,
-            az=az_rad,
-            location=self.observing_location,
-            obstime=self.observing_time,
-        )
-        skycoord = SkyCoord(aa)
-        return skycoord.icrs  # Return RA/Dec in ICRS frame
-
-    def convert_coordinates(self, grid_points):
-        """
-        Convert the grid points to RA/Dec coordinates if necessary.
-
-        Parameters
-        ----------
-        grid_points : list of dict
-            List of grid points, each represented as a dictionary with axis
-              names as keys and values.
-
-        Returns
-        -------
-        list of dict
-            The grid points with converted RA/Dec coordinates.
-        """
-        if self.coordinate_system == "ra_dec":
-            for point in grid_points:
-                if "zenith_angle" in point and "azimuth" in point:
-                    alt = (90.0 * u.deg) - point.pop("zenith_angle")
-                    az = point.pop("azimuth")
-                    radec = self.convert_altaz_to_radec(alt, az)
-                    point["ra"] = radec.ra.deg * u.deg
-                    point["dec"] = radec.dec.deg * u.deg
-        return grid_points
-
-    def serialize_grid_points(self, grid_points, output_file):
-        """Serialize the grid output and save to an ECSV table file."""
-        if Path(output_file).suffix.lower() != ".ecsv":
-            raise ValueError("Grid output file must use '.ecsv' extension.")
-
-        all_keys = self._collect_point_keys(grid_points)
-        rows, units = self._build_serialized_rows(grid_points, all_keys)
-
-        output_table = Table(rows=rows, names=all_keys)
-        for column_name, unit in units.items():
-            output_table[column_name].unit = u.Unit(unit)
-
-        output_table.meta = self._build_grid_metadata()
-
-        output_table.write(output_file, format="ascii.ecsv", overwrite=True)
-        self._logger.info(f"Output saved to {output_file}")
-
-    @staticmethod
-    def _collect_point_keys(grid_points):
-        """Collect all grid-point keys while preserving first-seen order."""
-        all_keys = []
-        for point in grid_points:
-            for key in point:
-                if key not in all_keys:
-                    all_keys.append(key)
-        return all_keys
-
-    def _serialize_grid_value(self, value):
-        """Serialize one grid value and return (value, unit)."""
-        if isinstance(value, u.Quantity):
-            serialized = round(float(value.value), self.serialization_round_decimals)
-            return serialized, str(value.unit)
-
-        if isinstance(value, dict) and "value" in value:
-            return value["value"], value.get("unit")
-
-        if value is None:
-            return np.nan, None
-
-        if isinstance(value, (np.floating, float)):
-            return round(float(value), self.serialization_round_decimals), None
-
-        if isinstance(value, (np.integer, int)):
-            return int(value), None
-
-        return value, None
-
-    def _build_serialized_rows(self, grid_points, all_keys):
-        """Build serialized row dictionaries and collect units."""
-        rows = []
-        units = {}
-
-        for point in grid_points:
-            row = {}
-            for key in all_keys:
-                serialized_value, unit = self._serialize_grid_value(point.get(key))
-                row[key] = serialized_value
-                if unit is not None:
-                    units.setdefault(key, unit)
-            rows.append(row)
-
-        return rows, units
-
-    def _build_grid_metadata(self):
-        """Build metadata for the output grid table."""
-        return {
-            "coordinate_system": self.coordinate_system,
-            "reference_frame": "ICRS (J2000)",
-            "observing_time_utc": self.observing_time.isot if self.observing_time else None,
-            "observing_time_scale": self.observing_time.scale if self.observing_time else None,
-            "telescope_ids": self.telescope_ids,
-            "lookup_table": str(Path(self.lookup_table)) if self.lookup_table else None,
-        }
-
-    def serialize_quantity(self, value):
-        """Serialize Quantity."""
-        if isinstance(value, u.Quantity):
-            serialized_value = float(value.value)
-            rounded_value = round(serialized_value, self.serialization_round_decimals)
-            return {"value": rounded_value, "unit": str(value.unit)}
-        if isinstance(value, float):
-            return round(value, self.serialization_round_decimals)
-        if isinstance(value, np.floating):
-            return round(float(value), self.serialization_round_decimals)
-        if isinstance(value, np.integer):
-            return int(value)
-        return value
diff --git a/src/simtools/production_configuration/job_grid_io.py b/src/simtools/production_configuration/job_grid_io.py
new file mode 100644
index 0000000000..00f8f28f1f
--- /dev/null
+++ b/src/simtools/production_configuration/job_grid_io.py
@@ -0,0 +1,143 @@
+"""Read and write executable job grids for production preparation."""
+
+from pathlib import Path
+
+from astropy import units as u
+from astropy.table import Table
+
+JOB_GRID_COLUMNS = [
+    "primary",
+    "azimuth_angle_value",
+    "azimuth_angle_unit",
+    "zenith_angle_value",
+    "zenith_angle_unit",
+    "energy_min_value",
+    "energy_min_unit",
+    "energy_max_value",
+    "energy_max_unit",
+    "core_scatter_number",
+    "core_scatter_max_value",
+    "core_scatter_max_unit",
+    "view_cone_min_value",
+    "view_cone_min_unit",
+    "view_cone_max_value",
+    "view_cone_max_unit",
+    "nshow",
+    "model_version",
+    "array_layout_name",
+    "corsika_le_interaction",
+    "corsika_he_interaction",
+    "run_number",
+]
+
+_QUANTITY_FIELDS = {
+    "azimuth_angle": ("azimuth_angle_value", "azimuth_angle_unit"),
+    "zenith_angle": ("zenith_angle_value", "zenith_angle_unit"),
+    "energy_min": ("energy_min_value", "energy_min_unit"),
+    "energy_max": ("energy_max_value", "energy_max_unit"),
+    "core_scatter_max": ("core_scatter_max_value", "core_scatter_max_unit"),
+    "view_cone_min": ("view_cone_min_value", "view_cone_min_unit"),
+    "view_cone_max": ("view_cone_max_value", "view_cone_max_unit"),
+}
+
+
+def _serialize_quantity(value):
+    """Serialize a Quantity to value/unit columns."""
+    return float(value.value), str(value.unit)
+
+
+def _deserialize_quantity(value, unit):
+    """Deserialize value/unit columns to a Quantity."""
+    return float(value) * u.Unit(unit)
+
+
+def _serialize_job_row(job_row):
+    """Serialize one job row to the on-disk schema."""
+    serialized_row = {
+        "primary": job_row["primary"],
+        "core_scatter_number": int(job_row["core_scatter_number"]),
+        "nshow": int(job_row["nshow"]),
+        "model_version": job_row["model_version"],
+        "array_layout_name": job_row["array_layout_name"],
+        "corsika_le_interaction": job_row["corsika_le_interaction"],
+        "corsika_he_interaction": job_row["corsika_he_interaction"],
+        "run_number": int(job_row["run_number"]),
+    }
+
+    for quantity_name, (value_key, unit_key) in _QUANTITY_FIELDS.items():
+        serialized_row[value_key], serialized_row[unit_key] = _serialize_quantity(
+            job_row[quantity_name]
+        )
+
+    return serialized_row
+
+
+def _deserialize_job_row(serialized_row):
+    """Deserialize one stored row to the in-memory job-row schema."""
+    job_row = {
+        "primary": serialized_row["primary"],
+        "core_scatter_number": int(serialized_row["core_scatter_number"]),
+        "nshow": int(serialized_row["nshow"]),
+        "model_version": serialized_row["model_version"],
+        "array_layout_name": serialized_row["array_layout_name"],
+        "corsika_le_interaction": serialized_row["corsika_le_interaction"],
+        "corsika_he_interaction": serialized_row["corsika_he_interaction"],
+        "run_number": int(serialized_row["run_number"]),
+    }
+
+    for quantity_name, (value_key, unit_key) in _QUANTITY_FIELDS.items():
+        job_row[quantity_name] = _deserialize_quantity(
+            serialized_row[value_key],
+            serialized_row[unit_key],
+        )
+
+    return job_row
+
+
+def serialize_job_grid(job_rows, output_file, metadata=None):
+    """
+    Serialize executable job rows to ECSV output.
+
+    Parameters
+    ----------
+    job_rows : list[dict]
+        Job rows in the in-memory schema.
+    output_file : str or Path
+        Output file path. Must use the ``.ecsv`` suffix.
+    metadata : dict, optional
+        Metadata to store alongside the rows.
+    """
+    output_path = Path(output_file)
+    serialized_rows = [_serialize_job_row(job_row) for job_row in job_rows]
+    metadata = metadata or {}
+
+    if output_path.suffix.lower() != ".ecsv":
+        raise ValueError("Job grid output file must use the '.ecsv' extension.")
+
+    output_table = Table(rows=serialized_rows, names=JOB_GRID_COLUMNS)
+    output_table.meta = metadata
+    output_table.write(output_path, format="ascii.ecsv", overwrite=True)
+
+
+def read_job_grid(input_file):
+    """
+    Read executable job rows from ECSV input.
+
+    Parameters
+    ----------
+    input_file : str or Path
+        Input file path.
+
+    Returns
+    -------
+    tuple[list[dict], dict]
+        Deserialized job rows and metadata.
+    """
+    input_path = Path(input_file)
+
+    if input_path.suffix.lower() != ".ecsv":
+        raise ValueError("Job grid input file must use the '.ecsv' extension.")
+
+    table = Table.read(input_path, format="ascii.ecsv")
+    rows = [{column_name: row[column_name] for column_name in table.colnames} for row in table]
+    return [_deserialize_job_row(row) for row in rows], dict(table.meta)
diff --git a/src/simtools/production_configuration/observation_grid.py b/src/simtools/production_configuration/observation_grid.py
new file mode 100644
index 0000000000..9c13ef7422
--- /dev/null
+++ b/src/simtools/production_configuration/observation_grid.py
@@ -0,0 +1,431 @@
+"""
+Generate observation sampling grids with direction-dependent CORSIKA physics limits.
+
+Samples observing positions (Alt/Az or RA/Dec) and interpolates physics limits
+(minimum energy, maximum scatter radius, maximum viewcone) from lookup tables.
+Supports linear/log/1/cos binning modes.
+"""
+
+import logging
+
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import AltAz, EarthLocation, SkyCoord
+from astropy.units import Quantity
+
+from simtools.production_configuration.corsika_limits_lookup import (
+    CorsikaLimitsLookup,
+    attach_lookup_limits_to_point,
+)
+
+
+class ProductionGridEngine:
+    """
+    Generate observation grids with direction-dependent CORSIKA physics limits.
+
+    Samples observing positions per configured axes, interpolates physics limits from
+    lookup tables, and optionally converts between coordinate systems (Alt/Az vs RA/Dec).
+    Results feed into :func:`production_configuration.simulation_jobs.build_simulation_jobs`
+    to expand into full job matrices.
+    """
+
+    def __init__(
+        self,
+        axes,
+        coordinate_system="horizontal",
+        observing_location=None,
+        time_of_observation=None,
+        lookup_table=None,
+        array_layout_name=None,
+    ):
+        """
+        Initialize the production-grid engine.
+
+        Parameters
+        ----------
+        axes : dict
+            Dictionary where each key is the axis name and the value is a dictionary
+            defining the axis properties (range, binning, scaling, etc.).
+        coordinate_system : str, optional
+            The coordinate system for the grid generation.
+        observing_location : EarthLocation, optional
+            The location of the observation (latitude, longitude, height).
+        time_of_observation : Time, optional
+            The time of observation required for RA/Dec transforms.
+        lookup_table : str, optional
+            Path to the lookup table file (ECSV format).
+        array_layout_name : str, optional
+            Array layout name used to select lookup-table rows.
+        """
+        self._logger = logging.getLogger(__name__)
+        self.axes = axes["axes"] if "axes" in axes else axes
+        self.coordinate_system = coordinate_system
+        self.observing_location = (
+            observing_location
+            if observing_location is not None
+            else EarthLocation(lat=0.0 * u.deg, lon=0.0 * u.deg, height=0 * u.m)
+        )
+        self.time_of_observation = time_of_observation
+        self.lookup_table = lookup_table
+        self.array_layout_name = array_layout_name
+        self.interpolated_limits = {}
+        self._limits_lookup = None
+        if self.lookup_table:
+            self._limits_lookup = CorsikaLimitsLookup(
+                lookup_table=lookup_table,
+                array_layout_name=array_layout_name,
+            )
+        self.target_values = self._generate_target_values()
+
+        if self.lookup_table:
+            if self.coordinate_system == "ra_dec":
+                self._prepare_lookup_table_limits_for_point_interpolation()
+            else:
+                self._apply_lookup_table_limits()
+
+    def _require_time_of_observation(self):
+        """Return observing time if available, else raise a clear error."""
+        if self.time_of_observation is None:
+            raise ValueError("Observing time is required for ra_dec grid generation.")
+        return self.time_of_observation
+
+    def _get_max_zenith_for_radec_mode(self):
+        """Read maximum zenith from axes for RA/Dec direction sampling."""
+        zenith_axis = self.axes.get("zenith_angle")
+        if not zenith_axis or "range" not in zenith_axis or len(zenith_axis["range"]) != 2:
+            raise ValueError(
+                "RA/Dec direction sampling requires 'zenith_angle' axis with a valid "
+                "two-element 'range' in the axes definition."
+            )
+        return float(zenith_axis["range"][1])
+
+    def _prepare_lookup_table_limits_for_point_interpolation(self):
+        """Prepare lookup arrays for per-point interpolation in RA/Dec grid mode."""
+        self._limits_lookup.prepare_point_interpolators()
+
+    def _generate_target_values(self):
+        """
+        Generate target axis values and store them as Quantities.
+
+        Returns
+        -------
+        dict
+            Dictionary of target values for each axis, stored as Quantity objects.
+        """
+        target_values = {}
+        for axis_name, axis in self.axes.items():
+            axis_range = axis["range"]
+            binning = axis["binning"]
+            scaling = axis.get("scaling", "linear")
+            units = axis.get("units", None)
+
+            if axis_name == "azimuth":
+                values = self.create_circular_binning(axis_range, binning)
+            elif scaling == "log":
+                values = np.logspace(np.log10(axis_range[0]), np.log10(axis_range[1]), binning)
+            elif scaling == "1/cos":
+                cos_min = np.cos(np.radians(axis_range[0]))
+                cos_max = np.cos(np.radians(axis_range[1]))
+                inv_cos_values = np.linspace(1 / cos_min, 1 / cos_max, binning)
+                values = np.degrees(np.arccos(1 / inv_cos_values))
+            else:
+                values = np.linspace(axis_range[0], axis_range[1], binning)
+
+            if units:
+                values = values * u.Unit(units)
+
+            target_values[axis_name] = values
+
+        return target_values
+
+    def _apply_lookup_table_limits(self):
+        """Apply limits from the lookup table and interpolate values."""
+        self.interpolated_limits = self._limits_lookup.interpolate_grid_limits(self.target_values)
+
+    def _generate_radec_grid_direction_points(self):
+        """Generate direction points from declination lines and hour-angle spacing."""
+        time_of_observation = self._require_time_of_observation()
+        max_zenith = self._get_max_zenith_for_radec_mode()
+        lst_deg = time_of_observation.sidereal_time(
+            "apparent", longitude=self.observing_location.lon
+        ).deg
+
+        direction_points = []
+        for declination in np.arange(-90.0, 91.0, 1.0):
+            cos_dec = np.cos(np.deg2rad(declination))
+            step_ha = 1.0 / cos_dec if cos_dec > 1e-6 else 360.0
+            n_ha = max(1, int(np.ceil(360.0 / step_ha)))
+            hour_angles = np.linspace(-180.0, 180.0, n_ha, endpoint=False)
+            ra_values = (lst_deg - hour_angles) % 360.0
+
+            skycoord = SkyCoord(
+                ra=ra_values * u.deg,
+                dec=np.full_like(ra_values, declination) * u.deg,
+                frame="icrs",
+            )
+            altaz = skycoord.transform_to(
+                AltAz(location=self.observing_location, obstime=time_of_observation)
+            )
+
+            zenith_values = (90.0 * u.deg - altaz.alt).to(u.deg).value
+            mask = (zenith_values >= 0.0) & (zenith_values <= max_zenith)
+
+            for idx in np.nonzero(mask)[0]:
+                direction_points.append(
+                    {
+                        "zenith_angle": zenith_values[idx] * u.deg,
+                        "azimuth": altaz.az.deg[idx] * u.deg,
+                    }
+                )
+        return direction_points
+
+    def _generate_extra_axis_combinations(self, excluded_keys):
+        """Generate combinations for all axes except the excluded ones."""
+        extra_axes = {
+            key: value for key, value in self.target_values.items() if key not in excluded_keys
+        }
+        if not extra_axes:
+            return list(extra_axes.keys()), [], [np.array([])]
+
+        extra_value_arrays = [value.value for value in extra_axes.values()]
+        extra_units = [value.unit for value in extra_axes.values()]
+        extra_grid = np.meshgrid(*extra_value_arrays, indexing="ij")
+        extra_combinations = np.vstack(list(map(np.ravel, extra_grid))).T
+        return list(extra_axes.keys()), extra_units, extra_combinations
+
+    def _interpolate_limits_for_point(self, zenith, azimuth, nsb):
+        """Interpolate lookup-table limits for a single point."""
+        return self._limits_lookup.interpolate_point(zenith, azimuth, nsb)
+
+    def _add_lookup_limits_to_point(self, point, zenith, azimuth):
+        """Interpolate and attach lookup-table limits to a grid point."""
+        if not self.lookup_table:
+            return
+
+        nsb_value = point.get("nsb_level", 1)
+        if isinstance(nsb_value, Quantity):
+            nsb_value = nsb_value.value
+        limits = self._interpolate_limits_for_point(
+            zenith=zenith,
+            azimuth=azimuth,
+            nsb=float(nsb_value),
+        )
+        attach_lookup_limits_to_point(point, limits)
+
+    def _generate_grid_from_radec_axes(self, include_horizontal_coordinates=False):
+        """Generate grid points from explicit RA/Dec axes definitions."""
+        time_of_observation = self._require_time_of_observation()
+
+        axis_keys = [key for key in self.target_values if key not in ("zenith_angle", "azimuth")]
+        value_arrays = [self.target_values[key].value for key in axis_keys]
+        units = [self.target_values[key].unit for key in axis_keys]
+        grid = np.meshgrid(*value_arrays, indexing="ij")
+        combinations = np.vstack(list(map(np.ravel, grid))).T
+
+        grid_points = []
+        for combination in combinations:
+            grid_point = {
+                key: Quantity(combination[i], units[i]) for i, key in enumerate(axis_keys)
+            }
+
+            skycoord = SkyCoord(
+                ra=grid_point["ra"].to(u.deg),
+                dec=grid_point["dec"].to(u.deg),
+                frame="icrs",
+            )
+            altaz = skycoord.transform_to(
+                AltAz(location=self.observing_location, obstime=time_of_observation)
+            )
+            zenith = (90.0 * u.deg - altaz.alt).to(u.deg)
+            azimuth = altaz.az.to(u.deg)
+
+            if include_horizontal_coordinates:
+                grid_point["zenith_angle"] = zenith
+                grid_point["azimuth"] = azimuth
+
+            self._add_lookup_limits_to_point(
+                grid_point,
+                zenith=zenith.value,
+                azimuth=azimuth.value,
+            )
+            grid_points.append(grid_point)
+
+        return grid_points
+
+    def _generate_grid_radec_mode(self, include_horizontal_coordinates=False):
+        """Generate grid points for RA/Dec mode."""
+        if "ra" in self.axes and "dec" in self.axes:
+            return self._generate_grid_from_radec_axes(
+                include_horizontal_coordinates=include_horizontal_coordinates
+            )
+
+        direction_points = self._generate_radec_grid_direction_points()
+        extra_keys, extra_units, extra_combinations = self._generate_extra_axis_combinations(
+            excluded_keys=("zenith_angle", "azimuth")
+        )
+
+        grid_points = []
+        for direction_point in direction_points:
+            for extra_combination in extra_combinations:
+                point = dict(direction_point)
+                for i, key in enumerate(extra_keys):
+                    point[key] = Quantity(extra_combination[i], extra_units[i])
+
+                self._add_lookup_limits_to_point(
+                    point,
+                    zenith=point["zenith_angle"].value,
+                    azimuth=point["azimuth"].value,
+                )
+                grid_points.append(point)
+
+        return self.convert_coordinates(
+            grid_points,
+            keep_horizontal_coordinates=include_horizontal_coordinates,
+        )
+
+    def create_circular_binning(self, azimuth_range, num_bins):
+        """
+        Create bin centers for azimuth angles, handling circular wraparound.
+
+        Parameters
+        ----------
+        azimuth_range : tuple
+            (min_azimuth, max_azimuth), can wrap around 0 deg.
+        num_bins : int
+            Number of bins.
+
+        Returns
+        -------
+        np.ndarray
+            Array of bin centers.
+        """
+        azimuth_min, azimuth_max = azimuth_range
+        azimuth_min %= 360
+        azimuth_max %= 360
+
+        clockwise_distance = (azimuth_max - azimuth_min) % 360
+        counterclockwise_distance = (azimuth_min - azimuth_max) % 360
+
+        if clockwise_distance <= counterclockwise_distance:
+            return (
+                np.linspace(azimuth_min, azimuth_min + clockwise_distance, num_bins, endpoint=True)
+                % 360
+            )
+        return (
+            np.linspace(
+                azimuth_min, azimuth_min - counterclockwise_distance, num_bins, endpoint=True
+            )
+            % 360
+        )
+
+    def _generate_horizontal_grid(self):
+        """Generate grid points for horizontal (zenith/azimuth) mode."""
+        value_arrays = [value.value for value in self.target_values.values()]
+        units = [value.unit for value in self.target_values.values()]
+        grid = np.meshgrid(*value_arrays, indexing="ij")
+        combinations = np.vstack(list(map(np.ravel, grid))).T
+        grid_points = []
+
+        for combination in combinations:
+            grid_point = {
+                key: Quantity(combination[i], units[i])
+                for i, key in enumerate(self.target_values.keys())
+            }
+
+            zenith_idx = np.searchsorted(
+                self.target_values["zenith_angle"].value, grid_point["zenith_angle"].value
+            )
+            azimuth_idx = np.searchsorted(
+                self.target_values["azimuth"].value, grid_point["azimuth"].value
+            )
+            nsb_idx = np.searchsorted(
+                self.target_values["nsb_level"].value,
+                grid_point["nsb_level"].value,
+            )
+
+            if "lower_energy_threshold" in self.interpolated_limits:
+                grid_point["lower_energy_threshold"] = (
+                    self.interpolated_limits["lower_energy_threshold"][
+                        zenith_idx, azimuth_idx, nsb_idx
+                    ]
+                    * u.TeV
+                )
+
+            if "upper_scatter_radius" in self.interpolated_limits:
+                grid_point["scatter_radius"] = (
+                    self.interpolated_limits["upper_scatter_radius"][
+                        zenith_idx, azimuth_idx, nsb_idx
+                    ]
+                    * u.m
+                )
+
+            if "viewcone_radius" in self.interpolated_limits:
+                grid_point["viewcone_radius"] = (
+                    self.interpolated_limits["viewcone_radius"][zenith_idx, azimuth_idx, nsb_idx]
+                    * u.deg
+                )
+
+            grid_points.append(grid_point)
+
+        return grid_points
+
+    def generate_simulation_grid(self):
+        """Generate observation grid with CORSIKA limits. Always includes Alt/Az for backends."""
+        if self.coordinate_system == "ra_dec":
+            return self._generate_grid_radec_mode(include_horizontal_coordinates=True)
+        return self._generate_horizontal_grid()
+
+    def convert_altaz_to_radec(self, alt, az):
+        """
+        Convert Altitude/Azimuth (AltAz) coordinates to RA/Dec.
+
+        Parameters
+        ----------
+        alt : float
+            Altitude angle in degrees.
+        az : float
+            Azimuth angle in degrees.
+
+        Returns
+        -------
+        SkyCoord
+            SkyCoord object containing the RA/Dec coordinates.
+        """
+        if self.time_of_observation is None:
+            raise ValueError("Conversion to RA/Dec requires time_of_observation to be set.")
+
+        aa = AltAz(
+            alt=alt.to(u.rad),
+            az=az.to(u.rad),
+            location=self.observing_location,
+            obstime=self.time_of_observation,
+        )
+        sky_coord = SkyCoord(aa)
+        return sky_coord.icrs  # Return RA/Dec in ICRS frame
+
+    def convert_coordinates(self, grid_points, keep_horizontal_coordinates=False):
+        """
+        Convert the grid points RA/Dec coordinates if necessary.
+
+        Parameters
+        ----------
+        grid_points : list of dict
+            List of grid points.
+
+        Returns
+        -------
+        list of dict
+            The grid points with converted RA/Dec coordinates.
+        """
+        if self.coordinate_system == "ra_dec":
+            for point in grid_points:
+                if "zenith_angle" in point and "azimuth" in point:
+                    alt = (90.0 * u.deg) - point["zenith_angle"]
+                    az = point["azimuth"]
+                    radec = self.convert_altaz_to_radec(alt, az)
+                    if not keep_horizontal_coordinates:
+                        point.pop("zenith_angle")
+                        point.pop("azimuth")
+                    point["ra"] = radec.ra.deg * u.deg
+                    point["dec"] = radec.dec.deg * u.deg
+        return grid_points
diff --git a/src/simtools/production_configuration/simulation_jobs.py b/src/simtools/production_configuration/simulation_jobs.py
new file mode 100644
index 0000000000..497dd65dda
--- /dev/null
+++ b/src/simtools/production_configuration/simulation_jobs.py
@@ -0,0 +1,622 @@
+"""Expand observation grids into full simulation job matrices.
+
+Combines grids (from ProductionGridEngine or explicit axes) with primaries, interactions,
+model versions, energy ranges, and run counts into a complete job parameter set.
+"""
+
+import itertools
+import shlex
+
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import EarthLocation
+from astropy.time import Time
+
+from simtools.configuration import defaults
+from simtools.configuration.commandline_parser import CommandLineParser
+from simtools.layout.array_layout_utils import resolve_array_layout_name
+from simtools.model.site_model import SiteModel
+from simtools.production_configuration.corsika_limits_lookup import (
+    CorsikaLimitsLookup,
+    attach_lookup_limits_to_point,
+)
+from simtools.production_configuration.observation_grid import ProductionGridEngine
+from simtools.utils.general import ensure_list
+
+_GRID_AXES = [
+    "primary",
+    "azimuth_angle",
+    "zenith_angle",
+    "model_version",
+    "corsika_le_interaction",
+    "corsika_he_interaction",
+]
+
+_GRID_AXIS_DEFAULTS = {
+    "corsika_le_interaction": defaults.CORSIKA_LE_INTERACTION,
+    "corsika_he_interaction": defaults.CORSIKA_HE_INTERACTION,
+}
+
+GRID_AXIS_ARGUMENTS = {
+    "azimuth": {
+        "engine_axis": "azimuth",
+        "unit": "deg",
+        "help": "Azimuth range (deg)",
+    },
+    "zenith": {
+        "engine_axis": "zenith_angle",
+        "unit": "deg",
+        "help": "Zenith angle range (deg)",
+    },
+    "ra": {
+        "engine_axis": "ra",
+        "unit": "deg",
+        "help": "Right ascension range (deg)",
+    },
+    "dec": {
+        "engine_axis": "dec",
+        "unit": "deg",
+        "help": "Declination range (deg)",
+    },
+    "nsb": {
+        "engine_axis": "nsb_level",
+        "unit": "MHz",
+        "help": "NSB level range (MHz)",
+    },
+    "offset": {
+        "engine_axis": "offset",
+        "unit": "deg",
+        "help": "Offset range (deg)",
+    },
+}
+
+_AXIS_SCALING_CHOICES = ("linear", "log", "1/cos")
+_HORIZONTAL_AXES = ("azimuth", "zenith")
+_RADEC_AXES = ("ra", "dec")
+_REQUIRED_AXES = ("nsb", "offset")
+
+
+def _parse_axis_range_tokens(range_tokens):
+    """Parse a quantity pair from CLI axis range tokens."""
+    if len(range_tokens) == 1:
+        return CommandLineParser.parse_quantity_pair(range_tokens[0])
+    if len(range_tokens) == 2:
+        return tuple(u.Quantity(value) for value in range_tokens)
+    if len(range_tokens) == 4:
+        return tuple(
+            u.Quantity(f"{range_tokens[index]} {range_tokens[index + 1]}")
+            for index in range(0, len(range_tokens), 2)
+        )
+    raise ValueError("Axis range must contain exactly two quantities.")
+
+
+def _normalize_axis_spec_tokens(axis_spec):
+    """Return axis specification tokens from CLI or configuration input."""
+    if isinstance(axis_spec, str):
+        return shlex.split(axis_spec)
+    if isinstance(axis_spec, (list, tuple)):
+        if len(axis_spec) == 1 and isinstance(axis_spec[0], str):
+            return shlex.split(axis_spec[0])
+        return [str(token) for token in axis_spec]
+    raise TypeError("Axis definitions must be strings or lists of CLI-style tokens.")
+
+
+def _parse_axis_spec(axis_spec):
+    """Parse one compact axis definition."""
+    tokens = _normalize_axis_spec_tokens(axis_spec)
+    if len(tokens) < 3:
+        raise ValueError(
+            "Axis definitions require at least an axis name, range, and binning value."
+        )
+
+    axis_name = tokens[0]
+    if axis_name not in GRID_AXIS_ARGUMENTS:
+        supported_axes = ", ".join(sorted(GRID_AXIS_ARGUMENTS))
+        raise ValueError(f"Unknown axis '{axis_name}'. Supported axes: {supported_axes}.")
+
+    scaling = "linear"
+    if tokens[-1] in _AXIS_SCALING_CHOICES:
+        scaling = tokens[-1]
+        binning_token = tokens[-2]
+        range_tokens = tokens[1:-2]
+    else:
+        binning_token = tokens[-1]
+        range_tokens = tokens[1:-1]
+
+    if not range_tokens:
+        raise ValueError(f"Axis '{axis_name}' is missing its range definition.")
+
+    try:
+        binning = int(binning_token)
+    except ValueError as exc:
+        raise ValueError(f"Axis '{axis_name}' binning must be an integer.") from exc
+
+    axis_range = _parse_axis_range_tokens(range_tokens)
+    axis_args = GRID_AXIS_ARGUMENTS[axis_name]
+    return axis_name, {
+        "range": [u.Quantity(value).to_value(axis_args["unit"]) for value in axis_range],
+        "binning": binning,
+        "scaling": scaling,
+        "units": axis_args["unit"],
+    }
+
+
+def _iter_compact_axis_specs(args_dict):
+    """Iterate over compact axis definitions from CLI or configuration."""
+    axis_specs = args_dict.get("axis") or []
+    if isinstance(axis_specs, str):
+        axis_specs = [axis_specs]
+
+    normalized_specs = []
+    for axis_spec in axis_specs:
+        if (
+            isinstance(axis_spec, (list, tuple))
+            and len(axis_spec) > 1
+            and all(isinstance(item, str) for item in axis_spec)
+            and str(axis_spec[0]).strip() not in GRID_AXIS_ARGUMENTS
+        ):
+            normalized_specs.extend(axis_spec)
+            continue
+        normalized_specs.append(axis_spec)
+    return normalized_specs
+
+
+def _resolve_axis_configs(args_dict):
+    """Resolve compact axis definitions into one normalized mapping."""
+    axis_configs = {}
+    for axis_spec in _iter_compact_axis_specs(args_dict):
+        axis_name, axis_config = _parse_axis_spec(axis_spec)
+        axis_configs[axis_name] = axis_config
+
+    return axis_configs
+
+
+def _resolve_coordinate_system(axis_configs):
+    """Resolve the coordinate system from axis definitions."""
+    has_horizontal_axes = all(axis_name in axis_configs for axis_name in _HORIZONTAL_AXES)
+    has_radec_axes = all(axis_name in axis_configs for axis_name in _RADEC_AXES)
+
+    if has_horizontal_axes and has_radec_axes:
+        raise ValueError("Cannot define both azimuth/zenith and ra/dec axes at the same time.")
+    if has_radec_axes:
+        return "ra_dec"
+    if has_horizontal_axes:
+        return "horizontal"
+    return None
+
+
+def _resolve_coordinate_system_from_args(args_dict):
+    """Resolve the coordinate system from raw CLI arguments."""
+    return _resolve_coordinate_system(_resolve_axis_configs(args_dict))
+
+
+def resolve_time_of_observation(time_of_observation, args_dict):
+    """Generate Time object for the observing time. Required if RA/Dec axes are present."""
+    coordinate_system = _resolve_coordinate_system_from_args(args_dict)
+    if coordinate_system == "ra_dec":
+        if not time_of_observation:
+            raise ValueError("time_of_observation is required when using RA/Dec axes.")
+        return Time(time_of_observation, scale="utc")
+    return None
+
+
+def resolve_single_model_version(model_version):
+    """Resolve one model version for helpers that require a scalar version."""
+    if isinstance(model_version, list):
+        return model_version[0]
+    return model_version
+
+
+def build_observing_location(site, model_version):
+    """Build observing location from the site model."""
+    site_model = SiteModel(model_version=resolve_single_model_version(model_version), site=site)
+    return EarthLocation(
+        lat=site_model.get_parameter_value_with_unit("reference_point_latitude"),
+        lon=site_model.get_parameter_value_with_unit("reference_point_longitude"),
+        height=site_model.get_parameter_value_with_unit("reference_point_altitude"),
+    )
+
+
+def build_axes_dict_from_cli_args(args_dict):
+    """Build ProductionGridEngine-compatible axes configuration from CLI arguments."""
+    axis_configs = _resolve_axis_configs(args_dict)
+    coordinate_system = _resolve_coordinate_system(axis_configs)
+
+    if coordinate_system is None:
+        raise ValueError("Must provide either both azimuth/zenith or both ra/dec axis definitions.")
+
+    missing_required_axes = [
+        axis_name for axis_name in _REQUIRED_AXES if axis_name not in axis_configs
+    ]
+    if missing_required_axes:
+        missing_axes = ", ".join(missing_required_axes)
+        raise ValueError(f"Missing required shared axis definition(s): {missing_axes}.")
+
+    direction_axes = _RADEC_AXES if coordinate_system == "ra_dec" else _HORIZONTAL_AXES
+    return {
+        GRID_AXIS_ARGUMENTS[axis_name]["engine_axis"]: axis_configs[axis_name]
+        for axis_name in (*_REQUIRED_AXES, *direction_axes)
+    }
+
+
+def build_production_grid_engine(args_dict, array_layout_name=None):
+    """Build a production-grid engine from application arguments."""
+    axes = build_axes_dict_from_cli_args(args_dict)
+    coordinate_system = _resolve_coordinate_system_from_args(args_dict)
+    if coordinate_system == "ra_dec":
+        observing_location = build_observing_location(
+            site=args_dict["site"],
+            model_version=args_dict["model_version"],
+        )
+    elif coordinate_system == "horizontal":
+        coordinate_system = "horizontal"
+        observing_location = None
+    else:
+        raise ValueError("Must provide either both azimuth/zenith or both ra/dec axis definitions.")
+    resolved_layout_name = array_layout_name or resolve_array_layout_name(
+        args_dict.get("array_layout_name"),
+        resolve_single_model_version(args_dict.get("model_version")),
+    )
+    return ProductionGridEngine(
+        axes=axes,
+        coordinate_system=coordinate_system,
+        observing_location=observing_location,
+        time_of_observation=resolve_time_of_observation(
+            args_dict.get("time_of_observation"),
+            args_dict,
+        ),
+        lookup_table=args_dict.get("corsika_limits"),
+        array_layout_name=resolved_layout_name,
+    )
+
+
+def build_job_grid_metadata(args_dict):
+    """Build metadata stored alongside serialized executable job grids."""
+    time_of_observation = resolve_time_of_observation(
+        args_dict.get("time_of_observation"),
+        args_dict,
+    )
+    coordinate_system = _resolve_coordinate_system_from_args(args_dict)
+    return {
+        "site": args_dict.get("site"),
+        "simulation_software": args_dict.get("simulation_software"),
+        "coordinate_system": coordinate_system,
+        "time_of_observation_utc": time_of_observation.isot if time_of_observation else None,
+        "time_of_observation_scale": time_of_observation.scale if time_of_observation else None,
+        "corsika_limits": (
+            str(args_dict["corsika_limits"]) if args_dict.get("corsika_limits") else None
+        ),
+    }
+
+
+def normalize_grid_axes(args_dict):
+    """Return normalized grid axes for cartesian product expansion."""
+    return {
+        axis: (
+            ensure_list(args_dict[axis])
+            if axis in args_dict and args_dict[axis] is not None
+            else [_GRID_AXIS_DEFAULTS[axis]]
+            if axis in _GRID_AXIS_DEFAULTS
+            else [None]
+        )
+        for axis in _GRID_AXES
+    }
+
+
+def normalize_energy_ranges(energy_range):
+    """Normalize energy range argument to a list of ``(e_min, e_max)`` pairs."""
+    if isinstance(energy_range, tuple) and len(energy_range) == 2:
+        return [energy_range]
+
+    if isinstance(energy_range, list):
+        if len(energy_range) == 2 and all(hasattr(item, "to") for item in energy_range):
+            return [(energy_range[0], energy_range[1])]
+        if all(isinstance(item, (list, tuple)) and len(item) == 2 for item in energy_range):
+            return [tuple(item) for item in energy_range]
+
+    raise ValueError(
+        "energy_range must be one pair (e_min, e_max) or a list of (e_min, e_max) pairs."
+    )
+
+
+def get_energy_range_for_zenith_angle(
+    zenith_angle, energy_range_pair, corsika_limits, azimuth_angle=None, nsb_level=1.0
+):
+    """
+    Return a zenith-dependent energy range pair or ``None`` to skip the step.
+
+    The lower energy bound is clipped to the lookup-table threshold for the
+    requested direction. If the threshold exceeds the configured upper bound,
+    the simulation step is skipped.
+    """
+    interpolated_limits = _interpolate_corsika_limits(
+        corsika_limits, zenith_angle, azimuth_angle, nsb_level
+    )
+    if interpolated_limits is None:
+        return energy_range_pair
+    return _clip_energy_range_from_threshold(
+        energy_range_pair,
+        interpolated_limits["lower_energy_threshold"] * u.TeV,
+    )
+
+
+def get_core_scatter_max_for_zenith_angle(
+    zenith_angle, core_scatter, corsika_limits, azimuth_angle=None, nsb_level=1.0
+):
+    """
+    Return zenith-dependent max core-scatter value.
+
+    The lookup-table scatter radius is treated as an upper limit and therefore
+    clipped against the user-provided maximum value.
+    """
+    interpolated_limits = _interpolate_corsika_limits(
+        corsika_limits, zenith_angle, azimuth_angle, nsb_level
+    )
+    if interpolated_limits is None:
+        return core_scatter[1]
+    return _clip_max_quantity(core_scatter[1], interpolated_limits["upper_scatter_radius"] * u.m)
+
+
+def get_viewcone_max_for_zenith_angle(
+    zenith_angle, view_cone, corsika_limits, azimuth_angle=None, nsb_level=1.0
+):
+    """Return zenith-dependent max viewcone value."""
+    interpolated_limits = _interpolate_corsika_limits(
+        corsika_limits, zenith_angle, azimuth_angle, nsb_level
+    )
+    if interpolated_limits is None:
+        return view_cone[1]
+    return _clip_max_quantity(view_cone[1], interpolated_limits["viewcone_radius"] * u.deg)
+
+
+def calculate_log_energy_midpoint(energy_range_pair):
+    """Return the geometric-mean energy for an energy range pair."""
+    energy_min, energy_max = energy_range_pair
+
+    if not isinstance(energy_min, u.Quantity) or not isinstance(energy_max, u.Quantity):
+        raise TypeError("energy_range_pair must contain astropy Quantity values.")
+
+    energy_min_tev = energy_min.to(u.TeV)
+    energy_max_tev = energy_max.to(u.TeV)
+
+    if energy_min_tev <= 0 * u.TeV or energy_max_tev <= 0 * u.TeV:
+        raise ValueError("Energy range values must be strictly positive.")
+
+    mean_log_energy = np.mean(
+        [
+            np.log10(energy_min_tev.value),
+            np.log10(energy_max_tev.value),
+        ]
+    )
+    return 10**mean_log_energy * u.TeV
+
+
+def calculate_scaled_nshow(
+    energy_range_pair,
+    baseline_nshow,
+    nshow_power_index=None,
+    reference_energy=None,
+):
+    """Return an energy-dependent nshow value."""
+    if baseline_nshow < 1:
+        raise ValueError("baseline_nshow must be a positive integer.")
+
+    if nshow_power_index is None:
+        return baseline_nshow
+
+    if reference_energy is None:
+        raise ValueError("reference_energy is required when nshow_power_index is configured.")
+
+    midpoint_energy = calculate_log_energy_midpoint(energy_range_pair)
+    scaling_factor = (midpoint_energy / reference_energy.to(midpoint_energy.unit)).to_value(
+        u.dimensionless_unscaled
+    ) ** nshow_power_index
+    scaled_nshow = int(np.ceil(baseline_nshow * scaling_factor))
+
+    if scaled_nshow < 1:
+        raise ValueError("Scaled nshow must be at least 1.")
+
+    return scaled_nshow
+
+
+def _clip_energy_range_from_threshold(energy_range_pair, lower_energy_threshold):
+    """Clip the lower energy bound of a configured energy range."""
+    if lower_energy_threshold is None:
+        return energy_range_pair
+
+    energy_min, energy_max = energy_range_pair
+    lower_energy_threshold = lower_energy_threshold.to(energy_min.unit)
+    if lower_energy_threshold > energy_max.to(lower_energy_threshold.unit):
+        return None
+    if lower_energy_threshold <= energy_min:
+        return energy_range_pair
+    return lower_energy_threshold, energy_max
+
+
+def _interpolate_corsika_limits(corsika_limits, zenith_angle, azimuth_angle=None, nsb_level=1.0):
+    """Return interpolated lookup limits for one pointing, or ``None`` if disabled."""
+    if corsika_limits is None:
+        return None
+    if not isinstance(corsika_limits, CorsikaLimitsLookup):
+        corsika_limits = CorsikaLimitsLookup(corsika_limits)
+    return corsika_limits.interpolate_point(
+        zenith_angle,
+        0.0 * u.deg if azimuth_angle is None else azimuth_angle,
+        nsb_level,
+    )
+
+
+def _clip_max_quantity(configured_max, lookup_max):
+    """Clip a configured maximum value against an interpolated lookup value."""
+    if lookup_max is None:
+        return configured_max
+    return min(configured_max, lookup_max.to(configured_max.unit))
+
+
+def _resolve_nshow_params(args_dict):
+    """Extract and convert nshow-scaling parameters from an args dict."""
+    nshow = args_dict["nshow"]
+    nshow_power_index = args_dict.get("nshow_power_index")
+    reference_energy = args_dict.get("nshow_reference_energy")
+    if nshow_power_index is not None and reference_energy is not None:
+        reference_energy = u.Quantity(reference_energy)
+    return nshow, nshow_power_index, reference_energy
+
+
+def _build_rows_for_point(
+    point_base,
+    energy_ranges,
+    lower_energy_threshold,
+    nshow,
+    nshow_power_index,
+    reference_energy,
+    number_of_runs,
+    run_number,
+):
+    """Build all simulation-run rows for a single grid point across all energy ranges."""
+    rows = []
+    for energy_range_pair in energy_ranges:
+        selected_energy_range = _clip_energy_range_from_threshold(
+            energy_range_pair, lower_energy_threshold
+        )
+        if selected_energy_range is None:
+            continue
+        selected_nshow = calculate_scaled_nshow(
+            selected_energy_range, nshow, nshow_power_index, reference_energy
+        )
+        for i in range(number_of_runs):
+            rows.append(
+                {
+                    **point_base,
+                    "energy_min": selected_energy_range[0],
+                    "energy_max": selected_energy_range[1],
+                    "nshow": selected_nshow,
+                    "run_number": run_number + i,
+                }
+            )
+    return rows
+
+
+def _generate_observation_points_from_axes(azimuth_values, zenith_values, corsika_limits):
+    """Sample azimuth * zenith grid with optional CORSIKA limits interpolation."""
+    points = []
+    for azimuth, zenith in itertools.product(azimuth_values, zenith_values):
+        point = {
+            "azimuth": azimuth,
+            "zenith_angle": zenith,
+        }
+        if corsika_limits is not None:
+            attach_lookup_limits_to_point(point, corsika_limits.interpolate_point(zenith, azimuth))
+        points.append(point)
+    return points
+
+
+def _generate_observation_grids_per_layout(args_dict, grid_axes):
+    """Generate observation grids per array layout.
+
+    Uses either ProductionGridEngine or explicit azimuth/zenith axes.
+    """
+    observation_grids_per_layout = {}
+    resolved_layout_names = {
+        model_version: resolve_array_layout_name(args_dict.get("array_layout_name"), model_version)
+        for model_version in grid_axes["model_version"]
+    }
+    use_shared_axes_definition = bool(args_dict.get("axis"))
+    corsika_limits_path = args_dict.get("corsika_limits")
+
+    for model_version in grid_axes["model_version"]:
+        resolved_layout_name = resolved_layout_names[model_version]
+        if resolved_layout_name in observation_grids_per_layout:
+            continue
+
+        if use_shared_axes_definition:
+            observation_grids_per_layout[resolved_layout_name] = build_production_grid_engine(
+                args_dict,
+                array_layout_name=resolved_layout_name,
+            ).generate_simulation_grid()
+            continue
+
+        corsika_limits = None
+        if corsika_limits_path is not None:
+            corsika_limits = CorsikaLimitsLookup(
+                corsika_limits_path,
+                array_layout_name=resolved_layout_name,
+            )
+        observation_grids_per_layout[resolved_layout_name] = _generate_observation_points_from_axes(
+            azimuth_values=grid_axes["azimuth_angle"],
+            zenith_values=grid_axes["zenith_angle"],
+            corsika_limits=corsika_limits,
+        )
+
+    return observation_grids_per_layout, resolved_layout_names
+
+
+def build_simulation_jobs(args_dict):
+    """
+    Expand production config into full simulation job matrix.
+
+    Cartesian product: primaries * model_versions * interactions * observation_directions
+    * energy_ranges * run_counts. Energy ranges clipped by direction-dependent CORSIKA
+    limits.
+    nshow optionally energy-scaled.
+
+    Activates ProductionGridEngine when axis-range CLI arguments are provided;
+    otherwise uses explicit azimuth * zenith axes.
+
+    Returns
+    -------
+    list[dict]
+        Each job: primary, model_version, interactions, directions (Alt/Az), energy_min/max
+        (clipped), nshow, run_number, scatter/viewcone values (clipped by physics limits).
+    """
+    grid_axes = normalize_grid_axes(args_dict)
+    energy_ranges = normalize_energy_ranges(args_dict["energy_range"])
+    nshow, nshow_power_index, reference_energy = _resolve_nshow_params(args_dict)
+    number_of_runs = int(args_dict.get("number_of_runs", 1))
+    run_number = int(args_dict.get("run_number") or 1)
+
+    core_scatter = args_dict["core_scatter"]
+    view_cone = args_dict["view_cone"]
+    core_scatter_number = int(core_scatter[0])
+    view_cone_min = view_cone[0]
+    observation_grids_per_layout, resolved_layout_names = _generate_observation_grids_per_layout(
+        args_dict, grid_axes
+    )
+    rows = []
+
+    for primary, model_version, corsika_le, corsika_he in itertools.product(
+        grid_axes["primary"],
+        grid_axes["model_version"],
+        grid_axes["corsika_le_interaction"],
+        grid_axes["corsika_he_interaction"],
+    ):
+        resolved_layout_name = resolved_layout_names[model_version]
+        for point in observation_grids_per_layout[resolved_layout_name]:
+            observation_params = {
+                "primary": primary,
+                "azimuth_angle": point["azimuth"],
+                "zenith_angle": point["zenith_angle"],
+                "model_version": model_version,
+                "array_layout_name": resolved_layout_name,
+                "corsika_le_interaction": corsika_le,
+                "corsika_he_interaction": corsika_he,
+                "core_scatter_number": core_scatter_number,
+                "core_scatter_max": _clip_max_quantity(
+                    core_scatter[1], point.get("scatter_radius")
+                ),
+                "view_cone_min": view_cone_min,
+                "view_cone_max": _clip_max_quantity(view_cone[1], point.get("viewcone_radius")),
+            }
+            rows.extend(
+                _build_rows_for_point(
+                    point_base=observation_params,
+                    energy_ranges=energy_ranges,
+                    lower_energy_threshold=point.get("lower_energy_threshold"),
+                    nshow=nshow,
+                    nshow_power_index=nshow_power_index,
+                    reference_energy=reference_energy,
+                    number_of_runs=number_of_runs,
+                    run_number=run_number,
+                )
+            )
+    return rows
diff --git a/src/simtools/utils/general.py b/src/simtools/utils/general.py
index 0700dc8160..ed51c6b0fa 100644
--- a/src/simtools/utils/general.py
+++ b/src/simtools/utils/general.py
@@ -2,6 +2,7 @@
 
 import datetime
 import glob
+import json
 import logging
 import os
 import tarfile
@@ -216,8 +217,16 @@ def ensure_list(value):
         return []
     if isinstance(value, list):
         return value
-    if isinstance(value, tuple):
+    if isinstance(value, (tuple, set)):
         return list(value)
+    if isinstance(value, str):
+        stripped = value.strip()
+        if stripped.startswith("["):
+            try:
+                return json.loads(stripped)
+            except json.JSONDecodeError:
+                pass
+        return [stripped]
     return [value]
 
 
diff --git a/src/simtools/utils/value_conversion.py b/src/simtools/utils/value_conversion.py
index e12bbedd49..506fd64f6b 100644
--- a/src/simtools/utils/value_conversion.py
+++ b/src/simtools/utils/value_conversion.py
@@ -231,3 +231,25 @@ def _unit_as_string(unit):
         unit = [unit]
     unit = [str(element) if element is not None else None for element in unit]
     return unit[0] if len(set(unit)) == 1 else unit
+
+
+def get_value_in_unit(value, unit=None):
+    """
+    Return the numerical value of a Quantity in the requested unit.
+
+    Parameters
+    ----------
+    value: astropy.units.Quantity or any
+        Value to be converted.
+    unit: astropy.units.Unit or str or None
+        Unit to convert the value to.
+
+    Returns
+    -------
+    float or any
+        Numerical value of the Quantity in the requested unit,
+        or the original value if it is not a Quantity.
+    """
+    if isinstance(value, u.Quantity):
+        return value.to_value(unit) if unit is not None else value.value
+    return value
diff --git a/src/simtools/visualization/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py
index 4696d5befa..a138ca77b0 100644
--- a/src/simtools/visualization/plot_production_grid.py
+++ b/src/simtools/visualization/plot_production_grid.py
@@ -67,7 +67,7 @@ def __init__(
         self.grid_metadata = {}
         self.grid_points = self._load_grid_points()
 
-        observation_time_utc = observation_time or self.grid_metadata.get("observing_time_utc")
+        observation_time_utc = observation_time or self.grid_metadata.get("time_of_observation_utc")
         if observation_time_utc is None:
             observation_time_utc = "2025-01-01 00:00:00"
         self.time = Time(observation_time_utc, scale="utc")
diff --git a/tests/integration_tests/config/plot_production_grid.yml b/tests/integration_tests/config/plot_production_grid.yml
index e4e8054f0c..a9b0096507 100644
--- a/tests/integration_tests/config/plot_production_grid.yml
+++ b/tests/integration_tests/config/plot_production_grid.yml
@@ -2,7 +2,7 @@
 applications:
 - application: simtools-plot-production-grid
   configuration:
-    grid_points_file: tests/resources/production_grid_points.ecsv
+    grid_points_file: tests/resources/production_grid_points_horizontal.ecsv
     output_path: simtools-output
     site: North
     model_version: 6.0.2
diff --git a/tests/integration_tests/config/production_derive_statistics.yml b/tests/integration_tests/config/production_derive_statistics.yml
index 9bdc2b6c0a..0b29f46f11 100644
--- a/tests/integration_tests/config/production_derive_statistics.yml
+++ b/tests/integration_tests/config/production_derive_statistics.yml
@@ -6,7 +6,7 @@ applications:
     - 180
     base_path: tests/resources/production_dl2_fits/
     file_name_template: prod6_LaPalma-{zenith}deg_gamma_cone.N.Am-4LSTs09MSTs_ID0_reduced.fits
-    grid_points_production_file: tests/resources/production_grid_points.ecsv
+    grid_points_production_file: tests/resources/production_grid_points_horizontal.ecsv
     metrics_file: tests/resources/production_simulation_config_metrics.yml
     nsb:
     - 1
diff --git a/tests/integration_tests/config/production_generate_grid.yml b/tests/integration_tests/config/production_generate_grid.yml
deleted file mode 100644
index 621f8d8404..0000000000
--- a/tests/integration_tests/config/production_generate_grid.yml
+++ /dev/null
@@ -1,18 +0,0 @@
----
-applications:
-- application: simtools-production-generate-grid
-  configuration:
-    axes: tests/resources/production_grid_generation_axes_definition.yml
-    coordinate_system: zenith_azimuth
-    lookup_table: tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv
-    model_version: 6.0.2
-    observing_time: '2017-09-16 00:00:00'
-    output_file: grid_output.ecsv
-    output_path: simtools-output
-    site: North
-    telescope_ids: [LSTN-01]
-  integration_tests:
-  - output_file: grid_output.ecsv
-  test_name: production_generate_grid
-schema_name: application_workflow.metaschema
-schema_version: 0.4.0
diff --git a/tests/integration_tests/config/production_generate_grid_horizontal.yml b/tests/integration_tests/config/production_generate_grid_horizontal.yml
new file mode 100644
index 0000000000..aca0cce3c2
--- /dev/null
+++ b/tests/integration_tests/config/production_generate_grid_horizontal.yml
@@ -0,0 +1,32 @@
+---
+applications:
+- application: simtools-production-generate-grid
+  configuration:
+    array_layout_name:
+      by_version:
+        "<7.0.0": 1lst
+        ">=7.0.0": 1lst
+    axis:
+    - azimuth 310 deg 20 deg 3 linear
+    - zenith 30 deg 40 deg 2 linear
+    - nsb 4 MHz 5 MHz 2 linear
+    - offset 0 deg 10 deg 2 linear
+    core_scatter: 10 500 m
+    corsika_he_interaction: epos
+    corsika_le_interaction: urqmd
+    energy_range: 30 GeV 300 GeV
+    corsika_limits: tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv
+    model_version: 6.0.2
+    nshow: 5
+    number_of_runs: 2
+    output_file: job_grid_horizontal.ecsv
+    output_path: simtools-output
+    primary: gamma
+    simulation_software: corsika_sim_telarray
+    site: North
+    view_cone: 0 deg 10 deg
+  integration_tests:
+  - output_file: job_grid_horizontal.ecsv
+  test_name: production_generate_grid_horizontal
+schema_name: application_workflow.metaschema
+schema_version: 0.4.0
diff --git a/tests/integration_tests/config/production_generate_grid_horizontal_explicit.yml b/tests/integration_tests/config/production_generate_grid_horizontal_explicit.yml
new file mode 100644
index 0000000000..46346d2c9e
--- /dev/null
+++ b/tests/integration_tests/config/production_generate_grid_horizontal_explicit.yml
@@ -0,0 +1,43 @@
+---
+applications:
+- application: simtools-production-generate-grid
+  configuration:
+    array_layout_name:
+      by_version:
+        "<7.0.0": alpha
+        ">=7.0.0": CTAO-North-Alpha
+    azimuth_angle:
+    - North
+    - South
+    core_scatter: 10 500 m
+    corsika_he_interaction:
+    - epos
+    - qgs3
+    corsika_le_interaction: urqmd
+    energy_range:
+    - 30 GeV 30 GeV
+    - 300 GeV 300 GeV
+    label: test-production-grid
+    log_level: DEBUG
+    model_version:
+    - 6.3.0
+    - 7.0.0
+    nshow: 5
+    number_of_runs: 2
+    output_file: job_grid_horizontal_explicit.ecsv
+    output_path: simtools-output
+    primary:
+    - gamma
+    - proton
+    run_number: 10
+    simulation_software: corsika_sim_telarray
+    site: North
+    view_cone: 0 deg 10 deg
+    zenith_angle:
+    - 20
+    - 40
+  integration_tests:
+  - output_file: job_grid_horizontal_explicit.ecsv
+  test_name: production_generate_grid_horizontal_explicit
+schema_name: application_workflow.metaschema
+schema_version: 0.4.0
diff --git a/tests/integration_tests/config/production_generate_grid_ra_dec.yml b/tests/integration_tests/config/production_generate_grid_ra_dec.yml
new file mode 100644
index 0000000000..dfa1eee1b6
--- /dev/null
+++ b/tests/integration_tests/config/production_generate_grid_ra_dec.yml
@@ -0,0 +1,34 @@
+---
+applications:
+- application: simtools-production-generate-grid
+  configuration:
+    array_layout_name:
+      by_version:
+        "<7.0.0": 1lst
+        ">=7.0.0": 1lst
+    axis:
+    - ra 0 deg 360 deg 36 linear
+    - dec -90 deg 90 deg 18 linear
+    - nsb 4 MHz 4 MHz 1 linear
+    - offset 0 deg 10 deg 2 linear
+    core_scatter: 10 500 m
+    corsika_he_interaction: epos
+    corsika_le_interaction: urqmd
+    energy_range: 30 GeV 300 GeV
+    corsika_limits: tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv
+    model_version: 6.0.2
+    nshow: 5
+    number_of_runs: 2
+    time_of_observation: '2017-09-16 00:00:00'
+    output_file: job_grid_ra_dec.ecsv
+    output_path: simtools-output
+    primary: gamma
+    simulation_software: corsika_sim_telarray
+    site: North
+    view_cone: 0 deg 10 deg
+  integration_tests:
+  - output_file: job_grid_ra_dec.ecsv
+  skip_integration_test: Requires model DB access for site-location lookup in ra_dec mode.
+  test_name: production_generate_grid_ra_dec
+schema_name: application_workflow.metaschema
+schema_version: 0.4.0
diff --git a/tests/integration_tests/config/production_generate_grid_radec.yml b/tests/integration_tests/config/production_generate_grid_radec.yml
deleted file mode 100644
index 52d14cf302..0000000000
--- a/tests/integration_tests/config/production_generate_grid_radec.yml
+++ /dev/null
@@ -1,18 +0,0 @@
----
-applications:
-- application: simtools-production-generate-grid
-  configuration:
-    axes: tests/resources/production_grid_generation_axes_definition_radec.yml
-    coordinate_system: ra_dec
-    lookup_table: tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv
-    model_version: 6.0.2
-    observing_time: '2017-09-16 00:00:00'
-    output_file: grid_output_radec.ecsv
-    output_path: simtools-output
-    site: North
-    telescope_ids: [LSTN-01]
-  integration_tests:
-  - output_file: grid_output_radec.ecsv
-  test_name: production_generate_grid_radec
-schema_name: application_workflow.metaschema
-schema_version: 0.4.0
diff --git a/tests/integration_tests/config/simulate_prod_htcondor_generator_gamma_20_deg_north.yml b/tests/integration_tests/config/simulate_prod_htcondor_generator_gamma_20_deg_north.yml
index 4c10c14b3d..7bd2e9a4e8 100644
--- a/tests/integration_tests/config/simulate_prod_htcondor_generator_gamma_20_deg_north.yml
+++ b/tests/integration_tests/config/simulate_prod_htcondor_generator_gamma_20_deg_north.yml
@@ -3,29 +3,11 @@ applications:
 - application: simtools-simulate-prod-htcondor-generator
   configuration:
     apptainer_image: tests/resources/dummy_apptainer_image.sif
-    array_layout_name:
-      by_version:
-        "<7.0.0": alpha
-        ">=7.0.0": CTAO-North-Alpha
-    azimuth_angle: North
-    core_scatter: 10 500 m
-    corsika_he_interaction: epos
-    corsika_le_interaction: urqmd
-    energy_range: 100 GeV 1 TeV
-    view_cone: 0 deg 10 deg
-    eslope: -2.0
+    job_grid_file: tests/resources/production_job_grid_horizontal_single.ecsv
     label: test-production-North
     log_level: DEBUG
-    model_version: 6.0.2
-    nshow: 5
-    number_of_runs: 1
     output_path: htcondor_submit
-    primary: gamma
     priority: 5
-    run_number: 10
-    simulation_software: corsika_sim_telarray
-    site: North
-    zenith_angle: 20
   integration_tests:
   - test_output_files:
     - file: simulate_prod.submit.condor
diff --git a/tests/integration_tests/config/simulate_prod_htcondor_generator_grid_example.yml b/tests/integration_tests/config/simulate_prod_htcondor_generator_grid_horizontal.yml
similarity index 58%
rename from tests/integration_tests/config/simulate_prod_htcondor_generator_grid_example.yml
rename to tests/integration_tests/config/simulate_prod_htcondor_generator_grid_horizontal.yml
index 037fe91037..4b5f12a267 100644
--- a/tests/integration_tests/config/simulate_prod_htcondor_generator_grid_example.yml
+++ b/tests/integration_tests/config/simulate_prod_htcondor_generator_grid_horizontal.yml
@@ -5,42 +5,12 @@ applications:
     apptainer_image:
       0.30.0: tests/resources/dummy_apptainer_image.sif
       0.29.0: tests/resources/dummy_apptainer_image.sif
-    array_layout_name:
-      by_version:
-        "<7.0.0": alpha
-        ">=7.0.0": CTAO-North-Alpha
-    azimuth_angle:
-    - North
-    - South
-    core_scatter: 10 500 m
-    corsika_he_interaction:
-    - epos
-    - qgs3
-    corsika_le_interaction: urqmd
-    energy_range:
-    - 30 GeV 30 GeV
-    - 300 GeV 300 GeV
-    view_cone: 0 deg 10 deg
-    eslope: -2.0
+    job_grid_file: tests/resources/production_job_grid_horizontal.ecsv
     label: test-production-grid
     log_level: DEBUG
-    model_version:
-    - 6.3.0
-    - 7.0.0
-    nshow: 5
-    number_of_runs: 2
     output_path: htcondor_submit
-    primary:
-    - gamma
-    - proton
     priority: 5
-    run_number: 10
     simulation_output: ./simulation_output_path
-    simulation_software: corsika_sim_telarray
-    site: North
-    zenith_angle:
-    - 20
-    - 40
   integration_tests:
   - test_output_files:
     - file: simulate_prod.submit.0.30.0.condor
@@ -53,6 +23,6 @@ applications:
       path_descriptor: output_path
     - file: simulate_prod.submit.sh
       path_descriptor: output_path
-  test_name: htcondor-grid-example
+  test_name: htcondor-grid-horizontal
 schema_name: application_workflow.metaschema
 schema_version: 0.4.0
diff --git a/tests/resources/production_grid_generation_axes_definition_radec.yml b/tests/resources/production_grid_generation_axes_definition_ra_dec.yml
similarity index 100%
rename from tests/resources/production_grid_generation_axes_definition_radec.yml
rename to tests/resources/production_grid_generation_axes_definition_ra_dec.yml
diff --git a/tests/resources/production_grid_points.ecsv b/tests/resources/production_grid_points_horizontal.ecsv
similarity index 74%
rename from tests/resources/production_grid_points.ecsv
rename to tests/resources/production_grid_points_horizontal.ecsv
index 8ed4192829..ea09d1cebc 100644
--- a/tests/resources/production_grid_points.ecsv
+++ b/tests/resources/production_grid_points_horizontal.ecsv
@@ -9,12 +9,12 @@
 # - {name: scatter_radius, unit: m, datatype: float64}
 # - {name: viewcone_radius, unit: deg, datatype: float64}
 # meta: !!omap
-# - {coordinate_system: zenith_azimuth}
+# - {coordinate_system: horizontal}
 # - {reference_frame: ICRS (J2000)}
-# - {observing_time_utc: '2017-09-16T00:00:00.000'}
-# - {observing_time_scale: utc}
+# - {time_of_observation_utc: '2017-09-16T00:00:00.000'}
+# - {time_of_observation_scale: utc}
 # - telescope_ids: [LSTN-01]
-# - {lookup_table: tests/resources/corsika_simulation_limits/corsika_simulation_limits_lookup_grid_test.ecsv}
+# - {corsika_limits: tests/resources/corsika_simulation_limits/corsika_simulation_limits_lookup_grid_test.ecsv}
 # schema: astropy-2.0
 azimuth zenith_angle nsb_level offset lower_energy_threshold scatter_radius viewcone_radius
 310.0 30.0 4.0 0.0 0.004588 2047.848197 9.986606
diff --git a/tests/resources/production_job_grid_horizontal.ecsv b/tests/resources/production_job_grid_horizontal.ecsv
new file mode 100644
index 0000000000..52cb266253
--- /dev/null
+++ b/tests/resources/production_job_grid_horizontal.ecsv
@@ -0,0 +1,33 @@
+# %ECSV 1.0
+# ---
+# datatype:
+# - {name: primary, datatype: string}
+# - {name: azimuth_angle_value, datatype: float64}
+# - {name: azimuth_angle_unit, datatype: string}
+# - {name: zenith_angle_value, datatype: float64}
+# - {name: zenith_angle_unit, datatype: string}
+# - {name: energy_min_value, datatype: float64}
+# - {name: energy_min_unit, datatype: string}
+# - {name: energy_max_value, datatype: float64}
+# - {name: energy_max_unit, datatype: string}
+# - {name: core_scatter_number, datatype: int64}
+# - {name: core_scatter_max_value, datatype: float64}
+# - {name: core_scatter_max_unit, datatype: string}
+# - {name: view_cone_min_value, datatype: float64}
+# - {name: view_cone_min_unit, datatype: string}
+# - {name: view_cone_max_value, datatype: float64}
+# - {name: view_cone_max_unit, datatype: string}
+# - {name: nshow, datatype: int64}
+# - {name: model_version, datatype: string}
+# - {name: array_layout_name, datatype: string}
+# - {name: corsika_le_interaction, datatype: string}
+# - {name: corsika_he_interaction, datatype: string}
+# - {name: run_number, datatype: int64}
+# meta: !!omap
+# - {site: North}
+# - {simulation_software: corsika_sim_telarray}
+# - {coordinate_system: horizontal}
+# schema: astropy-2.0
+primary azimuth_angle_value azimuth_angle_unit zenith_angle_value zenith_angle_unit energy_min_value energy_min_unit energy_max_value energy_max_unit core_scatter_number core_scatter_max_value core_scatter_max_unit view_cone_min_value view_cone_min_unit view_cone_max_value view_cone_max_unit nshow model_version array_layout_name corsika_le_interaction corsika_he_interaction run_number
+gamma 0.0 deg 20.0 deg 30.0 GeV 30.0 GeV 10 400.0 m 0.0 deg 6.0 deg 5 6.3.0 alpha urqmd epos 10
+proton 180.0 deg 40.0 deg 300.0 GeV 300.0 GeV 10 500.0 m 0.0 deg 10.0 deg 5 7.0.0 CTAO-North-Alpha urqmd qgs3 11
diff --git a/tests/resources/production_job_grid_horizontal_single.ecsv b/tests/resources/production_job_grid_horizontal_single.ecsv
new file mode 100644
index 0000000000..b0054e11b4
--- /dev/null
+++ b/tests/resources/production_job_grid_horizontal_single.ecsv
@@ -0,0 +1,32 @@
+# %ECSV 1.0
+# ---
+# datatype:
+# - {name: primary, datatype: string}
+# - {name: azimuth_angle_value, datatype: float64}
+# - {name: azimuth_angle_unit, datatype: string}
+# - {name: zenith_angle_value, datatype: float64}
+# - {name: zenith_angle_unit, datatype: string}
+# - {name: energy_min_value, datatype: float64}
+# - {name: energy_min_unit, datatype: string}
+# - {name: energy_max_value, datatype: float64}
+# - {name: energy_max_unit, datatype: string}
+# - {name: core_scatter_number, datatype: int64}
+# - {name: core_scatter_max_value, datatype: float64}
+# - {name: core_scatter_max_unit, datatype: string}
+# - {name: view_cone_min_value, datatype: float64}
+# - {name: view_cone_min_unit, datatype: string}
+# - {name: view_cone_max_value, datatype: float64}
+# - {name: view_cone_max_unit, datatype: string}
+# - {name: nshow, datatype: int64}
+# - {name: model_version, datatype: string}
+# - {name: array_layout_name, datatype: string}
+# - {name: corsika_le_interaction, datatype: string}
+# - {name: corsika_he_interaction, datatype: string}
+# - {name: run_number, datatype: int64}
+# meta: !!omap
+# - {site: North}
+# - {simulation_software: corsika_sim_telarray}
+# - {coordinate_system: horizontal}
+# schema: astropy-2.0
+primary azimuth_angle_value azimuth_angle_unit zenith_angle_value zenith_angle_unit energy_min_value energy_min_unit energy_max_value energy_max_unit core_scatter_number core_scatter_max_value core_scatter_max_unit view_cone_min_value view_cone_min_unit view_cone_max_value view_cone_max_unit nshow model_version array_layout_name corsika_le_interaction corsika_he_interaction run_number
+gamma 0.0 deg 20.0 deg 100.0 GeV 1.0 TeV 10 500.0 m 0.0 deg 10.0 deg 5 6.0.2 alpha urqmd epos 10
diff --git a/tests/unit_tests/applications/test_production_generate_grid.py b/tests/unit_tests/applications/test_production_generate_grid.py
new file mode 100644
index 0000000000..9cf075d4b8
--- /dev/null
+++ b/tests/unit_tests/applications/test_production_generate_grid.py
@@ -0,0 +1,66 @@
+from pathlib import Path
+from types import SimpleNamespace
+from unittest.mock import Mock, patch
+
+import simtools.applications.production_generate_grid as app
+from simtools.configuration.commandline_parser import CommandLineParser
+
+
+@patch("simtools.applications.production_generate_grid.serialize_job_grid")
+@patch("simtools.applications.production_generate_grid.build_job_grid_metadata")
+@patch("simtools.applications.production_generate_grid.build_simulation_jobs")
+@patch("simtools.applications.production_generate_grid.build_application")
+def test_main_serializes_job_grid(
+    mock_build_application,
+    mock_build_simulation_jobs,
+    mock_build_job_grid_metadata,
+    mock_serialize_job_grid,
+):
+    io_handler = Mock()
+    io_handler.get_output_file.return_value = Path("job_grid.ecsv")
+    args = {
+        "output_file": "job_grid.ecsv",
+    }
+    mock_build_application.return_value = SimpleNamespace(args=args, io_handler=io_handler)
+    mock_build_simulation_jobs.return_value = [{"primary": "gamma"}]
+    mock_build_job_grid_metadata.return_value = {"site": "North"}
+
+    app.main()
+
+    mock_build_simulation_jobs.assert_called_once_with(args)
+    mock_build_job_grid_metadata.assert_called_once_with(args)
+    mock_serialize_job_grid.assert_called_once_with(
+        job_rows=[{"primary": "gamma"}],
+        output_file=Path("job_grid.ecsv"),
+        metadata={"site": "North"},
+    )
+
+
+def test_add_arguments_accepts_compact_axis_definitions():
+    parser = CommandLineParser()
+    app._add_arguments(parser)
+
+    args = parser.parse_args(
+        [
+            "--axis",
+            "azimuth",
+            "310",
+            "deg",
+            "20",
+            "deg",
+            "3",
+            "linear",
+            "--axis",
+            "nsb",
+            "4",
+            "MHz",
+            "5",
+            "MHz",
+            "2",
+        ]
+    )
+
+    assert args.axis == [
+        ["azimuth", "310", "deg", "20", "deg", "3", "linear"],
+        ["nsb", "4", "MHz", "5", "MHz", "2"],
+    ]
diff --git a/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py b/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py
index bfa76c9bdb..7cb9f88060 100644
--- a/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py
+++ b/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py
@@ -1,5 +1,6 @@
 #!/usr/bin/python3
 
+import argparse
 from types import SimpleNamespace
 from unittest.mock import patch
 
@@ -17,10 +18,19 @@ def test_main_uses_standard_build_application(
 
     app.main()
 
-    assert mock_build_application.call_args.kwargs["initialization_kwargs"] == {
-        "db_config": False,
-        "preserve_by_version_keys": ["array_layout_name"],
-        "simulation_model": ["site", "layout", "telescope", "model_version"],
-        "simulation_configuration": {"software": None, "corsika_configuration": ["all"]},
-    }
+    assert mock_build_application.call_args.kwargs == {}
     mock_generate_submission_script.assert_called_once_with({"output_path": "htcondor_submit"})
+
+
+def test_add_arguments_registers_job_grid_argument():
+    parser = argparse.ArgumentParser()
+
+    app._add_arguments(parser)
+    args = parser.parse_args(
+        [
+            "--job_grid_file",
+            "job_grid.ecsv",
+        ]
+    )
+
+    assert args.job_grid_file == "job_grid.ecsv"
diff --git a/tests/unit_tests/job_execution/test_htcondor_script_generator.py b/tests/unit_tests/job_execution/test_htcondor_script_generator.py
index ffef52b577..c3ca50ef9a 100644
--- a/tests/unit_tests/job_execution/test_htcondor_script_generator.py
+++ b/tests/unit_tests/job_execution/test_htcondor_script_generator.py
@@ -5,64 +5,136 @@
 import pytest
 
 from simtools.job_execution.htcondor_script_generator import (
-    _build_job_specs,
+    _format_param_value,
     _format_quantity,
     _get_submit_file,
     _get_submit_script,
     _resolve_apptainer_images,
     _sanitize_label_for_filename,
     _write_params_file,
+    build_job_specs,
     generate_submission_script,
 )
 
 
 @pytest.fixture
-def args_dict():
+def args_dict(tmp_test_directory):
     return {
-        "output_path": "/test_output",
+        "output_path": str(Path(tmp_test_directory) / "output"),
         "apptainer_image": "/path/to/apptainer_image.sif",
         "priority": 5,
-        "azimuth_angle": 45 * u.deg,
-        "zenith_angle": 20 * u.deg,
-        "energy_range": [1 * u.GeV, 10 * u.GeV],
-        "core_scatter": [0, 100 * u.m],
-        "view_cone": [0 * u.deg, 10 * u.deg],
+        "job_grid_file": str(Path(tmp_test_directory) / "output" / "job_grid.ecsv"),
         "label": "test_label",
         "simulation_software": "simtools",
-        "model_version": "v1.0",
         "site": "test_site",
-        "array_layout_name": "test_layout",
-        "primary": "gamma",
-        "nshow": 1000,
-        "run_number_offset": 1,
-        "run_number": 1,
-        "number_of_runs": 10,
         "log_level": "INFO",
-        "pack_for_grid_register": "simtools-output",
-        "corsika_le_interaction": "urqmd",
-        "corsika_he_interaction": "epos",
+        "simulation_output": "simtools-output",
     }
 
 
+@pytest.fixture
+def job_rows():
+    return [
+        {
+            "primary": "gamma",
+            "azimuth_angle": 45 * u.deg,
+            "zenith_angle": 20 * u.deg,
+            "energy_min": 1 * u.GeV,
+            "energy_max": 10 * u.GeV,
+            "core_scatter_number": 10,
+            "core_scatter_max": 100 * u.m,
+            "view_cone_min": 0 * u.deg,
+            "view_cone_max": 5 * u.deg,
+            "nshow": 1000,
+            "model_version": "7.0.0",
+            "array_layout_name": "CTAO-North-Alpha",
+            "corsika_le_interaction": "urqmd",
+            "corsika_he_interaction": "epos",
+            "run_number": 1,
+        }
+    ]
+
+
+@pytest.fixture
+def job_grid_metadata():
+    return {
+        "site": "test_site",
+        "simulation_software": "simtools",
+    }
+
+
+@mock.patch("simtools.job_execution.htcondor_script_generator.read_job_grid")
 @mock.patch("simtools.job_execution.htcondor_script_generator.Path.mkdir")
 @mock.patch("simtools.job_execution.htcondor_script_generator.open", new_callable=mock.mock_open)
 @mock.patch("simtools.job_execution.htcondor_script_generator.Path.chmod")
 @mock.patch("simtools.job_execution.htcondor_script_generator.Path.is_file", return_value=True)
-def test_generate_submission_script(mock_is_file, mock_chmod, mock_open, mock_mkdir, args_dict):
+def test_generate_submission_script(
+    mock_is_file,
+    mock_chmod,
+    mock_open,
+    mock_mkdir,
+    mock_read_job_grid,
+    args_dict,
+    job_rows,
+    job_grid_metadata,
+):
+    mock_read_job_grid.return_value = (job_rows, job_grid_metadata)
+
     generate_submission_script(args_dict)
 
     work_dir = Path(args_dict["output_path"])
     submit_file_name = "simulate_prod.submit"
 
     mock_is_file.assert_called_once()
+    mock_read_job_grid.assert_called_once_with(args_dict["job_grid_file"])
     mock_mkdir.assert_any_call(parents=True, exist_ok=True)
-
-    # Check if files are created and written
     mock_open.assert_any_call(work_dir / f"{submit_file_name}.condor", "w", encoding="utf-8")
     mock_open.assert_any_call(work_dir / f"{submit_file_name}.params.txt", "w", encoding="utf-8")
     mock_open.assert_any_call(work_dir / f"{submit_file_name}.sh", "w", encoding="utf-8")
+    mock_chmod.assert_called_once_with(0o755)
 
-    # Check if chmod is called
+
+@mock.patch("simtools.job_execution.htcondor_script_generator.read_job_grid")
+@mock.patch("simtools.job_execution.htcondor_script_generator.Path.mkdir")
+@mock.patch("simtools.job_execution.htcondor_script_generator.open", new_callable=mock.mock_open)
+@mock.patch("simtools.job_execution.htcondor_script_generator.Path.chmod")
+@mock.patch("simtools.job_execution.htcondor_script_generator.Path.is_file", return_value=True)
+def test_generate_submission_script_writes_label_specific_files(
+    mock_is_file,
+    mock_chmod,
+    mock_open,
+    mock_mkdir,
+    mock_read_job_grid,
+    args_dict,
+    tmp_test_directory,
+    job_rows,
+    job_grid_metadata,
+):
+    args_dict["output_path"] = str(Path(tmp_test_directory) / "output")
+    args_dict["htcondor_log_path"] = str(Path(tmp_test_directory) / "custom_logs")
+    args_dict["apptainer_image"] = {
+        "prod 7.0.0": "/path/to/prod.sif",
+        "beta": "/path/to/beta.sif",
+    }
+    mock_read_job_grid.return_value = (job_rows, job_grid_metadata)
+
+    generate_submission_script(args_dict)
+
+    work_dir = Path(args_dict["output_path"])
+
+    mock_is_file.assert_has_calls([mock.call(), mock.call()], any_order=False)
+    mock_open.assert_any_call(
+        work_dir / "simulate_prod.submit.prod_7.0.0.condor", "w", encoding="utf-8"
+    )
+    mock_open.assert_any_call(
+        work_dir / "simulate_prod.submit.prod_7.0.0.params.txt", "w", encoding="utf-8"
+    )
+    mock_open.assert_any_call(work_dir / "simulate_prod.submit.beta.condor", "w", encoding="utf-8")
+    mock_open.assert_any_call(
+        work_dir / "simulate_prod.submit.beta.params.txt", "w", encoding="utf-8"
+    )
+    mock_open.assert_any_call(work_dir / "simulate_prod.submit.sh", "w", encoding="utf-8")
+    mock_mkdir.assert_any_call(parents=True, exist_ok=True)
     mock_chmod.assert_called_once_with(0o755)
 
 
@@ -79,11 +151,6 @@ def test_generate_submission_script_raises_for_missing_apptainer_image(
 
 
 def test_get_submit_script(args_dict):
-    n_core_scatter = args_dict["core_scatter"][0]
-    view_cone_low = args_dict["view_cone"][0].to(u.deg).value
-    view_cone_high = args_dict["view_cone"][1].to(u.deg).value
-    container_output_path = str(Path("/", "tmp", "simtools-output"))
-
     expected_script = f"""#!/usr/bin/env bash
 
 # Process ID used to generate run number
@@ -92,12 +159,12 @@ def test_get_submit_script(args_dict):
 set -a; source "$2"
 apptainer_label="${{3}}"
 primary="${{4}}"
-model_version="${{14}}"
-array_layout_name="${{15}}"
-corsika_le_interaction="${{16}}"
-corsika_he_interaction="${{17}}"
-run_number="${{18}}"
-pack_for_grid_register="${{19}}"
+model_version="${{19}}"
+array_layout_name="${{20}}"
+corsika_le_interaction="${{21}}"
+corsika_he_interaction="${{22}}"
+run_number="${{23}}"
+pack_for_grid_register="${{24}}"
 energy_range_tag="erange-${{7}}${{8}}-${{9}}${{10}}"
 job_label="{args_dict["label"]}_${{corsika_he_interaction}}-${{corsika_le_interaction}}_${{energy_range_tag}}"
 
@@ -110,16 +177,16 @@ def test_get_submit_script(args_dict):
     --primary "$primary" \\
     --azimuth_angle "${{5}}" \\
     --zenith_angle "${{6}}" \\
-    --nshow "${{13}}" \\
+    --nshow "${{18}}" \\
     --energy_range "${{7}} ${{8}} ${{9}} ${{10}}" \\
-    --core_scatter "{n_core_scatter} ${{11}} ${{12}}" \\
-    --view_cone "{view_cone_low} deg {view_cone_high} deg" \\
+    --core_scatter "${{11}} ${{12}} ${{13}}" \\
+    --view_cone "${{14}} ${{15}} ${{16}} ${{17}}" \\
     --corsika_le_interaction "$corsika_le_interaction" \\
     --corsika_he_interaction "$corsika_he_interaction" \\
     --run_number "$run_number" \\
-    --run_number_offset {args_dict["run_number_offset"]} \\
+    --run_number_offset 0 \\
     --save_reduced_event_lists \\
-    --output_path {container_output_path} \\
+    --output_path /tmp/simtools-output \\
     --log_level {args_dict["log_level"]} \\
     --pack_for_grid_register "$pack_for_grid_register"
 """
@@ -132,21 +199,21 @@ def test_get_submit_file_uses_queue_from_params(tmp_test_directory):
     log_dir = Path(tmp_test_directory) / "htcondor_logs" / "log"
     error_dir = Path(tmp_test_directory) / "htcondor_logs" / "error"
     output_dir = Path(tmp_test_directory) / "htcondor_logs" / "output"
+    htcondor_dirs = {"log": log_dir, "error": error_dir, "output": output_dir}
     content = _get_submit_file(
         executable="simulate_prod.submit.sh",
         apptainer_image=apptainer_image,
         priority=1,
         params_file_name="simulate_prod.submit.params.txt",
-        log_dir=log_dir,
-        error_dir=error_dir,
-        output_dir=output_dir,
+        htcondor_dirs=htcondor_dirs,
     )
 
     assert "queue apptainer_label,primary" in content
-    assert "energy_min_value,energy_min_unit,energy_max_value,energy_max_unit" in content
-    assert "core_scatter_max_value,core_scatter_max_unit" in content
+    assert "core_scatter_number,core_scatter_max_value,core_scatter_max_unit" in content
+    assert (
+        "view_cone_min_value,view_cone_min_unit,view_cone_max_value,view_cone_max_unit" in content
+    )
     assert "nshow,model_version,array_layout_name" in content
-    assert "model_version,array_layout_name,corsika_le_interaction" in content
     assert "from simulate_prod.submit.params.txt" in content
     assert 'arguments = "$(process) env.txt' in content
     assert str(log_dir) in content
@@ -191,6 +258,17 @@ def test_resolve_apptainer_images_empty_dict():
         _resolve_apptainer_images({})
 
 
+def test_resolve_apptainer_images_raises_for_truthy_empty_mapping():
+    class TruthyEmptyDict(dict):
+        def __bool__(self):
+            return True
+
+    with pytest.raises(
+        ValueError, match="At least one apptainer image label/path must be configured"
+    ):
+        _resolve_apptainer_images(TruthyEmptyDict())
+
+
 def test_resolve_apptainer_images_invalid_type(tmp_test_directory):
     with pytest.raises(
         TypeError, match="apptainer_image must be a string path or a label-to-path dictionary"
@@ -223,125 +301,59 @@ def test_format_quantity_full_coverage():
     assert unit is None
 
 
+def test_format_param_value_raises_for_missing_required_value():
+    with pytest.raises(ValueError, match="Missing required value for field 'primary'"):
+        _format_param_value(None, "primary")
+
+
 def test_sanitize_label_for_filename():
     assert _sanitize_label_for_filename("  my label:7/0*0?  ") == "my_label_7_0_0_"
     assert _sanitize_label_for_filename("v7.0.0-beta_1") == "v7.0.0-beta_1"
     assert _sanitize_label_for_filename(42) == "42"
 
 
-def test_build_job_specs_expands_model_version_list(args_dict):
-    args_dict["model_version"] = ["6.3.0", "7.0.0"]
-
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-    model_versions = {job_spec["model_version"] for job_spec in job_specs}
-
-    assert model_versions == {"6.3.0", "7.0.0"}
-    assert len(job_specs) == 2 * args_dict["number_of_runs"]
-
-
-def test_build_job_specs_expands_energy_range_list_of_pairs(args_dict):
-    args_dict["number_of_runs"] = 1
-    args_dict["energy_range"] = [
-        (30 * u.GeV, 30 * u.GeV),
-        (300 * u.GeV, 300 * u.GeV),
-    ]
-
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-    energy_pairs = {(job_spec["energy_min"], job_spec["energy_max"]) for job_spec in job_specs}
-
-    assert len(job_specs) == 2
-    assert energy_pairs == {
-        (30 * u.GeV, 30 * u.GeV),
-        (300 * u.GeV, 300 * u.GeV),
-    }
-
-
-def test_build_job_specs_uses_default_interaction_models_when_missing(args_dict):
-    args_dict.pop("corsika_le_interaction")
-    args_dict.pop("corsika_he_interaction")
+def test_build_job_specs_reads_grid_file(args_dict, job_rows, job_grid_metadata):
+    with mock.patch(
+        "simtools.job_execution.htcondor_script_generator.read_job_grid",
+        return_value=(job_rows, job_grid_metadata),
+    ) as mock_read_job_grid:
+        job_specs, metadata = build_job_specs(args_dict, ["7.0.0"])
 
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-
-    assert {job_spec["corsika_le_interaction"] for job_spec in job_specs} == {"urqmd"}
-    assert {job_spec["corsika_he_interaction"] for job_spec in job_specs} == {"epos"}
-
-
-def test_build_job_specs_increments_run_number(args_dict):
-    args_dict["number_of_runs"] = 2
-    args_dict["run_number"] = 10
-    args_dict["model_version"] = ["6.3.0", "7.0.0"]
-
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-    run_numbers = [job_spec["run_number"] for job_spec in job_specs]
-
-    assert run_numbers == [10, 11, 12, 13]
-
-
-def test_write_params_file_resolves_array_layout_name_by_model_version(
-    args_dict, tmp_test_directory
-):
-    args_dict["number_of_runs"] = 1
-    args_dict["model_version"] = ["6.3.0", "7.0.0"]
-    args_dict["array_layout_name"] = {
-        "by_version": {
-            "<7.0.0": "alpha",
-            ">=7.0.0": "CTAO-North-Alpha",
-        }
-    }
-
-    params_file_path = Path(tmp_test_directory) / "params.txt"
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
+    mock_read_job_grid.assert_called_once_with(args_dict["job_grid_file"])
+    assert metadata == job_grid_metadata
+    assert job_specs[0]["image_label"] == "7.0.0"
+    assert job_specs[0]["pack_for_grid_register"] == "simtools-output/7.0.0"
+    assert job_specs[0]["array_layout_name"] == "CTAO-North-Alpha"
 
-    _write_params_file(params_file_path, job_specs)
 
-    params_lines = params_file_path.read_text(encoding="utf-8").splitlines()
-
-    assert "6.3.0 alpha urqmd epos 1 simtools-output/7.0.0" in params_lines[0]
-    assert "7.0.0 CTAO-North-Alpha urqmd epos 2 simtools-output/7.0.0" in params_lines[1]
-
-
-def test_write_params_file_resolves_stringified_by_version_layout(args_dict, tmp_test_directory):
-    args_dict["number_of_runs"] = 1
-    args_dict["model_version"] = ["6.3.0", "7.0.0"]
-    args_dict["array_layout_name"] = str(
-        {
-            "by_version": {
-                "<7.0.0": "alpha",
-                ">=7.0.0": "CTAO-North-Alpha",
-            }
-        }
-    )
-
-    params_file_path = Path(tmp_test_directory) / "params.txt"
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-
-    _write_params_file(params_file_path, job_specs)
-
-    params_lines = params_file_path.read_text(encoding="utf-8").splitlines()
-
-    assert "6.3.0 alpha urqmd epos 1 simtools-output/7.0.0" in params_lines[0]
-    assert "7.0.0 CTAO-North-Alpha urqmd epos 2 simtools-output/7.0.0" in params_lines[1]
+def test_build_job_specs_raises_for_missing_required_metadata(args_dict, job_rows):
+    with mock.patch(
+        "simtools.job_execution.htcondor_script_generator.read_job_grid",
+        return_value=(job_rows, {"coordinate_system": "horizontal"}),
+    ):
+        with pytest.raises(
+            ValueError, match=r"missing required field\(s\): site, simulation_software"
+        ):
+            build_job_specs(args_dict, ["7.0.0"])
 
 
 def test_write_params_file_keeps_energy_units(tmp_test_directory):
     params_file_path = Path(tmp_test_directory) / "params.txt"
     label_job_specs = [
         {
-            "apptainer_label": "7.0.0",
+            "image_label": "7.0.0",
             "primary": "gamma",
             "azimuth_angle": 0 * u.deg,
             "zenith_angle": 20 * u.deg,
             "energy_min": 30 * u.GeV,
             "energy_max": 10 * u.TeV,
+            "core_scatter_number": 10,
             "core_scatter_max": 200 * u.m,
+            "view_cone_min": 0 * u.deg,
+            "view_cone_max": 5 * u.deg,
             "nshow": 1000,
             "model_version": "7.0.0",
-            "array_layout_name": {
-                "by_version": {
-                    "<7.0.0": "alpha",
-                    ">=7.0.0": "CTAO-North-Alpha",
-                }
-            },
+            "array_layout_name": "CTAO-North-Alpha",
             "corsika_le_interaction": "urqmd",
             "corsika_he_interaction": "epos",
             "run_number": 10,
@@ -352,7 +364,7 @@ def test_write_params_file_keeps_energy_units(tmp_test_directory):
     _write_params_file(params_file_path, label_job_specs)
 
     assert params_file_path.read_text(encoding="utf-8") == (
-        "7.0.0 gamma 0.0 20.0 30.0 GeV 10.0 TeV 200.0 m "
+        "7.0.0 gamma 0.0 20.0 30.0 GeV 10.0 TeV 10 200.0 m 0.0 deg 5.0 deg "
         "1000 7.0.0 CTAO-North-Alpha urqmd epos 10 simtools-output/7.0.0\n"
     )
 
@@ -361,13 +373,16 @@ def test_write_params_file_replaces_whitespace_in_apptainer_label(tmp_test_direc
     params_file_path = Path(tmp_test_directory) / "params.txt"
     label_job_specs = [
         {
-            "apptainer_label": "grid label 7.0.0",
+            "image_label": "grid label 7.0.0",
             "primary": "gamma",
             "azimuth_angle": 0 * u.deg,
             "zenith_angle": 20 * u.deg,
             "energy_min": 30 * u.GeV,
             "energy_max": 10 * u.TeV,
+            "core_scatter_number": 10,
             "core_scatter_max": 200 * u.m,
+            "view_cone_min": 0 * u.deg,
+            "view_cone_max": 5 * u.deg,
             "nshow": 1000,
             "model_version": "7.0.0",
             "array_layout_name": "CTAO-North-Alpha",
@@ -381,35 +396,6 @@ def test_write_params_file_replaces_whitespace_in_apptainer_label(tmp_test_direc
     _write_params_file(params_file_path, label_job_specs)
 
     assert params_file_path.read_text(encoding="utf-8") == (
-        "grid_label_7.0.0 gamma 0.0 20.0 30.0 GeV 10.0 TeV 200.0 m "
+        "grid_label_7.0.0 gamma 0.0 20.0 30.0 GeV 10.0 TeV 10 200.0 m 0.0 deg 5.0 deg "
         "1000 7.0.0 CTAO-North-Alpha urqmd epos 10 simtools-output/grid_label_7.0.0\n"
     )
-
-
-@mock.patch(
-    "simtools.job_execution.htcondor_script_generator._get_energy_range_for_zenith_angle",
-    return_value=None,
-)
-def test_build_job_specs_skips_entries_when_energy_range_is_none(mock_energy_range, args_dict):
-    args_dict["corsika_limits"] = "limits.ecsv"
-    args_dict["number_of_runs"] = 1
-
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-
-    assert job_specs == []
-    mock_energy_range.assert_called_once()
-
-
-@mock.patch(
-    "simtools.job_execution.htcondor_script_generator._get_nshow_for_energy_range_and_zenith_angle",
-    return_value=777,
-)
-def test_build_job_specs_uses_dummy_nshow_when_corsika_limits_set(mock_nshow, args_dict):
-    args_dict["corsika_limits"] = "limits.ecsv"
-    args_dict["number_of_runs"] = 1
-
-    job_specs = _build_job_specs(args_dict, ["7.0.0"])
-
-    assert len(job_specs) == 1
-    assert job_specs[0]["nshow"] == 777
-    mock_nshow.assert_called_once()
diff --git a/tests/unit_tests/layout/test_array_layout_utils.py b/tests/unit_tests/layout/test_array_layout_utils.py
index e454cf55fa..cc67ffddeb 100644
--- a/tests/unit_tests/layout/test_array_layout_utils.py
+++ b/tests/unit_tests/layout/test_array_layout_utils.py
@@ -1008,19 +1008,6 @@ def test_create_regular_array_metadata():
     assert list(table["telescope_name"]) == sorted(table["telescope_name"])
 
 
-def test_get_array_name_valid():
-    assert array_layout_utils._get_array_name("4MST") == ("MST", 4)
-    assert array_layout_utils._get_array_name("1LST") == ("LST", 1)
-    assert array_layout_utils._get_array_name("2SST") == ("SST", 2)
-
-
-def test_get_array_name_invalid():
-    with pytest.raises(ValueError, match="Invalid array_name: 'MST'"):
-        array_layout_utils._get_array_name("MST")
-    with pytest.raises(ValueError, match="Invalid array_name: 'A4MST'"):
-        array_layout_utils._get_array_name("A4MST")
-
-
 def test_write_array_elements_from_file_to_repository_utm(tmp_test_directory):
     array_layout_utils.write_array_elements_from_file_to_repository(
         coordinate_system="utm",
@@ -1325,3 +1312,32 @@ def test_get_array_elements_from_db_for_layouts_all_with_no_layouts(mocker):
     assert result == {}
     instance.get_list_of_array_layouts.assert_called_once()
     instance.get_array_elements_for_layout.assert_not_called()
+
+
+def test_resolve_array_layout_name_resolves_stringified_by_version_layout():
+    array_layout_name = str(
+        {
+            "by_version": {
+                "<7.0.0": "alpha",
+                ">=7.0.0": "CTAO-North-Alpha",
+            }
+        }
+    )
+
+    assert (
+        array_layout_utils.resolve_array_layout_name(array_layout_name, "7.0.0")
+        == "CTAO-North-Alpha"
+    )
+
+
+def test_resolve_array_layout_name_unwraps_single_item_list():
+    assert (
+        array_layout_utils.resolve_array_layout_name(["CTAO-North-Alpha"], "7.0.0")
+        == "CTAO-North-Alpha"
+    )
+
+
+def test_resolve_array_layout_name_keeps_invalid_stringified_dict():
+    invalid_layout = "{not valid"
+
+    assert array_layout_utils.resolve_array_layout_name(invalid_layout, "7.0.0") == invalid_layout
diff --git a/tests/unit_tests/production_configuration/test_corsika_limits_lookup.py b/tests/unit_tests/production_configuration/test_corsika_limits_lookup.py
new file mode 100644
index 0000000000..79370e9995
--- /dev/null
+++ b/tests/unit_tests/production_configuration/test_corsika_limits_lookup.py
@@ -0,0 +1,84 @@
+import numpy as np
+import pytest
+from astropy import units as u
+
+from simtools.production_configuration.corsika_limits_lookup import CorsikaLimitsLookup
+
+
+def test_load_matching_lookup_arrays_filters_by_array_layout_name():
+    lookup = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv",
+        array_layout_name="1mst",
+    )
+
+    arrays = lookup.load_matching_lookup_arrays()
+
+    assert len(arrays["points"]) == 7
+    assert arrays["lower_energy_threshold"][0] == pytest.approx(0.022)
+
+
+def test_load_matching_lookup_arrays_raises_for_unknown_array_layout():
+    lookup = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv",
+        array_layout_name="does-not-exist",
+    )
+
+    with pytest.raises(ValueError, match="array_layout_name"):
+        lookup.load_matching_lookup_arrays()
+
+
+def test_load_matching_lookup_arrays_without_layout_returns_all_rows():
+    lookup = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv"
+    )
+
+    arrays = lookup.load_matching_lookup_arrays()
+
+    assert len(arrays["points"]) == 14
+
+
+def test_prepare_point_interpolators_builds_interpolator_state():
+    lookup = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv",
+        array_layout_name="1mst",
+    )
+
+    interpolators = lookup.prepare_point_interpolators()
+
+    assert lookup.lookup_points_for_interpolation.shape[1] == 3
+    assert set(interpolators) == {
+        "lower_energy_threshold",
+        "upper_scatter_radius",
+        "viewcone_radius",
+    }
+
+
+def test_interpolate_grid_limits_returns_requested_grid_shape():
+    lookup = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv",
+        array_layout_name="1mst",
+    )
+
+    interpolated = lookup.interpolate_grid_limits(
+        {
+            "zenith_angle": [20, 40] * u.deg,
+            "azimuth": [0, 180] * u.deg,
+            "nsb_level": [1, 5] * u.dimensionless_unscaled,
+        }
+    )
+
+    assert interpolated["lower_energy_threshold"].shape == (2, 2, 2)
+    assert np.isfinite(interpolated["upper_scatter_radius"][0, 0, 0])
+
+
+def test_interpolate_point_returns_interpolated_values():
+    lookup = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv",
+        array_layout_name="1mst",
+    )
+
+    interpolated = lookup.interpolate_point(20 * u.deg, 0 * u.deg, nsb=1)
+
+    assert interpolated["lower_energy_threshold"] == pytest.approx(0.022)
+    assert interpolated["upper_scatter_radius"] == pytest.approx(1000.0)
+    assert interpolated["viewcone_radius"] == pytest.approx(9.25)
diff --git a/tests/unit_tests/production_configuration/test_generate_production_grid.py b/tests/unit_tests/production_configuration/test_generate_production_grid.py
deleted file mode 100644
index fd5ea7ca55..0000000000
--- a/tests/unit_tests/production_configuration/test_generate_production_grid.py
+++ /dev/null
@@ -1,667 +0,0 @@
-import logging
-import warnings
-from pathlib import Path
-
-import numpy as np
-import pytest
-import yaml
-from astropy import units as u
-from astropy.coordinates import EarthLocation, SkyCoord
-from astropy.table import Table
-from astropy.time import Time
-from astropy.units import Quantity
-from astropy.utils import iers
-from astropy.utils.iers import IERSWarning
-from scipy.spatial import QhullError
-
-from simtools.production_configuration.generate_production_grid import GridGeneration
-
-
-@pytest.fixture(autouse=True, scope="module")
-def disable_iers_auto_download():
-    """Disable IERS auto-download during tests to avoid network dependency."""
-    previous_auto_download = iers.conf.auto_download
-    iers.conf.auto_download = False
-    try:
-        yield
-    finally:
-        iers.conf.auto_download = previous_auto_download
-
-
-def _create_grid_generation(
-    axes,
-    coordinate_system,
-    observing_location,
-    observing_time,
-    lookup_table,
-    simtel_file=None,
-):
-    """Create a GridGeneration instance with a standard telescope selection."""
-    return GridGeneration(
-        axes=axes,
-        coordinate_system=coordinate_system,
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-        telescope_ids=["LSTN-01"],
-        simtel_file=simtel_file,
-    )
-
-
-def _build_single_point_radec_axes_definition(source_radec):
-    """Build one-bin RA/Dec axes around a source coordinate."""
-    return {
-        "axes": {
-            "ra": {
-                "range": [source_radec.ra.deg, source_radec.ra.deg],
-                "binning": 1,
-                "scaling": "linear",
-                "units": "deg",
-            },
-            "dec": {
-                "range": [source_radec.dec.deg, source_radec.dec.deg],
-                "binning": 1,
-                "scaling": "linear",
-                "units": "deg",
-            },
-            "nsb_level": {
-                "range": [4, 4],
-                "binning": 1,
-                "scaling": "linear",
-                "units": "MHz",
-            },
-        }
-    }
-
-
-@pytest.fixture
-def axes_definition():
-    """Load the axes definition from the YAML file."""
-    axes_file = Path("tests/resources/production_grid_generation_axes_definition.yml")
-    with open(axes_file) as f:
-        return yaml.safe_load(f)
-
-
-@pytest.fixture
-def lookup_table():
-    """Load the lookup table from the resources directory."""
-    return str(
-        Path("tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv")
-    )
-
-
-@pytest.fixture
-def observing_location():
-    """Return a mock observing location."""
-    latitude = 28.7622  # degrees
-    longitude = -17.8920  # degrees
-    return EarthLocation(lon=longitude * u.deg, lat=latitude * u.deg, height=2000 * u.m)
-
-
-@pytest.fixture
-def observing_time():
-    """Return a mock observing time."""
-    return Time("2017-09-16 00:00:00")
-
-
-@pytest.fixture
-def grid_gen(axes_definition, lookup_table, observing_location, observing_time):
-    """Create a GridGeneration object with the provided fixtures."""
-    return _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="zenith_azimuth",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-    )
-
-
-def test_generate_grid(grid_gen):
-    grid_points = grid_gen.generate_grid()
-    assert isinstance(grid_points, list)
-    assert len(grid_points) > 0
-    assert all(isinstance(point, dict) for point in grid_points)
-    assert all("zenith_angle" in point for point in grid_points)
-    assert all("azimuth" in point for point in grid_points)
-    assert all("nsb_level" in point for point in grid_points)
-
-
-def test_generate_grid_log_scaling(
-    axes_definition, lookup_table, observing_location, observing_time
-):
-    """Test grid generation with logarithmic scaling for nsb_level axis."""
-    axes_definition["axes"]["nsb_level"] = {
-        "range": [2, 5],
-        "binning": 4,
-        "scaling": "log",
-        "units": "MHz",
-    }
-
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="zenith_azimuth",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-    )
-
-    grid_points = grid_gen.generate_grid()
-    nsb_values = [point["nsb_level"].value for point in grid_points]
-    unique_nsb_values = np.unique(nsb_values)
-
-    expected_values = np.logspace(
-        np.log10(axes_definition["axes"]["nsb_level"]["range"][0]),
-        np.log10(axes_definition["axes"]["nsb_level"]["range"][1]),
-        axes_definition["axes"]["nsb_level"]["binning"],
-    )
-
-    assert len(unique_nsb_values) == len(expected_values)
-    assert np.allclose(unique_nsb_values, expected_values, rtol=1e-4)
-
-
-def test_generate_grid_1_over_cos_scaling(
-    axes_definition, lookup_table, observing_location, observing_time
-):
-    """Test grid generation with 1/cos scaling for zenith_angle axis."""
-    axes_definition["axes"]["zenith_angle"] = {
-        "range": [30, 60],
-        "binning": 5,
-        "scaling": "1/cos",
-        "units": "deg",
-    }
-
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="zenith_azimuth",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-    )
-
-    grid_points = grid_gen.generate_grid()
-    zenith_values = [point["zenith_angle"].value for point in grid_points]
-    unique_zenith_values = np.unique(zenith_values)
-
-    cos_min = np.cos(np.radians(axes_definition["axes"]["zenith_angle"]["range"][0]))
-    cos_max = np.cos(np.radians(axes_definition["axes"]["zenith_angle"]["range"][1]))
-    cos_values = np.linspace(
-        1 / cos_min, 1 / cos_max, axes_definition["axes"]["zenith_angle"]["binning"]
-    )
-    expected_values = np.degrees(np.arccos(1 / cos_values))
-
-    assert len(unique_zenith_values) == len(expected_values)
-    assert np.allclose(unique_zenith_values, expected_values, rtol=1e-4)
-
-
-def test_generate_grid_radec_mode_minimal(observing_location, observing_time):
-    """Generate a minimal RA/Dec-native grid and apply zenith cut."""
-    axes_definition = {
-        "axes": {
-            "zenith_angle": {
-                "range": [0, 10],
-                "binning": 2,
-                "scaling": "linear",
-                "units": "deg",
-            },
-            "azimuth": {
-                "range": [0, 10],
-                "binning": 2,
-                "scaling": "linear",
-                "units": "deg",
-            },
-            "nsb_level": {
-                "range": [4, 4],
-                "binning": 1,
-                "scaling": "linear",
-                "units": "MHz",
-            },
-        }
-    }
-
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="ra_dec",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=None,
-    )
-
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore", IERSWarning)
-        grid_points = grid_gen.generate_grid()
-
-    assert len(grid_points) > 0
-    assert all("zenith_angle" in point and "azimuth" in point for point in grid_points)
-    assert all(point["zenith_angle"].value <= 10 for point in grid_points)
-
-
-def test_generate_grid_radec_axes_mode(observing_location, observing_time):
-    """Generate a grid directly from explicit RA/Dec axes."""
-    source_altaz = SkyCoord(
-        alt=65.0 * u.deg,
-        az=210.0 * u.deg,
-        frame="altaz",
-        obstime=observing_time,
-        location=observing_location,
-    )
-    source_radec = source_altaz.icrs
-
-    axes_definition = _build_single_point_radec_axes_definition(source_radec)
-
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="ra_dec",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=None,
-    )
-
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore", IERSWarning)
-        grid_points = grid_gen.generate_grid()
-
-    assert len(grid_points) == 1
-    assert "ra" in grid_points[0]
-    assert "dec" in grid_points[0]
-    assert "zenith_angle" not in grid_points[0]
-    assert "azimuth" not in grid_points[0]
-    assert grid_points[0]["ra"].value == pytest.approx(source_radec.ra.deg, abs=0.05)
-    assert grid_points[0]["dec"].value == pytest.approx(source_radec.dec.deg, abs=0.05)
-
-
-def test_generate_grid_radec_axes_mode_keeps_below_horizon_points(
-    observing_location, observing_time
-):
-    """Keep explicit RA/Dec YAML points even when they are below the horizon."""
-    source_altaz = SkyCoord(
-        alt=-20.0 * u.deg,
-        az=45.0 * u.deg,
-        frame="altaz",
-        obstime=observing_time,
-        location=observing_location,
-    )
-    source_radec = source_altaz.icrs
-
-    axes_definition = _build_single_point_radec_axes_definition(source_radec)
-
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="ra_dec",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=None,
-    )
-
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore", IERSWarning)
-        grid_points = grid_gen.generate_grid()
-
-    assert len(grid_points) == 1
-    assert grid_points[0]["ra"].value == pytest.approx(source_radec.ra.deg, abs=0.05)
-    assert grid_points[0]["dec"].value == pytest.approx(source_radec.dec.deg, abs=0.05)
-
-
-def test_generate_grid_radec_axes_mode_with_lookup(
-    lookup_table, observing_location, observing_time
-):
-    """Interpolate production limits for explicit RA/Dec axes points."""
-    source_altaz = SkyCoord(
-        alt=50.0 * u.deg,
-        az=0.0 * u.deg,
-        frame="altaz",
-        obstime=observing_time,
-        location=observing_location,
-    )
-    source_radec = source_altaz.icrs
-    axes_definition = _build_single_point_radec_axes_definition(source_radec)
-    axes_definition["axes"]["nsb_level"]["range"] = [1, 1]
-
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="ra_dec",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-    )
-
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore", IERSWarning)
-        grid_points = grid_gen.generate_grid()
-
-    assert len(grid_points) == 1
-    assert "lower_energy_threshold" in grid_points[0]
-    assert "scatter_radius" in grid_points[0]
-    assert "viewcone_radius" in grid_points[0]
-    assert grid_points[0]["lower_energy_threshold"].value == pytest.approx(0.007, rel=1e-2)
-    assert grid_points[0]["scatter_radius"].value == pytest.approx(1100.0, rel=1e-2)
-    assert grid_points[0]["viewcone_radius"].value == pytest.approx(7.0, rel=1e-2)
-
-
-def test_generate_grid_radec_axes_mode_with_sparse_lookup_raises_clear_error(
-    lookup_table, observing_location, observing_time, monkeypatch
-):
-    """Raise a user-facing error when lookup points are too sparse for 3D interpolation."""
-    source_altaz = SkyCoord(
-        alt=50.0 * u.deg,
-        az=0.0 * u.deg,
-        frame="altaz",
-        obstime=observing_time,
-        location=observing_location,
-    )
-    source_radec = source_altaz.icrs
-    axes_definition = _build_single_point_radec_axes_definition(source_radec)
-
-    def _raise_qhull(*args, **kwargs):
-        raise QhullError("mocked sparse triangulation failure")
-
-    monkeypatch.setattr(
-        "simtools.production_configuration.generate_production_grid.LinearNDInterpolator",
-        _raise_qhull,
-    )
-
-    with pytest.raises(ValueError, match="does not contain enough unique points"):
-        _create_grid_generation(
-            axes=axes_definition,
-            coordinate_system="ra_dec",
-            observing_location=observing_location,
-            observing_time=observing_time,
-            lookup_table=lookup_table,
-        )
-
-
-def test_interpolated_limits(grid_gen):
-    grid_gen.interpolated_limits = {
-        "lower_energy_threshold": np.array(
-            [[[0.1, 0.2], [0.3, 0.4], [0.5, 0.6]], [[0.7, 0.8], [0.9, 1.0], [1.1, 1.2]]]
-        ),
-        "upper_scatter_radius": np.array(
-            [[[100, 200], [300, 400], [500, 600]], [[700, 800], [900, 1000], [1100, 1200]]]
-        ),
-        "viewcone_radius": np.array([[[1, 2], [3, 4], [5, 6]], [[7, 8], [9, 10], [11, 12]]]),
-        "target_zeniths": np.array([30, 40]),
-        "target_azimuths": np.array([310, 345, 20]),
-        "target_nsb": np.array([4, 5]),
-    }
-
-    grid_points = grid_gen.generate_grid()
-
-    # Check that interpolated values are correctly assigned
-    for point in grid_points:
-        assert "lower_energy_threshold" in point
-        assert "scatter_radius" in point
-        assert "viewcone_radius" in point
-        assert isinstance(point["lower_energy_threshold"], Quantity)
-        assert isinstance(point["scatter_radius"], Quantity)
-        assert isinstance(point["viewcone_radius"], Quantity)
-
-
-def test_load_matching_lookup_arrays_with_simtel_id_mapping(grid_gen, tmp_test_directory):
-    """Match lookup rows when telescope IDs are stored as sim_telarray numeric IDs."""
-    lookup_table_path = tmp_test_directory / "lookup_simtel_ids.ecsv"
-    Table(
-        {
-            "telescope_ids": ["[1]"],
-            "zenith": [20.0],
-            "azimuth": [0.0],
-            "nsb_level": [4.0],
-            "lower_energy_limit": [0.01],
-            "upper_radius_limit": [1200.0],
-            "viewcone_radius": [5.0],
-        }
-    ).write(lookup_table_path, format="ascii.ecsv", overwrite=True)
-
-    grid_gen.lookup_table = str(lookup_table_path)
-    grid_gen.telescope_ids = ["LSTN-01"]
-    grid_gen._simtel_id_to_name = {1: "LSTN-01"}
-
-    lookup_arrays = grid_gen._load_matching_lookup_arrays()
-
-    assert lookup_arrays["points"].shape[0] == 1
-    assert lookup_arrays["lower_energy_threshold"][0] == pytest.approx(0.01)
-
-
-def test_load_matching_lookup_arrays_numeric_ids_require_simtel_file(grid_gen, tmp_test_directory):
-    """Raise a clear error when numeric lookup telescope IDs are present without mapping input."""
-    lookup_table_path = tmp_test_directory / "lookup_simtel_ids_no_mapping.ecsv"
-    Table(
-        {
-            "telescope_ids": ["[1]"],
-            "zenith": [20.0],
-            "azimuth": [0.0],
-            "nsb_level": [4.0],
-            "lower_energy_limit": [0.01],
-            "upper_radius_limit": [1200.0],
-            "viewcone_radius": [5.0],
-        }
-    ).write(lookup_table_path, format="ascii.ecsv", overwrite=True)
-
-    grid_gen.lookup_table = str(lookup_table_path)
-    grid_gen.telescope_ids = ["LSTN-01"]
-    grid_gen.simtel_file = None
-    grid_gen._simtel_id_to_name = {}
-
-    with pytest.raises(ValueError, match="Provide --simtel_file"):
-        grid_gen._load_matching_lookup_arrays()
-
-
-def test_serialize_grid_points_with_output_file(grid_gen, tmp_test_directory, caplog):
-    """Test serialize_grid_points when an output file is provided."""
-    grid_points = [
-        {
-            "zenith_angle": 30 * u.deg,
-            "azimuth": 310 * u.deg,
-            "nsb_level": 4,
-            "lower_energy_threshold": 0.1 * u.TeV,
-            "scatter_radius": 100 * u.m,
-            "viewcone_radius": 1 * u.deg,
-        },
-        {
-            "zenith_angle": 40 * u.deg,
-            "azimuth": 345 * u.deg,
-            "nsb_level": 5,
-            "lower_energy_threshold": 0.2 * u.TeV,
-            "scatter_radius": 200 * u.m,
-            "viewcone_radius": 2 * u.deg,
-        },
-    ]
-
-    output_file = tmp_test_directory / "grid_output.ecsv"
-    with caplog.at_level(logging.INFO):
-        grid_gen.serialize_grid_points(grid_points, output_file=output_file)
-    assert output_file.exists()
-
-    output_data = Table.read(output_file, format="ascii.ecsv")
-    assert "zenith_angle" in output_data.colnames
-    assert "lower_energy_threshold" in output_data.colnames
-    assert output_data.meta["coordinate_system"] == grid_gen.coordinate_system
-    assert output_data.meta["reference_frame"] == "ICRS (J2000)"
-
-    assert f"Output saved to {output_file}" in caplog.text
-
-
-def test_serialize_quantity(grid_gen):
-    # Case 1: Value is a Quantity (single value)
-    quantity = 5 * u.m
-    serialized = grid_gen.serialize_quantity(quantity)
-    assert serialized == {"value": 5, "unit": "m"}
-
-    # Case 2: Value is not a Quantity (single value)
-    value = 5
-    serialized_value = grid_gen.serialize_quantity(value)
-    assert serialized_value == value
-
-
-@pytest.mark.xfail(reason="May fail due to IERS data download timeout", strict=False)
-def test_convert_altaz_to_radec_and_coordinates(grid_gen):
-    warnings.simplefilter("error", IERSWarning)
-    # Case 1: Valid AltAz to RA/Dec conversion
-    alt, az = 45.0 * u.deg, 30.0 * u.deg
-    radec = grid_gen.convert_altaz_to_radec(alt, az)
-    assert isinstance(radec, SkyCoord)
-    assert radec.ra.unit == u.deg
-    assert radec.dec.unit == u.deg
-
-    # Case 2: Valid coordinate conversion
-    grid_gen.coordinate_system = "ra_dec"
-    grid_points = [
-        {"zenith_angle": 30 * u.deg, "azimuth": 45 * u.deg},
-        {"zenith_angle": 20 * u.deg, "azimuth": 60 * u.deg},
-    ]
-    converted_points = grid_gen.convert_coordinates(grid_points)
-    assert "ra" in converted_points[0]
-    assert "dec" in converted_points[0]
-
-    # Case 3: Missing zenith_angle or azimuth
-    grid_points = [{"azimuth": 45 * u.deg}]
-    converted_points = grid_gen.convert_coordinates(grid_points)
-    assert "ra" not in converted_points[0]
-    assert "dec" not in converted_points[0]
-
-
-def test_create_circular_binning(grid_gen):
-    # Case 1: No wraparound
-    bins = grid_gen.create_circular_binning((0, 150), 6)
-    assert len(bins) == 6
-    assert bins[0] == 0
-    assert bins[-1] == 150
-    assert bins[1] == 30
-    assert bins[2] == 60
-
-    # Case 2: Wraparound
-    bins = grid_gen.create_circular_binning((300, 20), 5)
-    assert len(bins) == 5
-    assert bins[0] == 300
-    assert bins[-1] == 20
-    assert bins[1] == 320
-
-    # Case 3: Single bin
-    bins = grid_gen.create_circular_binning((0, 360), 1)
-    assert len(bins) == 1
-    assert bins[0] == 0
-
-
-def test_create_circular_binning_with_shortest_path(grid_gen):
-    # Case 1: Clockwise path (shortest distance)
-    bins = grid_gen.create_circular_binning((350, 10), 3)
-    expected_bins = [350, 0, 10]
-    assert np.allclose(bins, expected_bins)
-
-    # Case 2: Counterclockwise path (shortest distance)
-    bins = grid_gen.create_circular_binning((10, 350), 3)
-    expected_bins = [10, 0, 350]
-    assert np.allclose(bins, expected_bins)
-
-
-def test_apply_lookup_table_limits(lookup_table, observing_location, observing_time):
-    axes_definition = {
-        "axes": {
-            "azimuth": {
-                "range": [0, 180],
-                "binning": 2,
-                "scaling": "linear",
-                "units": "deg",
-            },
-            "zenith_angle": {
-                "range": [20, 40],
-                "binning": 2,
-                "scaling": "linear",
-                "units": "deg",
-            },
-            "nsb_level": {
-                "range": [1, 1],
-                "binning": 1,
-                "scaling": "linear",
-                "units": "MHz",
-            },
-        }
-    }
-    grid_gen = _create_grid_generation(
-        axes=axes_definition,
-        coordinate_system="zenith_azimuth",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-    )
-
-    assert "lower_energy_threshold" in grid_gen.interpolated_limits
-    assert "upper_scatter_radius" in grid_gen.interpolated_limits
-    assert "viewcone_radius" in grid_gen.interpolated_limits
-    assert np.isclose(
-        grid_gen.interpolated_limits["lower_energy_threshold"][0][0][0], 0.007, rtol=1e-2
-    )
-    assert np.isclose(
-        grid_gen.interpolated_limits["upper_scatter_radius"][0][0][0], 925.0, rtol=1e-2
-    )
-    assert np.isclose(grid_gen.interpolated_limits["viewcone_radius"][0][0][0], 9.25, rtol=1e-2)
-
-    assert np.shape(grid_gen.interpolated_limits["lower_energy_threshold"]) == (2, 2, 1)
-    assert np.shape(grid_gen.interpolated_limits["upper_scatter_radius"]) == (2, 2, 1)
-    assert np.shape(grid_gen.interpolated_limits["viewcone_radius"]) == (2, 2, 1)
-
-
-def test_no_matching_rows_in_lookup_table(axes_definition, observing_location, observing_time):
-    """Test behavior when no matching rows are found in the lookup table."""
-    with pytest.raises(
-        ValueError, match=r"No matching rows in the lookup table for telescope_ids: \[999\]"
-    ):
-        GridGeneration(
-            axes=axes_definition,
-            coordinate_system="zenith_azimuth",
-            observing_location=observing_location,
-            observing_time=observing_time,
-            lookup_table="tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv",
-            telescope_ids=[999],
-        )
-
-
-def test_matching_rows_with_string_telescope_id(
-    axes_definition, lookup_table, observing_location, observing_time
-):
-    """Match legacy numeric lookup-table rows using string telescope IDs."""
-    grid_gen = GridGeneration(
-        axes=axes_definition,
-        coordinate_system="zenith_azimuth",
-        observing_location=observing_location,
-        observing_time=observing_time,
-        lookup_table=lookup_table,
-        telescope_ids=["LSTN-01"],
-    )
-
-    lookup_arrays = grid_gen._load_matching_lookup_arrays()
-    assert lookup_arrays["points"].shape[0] > 0
-
-
-def test_missing_observing_time(grid_gen):
-    """Test behavior when observing_time is not set."""
-    grid_gen.observing_time = None
-
-    with pytest.raises(ValueError, match="Observing time is not set"):
-        grid_gen.convert_altaz_to_radec(45 * u.deg, 30 * u.deg)
-
-
-def test_iers_not_modified_without_env(monkeypatch):
-    from simtools.application_control import _configure_iers_from_env
-
-    iers.conf.auto_download = True
-
-    monkeypatch.delenv("SIMTOOLS_OFFLINE_IERS", raising=False)
-
-    _configure_iers_from_env()
-
-    GridGeneration(axes={"axes": {}})
-
-    assert iers.conf.auto_download is True
-
-
-def test_iers_disabled_with_env(monkeypatch):
-    from simtools.application_control import _configure_iers_from_env
-
-    iers.conf.auto_download = True
-
-    monkeypatch.setenv("SIMTOOLS_OFFLINE_IERS", "1")
-
-    _configure_iers_from_env()
-
-    GridGeneration(axes={"axes": {}})
-
-    assert iers.conf.auto_download is False
diff --git a/tests/unit_tests/production_configuration/test_job_grid_io.py b/tests/unit_tests/production_configuration/test_job_grid_io.py
new file mode 100644
index 0000000000..3d6f3333a9
--- /dev/null
+++ b/tests/unit_tests/production_configuration/test_job_grid_io.py
@@ -0,0 +1,63 @@
+from pathlib import Path
+
+import astropy.units as u
+import pytest
+
+from simtools.production_configuration.job_grid_io import read_job_grid, serialize_job_grid
+
+
+def _job_rows():
+    return [
+        {
+            "primary": "gamma",
+            "azimuth_angle": 45 * u.deg,
+            "zenith_angle": 20 * u.deg,
+            "energy_min": 30 * u.GeV,
+            "energy_max": 10 * u.TeV,
+            "core_scatter_number": 10,
+            "core_scatter_max": 200 * u.m,
+            "view_cone_min": 0 * u.deg,
+            "view_cone_max": 5 * u.deg,
+            "nshow": 1000,
+            "model_version": "7.0.0",
+            "array_layout_name": "CTAO-North-Alpha",
+            "corsika_le_interaction": "urqmd",
+            "corsika_he_interaction": "epos",
+            "run_number": 10,
+        }
+    ]
+
+
+def _metadata():
+    return {
+        "site": "North",
+        "simulation_software": "corsika_sim_telarray",
+        "coordinate_system": "ra_dec",
+    }
+
+
+def test_serialize_and_read_job_grid_ecsv(tmp_test_directory):
+    output_file = Path(tmp_test_directory) / "job_grid.ecsv"
+
+    serialize_job_grid(_job_rows(), output_file, metadata=_metadata())
+    rows, metadata = read_job_grid(output_file)
+
+    assert metadata["site"] == "North"
+    assert rows[0]["energy_min"] == 30 * u.GeV
+    assert rows[0]["core_scatter_number"] == 10
+    assert rows[0]["array_layout_name"] == "CTAO-North-Alpha"
+
+
+def test_serialize_job_grid_rejects_non_ecsv_output(tmp_test_directory):
+    output_file = Path(tmp_test_directory) / "job_grid.txt"
+
+    with pytest.raises(ValueError, match="\\.ecsv"):
+        serialize_job_grid(_job_rows(), output_file, metadata=_metadata())
+
+
+def test_read_job_grid_rejects_non_ecsv_input(tmp_test_directory):
+    input_file = Path(tmp_test_directory) / "job_grid.txt"
+    input_file.write_text("dummy", encoding="utf-8")
+
+    with pytest.raises(ValueError, match="\\.ecsv"):
+        read_job_grid(input_file)
diff --git a/tests/unit_tests/production_configuration/test_observation_grid.py b/tests/unit_tests/production_configuration/test_observation_grid.py
new file mode 100644
index 0000000000..90b990dd6a
--- /dev/null
+++ b/tests/unit_tests/production_configuration/test_observation_grid.py
@@ -0,0 +1,360 @@
+from unittest.mock import Mock, patch
+
+import numpy as np
+import pytest
+from astropy import units as u
+from astropy.coordinates import EarthLocation, SkyCoord
+from astropy.tests.helper import assert_quantity_allclose
+from astropy.time import Time
+
+from simtools.production_configuration.observation_grid import ProductionGridEngine
+
+
+def test_generate_simulation_grid_keeps_horizontal_coordinates_for_radec_axes():
+    axes = {
+        "axes": {
+            "ra": {"range": [0, 0], "binning": 1, "scaling": "linear", "units": "deg"},
+            "dec": {"range": [0, 0], "binning": 1, "scaling": "linear", "units": "deg"},
+        }
+    }
+    engine = ProductionGridEngine(
+        axes=axes,
+        coordinate_system="ra_dec",
+        observing_location=EarthLocation(lat=28.76 * u.deg, lon=-17.89 * u.deg, height=2200 * u.m),
+        time_of_observation=Time("2017-09-16 00:00:00", scale="utc"),
+        lookup_table=None,
+    )
+
+    simulation_grid = engine.generate_simulation_grid()
+
+    assert "zenith_angle" in simulation_grid[0]
+    assert "azimuth" in simulation_grid[0]
+    assert "ra" in simulation_grid[0]
+    assert "dec" in simulation_grid[0]
+
+
+def test_require_time_of_observation_raises_without_time():
+    engine = ProductionGridEngine(axes={"axes": {}}, time_of_observation=None)
+
+    with pytest.raises(ValueError, match="Observing time"):
+        engine._require_time_of_observation()
+
+
+def test_get_max_zenith_for_radec_mode_reads_axis_range():
+    engine = ProductionGridEngine(
+        axes={"axes": {"zenith_angle": {"range": [10, 60], "binning": 2, "units": "deg"}}}
+    )
+
+    assert engine._get_max_zenith_for_radec_mode() == 60
+
+
+def test_get_max_zenith_for_radec_mode_raises_for_invalid_axis_definition():
+    engine = ProductionGridEngine(axes={"axes": {}})
+
+    with pytest.raises(ValueError, match="valid two-element 'range'"):
+        engine._get_max_zenith_for_radec_mode()
+
+
+def test_generate_target_values_supports_log_and_inverse_cos_scaling():
+    engine = ProductionGridEngine(
+        axes={
+            "axes": {
+                "energy": {"range": [1, 100], "binning": 3, "scaling": "log", "units": "TeV"},
+                "zenith_angle": {
+                    "range": [0, 60],
+                    "binning": 3,
+                    "scaling": "1/cos",
+                    "units": "deg",
+                },
+            }
+        }
+    )
+
+    assert_quantity_allclose(engine.target_values["energy"], [1, 10, 100] * u.TeV)
+    assert len(engine.target_values["zenith_angle"]) == 3
+
+
+def test_add_lookup_limits_to_point_uses_quantity_nsb_level():
+    engine = ProductionGridEngine(axes={"axes": {}}, lookup_table=None)
+    engine.lookup_table = "limits.ecsv"
+    engine._interpolate_limits_for_point = Mock(
+        return_value={
+            "lower_energy_threshold": 0.02,
+            "upper_scatter_radius": 120.0,
+            "viewcone_radius": 3.5,
+        }
+    )
+    point = {"nsb_level": 5 * u.dimensionless_unscaled}
+
+    engine._add_lookup_limits_to_point(point, zenith=20.0, azimuth=180.0)
+
+    engine._interpolate_limits_for_point.assert_called_once_with(
+        zenith=20.0, azimuth=180.0, nsb=5.0
+    )
+    assert_quantity_allclose(point["lower_energy_threshold"], 0.02 * u.TeV)
+    assert_quantity_allclose(point["scatter_radius"], 120 * u.m)
+    assert_quantity_allclose(point["viewcone_radius"], 3.5 * u.deg)
+
+
+def test_interpolate_limits_for_point_delegates_to_lookup_helper():
+    engine = ProductionGridEngine(axes={"axes": {}}, lookup_table=None)
+    engine._limits_lookup = Mock()
+    engine._limits_lookup.interpolate_point.return_value = {"lower_energy_threshold": 0.02}
+
+    limits = engine._interpolate_limits_for_point(20.0, 180.0, 1.0)
+
+    assert limits == {"lower_energy_threshold": 0.02}
+    engine._limits_lookup.interpolate_point.assert_called_once_with(20.0, 180.0, 1.0)
+
+
+def test_prepare_lookup_table_limits_for_point_interpolation_prepares():
+    engine = ProductionGridEngine(axes={"axes": {}}, lookup_table=None)
+    engine._limits_lookup = Mock()
+
+    engine._prepare_lookup_table_limits_for_point_interpolation()
+
+    engine._limits_lookup.prepare_point_interpolators.assert_called_once_with()
+
+
+def test_generate_extra_axis_combinations_returns_empty_placeholder_for_no_extra_axes():
+    engine = ProductionGridEngine(
+        axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}}
+    )
+
+    keys, units, combinations = engine._generate_extra_axis_combinations(("zenith_angle",))
+
+    assert keys == []
+    assert units == []
+    assert len(combinations) == 1
+
+
+def test_generate_extra_axis_combinations_returns_mesh_for_remaining_axes():
+    engine = ProductionGridEngine(
+        axes={
+            "axes": {
+                "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"},
+                "nsb_level": {"range": [1, 2], "binning": 2, "units": "1"},
+                "offset": {"range": [0, 1], "binning": 2, "units": "deg"},
+            }
+        }
+    )
+
+    keys, units, combinations = engine._generate_extra_axis_combinations(("zenith_angle",))
+
+    assert keys == ["nsb_level", "offset"]
+    assert units == [u.dimensionless_unscaled, u.deg]
+    assert combinations.shape == (4, 2)
+
+
+def test_create_circular_binning_uses_shorter_counterclockwise_path():
+    engine = ProductionGridEngine(axes={"axes": {}})
+
+    binning = engine.create_circular_binning((10, 350), 3)
+
+    assert np.allclose(binning, [10, 0, 350])
+
+
+def test_create_circular_binning_uses_clockwise_path_when_shorter():
+    engine = ProductionGridEngine(axes={"axes": {}})
+
+    binning = engine.create_circular_binning((350, 10), 3)
+
+    assert np.allclose(binning, [350, 0, 10])
+
+
+def test_convert_altaz_to_radec_raises_without_time_of_observation():
+    engine = ProductionGridEngine(axes={"axes": {}}, time_of_observation=None)
+
+    with pytest.raises(ValueError, match="time_of_observation"):
+        engine.convert_altaz_to_radec(45 * u.deg, 180 * u.deg)
+
+
+def test_convert_altaz_to_radec_returns_icrs_coordinates():
+    engine = ProductionGridEngine(
+        axes={"axes": {}},
+        observing_location=EarthLocation(lat=28.76 * u.deg, lon=-17.89 * u.deg, height=2200 * u.m),
+        time_of_observation=Time("2017-09-16 00:00:00", scale="utc"),
+    )
+
+    radec = engine.convert_altaz_to_radec(45 * u.deg, 180 * u.deg)
+
+    assert isinstance(radec, SkyCoord)
+    assert radec.frame.name == "icrs"
+
+
+def test_convert_coordinates_drops_horizontal_coordinates_when_requested():
+    engine = ProductionGridEngine(axes={"axes": {}}, coordinate_system="ra_dec")
+    engine.convert_altaz_to_radec = Mock(return_value=Mock(ra=Mock(deg=10), dec=Mock(deg=-20)))
+    points = [{"zenith_angle": 20 * u.deg, "azimuth": 180 * u.deg}]
+
+    converted = engine.convert_coordinates(points, keep_horizontal_coordinates=False)
+
+    assert "zenith_angle" not in converted[0]
+    assert "azimuth" not in converted[0]
+    assert_quantity_allclose(converted[0]["ra"], 10 * u.deg)
+    assert_quantity_allclose(converted[0]["dec"], -20 * u.deg)
+
+
+def test_convert_coordinates_keeps_points_without_horizontal_values():
+    engine = ProductionGridEngine(axes={"axes": {}}, coordinate_system="ra_dec")
+    points = [{"ra": 10 * u.deg, "dec": -20 * u.deg}]
+
+    converted = engine.convert_coordinates(points, keep_horizontal_coordinates=False)
+
+    assert converted == points
+
+
+@patch.object(ProductionGridEngine, "_prepare_lookup_table_limits_for_point_interpolation")
+@patch("simtools.production_configuration.observation_grid.CorsikaLimitsLookup")
+def test_init_with_radec_lookup_prepares_point_interpolation(
+    mock_corsika_limits_lookup,
+    mock_prepare_lookup_table_limits,
+):
+    ProductionGridEngine(
+        axes={"axes": {"ra": {"range": [0, 0], "binning": 1, "units": "deg"}}},
+        coordinate_system="ra_dec",
+        time_of_observation=Time("2017-09-16 00:00:00", scale="utc"),
+        lookup_table="limits.ecsv",
+        array_layout_name="alpha",
+    )
+
+    mock_corsika_limits_lookup.assert_called_once_with(
+        lookup_table="limits.ecsv",
+        array_layout_name="alpha",
+    )
+    mock_prepare_lookup_table_limits.assert_called_once_with()
+
+
+@patch("simtools.production_configuration.observation_grid.CorsikaLimitsLookup")
+def test_init_with_horizontal_lookup_applies_interpolated_limits(mock_corsika_limits_lookup):
+    mock_lookup = mock_corsika_limits_lookup.return_value
+    mock_lookup.interpolate_grid_limits.return_value = {
+        "lower_energy_threshold": np.array([[[0.02]]])
+    }
+
+    engine = ProductionGridEngine(
+        axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}},
+        lookup_table="limits.ecsv",
+        array_layout_name="alpha",
+    )
+
+    assert "lower_energy_threshold" in engine.interpolated_limits
+
+
+def test_generate_horizontal_grid_uses_interpolated_limit_arrays():
+    engine = ProductionGridEngine(
+        axes={
+            "axes": {
+                "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"},
+                "azimuth": {"range": [0, 0], "binning": 1, "units": "deg"},
+                "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"},
+            }
+        }
+    )
+    engine.interpolated_limits = {
+        "lower_energy_threshold": np.array([[[0.02]]]),
+        "upper_scatter_radius": np.array([[[120.0]]]),
+        "viewcone_radius": np.array([[[3.0]]]),
+    }
+
+    grid = engine._generate_horizontal_grid()
+
+    assert_quantity_allclose(grid[0]["lower_energy_threshold"], 0.02 * u.TeV)
+    assert_quantity_allclose(grid[0]["scatter_radius"], 120 * u.m)
+    assert_quantity_allclose(grid[0]["viewcone_radius"], 3 * u.deg)
+
+
+@patch("simtools.production_configuration.observation_grid.AltAz")
+@patch("simtools.production_configuration.observation_grid.SkyCoord")
+@patch("simtools.production_configuration.observation_grid.np.arange", return_value=np.array([0.0]))
+def test_generate_radec_grid_direction_points_filters_by_max_zenith(
+    mock_arange, mock_skycoord, mock_altaz
+):
+    engine = ProductionGridEngine(
+        axes={"axes": {"zenith_angle": {"range": [0, 20], "binning": 2, "units": "deg"}}},
+        observing_location=EarthLocation(lat=28.76 * u.deg, lon=-17.89 * u.deg, height=2200 * u.m),
+        time_of_observation=Mock(),
+    )
+    engine.time_of_observation.sidereal_time.return_value = Mock(deg=30.0)
+    mock_skycoord.return_value.transform_to.return_value = Mock(
+        alt=np.array([80.0, 10.0]) * u.deg,
+        az=Mock(deg=np.array([100.0, 200.0])),
+    )
+
+    direction_points = engine._generate_radec_grid_direction_points()
+
+    assert len(direction_points) == 1
+    assert_quantity_allclose(direction_points[0]["zenith_angle"], 10 * u.deg)
+    assert_quantity_allclose(direction_points[0]["azimuth"], 100 * u.deg)
+    mock_arange.assert_called_once()
+    mock_altaz.assert_called_once()
+
+
+def test_generate_horizontal_grid_handles_partial_interpolated_limit_arrays():
+    engine = ProductionGridEngine(
+        axes={
+            "axes": {
+                "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"},
+                "azimuth": {"range": [0, 0], "binning": 1, "units": "deg"},
+                "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"},
+            }
+        }
+    )
+    engine.interpolated_limits = {
+        "upper_scatter_radius": np.array([[[120.0]]]),
+        "viewcone_radius": np.array([[[3.0]]]),
+    }
+
+    grid = engine._generate_horizontal_grid()
+
+    assert "lower_energy_threshold" not in grid[0]
+    assert_quantity_allclose(grid[0]["scatter_radius"], 120 * u.m)
+    assert_quantity_allclose(grid[0]["viewcone_radius"], 3 * u.deg)
+
+
+def test_generate_grid_radec_mode_adds_extra_axis_quantities():
+    engine = ProductionGridEngine(
+        axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}}
+    )
+    engine.coordinate_system = "ra_dec"
+    engine._generate_radec_grid_direction_points = Mock(
+        return_value=[{"zenith_angle": 20 * u.deg, "azimuth": 180 * u.deg}]
+    )
+    engine._generate_extra_axis_combinations = Mock(
+        return_value=(["nsb_level"], [u.dimensionless_unscaled], [np.array([5.0])])
+    )
+    engine._add_lookup_limits_to_point = Mock()
+    engine.convert_coordinates = Mock(return_value=[{"ra": 10 * u.deg, "dec": -20 * u.deg}])
+
+    grid = engine._generate_grid_radec_mode()
+
+    assert grid == [{"ra": 10 * u.deg, "dec": -20 * u.deg}]
+    added_point = engine.convert_coordinates.call_args.args[0][0]
+    assert_quantity_allclose(added_point["nsb_level"], 5 * u.dimensionless_unscaled)
+
+
+def test_generate_grid_radec_mode_uses_direction_points_when_ra_dec_axes_missing():
+    engine = ProductionGridEngine(
+        axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1}}}
+    )
+    engine.coordinate_system = "ra_dec"
+    engine._generate_radec_grid_direction_points = Mock(
+        return_value=[{"zenith_angle": 20 * u.deg, "azimuth": 180 * u.deg}]
+    )
+    engine._generate_extra_axis_combinations = Mock(return_value=([], [], [np.array([])]))
+    engine.convert_coordinates = Mock(return_value=[{"ra": 10 * u.deg, "dec": -20 * u.deg}])
+
+    grid = engine._generate_grid_radec_mode()
+
+    assert grid == [{"ra": 10 * u.deg, "dec": -20 * u.deg}]
+    engine.convert_coordinates.assert_called_once()
+
+
+def test_generate_simulation_grid_uses_horizontal_grid_for_horizontal_mode():
+    engine = ProductionGridEngine(axes={"axes": {}})
+    engine.coordinate_system = "horizontal"
+    engine._generate_horizontal_grid = Mock(return_value=[{"azimuth": 0 * u.deg}])
+
+    grid = engine.generate_simulation_grid()
+
+    assert grid == [{"azimuth": 0 * u.deg}]
diff --git a/tests/unit_tests/production_configuration/test_simulation_jobs.py b/tests/unit_tests/production_configuration/test_simulation_jobs.py
new file mode 100644
index 0000000000..e21d4366ff
--- /dev/null
+++ b/tests/unit_tests/production_configuration/test_simulation_jobs.py
@@ -0,0 +1,655 @@
+from unittest.mock import Mock, patch
+
+import pytest
+from astropy import units as u
+from astropy.coordinates import EarthLocation
+from astropy.tests.helper import assert_quantity_allclose
+
+from simtools.production_configuration.corsika_limits_lookup import CorsikaLimitsLookup
+from simtools.production_configuration.simulation_jobs import (
+    _build_rows_for_point,
+    _clip_energy_range_from_threshold,
+    _clip_max_quantity,
+    _generate_observation_grids_per_layout,
+    _generate_observation_points_from_axes,
+    _resolve_nshow_params,
+    build_axes_dict_from_cli_args,
+    build_job_grid_metadata,
+    build_observing_location,
+    build_production_grid_engine,
+    build_simulation_jobs,
+    calculate_log_energy_midpoint,
+    calculate_scaled_nshow,
+    get_core_scatter_max_for_zenith_angle,
+    get_energy_range_for_zenith_angle,
+    get_viewcone_max_for_zenith_angle,
+    normalize_energy_ranges,
+    normalize_grid_axes,
+    resolve_single_model_version,
+    resolve_time_of_observation,
+)
+
+
+def test_resolve_single_model_version_uses_first_list_entry():
+    assert resolve_single_model_version(["7.0.0", "7.1.0"]) == "7.0.0"
+    assert resolve_single_model_version("7.0.0") == "7.0.0"
+
+
+def test_resolve_time_of_observation_returns_none_for_horizontal_without_input():
+    args_dict = {
+        "axis": [["azimuth", "0", "deg", "1", "deg", "2"], ["zenith", "0", "deg", "1", "deg", "2"]]
+    }
+    assert resolve_time_of_observation(None, args_dict) is None
+
+
+def test_build_axes_dict_from_cli_args_builds_horizontal_grid():
+    axes = build_axes_dict_from_cli_args(
+        {
+            "axis": [
+                ["azimuth", "310", "deg", "20", "deg", "3", "linear"],
+                ["zenith", "30", "deg", "40", "deg", "2", "linear"],
+                ["nsb", "4", "MHz", "5", "MHz", "2", "linear"],
+                ["offset", "0", "deg", "10", "deg", "2", "linear"],
+            ]
+        }
+    )
+
+    assert "azimuth" in axes
+    assert "zenith_angle" in axes
+    assert "ra" not in axes
+    assert "dec" not in axes
+
+
+def test_build_axes_dict_from_cli_args_accepts_compact_axis_definitions():
+    axes = build_axes_dict_from_cli_args(
+        {
+            "axis": [
+                ["azimuth", "310", "deg", "20", "deg", "3", "linear"],
+                ["zenith", "30", "deg", "40", "deg", "2"],
+                ["nsb", "4", "MHz", "5", "MHz", "2"],
+                ["offset", "0", "deg", "10", "deg", "2"],
+            ]
+        }
+    )
+
+    assert axes == {
+        "nsb_level": {"range": [4.0, 5.0], "binning": 2, "scaling": "linear", "units": "MHz"},
+        "offset": {"range": [0.0, 10.0], "binning": 2, "scaling": "linear", "units": "deg"},
+        "azimuth": {"range": [310.0, 20.0], "binning": 3, "scaling": "linear", "units": "deg"},
+        "zenith_angle": {
+            "range": [30.0, 40.0],
+            "binning": 2,
+            "scaling": "linear",
+            "units": "deg",
+        },
+    }
+
+
+def test_build_axes_dict_from_cli_args_accepts_config_wrapped_axis_strings():
+    axes = build_axes_dict_from_cli_args(
+        {
+            "axis": [
+                ["azimuth 310 deg 20 deg 3 linear"],
+                ["zenith 30 deg 40 deg 2"],
+                ["nsb 4 MHz 5 MHz 2"],
+                ["offset 0 deg 10 deg 2"],
+            ]
+        }
+    )
+
+    assert axes["azimuth"]["range"] == [310.0, 20.0]
+    assert axes["zenith_angle"]["binning"] == 2
+    assert axes["nsb_level"]["units"] == "MHz"
+
+
+def test_build_axes_dict_from_cli_args_accepts_merged_config_axis_list():
+    axes = build_axes_dict_from_cli_args(
+        {
+            "axis": [
+                [
+                    "azimuth 310 deg 20 deg 3 linear",
+                    "zenith 30 deg 40 deg 2",
+                    "nsb 4 MHz 5 MHz 2",
+                    "offset 0 deg 10 deg 2",
+                ]
+            ]
+        }
+    )
+
+    assert axes["azimuth"]["range"] == [310.0, 20.0]
+    assert axes["zenith_angle"]["binning"] == 2
+    assert axes["offset"]["range"] == [0.0, 10.0]
+
+
+def test_build_axes_dict_from_cli_args_rejects_multiple_direction_coordinate_systems():
+    with pytest.raises(ValueError, match="Cannot define both azimuth/zenith and ra/dec axes"):
+        build_axes_dict_from_cli_args(
+            {
+                "axis": [
+                    ["azimuth", "310", "deg", "20", "deg", "3"],
+                    ["zenith", "30", "deg", "40", "deg", "2"],
+                    ["ra", "0", "deg", "360", "deg", "36"],
+                    ["dec", "-90", "deg", "90", "deg", "18"],
+                    ["nsb", "4", "MHz", "5", "MHz", "2"],
+                    ["offset", "0", "deg", "10", "deg", "2"],
+                ]
+            }
+        )
+
+
+def test_resolve_time_of_observation_raises_for_radec_without_input():
+    args_dict = {
+        "axis": [["ra", "0", "deg", "1", "deg", "2"], ["dec", "0", "deg", "1", "deg", "2"]]
+    }
+    with pytest.raises(ValueError, match="time_of_observation"):
+        resolve_time_of_observation(None, args_dict)
+
+
+def test_build_job_grid_metadata_includes_job_context():
+    metadata = build_job_grid_metadata(
+        {
+            "site": "North",
+            "simulation_software": "corsika_sim_telarray",
+            "axis": [
+                ["ra", "0", "deg", "1", "deg", "2"],
+                ["dec", "0", "deg", "1", "deg", "2"],
+            ],
+            "time_of_observation": "2017-09-16 00:00:00",
+            "corsika_limits": "limits.ecsv",
+        }
+    )
+
+    assert metadata["site"] == "North"
+    assert metadata["simulation_software"] == "corsika_sim_telarray"
+    assert metadata["coordinate_system"] == "ra_dec"
+    assert metadata["time_of_observation_utc"].startswith("2017-09-16T00:00:00")
+    assert metadata["corsika_limits"] == "limits.ecsv"
+
+
+@patch("simtools.production_configuration.simulation_jobs.SiteModel")
+def test_build_observing_location_uses_site_model(mock_site_model):
+    model = mock_site_model.return_value
+    model.get_parameter_value_with_unit.side_effect = [1 * u.deg, 2 * u.deg, 3 * u.m]
+
+    location = build_observing_location("North", ["7.0.0"])
+
+    assert isinstance(location, EarthLocation)
+    mock_site_model.assert_called_once_with(model_version="7.0.0", site="North")
+
+
+@patch("simtools.production_configuration.simulation_jobs.ProductionGridEngine")
+def test_build_production_grid_engine_resolves_layout_name(
+    mock_production_grid_engine,
+):
+    args_dict = {
+        "array_layout_name": {"by_version": {"<7.0.0": "alpha", ">=7.0.0": "beta"}},
+        "model_version": ["7.0.0"],
+        "time_of_observation": None,
+        "corsika_limits": "limits.ecsv",
+        "axis": [
+            ["azimuth", "310", "deg", "20", "deg", "3", "linear"],
+            ["zenith", "30", "deg", "40", "deg", "2", "linear"],
+            ["nsb", "4", "MHz", "5", "MHz", "2", "linear"],
+            ["offset", "0", "deg", "10", "deg", "2", "linear"],
+        ],
+    }
+
+    build_production_grid_engine(args_dict)
+
+    mock_production_grid_engine.assert_called_once_with(
+        axes={
+            "nsb_level": {"range": [4.0, 5.0], "binning": 2, "scaling": "linear", "units": "MHz"},
+            "offset": {"range": [0.0, 10.0], "binning": 2, "scaling": "linear", "units": "deg"},
+            "azimuth": {
+                "range": [310.0, 20.0],
+                "binning": 3,
+                "scaling": "linear",
+                "units": "deg",
+            },
+            "zenith_angle": {
+                "range": [30.0, 40.0],
+                "binning": 2,
+                "scaling": "linear",
+                "units": "deg",
+            },
+        },
+        coordinate_system="horizontal",
+        observing_location=None,
+        time_of_observation=None,
+        lookup_table="limits.ecsv",
+        array_layout_name="beta",
+    )
+
+
+@patch("simtools.production_configuration.simulation_jobs.build_observing_location")
+@patch("simtools.production_configuration.simulation_jobs.ProductionGridEngine")
+def test_build_production_grid_engine_builds_observing_location_for_radec(
+    mock_production_grid_engine,
+    mock_build_observing_location,
+):
+    location = EarthLocation(lat=1 * u.deg, lon=2 * u.deg, height=3 * u.m)
+    mock_build_observing_location.return_value = location
+
+    build_production_grid_engine(
+        {
+            "site": "North",
+            "array_layout_name": "alpha",
+            "model_version": ["7.0.0"],
+            "time_of_observation": "2017-09-16 00:00:00",
+            "corsika_limits": None,
+            "axis": [
+                ["ra", "0", "deg", "360", "deg", "36", "linear"],
+                ["dec", "-90", "deg", "90", "deg", "18", "linear"],
+                ["nsb", "4", "MHz", "4", "MHz", "1", "linear"],
+                ["offset", "0", "deg", "10", "deg", "2", "linear"],
+            ],
+        }
+    )
+
+    mock_build_observing_location.assert_called_once_with(site="North", model_version=["7.0.0"])
+    assert mock_production_grid_engine.call_args.kwargs["observing_location"] == location
+
+
+def test_normalize_grid_axes_applies_defaults():
+    normalized = normalize_grid_axes(
+        {
+            "primary": "gamma",
+            "azimuth_angle": 180 * u.deg,
+            "zenith_angle": 20 * u.deg,
+            "model_version": "7.0.0",
+        }
+    )
+
+    assert normalized["primary"] == ["gamma"]
+    assert normalized["corsika_le_interaction"] == ["urqmd"]
+    assert normalized["corsika_he_interaction"] == ["epos"]
+
+
+@pytest.mark.parametrize(
+    ("energy_range", "expected_length"),
+    [
+        ((30 * u.GeV, 100 * u.GeV), 1),
+        ([[30 * u.GeV, 100 * u.GeV], [100 * u.GeV, 1 * u.TeV]], 2),
+    ],
+)
+def test_normalize_energy_ranges_accepts_supported_shapes(energy_range, expected_length):
+    assert len(normalize_energy_ranges(energy_range)) == expected_length
+
+
+def test_normalize_energy_ranges_rejects_invalid_shape():
+    with pytest.raises(ValueError, match="energy_range"):
+        normalize_energy_ranges([30 * u.GeV])
+
+
+def test_normalize_energy_ranges_accepts_quantity_pair_list():
+    normalized = normalize_energy_ranges([30 * u.GeV, 100 * u.GeV])
+
+    assert normalized == [(30 * u.GeV, 100 * u.GeV)]
+
+
+def test_get_energy_range_for_zenith_angle_without_lookup_returns_configured_pair():
+    energy_range = (30 * u.GeV, 100 * u.GeV)
+
+    assert get_energy_range_for_zenith_angle(20 * u.deg, energy_range, None) == energy_range
+
+
+@patch.object(CorsikaLimitsLookup, "__init__", return_value=None)
+@patch.object(
+    CorsikaLimitsLookup,
+    "interpolate_point",
+    return_value={
+        "lower_energy_threshold": 0.2,
+        "upper_scatter_radius": 100.0,
+        "viewcone_radius": 2.0,
+    },
+)
+def test_get_energy_range_for_zenith_angle_wraps_lookup_path_and_skips_step(
+    mock_interpolate_point,
+    mock_init,
+):
+
+    energy_range = get_energy_range_for_zenith_angle(
+        20 * u.deg,
+        (30 * u.GeV, 100 * u.GeV),
+        "limits.ecsv",
+    )
+
+    assert energy_range is None
+    mock_init.assert_called_once_with("limits.ecsv")
+    mock_interpolate_point.assert_called_once()
+
+
+@patch.object(CorsikaLimitsLookup, "__init__", return_value=None)
+@patch.object(
+    CorsikaLimitsLookup,
+    "interpolate_point",
+    return_value={
+        "lower_energy_threshold": 0.01,
+        "upper_scatter_radius": 100.0,
+        "viewcone_radius": 2.0,
+    },
+)
+def test_get_energy_range_for_zenith_angle_keeps_range_below_threshold(
+    mock_interpolate_point,
+    mock_init,
+):
+
+    energy_range = get_energy_range_for_zenith_angle(
+        20 * u.deg,
+        (30 * u.GeV, 100 * u.GeV),
+        "limits.ecsv",
+    )
+
+    assert energy_range == (30 * u.GeV, 100 * u.GeV)
+    mock_init.assert_called_once_with("limits.ecsv")
+    mock_interpolate_point.assert_called_once()
+
+
+def test_get_energy_range_for_zenith_angle_clips_threshold():
+    corsika_limits = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv"
+    )
+    corsika_limits.interpolate_point = Mock(
+        return_value={
+            "lower_energy_threshold": 0.05,
+            "upper_scatter_radius": 100.0,
+            "viewcone_radius": 2.0,
+        }
+    )
+
+    energy_range = get_energy_range_for_zenith_angle(
+        20 * u.deg,
+        (30 * u.GeV, 100 * u.GeV),
+        corsika_limits,
+    )
+
+    assert_quantity_allclose(energy_range[0], 50 * u.GeV)
+    assert_quantity_allclose(energy_range[1], 100 * u.GeV)
+
+
+def test_get_core_scatter_max_for_zenith_angle_clips_value():
+    corsika_limits = CorsikaLimitsLookup(
+        "tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv"
+    )
+    corsika_limits.interpolate_point = Mock(
+        return_value={
+            "lower_energy_threshold": 0.01,
+            "upper_scatter_radius": 100.0,
+            "viewcone_radius": 2.0,
+        }
+    )
+
+    scatter_max = get_core_scatter_max_for_zenith_angle(
+        20 * u.deg,
+        [10, 150 * u.m],
+        corsika_limits,
+    )
+
+    assert scatter_max == 100 * u.m
+
+
+def test_get_core_scatter_max_for_zenith_angle_without_lookup_uses_configured_max():
+    assert get_core_scatter_max_for_zenith_angle(20 * u.deg, [10, 150 * u.m], None) == 150 * u.m
+
+
+@patch.object(CorsikaLimitsLookup, "__init__", return_value=None)
+@patch.object(
+    CorsikaLimitsLookup,
+    "interpolate_point",
+    return_value={
+        "lower_energy_threshold": 0.01,
+        "upper_scatter_radius": 100.0,
+        "viewcone_radius": 2.0,
+    },
+)
+def test_get_core_scatter_max_for_zenith_angle_uses_lookup_path(mock_interpolate_point, mock_init):
+    scatter_max = get_core_scatter_max_for_zenith_angle(
+        20 * u.deg,
+        [10, 150 * u.m],
+        "limits.ecsv",
+    )
+
+    assert scatter_max == 100 * u.m
+    mock_init.assert_called_once_with("limits.ecsv")
+    mock_interpolate_point.assert_called_once()
+
+
+def test_get_viewcone_max_for_zenith_angle_without_lookup_uses_configured_max():
+    assert get_viewcone_max_for_zenith_angle(20 * u.deg, [0 * u.deg, 5 * u.deg], None) == 5 * u.deg
+
+
+@patch.object(CorsikaLimitsLookup, "__init__", return_value=None)
+@patch.object(
+    CorsikaLimitsLookup,
+    "interpolate_point",
+    return_value={
+        "lower_energy_threshold": 0.01,
+        "upper_scatter_radius": 100.0,
+        "viewcone_radius": 2.0,
+    },
+)
+def test_get_viewcone_max_for_zenith_angle_uses_lookup_path(mock_interpolate_point, mock_init):
+
+    viewcone_max = get_viewcone_max_for_zenith_angle(
+        20 * u.deg,
+        [0 * u.deg, 5 * u.deg],
+        "limits.ecsv",
+    )
+
+    assert viewcone_max == 2 * u.deg
+    mock_init.assert_called_once_with("limits.ecsv")
+    mock_interpolate_point.assert_called_once()
+
+
+def test_calculate_log_energy_midpoint_returns_geometric_mean():
+    midpoint = calculate_log_energy_midpoint((10 * u.GeV, 1 * u.TeV))
+
+    assert_quantity_allclose(midpoint, 100 * u.GeV)
+
+
+def test_calculate_log_energy_midpoint_raises_for_non_quantities():
+    with pytest.raises(TypeError, match="Quantity"):
+        calculate_log_energy_midpoint((10, 100 * u.GeV))
+
+
+def test_calculate_log_energy_midpoint_raises_for_non_positive_energy():
+    with pytest.raises(ValueError, match="strictly positive"):
+        calculate_log_energy_midpoint((0 * u.GeV, 100 * u.GeV))
+
+
+def test_calculate_scaled_nshow_returns_baseline_without_power_index():
+    assert calculate_scaled_nshow((30 * u.GeV, 100 * u.GeV), 5) == 5
+
+
+def test_calculate_scaled_nshow_raises_for_non_positive_baseline():
+    with pytest.raises(ValueError, match="positive integer"):
+        calculate_scaled_nshow((30 * u.GeV, 100 * u.GeV), 0)
+
+
+def test_calculate_scaled_nshow_raises_for_missing_reference_energy():
+    with pytest.raises(ValueError, match="reference_energy"):
+        calculate_scaled_nshow((30 * u.GeV, 100 * u.GeV), 5, nshow_power_index=1.0)
+
+
+def test_calculate_scaled_nshow_scales_from_midpoint_energy():
+    scaled_nshow = calculate_scaled_nshow(
+        (10 * u.GeV, 1 * u.TeV),
+        5,
+        nshow_power_index=1.0,
+        reference_energy=10 * u.GeV,
+    )
+
+    assert scaled_nshow == 50
+
+
+@patch("simtools.production_configuration.simulation_jobs.np.ceil", return_value=0)
+def test_calculate_scaled_nshow_raises_when_scaled_value_is_below_one(mock_ceil):
+    with pytest.raises(ValueError, match="at least 1"):
+        calculate_scaled_nshow(
+            (10 * u.GeV, 1 * u.TeV),
+            5,
+            nshow_power_index=1.0,
+            reference_energy=10 * u.GeV,
+        )
+    mock_ceil.assert_called_once()
+
+
+def test_clip_energy_range_from_threshold_returns_none_above_max():
+    assert (
+        _clip_energy_range_from_threshold(
+            (30 * u.GeV, 100 * u.GeV),
+            200 * u.GeV,
+        )
+        is None
+    )
+
+
+def test_clip_energy_range_from_threshold_returns_original_without_threshold():
+    energy_range = (30 * u.GeV, 100 * u.GeV)
+
+    assert _clip_energy_range_from_threshold(energy_range, None) == energy_range
+
+
+def test_clip_max_quantity_returns_configured_value_without_lookup():
+    assert _clip_max_quantity(5 * u.deg, None) == 5 * u.deg
+
+
+def test_resolve_nshow_params_converts_reference_energy():
+    nshow, power_index, reference_energy = _resolve_nshow_params(
+        {"nshow": 5, "nshow_power_index": 1.0, "nshow_reference_energy": "100 GeV"}
+    )
+
+    assert nshow == 5
+    assert power_index == pytest.approx(1.0)
+    assert reference_energy == 100 * u.GeV
+
+
+def test_build_rows_for_point_skips_energy_ranges_below_threshold():
+    rows = _build_rows_for_point(
+        point_base={"primary": "gamma"},
+        energy_ranges=[(30 * u.GeV, 40 * u.GeV), (50 * u.GeV, 100 * u.GeV)],
+        lower_energy_threshold=45 * u.GeV,
+        nshow=5,
+        nshow_power_index=None,
+        reference_energy=None,
+        number_of_runs=2,
+        run_number=10,
+    )
+
+    assert [row["run_number"] for row in rows] == [10, 11]
+    assert all(row["energy_min"] == 50 * u.GeV for row in rows)
+
+
+def test_generate_observation_points_from_axes_adds_lookup_limits():
+    corsika_limits = Mock()
+    corsika_limits.interpolate_point.return_value = {
+        "lower_energy_threshold": 0.05,
+        "upper_scatter_radius": 150.0,
+        "viewcone_radius": 3.0,
+    }
+
+    points = _generate_observation_points_from_axes([180 * u.deg], [20 * u.deg], corsika_limits)
+
+    assert len(points) == 1
+    assert_quantity_allclose(points[0]["lower_energy_threshold"], 0.05 * u.TeV)
+    assert_quantity_allclose(points[0]["scatter_radius"], 150 * u.m)
+    assert_quantity_allclose(points[0]["viewcone_radius"], 3 * u.deg)
+
+
+@patch("simtools.production_configuration.simulation_jobs._generate_observation_points_from_axes")
+@patch("simtools.production_configuration.simulation_jobs.CorsikaLimitsLookup")
+def test_generate_observation_grids_per_layout_uses_layout_specific_lookup(
+    mock_corsika_limits_lookup,
+    mock_generate_observation_points_from_axes,
+):
+    mock_generate_observation_points_from_axes.return_value = [{"azimuth": 0 * u.deg}]
+    args_dict = {
+        "array_layout_name": {"by_version": {"<7.0.0": "alpha", ">=7.0.0": "beta"}},
+        "corsika_limits": "limits.ecsv",
+    }
+    grid_axes = {
+        "model_version": ["6.3.0", "7.0.0"],
+        "azimuth_angle": [0 * u.deg],
+        "zenith_angle": [20 * u.deg],
+    }
+
+    observation_grids, resolved_layout_names = _generate_observation_grids_per_layout(
+        args_dict, grid_axes
+    )
+
+    assert set(observation_grids) == {"alpha", "beta"}
+    assert resolved_layout_names == {"6.3.0": "alpha", "7.0.0": "beta"}
+    mock_corsika_limits_lookup.assert_any_call("limits.ecsv", array_layout_name="alpha")
+    mock_corsika_limits_lookup.assert_any_call("limits.ecsv", array_layout_name="beta")
+
+
+@patch("simtools.production_configuration.simulation_jobs.build_production_grid_engine")
+def test_generate_observation_grids_per_layout_uses_shared_axes_and_skips_duplicate_layouts(
+    mock_build_production_grid_engine,
+):
+    mock_build_production_grid_engine.return_value.generate_simulation_grid.return_value = [
+        {"azimuth": 0 * u.deg}
+    ]
+
+    observation_grids, resolved_layout_names = _generate_observation_grids_per_layout(
+        {
+            "axis": [
+                ["azimuth", "310", "deg", "20", "deg", "3"],
+                ["zenith", "20", "deg", "40", "deg", "2"],
+                ["nsb", "4", "MHz", "5", "MHz", "2"],
+                ["offset", "0", "deg", "10", "deg", "2"],
+            ],
+            "array_layout_name": {"by_version": {"<7.0.0": "alpha", ">=7.0.0": "alpha"}},
+        },
+        {
+            "model_version": ["6.3.0", "7.0.0"],
+            "azimuth_angle": [0 * u.deg],
+            "zenith_angle": [20 * u.deg],
+        },
+    )
+
+    assert observation_grids == {"alpha": [{"azimuth": 0 * u.deg}]}
+    assert resolved_layout_names == {"6.3.0": "alpha", "7.0.0": "alpha"}
+    mock_build_production_grid_engine.assert_called_once()
+
+
+@patch("simtools.production_configuration.simulation_jobs._generate_observation_grids_per_layout")
+def test_build_simulation_jobs_expands_runs_from_observation_grid(
+    mock_generate_observation_grids_per_layout,
+):
+    mock_generate_observation_grids_per_layout.return_value = (
+        {
+            "alpha": [
+                {
+                    "azimuth": 180 * u.deg,
+                    "zenith_angle": 20 * u.deg,
+                    "lower_energy_threshold": 40 * u.GeV,
+                    "scatter_radius": 100 * u.m,
+                    "viewcone_radius": 2 * u.deg,
+                }
+            ]
+        },
+        {"6.3.0": "alpha"},
+    )
+    rows = build_simulation_jobs(
+        {
+            "primary": ["gamma"],
+            "azimuth_angle": [180 * u.deg],
+            "zenith_angle": [20 * u.deg],
+            "model_version": ["6.3.0"],
+            "corsika_le_interaction": "urqmd",
+            "corsika_he_interaction": "epos",
+            "energy_range": [(30 * u.GeV, 100 * u.GeV)],
+            "core_scatter": [10, 200 * u.m],
+            "view_cone": [0 * u.deg, 5 * u.deg],
+            "nshow": 5,
+            "number_of_runs": 2,
+            "run_number": 11,
+            "array_layout_name": "alpha",
+        }
+    )
+
+    assert [row["run_number"] for row in rows] == [11, 12]
+    assert rows[0]["array_layout_name"] == "alpha"
+    assert rows[0]["core_scatter_max"] == 100 * u.m
+    assert rows[0]["view_cone_max"] == 2 * u.deg
diff --git a/tests/unit_tests/visualization/test_plot_production_grid.py b/tests/unit_tests/visualization/test_plot_production_grid.py
index f89fb8f663..2591753870 100644
--- a/tests/unit_tests/visualization/test_plot_production_grid.py
+++ b/tests/unit_tests/visualization/test_plot_production_grid.py
@@ -60,7 +60,7 @@ def _write_grid_file(tmp_test_directory, file_name, grid_points):
         table["nsb_level"].unit = "MHz"
     if "offset" in table.colnames:
         table["offset"].unit = "deg"
-    table.meta["observing_time_utc"] = "2025-01-01T00:00:00.000"
+    table.meta["time_of_observation_utc"] = "2025-01-01T00:00:00.000"
     table.write(file_path, format="ascii.ecsv", overwrite=True)
     return file_path
 
@@ -250,7 +250,7 @@ def test_load_grid_points_from_ecsv(tmp_test_directory):
     table["azimuth"].unit = "deg"
     table["zenith_angle"].unit = "deg"
     table["nsb_level"].unit = "MHz"
-    table.meta["observing_time_utc"] = "2025-01-01T00:00:00.000"
+    table.meta["time_of_observation_utc"] = "2025-01-01T00:00:00.000"
     table.write(grid_file, format="ascii.ecsv", overwrite=True)
 
     plotter = _create_plotter(