From 26a47c2584b384f2d966c5970daa27732d399043 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Luis=20Eduardo=20Ram=C3=ADrez=20Montoya?= Date: Sat, 12 Aug 2023 11:33:32 -0600 Subject: [PATCH] Change reduce from Date to TimeOfDay --- ext/DataFramesExt.jl | 17 +++---- src/observations/topocentric.jl | 80 +++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 8 deletions(-) diff --git a/ext/DataFramesExt.jl b/ext/DataFramesExt.jl index d1f07cd2..c35cd876 100644 --- a/ext/DataFramesExt.jl +++ b/ext/DataFramesExt.jl @@ -1,10 +1,11 @@ module DataFramesExt -using Dates: Date, Period, Hour, Day, datetime2julian, julian2datetime +using Dates: Period, Hour, Day, datetime2julian, julian2datetime using TaylorSeries: get_numvars using PlanetaryEphemeris: J2000, selecteph, su, ea, yr, daysec, auday2kmsec using NEOs: RadecMPC, date, gauss_triplets, propagate, RNp1BP_pN_A_J23E_J2S_eph_threads!, order, abstol, sseph, - scaled_variables, gauss_method, residuals, bwdfwdeph, newtonls, diffcorr, nrms, hascoord, tryls + scaled_variables, gauss_method, residuals, bwdfwdeph, newtonls, diffcorr, nrms, hascoord, tryls, + TimeOfDay import Base: convert import NEOs: AbstractAstrometry, extrapolation, reduce_nights, gaussinitcond, relax_factor @@ -72,9 +73,9 @@ via polynomial interpolation. function reduce_nights(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} # Convert to DataFrame df = DataFrame(radec) - # Group by observatory and Date - df.Date = Date.(df.date) - gdf = groupby(df, [:observatory, :Date]) + # Group by observatory and TimeOfDay + df.TimeOfDay = TimeOfDay.(df.date, df.observatory) + gdf = groupby(df, [:observatory, :TimeOfDay]) # Interpolate observation nights cdf = combine(gdf, extrapolation, keepkeys = false) # Eliminate unsuccesful interpolations @@ -89,9 +90,9 @@ end function relax_factor(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} # Convert to DataFrame df = DataFrame(radec) - # Group by observatory and Date - df.Date = Date.(df.date) - gdf = groupby(df, [:observatory, :Date]) + # Group by observatory and TimeOfDay + df.TimeOfDay = TimeOfDay.(df.date, df.observatory) + gdf = groupby(df, [:observatory, :TimeOfDay]) # Interpolate observation nights cdf = combine(gdf, nrow) # Count observations in each group diff --git a/src/observations/topocentric.jl b/src/observations/topocentric.jl index b9d5bb3e..e3d193cd 100644 --- a/src/observations/topocentric.jl +++ b/src/observations/topocentric.jl @@ -1,6 +1,86 @@ # Earth orientation parameters (eop) 2000 const eop_IAU2000A::EopIau2000A = fetch_iers_eop(Val(:IAU2000A)) +@auto_hash_equals struct TimeOfDay + light::Symbol + start::Date + stop::Date + utc::Int + function TimeOfDay(date::DateTime, observatory::ObservatoryMPC{T}) where {T <: AbstractFloat} + # Hours from UTC + utc = hours_from_UTC(observatory) + # Today's sunrise / sunset + today = sunriseset(date, observatory) + # Yesterday's sunrise / sunset + yesterday = sunriseset(date - Day(1), observatory) + # Tomorrow's sunrise / sunset + tomorrow = sunriseset(date + Day(1), observatory) + # Selection + if yesterday[2] <= date <= today[1] + return new(:night, yesterday[2], today[1], utc) + elseif today[1] <= date <= today[2] + return new(:day, today[1], today[2], utc) + elseif today[2] <= date <= tomorrow[1] + return new(:night, today[2], tomorrow[1], utc) + end + + end +end + +isday(x::TimeOfDay) = x.light == :day +isnight(x::TimeOfDay) = x.light == :night + +# Print method for RadecMPC +# Examples: +# N00hp15 α: 608995.65 δ: -25653.3 t: 2020-12-04T10:41:43.209 obs: 703 +# 99942 α: 422475.3 δ: 97289.49 t: 2021-05-12T06:28:35.904 obs: F51 +function show(io::IO, m::TimeOfDay) + print(io, uppercasefirst(string(m.light)), " from ", m.start, " to ", m.stop, " at UTC", @sprintf("%+d", m.utc)) +end + +hours_from_UTC(lon::T) where {T <: Number} = ceil(Int, 12*lon/π - 0.5) +function hours_from_UTC(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat} + lon, _ = lonlat(observatory) + return hours_from_UTC(lon) +end + +function lonlat(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat} + # ECEF [km] + p_ECEF = obs_pos_ECEF(observatory) + # ECEF [m] -> Geodetic [m] + lat_geodetic, lon, altitude = ecef_to_geodetic(1_000 * p_ECEF) + # Geodetic [m] -> Geocentric [m] + lat_geocentric, _ = geodetic_to_geocentric(lat_geodetic, altitude) + + return lon, lat_geocentric +end + +function sunriseset(date::DateTime, observatory::ObservatoryMPC{T}) where {T <: AbstractFloat} + # Fractional year [rad] + γ = 2π * (dayofyear(date) - 1 + (hour(date)-12)/24) / daysinyear(date) + # Equation of time [min] + eqtime = 229.18 * (0.000075 + 0.001868*cos(γ) - 0.032077*sin(γ) - 0.014615*cos(2γ) - 0.040849*sin(2γ) ) + # Solar declination angle [rad] + δ = 0.006918 - 0.399912*cos(γ) + 0.070257*sin(γ) - 0.006758*cos(2γ) + 0.000907*sin(2γ) - 0.002697*cos(3γ) + 0.00148*sin(3γ) + # Longitude and latitude [rad] + lon, lat = lonlat(observatory) + # Solar hour angle [deg] + ha_sunrise = acosd(cosd(90.833)/cos(lat)/cos(δ) - tan(lat)*tan(δ)) + ha_sunset = -ha_sunrise + # Day [DateTime] + day = DateTime(Date(date)) + # UTC time of sunrise [min] + sunrise_min = 720 - 4*(rad2deg(lon) + ha_sunrise) - eqtime + # UTC time of sunrise [DateTime] + sunrise = day + Microsecond(round(Int, sunrise_min*6e7)) + # UTC time of sunset [min] + sunset_min = 720 - 4*(rad2deg(lon) + ha_sunset) - eqtime + # UTC time of sunset [DateTime] + sunset = day + Microsecond(round(Int, sunset_min*6e7)) + + return sunrise, sunset +end + @doc raw""" obs_pos_ECEF(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat} obs_pos_ECEF(x::RadecMPC{T}) where {T <: AbstractFloat}