Commit 11f7f797 authored by yuedong0607's avatar yuedong0607
Browse files

fix the PST interpolation in bright star module refactorization

parent c6dcef7d

observation_sim/sky_background/native/README.md

+4 −0
Original line number Diff line number Diff line
@@ -10,6 +10,10 @@ All inputs are explicit arrays and all outputs are dimensionless weights. 
Catalog parsing, reference-file loading, spectral normalization, throughput

composition, and detector-rate integration remain in Python.



The versioned PST table is stored as base-10 logarithmic transmission. Both

the native and NumPy paths interpolate that logarithmic grid and exponentiate

the interpolated value, rather than interpolating linearized samples.



The setuptools build creates libstraylight_geometry.so. A standalone build is

also available with:

observation_sim/sky_background/native/geometry.py

+6 −3
Original line number Diff line number Diff line
@@ -11,12 +11,14 @@ import numpy as np 




def load_pst(path) -> np.ndarray:

    """Load the log10 PST grid and return linear transmission fractions."""

    """Load and retain the base-10 logarithmic PST grid."""



    values = np.loadtxt(path, dtype=np.float64)

    if values.ndim != 2 or values.shape[0] < 2 or values.shape[1] < 2:

        raise ValueError("PST must be a two-dimensional grid")

    return np.ascontiguousarray(np.power(10.0, values), dtype=np.float64)

    if not np.all(np.isfinite(values)):

        raise ValueError("PST grid contains non-finite values")

    return np.ascontiguousarray(values, dtype=np.float64)





