From b77bf4b31e9aefb101e112e3fb28981268778e1a Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Sun, 6 Sep 2026 23:05:32 -0700 Subject: [PATCH] Add survey area and depth-coverage API via Fortran and Python. Expose polygon area, fill factor, and efficiency evaluation so callers can build effective area vs magnitude from characterization directories. Co-authored-by: Cursor --- .gitignore | 4 + F95/poly_lib.f95 | 94 ++++++++++++ F95/surveysub.f95 | 168 ++++++++++++++++++++++ F95/tests/Makefile | 22 ++- F95/tests/test_pointing_geom.f95 | 56 ++++++++ F95/tests/test_polygon_area.f95 | 70 +++++++++ src/ossssim/__init__.py | 2 + src/ossssim/survey.py | 231 ++++++++++++++++++++++++++++++ src/tests/test_survey_coverage.py | 79 ++++++++++ 9 files changed, 723 insertions(+), 3 deletions(-) create mode 100644 F95/tests/test_pointing_geom.f95 create mode 100644 F95/tests/test_polygon_area.f95 create mode 100644 src/ossssim/survey.py create mode 100644 src/tests/test_survey_coverage.py diff --git a/.gitignore b/.gitignore index a72300f..5dd06d9 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,7 @@ uv.lock bdir .vscode *.egg-info + +# Fortran local build products +F95/libsrc.a +F95/ossssimlib/ diff --git a/F95/poly_lib.f95 b/F95/poly_lib.f95 index 0ec0096..4355644 100644 --- a/F95/poly_lib.f95 +++ b/F95/poly_lib.f95 @@ -1,6 +1,7 @@ module poly_lib use poly_dec + use parameters contains !-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* @@ -16,6 +17,10 @@ module poly_lib ! which tells if the point "p" is inside, outside or touching the ! polygon "poly". ! +! Also provides: +! polygon_area (poly) +! spherical area of a simple RA/Dec footprint, in square degrees. +! !-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* integer function point_in_polygon(p, poly) @@ -179,6 +184,95 @@ integer function calc_walk_summand(p1, p2) end function calc_walk_summand + real (kind=8) function polygon_area(poly) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Spherical area of a simple closed footprint on the sky. +! +! Vertices are RA (x) and Dec (y) in radians, as stored after read_sur. +! The polygon is fan-triangulated from the first vertex; each spherical +! triangle area uses the van Oosterom & Strackee formula on unit vectors. +! Result is absolute area in square degrees. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! Version 1 : September 2026 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! poly : Polygon structure (polygon), RA/Dec [rad] +! +! OUTPUT +! polygon_area : area in square degrees (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) poly + implicit none + + type(t_polygon), intent(in) :: poly + integer :: i, n + real (kind=8) :: area_sr, tri + real (kind=8) :: v1(3), v2(3), v3(3) + + polygon_area = 0.d0 + n = poly%n + if (n .lt. 3) return + + call ra_dec_to_xyz(poly%x(1), poly%y(1), v1) + area_sr = 0.d0 + do i = 2, n - 1 + call ra_dec_to_xyz(poly%x(i), poly%y(i), v2) + call ra_dec_to_xyz(poly%x(i+1), poly%y(i+1), v3) + tri = spherical_triangle_area(v1, v2, v3) + area_sr = area_sr + tri + end do + polygon_area = dabs(area_sr)*(180.d0/Pi)**2 + return + end function polygon_area + + subroutine ra_dec_to_xyz(ra, dec, v) +! Convert RA/Dec [rad] to a unit vector. +!f2py intent(in) ra +!f2py intent(in) dec +!f2py intent(out) v + implicit none + real (kind=8), intent(in) :: ra, dec + real (kind=8), intent(out) :: v(3) + real (kind=8) :: cd + + cd = dcos(dec) + v(1) = cd*dcos(ra) + v(2) = cd*dsin(ra) + v(3) = dsin(dec) + return + end subroutine ra_dec_to_xyz + + real (kind=8) function spherical_triangle_area(a, b, c) +! Signed spherical triangle area [sr] for unit vectors a,b,c. +! van Oosterom & Strackee: 2*atan2(det, 1 + a·b + b·c + c·a) +!f2py intent(in) a +!f2py intent(in) b +!f2py intent(in) c + implicit none + real (kind=8), intent(in) :: a(3), b(3), c(3) + real (kind=8) :: det, denom + + det = a(1)*(b(2)*c(3) - b(3)*c(2)) & + + a(2)*(b(3)*c(1) - b(1)*c(3)) & + + a(3)*(b(1)*c(2) - b(2)*c(1)) + denom = 1.d0 + dot3(a, b) + dot3(b, c) + dot3(c, a) + if ((det .eq. 0.d0) .and. (denom .eq. 0.d0)) then + spherical_triangle_area = 0.d0 + else + spherical_triangle_area = 2.d0*datan2(det, denom) + end if + return + end function spherical_triangle_area + + real (kind=8) function dot3(u, v) + implicit none + real (kind=8), intent(in) :: u(3), v(3) + dot3 = u(1)*v(1) + u(2)*v(2) + u(3)*v(3) + return + end function dot3 + subroutine check_polygon(poly) !-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- ! NOTE: user must make sure the input array is at least (2,n+1) long. diff --git a/F95/surveysub.f95 b/F95/surveysub.f95 index f73c29a..36bce3a 100644 --- a/F95/surveysub.f95 +++ b/F95/surveysub.f95 @@ -10,6 +10,12 @@ module surveysub use getsur use ioutils + ! Saved survey characterization for library / Python access + logical, save :: survey_loaded = .false. + integer, save :: n_sur_loaded = 0 + type(t_pointing), save, private :: points_loaded(n_sur_max) + real (kind=8), save, private :: sur_mm_loaded(n_sur_max) + contains @@ -501,6 +507,168 @@ end subroutine Detos1 subroutine reset_simulator() first = .true. iff = 0 + survey_loaded = .false. + n_sur_loaded = 0 end subroutine reset_simulator + subroutine survey_load(survey, lun_s, n_sur, ierr) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Load a survey directory into module storage for index-based queries. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) survey +!f2py intent(in) lun_s +!f2py intent(out) n_sur +!f2py intent(out) ierr + implicit none + character(*), intent(in) :: survey + integer, intent(in) :: lun_s + integer, intent(out) :: n_sur, ierr + + call GetSurvey(survey, lun_s, n_sur_loaded, points_loaded, & + sur_mm_loaded, ierr) + if (ierr .eq. 0) then + survey_loaded = .true. + else + survey_loaded = .false. + n_sur_loaded = 0 + end if + n_sur = n_sur_loaded + return + end subroutine survey_load + + subroutine pointing_geom(idx, area_deg2, fill) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Geometric spherical area [deg^2] and fill factor for pointing idx (1-based). +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) idx +!f2py intent(out) area_deg2 +!f2py intent(out) fill + implicit none + integer, intent(in) :: idx + real (kind=8), intent(out) :: area_deg2, fill + + area_deg2 = 0.d0 + fill = 0.d0 + if ((.not. survey_loaded) .or. (idx .lt. 1) .or. (idx .gt. n_sur_loaded)) & + return + area_deg2 = polygon_area(points_loaded(idx)%poly) + fill = points_loaded(idx)%ff + return + end subroutine pointing_geom + + subroutine pointing_center(idx, ra, dec) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Approximate sky centre as the mean of polygon vertices [rad]. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) idx +!f2py intent(out) ra +!f2py intent(out) dec + implicit none + integer, intent(in) :: idx + real (kind=8), intent(out) :: ra, dec + integer :: i, n + + ra = 0.d0 + dec = 0.d0 + if ((.not. survey_loaded) .or. (idx .lt. 1) .or. (idx .gt. n_sur_loaded)) & + return + n = points_loaded(idx)%poly%n + if (n .le. 0) return + do i = 1, n + ra = ra + points_loaded(idx)%poly%x(i) + dec = dec + points_loaded(idx)%poly%y(i) + end do + ra = ra/dble(n) + dec = dec/dble(n) + return + end subroutine pointing_center + + subroutine pointing_meta(idx, efnam, epoch, code, mag_lim, rate_mid_asphr) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Metadata for pointing idx: efficiency filename, epoch [JD], obs code, +! limiting magnitude, and midpoint rate_cut ["/hr]. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) idx +!f2py intent(out) efnam +!f2py intent(out) epoch +!f2py intent(out) code +!f2py intent(out) mag_lim +!f2py intent(out) rate_mid_asphr + implicit none + integer, intent(in) :: idx + character(80), intent(out) :: efnam + real (kind=8), intent(out) :: epoch, mag_lim, rate_mid_asphr + integer, intent(out) :: code + real (kind=8) :: rmid + + efnam = ' ' + epoch = 0.d0 + code = 0 + mag_lim = 0.d0 + rate_mid_asphr = 0.d0 + if ((.not. survey_loaded) .or. (idx .lt. 1) .or. (idx .gt. n_sur_loaded)) & + return + efnam = points_loaded(idx)%efnam + epoch = points_loaded(idx)%o_pos(1)%jday + code = points_loaded(idx)%code + mag_lim = sur_mm_loaded(idx) + rmid = 0.5d0*(points_loaded(idx)%c%r_cut%min + & + points_loaded(idx)%c%r_cut%max) + rate_mid_asphr = rmid/drad*3600.d0/24.d0 + return + end subroutine pointing_meta + + real (kind=8) function pointing_eta(idx, mag, rate_asphr, maglim) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Detection efficiency at magnitude and on-sky rate ["/hr]. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) idx +!f2py intent(in) mag +!f2py intent(in) rate_asphr +!f2py intent(out) maglim + implicit none + integer, intent(in) :: idx + real (kind=8), intent(in) :: mag, rate_asphr + real (kind=8), intent(out) :: maglim + real (kind=8) :: rate + + pointing_eta = 0.d0 + maglim = 0.d0 + if ((.not. survey_loaded) .or. (idx .lt. 1) .or. (idx .gt. n_sur_loaded)) & + return + rate = rate_asphr*24.d0/3600.d0*drad + pointing_eta = eta(points_loaded(idx)%c%eff_p, points_loaded(idx)%c%nr, & + mag, rate, maglim) + return + end function pointing_eta + + subroutine pointing_eta_grid(idx, mags, n, rate_asphr, etas) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Vectorized efficiency evaluation on a magnitude grid. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) idx +!f2py intent(in) mags +!f2py intent(hide) n +!f2py intent(in) rate_asphr +!f2py intent(out) etas +!f2py depend(n) mags +!f2py depend(n) etas + implicit none + integer, intent(in) :: idx, n + real (kind=8), intent(in) :: mags(n), rate_asphr + real (kind=8), intent(out) :: etas(n) + integer :: i + real (kind=8) :: maglim, rate + + etas = 0.d0 + if ((.not. survey_loaded) .or. (idx .lt. 1) .or. (idx .gt. n_sur_loaded)) & + return + rate = rate_asphr*24.d0/3600.d0*drad + do i = 1, n + etas(i) = eta(points_loaded(idx)%c%eff_p, points_loaded(idx)%c%nr, & + mags(i), rate, maglim) + end do + return + end subroutine pointing_eta_grid + end module surveysub diff --git a/F95/tests/Makefile b/F95/tests/Makefile index f25ccdf..d1d5a83 100644 --- a/F95/tests/Makefile +++ b/F95/tests/Makefile @@ -4,7 +4,7 @@ LANGUAGE = F95 %: %.sh -test: clean $(MODELS) HST test_eta_rate +test: clean $(MODELS) HST test_eta_rate test_polygon_area test_pointing_geom $(MODELS): clean $(MAKE) -s -C ../ Driver GIMEOBJ=$@ @@ -24,11 +24,27 @@ test_eta_rate: \rm -f test_eta_rate $(MAKE) -s -C ../ clean -.PHONY: $(MODELS) test_eta_rate +# Unit test for spherical polygon_area +test_polygon_area: + $(MAKE) -s -C ../ libsrc.a + gfortran -O0 -I../ -o test_polygon_area test_polygon_area.f95 -L../ -lsrc + ./test_polygon_area + \rm -f test_polygon_area + $(MAKE) -s -C ../ rmtmp + +# Unit test for survey_load / pointing_geom against CFEPS +test_pointing_geom: + $(MAKE) -s -C ../ libsrc.a + gfortran -O0 -I../ -o test_pointing_geom test_pointing_geom.f95 -L../ -lsrc + ./test_pointing_geom + \rm -f test_pointing_geom + $(MAKE) -s -C ../ rmtmp + +.PHONY: $(MODELS) test_eta_rate test_polygon_area test_pointing_geom .DEFAULT: .PHONY: clean clean: - \rm -f ModelUsed.dat SimulDetect.dat SimulTrack.dat LOG test test_eta_rate + \rm -f ModelUsed.dat SimulDetect.dat SimulTrack.dat LOG test test_polygon_area test_pointing_geom diff --git a/F95/tests/test_pointing_geom.f95 b/F95/tests/test_pointing_geom.f95 new file mode 100644 index 0000000..25a1f34 --- /dev/null +++ b/F95/tests/test_pointing_geom.f95 @@ -0,0 +1,56 @@ +program test_pointing_geom +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Load CFEPS characterization and check pointing_geom fill factor for +! the first pointing (L3f-smooth.eff, fill=0.80). +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + + use surveysub + + implicit none + + integer :: n_sur, ierr, nfail + real (kind=8) :: area, fill + character(80) :: efnam + real (kind=8) :: epoch, mag_lim, rate_mid + integer :: code + + nfail = 0 + call reset_simulator() + call survey_load('../../src/ossssim/Characterizations/CFEPS', 21, n_sur, ierr) + if ((ierr .ne. 0) .or. (n_sur .le. 0)) then + write (6, *) 'FAIL: survey_load ierr=', ierr, ' n_sur=', n_sur + stop 1 + end if + write (6, *) 'PASS: loaded n_sur=', n_sur + + call pointing_geom(1, area, fill) + if (area .le. 0.d0) then + write (6, *) 'FAIL: area <= 0:', area + nfail = nfail + 1 + else + write (6, *) 'PASS: area=', area + end if + if (dabs(fill - 0.80d0) .gt. 1.d-6) then + write (6, *) 'FAIL: fill expected 0.80 got', fill + nfail = nfail + 1 + else + write (6, *) 'PASS: fill=', fill + end if + + call pointing_meta(1, efnam, epoch, code, mag_lim, rate_mid) + if (index(efnam, 'L3f-smooth') .le. 0) then + write (6, *) 'FAIL: unexpected efnam ', efnam + nfail = nfail + 1 + else + write (6, *) 'PASS: efnam=', efnam + end if + + if (nfail .eq. 0) then + write (6, *) 'ALL TESTS PASSED' + stop 0 + else + write (6, *) 'FAILURES:', nfail + stop 1 + end if + +end program test_pointing_geom diff --git a/F95/tests/test_polygon_area.f95 b/F95/tests/test_polygon_area.f95 new file mode 100644 index 0000000..4425a0a --- /dev/null +++ b/F95/tests/test_polygon_area.f95 @@ -0,0 +1,70 @@ +program test_polygon_area +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! Unit tests for polygon_area: a small equatorial rectangle should have +! area approximately width_deg * height_deg. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + + use parameters + use datadec + use poly_lib + use getsur + + implicit none + + type(t_polygon) :: poly + real (kind=8) :: w, h, ra, dec, area, expected, rel + integer :: nfail + + nfail = 0 + + ! 1 deg x 1 deg rectangle at Dec = 0 (half-width/height in radians) + w = 0.5d0*drad + h = 0.5d0*drad + ra = 0.d0 + dec = 0.d0 + call create_rectangle(w, h, ra, dec, poly) + call check_polygon(poly) + area = polygon_area(poly) + expected = 1.d0 + rel = dabs(area - expected)/expected + if (rel .gt. 1.d-3) then + write (6, *) 'FAIL: 1x1 deg rectangle area=', area, ' expected~', expected + nfail = nfail + 1 + else + write (6, *) 'PASS: 1x1 deg rectangle area=', area + end if + + ! 2 deg x 3 deg rectangle + w = 1.0d0*drad + h = 1.5d0*drad + call create_rectangle(w, h, ra, dec, poly) + call check_polygon(poly) + area = polygon_area(poly) + expected = 6.d0 + rel = dabs(area - expected)/expected + if (rel .gt. 1.d-3) then + write (6, *) 'FAIL: 2x3 deg rectangle area=', area, ' expected~', expected + nfail = nfail + 1 + else + write (6, *) 'PASS: 2x3 deg rectangle area=', area + end if + + ! Degenerate polygon + poly%n = 2 + area = polygon_area(poly) + if (area .ne. 0.d0) then + write (6, *) 'FAIL: n<3 should give 0, got', area + nfail = nfail + 1 + else + write (6, *) 'PASS: n<3 returns 0' + end if + + if (nfail .eq. 0) then + write (6, *) 'ALL TESTS PASSED' + stop 0 + else + write (6, *) 'FAILURES: ', nfail + stop 1 + end if + +end program test_polygon_area diff --git a/src/ossssim/__init__.py b/src/ossssim/__init__.py index 0bacb95..77da66e 100644 --- a/src/ossssim/__init__.py +++ b/src/ossssim/__init__.py @@ -14,3 +14,5 @@ from .models import DetectFile, ModelFile, ModelOutputFile, ModelFileEcsv, ModelFileOld from .ephem import Ephem from .characterization import * +from .survey import SurveyCharacterization, Pointing +from . import Characterizations diff --git a/src/ossssim/survey.py b/src/ossssim/survey.py new file mode 100644 index 0000000..af86462 --- /dev/null +++ b/src/ossssim/survey.py @@ -0,0 +1,231 @@ +""" +Survey characterization access: areas, fill factors, and efficiency vs depth. + +Example +------- +>>> from ossssim import SurveyCharacterization, Characterizations +>>> import numpy as np +>>> survey = SurveyCharacterization.from_directory( +... Characterizations.surveys['CFEPS']) +>>> mags = np.arange(21.0, 26.0, 0.1) +>>> m, area = survey.coverage_vs_magnitude(mags) +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path +from typing import Dict, Iterable, List, Optional, Sequence, Tuple, Union +import numpy as np + +import ossssimlib + +MagLike = Union[float, np.ndarray, Sequence[float]] + + +def _as_str(value) -> str: + if isinstance(value, bytes): + return value.decode('utf-8', errors='ignore').strip() + return str(value).strip() + + +@dataclass +class Pointing: + """One survey pointing with geometric area, fill factor, and efficiency.""" + + id: str + index: int # 1-based Fortran index into the loaded survey + efnam: str + ra: float # radians + dec: float # radians + epoch: float # JD + area_deg2: float + fill_factor: float + mag_lim: float + obs_code: int + rate_mid_asphr: float + _lib: object = field(repr=False, default=None) + + def efficiency( + self, + mag: MagLike, + rate_asphr: Optional[float] = None, + ) -> Union[float, np.ndarray]: + """Evaluate η at magnitude(s). Default rate is the pointing rate_cut midpoint.""" + rate = self.rate_mid_asphr if rate_asphr is None else float(rate_asphr) + lib = self._lib or ossssimlib.surveysub + arr = np.asarray(mag, dtype=float) + if arr.ndim == 0: + eta, _maglim = lib.pointing_eta(self.index, float(arr), rate) + return float(eta) + etas = np.empty(arr.size, dtype=float) + lib.pointing_eta_grid( + self.index, arr.ravel(), arr.size, rate, etas + ) + return etas.reshape(arr.shape) + + def effective_area( + self, + mag: MagLike, + rate_asphr: Optional[float] = None, + include_fill: bool = True, + ) -> Union[float, np.ndarray]: + """Return area * fill * η(m) (fill omitted when include_fill is False).""" + eta = self.efficiency(mag, rate_asphr=rate_asphr) + scale = self.area_deg2 * (self.fill_factor if include_fill else 1.0) + return scale * eta + + +class SurveyCharacterization: + """ + Loaded survey characterization directory with keyed pointing access. + + Pointing keys use the efficiency-file basename without ``.eff``. When + several pointings share the same efficiency file they are disambiguated + as ``name#0``, ``name#1``, ... in ``pointings.list`` order. Use + ``by_index`` for unambiguous iteration. + """ + + def __init__( + self, + directory: Path, + pointings: Dict[str, Pointing], + by_index: List[Pointing], + ): + self.directory = Path(directory) + self.pointings = pointings + self.by_index = by_index + + @classmethod + def from_directory( + cls, + directory: Union[str, Path], + lun: int = 21, + ) -> "SurveyCharacterization": + """Load ``pointings.list`` and associated ``.eff`` files via Fortran.""" + directory = Path(directory).resolve() + if not (directory / 'pointings.list').is_file(): + raise FileNotFoundError(f"No pointings.list in {directory}") + + lib = ossssimlib.surveysub + lib.reset_simulator() + n_sur, ierr = lib.survey_load(str(directory), lun) + if ierr != 0 or n_sur <= 0: + raise RuntimeError( + f"survey_load failed for {directory}: ierr={ierr}, n_sur={n_sur}" + ) + + # First pass: collect raw metadata and count duplicate efnam stems + raw = [] + counts: Dict[str, int] = {} + for i in range(1, n_sur + 1): + area, fill = lib.pointing_geom(i) + ra, dec = lib.pointing_center(i) + efnam, epoch, code, mag_lim, rate_mid = lib.pointing_meta(i) + efnam_s = _as_str(efnam) + stem = Path(efnam_s).name + if stem.lower().endswith('.eff'): + stem = stem[:-4] + counts[stem] = counts.get(stem, 0) + 1 + raw.append( + dict( + index=i, + efnam=efnam_s, + stem=stem, + ra=float(ra), + dec=float(dec), + epoch=float(epoch), + area_deg2=float(area), + fill_factor=float(fill), + mag_lim=float(mag_lim), + obs_code=int(code), + rate_mid_asphr=float(rate_mid), + ) + ) + + # Second pass: assign ids (plain stem if unique, else stem#k) + seen: Dict[str, int] = {} + by_index: List[Pointing] = [] + pointings: Dict[str, Pointing] = {} + for item in raw: + stem = item['stem'] + if counts[stem] == 1: + pid = stem + else: + k = seen.get(stem, 0) + seen[stem] = k + 1 + pid = f"{stem}#{k}" + p = Pointing( + id=pid, + index=item['index'], + efnam=item['efnam'], + ra=item['ra'], + dec=item['dec'], + epoch=item['epoch'], + area_deg2=item['area_deg2'], + fill_factor=item['fill_factor'], + mag_lim=item['mag_lim'], + obs_code=item['obs_code'], + rate_mid_asphr=item['rate_mid_asphr'], + _lib=lib, + ) + by_index.append(p) + pointings[pid] = p + + return cls(directory, pointings, by_index) + + def keys(self) -> List[str]: + return list(self.pointings.keys()) + + def __getitem__(self, key: Union[str, int]) -> Pointing: + if isinstance(key, int): + return self.by_index[key] + return self.pointings[key] + + def __iter__(self): + return iter(self.by_index) + + def __len__(self) -> int: + return len(self.by_index) + + def coverage_vs_magnitude( + self, + mags: MagLike, + rate_asphr: Optional[float] = None, + pointing_ids: Optional[Iterable[str]] = None, + include_fill: bool = True, + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Effective survey area as a function of magnitude. + + A_eff(m) = sum_i area_i * fill_i * eta_i(m, rate_i) + + Parameters + ---------- + mags : + Magnitude grid. + rate_asphr : + On-sky rate ["/hr] applied to every pointing. When None, each + pointing uses its own rate_cut midpoint. + pointing_ids : + Subset of pointing keys; default is all pointings. + include_fill : + If False, omit fill factors from the sum. + + Returns + ------- + mags, A_eff : ndarray + Magnitude grid and effective area [deg^2]. + """ + m = np.asarray(mags, dtype=float) + if pointing_ids is None: + selected = self.by_index + else: + selected = [self.pointings[k] for k in pointing_ids] + + total = np.zeros(m.shape, dtype=float) + for p in selected: + total = total + np.asarray( + p.effective_area(m, rate_asphr=rate_asphr, include_fill=include_fill), + dtype=float, + ) + return m, total diff --git a/src/tests/test_survey_coverage.py b/src/tests/test_survey_coverage.py new file mode 100644 index 0000000..0bc61f6 --- /dev/null +++ b/src/tests/test_survey_coverage.py @@ -0,0 +1,79 @@ +"""Tests for SurveyCharacterization area / fill / depth coverage API.""" +import pathlib +import unittest + +import numpy as np + +from ossssim import Characterizations, SurveyCharacterization + +script_directory = pathlib.Path(__file__).parent.resolve() +CFEPS = Characterizations.surveys['CFEPS'] + + +class SurveyCharacterizationTest(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.survey = SurveyCharacterization.from_directory(CFEPS) + + def test_loads_pointings(self): + self.assertGreater(len(self.survey), 0) + self.assertIn('L3f-smooth', self.survey.keys()) + self.assertIn('L3h-smooth', self.survey.pointings) + + def test_duplicate_eff_keys_disambiguated(self): + # L3w has two pointings sharing L3w-smooth.eff + keys = self.survey.keys() + self.assertTrue( + any(k.startswith('L3w-smooth') for k in keys), + msg=f'expected L3w-smooth keys among {keys[:10]}...', + ) + l3w = [k for k in keys if k.startswith('L3w-smooth')] + self.assertGreaterEqual(len(l3w), 2) + self.assertIn('L3w-smooth#0', l3w) + self.assertIn('L3w-smooth#1', l3w) + + def test_area_and_fill(self): + p = self.survey['L3f-smooth'] + self.assertGreater(p.area_deg2, 0.0) + self.assertGreater(p.fill_factor, 0.0) + self.assertLessEqual(p.fill_factor, 1.0) + self.assertAlmostEqual(p.fill_factor, 0.80, places=2) + + def test_efficiency_and_effective_area(self): + p = self.survey['L3h-smooth'] + eta = p.efficiency(23.0) + self.assertGreater(eta, 0.0) + self.assertLessEqual(eta, 1.0) + a_eff = p.effective_area(23.0) + self.assertAlmostEqual( + a_eff, p.area_deg2 * p.fill_factor * eta, places=10 + ) + + def test_coverage_vs_magnitude_decreases(self): + mags = np.arange(21.0, 26.01, 0.5) + m, area = self.survey.coverage_vs_magnitude(mags) + self.assertEqual(list(m), list(mags)) + self.assertTrue(np.all(area >= 0.0)) + # Fainter end should not be larger than bright end for survey coverage + self.assertLessEqual(area[-1], area[0] + 1e-9) + # Check consistency with single-pointing sum at one mag + mid = 23.0 + _, a_mid = self.survey.coverage_vs_magnitude([mid]) + expected = sum(p.effective_area(mid) for p in self.survey) + self.assertAlmostEqual(float(a_mid[0]), float(expected), places=8) + + def test_subset_pointing_ids(self): + mags = np.array([22.0, 24.0]) + _, a_all = self.survey.coverage_vs_magnitude(mags) + _, a_one = self.survey.coverage_vs_magnitude( + mags, pointing_ids=['L3f-smooth'] + ) + self.assertTrue(np.all(a_one <= a_all + 1e-9)) + p = self.survey['L3f-smooth'] + np.testing.assert_allclose( + a_one, p.effective_area(mags), rtol=1e-10 + ) + + +if __name__ == '__main__': + unittest.main()