Giant Meteor 2025

library(tidyr)
library(data.table)
library(ggplot2)
library(lubridate)

# Regarding polymarket
# https://polymarket.com/event/5kt-meteor-strike-in-2025/5kt-meteor-strike-in-2025?

#  https://cneos.jpl.nasa.gov/fireballs/
dt = fread("cneos_fireball_data.csv")

cols = c(date = "Peak Brightness Date/Time (UT)",
  kt = "Calculated Total Impact Energy (kt)")

setnames(dt, cols, names(cols))
dt[lubridate::yday(date)>163]
dt |>
  mutate(yday = lubridate::yday(date)) |>
  mutate(month = lubridate::month(date)) |>
  mutate(year = lubridate::year(date)) -> dt

# obs even over months
table(dt$month)
#  1  2  3  4  5  6  7  8  9 10 11 12
# 89 92 88 90 83 77 86 75 77 96 87 83

ggplot(data=dt, aes(x=month)) +
  geom_histogram(bins=12,fill=NA, color='black') +
  theme_minimal() + labs(title="Counts of any meteor by month")

# Obs sparse before ~1996 to 2000 ish... call it good in 2000
table(dt$year)
# 1988 1990 1991 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010
#    1    1    1    2   13   21   33   21   13   30   37   22   29   37   39   45   37   23   37   34   36
# 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025
#   31   37   24   33   45   33   29   39   42   41   32   41   34   32   18

today = as.POSIXct("2025-06-12")
today_day_of_year = lubridate::yday(today)
first_good_day = as.POSIXct("2000-01-01")
last_good_day = as.POSIXct("2025-06-01")

# Arrival fraction of 5kt meteors by year
ggplot(dt[, .(fract = .SD[kt>=5, .N]/.N), year],
       aes(x=year, y = fract)) +
  geom_point() + theme_minimal() +
  labs(title= "Arrival fraction of 5kt meteors by year")

# Counts by month seem a little mal distributed
ggplot(data=dt[kt>=5 & date >=first_good_day], aes(x=month)) +
  geom_histogram(bins=12,fill=NA, color='black') +
  theme_minimal() + labs(title="Counts of >=5kt by month")

# Are the monthly counts poisson? Sure. Why not
table(table(dt[kt>=5 & date >=first_good_day, month]))
table(dt[kt>=5 & date >=first_good_day, month])
# 1 2 4 5
# 5 3 2 1

arrival_rate_per_any_kt = dt[date >= first_good_day & date <= last_good_day, .SD[kt>=5, .N]/.N]

# 0.00258565
arrival_rate_per_day_5kt = dt[kt>=5 &
                                date >= first_good_day &
                                date <= last_good_day,
                              .N/as.integer(last_good_day - first_good_day)]
######################################
# At least one arrival left: 40.7%!
arrive_at_least_once = 1-(1-arrival_rate_per_day_5kt )^(365-today_day_of_year)
arrive_at_least_once