def _library_path() -> Path | None:
@@ -97,7 +99,7 @@ def _numpy_directional_weights( 
    lower_angle = np.maximum(upper_angle - 1, 0)

    azimuth_fraction = azimuth_coordinate - lower_azimuth

    angle_fraction = angle_coordinate - upper_angle

    weights = (

    interpolated_log10 = (

        pst[lower_azimuth, lower_angle]

        * (1.0 - azimuth_fraction)

        * (1.0 - angle_fraction)
@@ -111,6 +113,7 @@ def _numpy_directional_weights( 
        * azimuth_fraction

        * angle_fraction

    )

    weights = np.power(10.0, interpolated_log10)

    weights[~valid] = 0.0

    return weights

observation_sim/sky_background/native/include/geometry.hpp

+6 −6
Original line number Diff line number Diff line
@@ -3,16 +3,16 @@ 


#include <cstddef>



// Evaluate the point-source-transmittance (PST) table for arbitrary source

// directions. The returned array contains one dimensionless PST weight per

// source. Directions outside the modeled field, inside the critical exclusion

// angle, or with invalid length receive zero weight.

// Interpolate the base-10 logarithmic point-source-transmittance (PST) table

// for arbitrary source directions and return one linear, dimensionless PST

// weight per source. Directions outside the modeled field, inside the critical

// exclusion angle, or with invalid length receive zero weight.

extern "C" std::size_t csst_directional_pst_weights(

    const double* source_directions,

    std::size_t source_count,

    const double* pointing_direction,

    const double* focal_plane_y_axis,

    const double* pst_grid,

    const double* log10_pst_grid,

    std::size_t pst_azimuth_count,

    std::size_t pst_angle_count,

    double critical_angle_rad,
@@ -26,7 +26,7 @@ extern "C" std::size_t csst_pst_weights( 
    const double* satellite_position_km,

    const double* pointing_direction,

    const double* focal_plane_y_axis,

    const double* pst_grid,

    const double* log10_pst_grid,

    std::size_t pst_azimuth_count,

    std::size_t pst_angle_count,

    double critical_angle_rad,

observation_sim/sky_background/native/src/geometry.cpp

+23 −19
Original line number Diff line number Diff line
@@ -21,9 +21,11 @@ double clamp_unit(double value) { 
    return std::max(-1.0, std::min(1.0, value));

}



// Interpolate the historical PST grid. Its azimuth samples are spaced by five

// degrees and its off-axis samples by two degrees, beginning at two degrees.

double interpolate_pst(const double* grid,

// Interpolate the historical base-10 logarithmic PST grid, then convert the

// interpolated value to a linear transmission. The azimuth samples are spaced

// by five degrees and the off-axis samples by two degrees, beginning at two

// degrees.

double interpolate_pst(const double* log10_grid,

                       std::size_t azimuth_count,

                       std::size_t angle_count,

                       double azimuth,
@@ -50,9 +52,10 @@ double interpolate_pst(const double* grid, 


    const auto value = [&](std::size_t azimuth_index,

                           std::size_t angle_index) {

        return grid[azimuth_index * angle_count + angle_index];

        return log10_grid[azimuth_index * angle_count + angle_index];

    };

    return value(lower_azimuth, lower_angle) * (1.0 - azimuth_fraction) *

    const double interpolated_log10 =

        value(lower_azimuth, lower_angle) * (1.0 - azimuth_fraction) *

            (1.0 - angle_fraction) +

        value(lower_azimuth, upper_angle) * (1.0 - azimuth_fraction) *

            angle_fraction +
@@ -60,19 +63,20 @@ double interpolate_pst(const double* grid, 
            (1.0 - angle_fraction) +

        value(upper_azimuth, upper_angle) * azimuth_fraction *

            angle_fraction;

    return std::pow(10.0, interpolated_log10);

}



std::size_t directional_weights(const double* source_directions,

                                std::size_t source_count,

                                const double* pointing_direction,

                                const double* focal_plane_y_axis,

                                const double* pst_grid,

                                const double* log10_pst_grid,

                                std::size_t pst_azimuth_count,

                                std::size_t pst_angle_count,

                                double critical_angle_rad,

                                double* weights) {

    if (source_directions == nullptr || pointing_direction == nullptr ||

        focal_plane_y_axis == nullptr || pst_grid == nullptr ||

        focal_plane_y_axis == nullptr || log10_pst_grid == nullptr ||

        weights == nullptr || pst_azimuth_count < 2 || pst_angle_count < 2) {

        return 0;

    }
@@ -119,7 +123,7 @@ std::size_t directional_weights(const double* source_directions, 
            dot3(projection, focal_plane_y_axis) /

            (projection_norm * y_axis_norm)));

        weights[index] = interpolate_pst(

            pst_grid,

            log10_pst_grid,

            pst_azimuth_count,

            pst_angle_count,

            azimuth,
@@ -136,7 +140,7 @@ extern "C" std::size_t csst_directional_pst_weights( 
    std::size_t source_count,

    const double* pointing_direction,

    const double* focal_plane_y_axis,

    const double* pst_grid,

    const double* log10_pst_grid,

    std::size_t pst_azimuth_count,

    std::size_t pst_angle_count,

    double critical_angle_rad,
@@ -145,7 +149,7 @@ extern "C" std::size_t csst_directional_pst_weights( 
                               source_count,

                               pointing_direction,

                               focal_plane_y_axis,

                               pst_grid,

                               log10_pst_grid,

                               pst_azimuth_count,

                               pst_angle_count,

                               critical_angle_rad,
@@ -158,7 +162,7 @@ extern "C" std::size_t csst_pst_weights( 
    const double* satellite_position_km,

    const double* pointing_direction,

    const double* focal_plane_y_axis,

    const double* pst_grid,

    const double* log10_pst_grid,

    std::size_t pst_azimuth_count,

    std::size_t pst_angle_count,

    double critical_angle_rad,
@@ -171,7 +175,7 @@ extern "C" std::size_t csst_pst_weights( 
                                                   source_count,

                                                   pointing_direction,

                                                   focal_plane_y_axis,

                                                   pst_grid,

                                                   log10_pst_grid,

                                                   pst_azimuth_count,

                                                   pst_angle_count,

                                                   critical_angle_rad,

tests/test_sky_background_refactor.py

+31 −0
Original line number Diff line number Diff line
@@ -97,6 +97,37 @@ def test_compiled_geometry_matches_numpy_fallback() -> None: 
    )





def test_pst_is_interpolated_in_log_space() -> None:

    log10_pst = np.full((73, 45), -8.0)

    log10_pst[0:2, 0] = -4.0

    pointing = np.asarray([1.0, 0.0, 0.0])

    y_axis = np.asarray([0.0, 1.0, 0.0])

    off_axis_angle = np.deg2rad(3.0)

    source = np.asarray(

        [[np.cos(off_axis_angle), np.sin(off_axis_angle), 0.0]]

    )



    weight = _numpy_directional_weights(

        source,

        pointing,

        y_axis,

        log10_pst,

        critical_angle_rad=0.0,

    )

    public_weight = compute_directional_pst_weights(

        source,

        pointing,

        y_axis,

        log10_pst,

        critical_angle_rad=0.0,

    )



    # Three degrees lies halfway between the two- and four-degree PST samples.

    # Interpolating -4 and -8 in log10 space gives a transmission of 1e-6.

    assert weight[0] == pytest.approx(1e-6, rel=1e-12)

    assert public_weight[0] == pytest.approx(1e-6, rel=1e-12)





def test_runtime_catalog_and_sed_backend_are_consistent() -> None:

    repository = InstrumentRepository.load_builtin("next-v1")

    backend = BrightStarSEDBackend.from_repository(repository)