Add survey area and depth-coverage API

.gitignore

@@ -24,3 +24,7 @@ ossssimlib
 bdir
 .vscode
 *.egg-info
+
+# Fortran local build products
+F95/libsrc.a
+F95/ossssimlib/

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.

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

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