use crate::coords::normalize_degrees;
use crate::error::{JyotishError, Result};
use crate::planet::Planet;
use crate::sun::solar_longitude;
use serde::{Deserialize, Serialize};
use std::fmt;
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
#[non_exhaustive]
pub enum Season {
VernalEquinox,
SummerSolstice,
AutumnalEquinox,
WinterSolstice,
}
impl Season {
pub fn target_longitude(self) -> f64 {
match self {
Self::VernalEquinox => 0.0,
Self::SummerSolstice => 90.0,
Self::AutumnalEquinox => 180.0,
Self::WinterSolstice => 270.0,
}
}
}
impl fmt::Display for Season {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
Self::VernalEquinox => write!(f, "Vernal Equinox"),
Self::SummerSolstice => write!(f, "Summer Solstice"),
Self::AutumnalEquinox => write!(f, "Autumnal Equinox"),
Self::WinterSolstice => write!(f, "Winter Solstice"),
}
}
}
#[derive(Debug, Clone, Copy)]
pub struct SynodicEvent {
pub body1: Planet,
pub body2: Planet,
pub jd: f64,
pub separation: f64,
pub is_opposition: bool,
}
pub fn next_season(season: Season, jd: f64) -> Result<f64> {
let target = season.target_longitude();
let current_lon = solar_longitude(jd);
let mut angle_ahead = target - current_lon;
if angle_ahead <= 0.0 {
angle_ahead += 360.0;
}
let estimate = jd + angle_ahead;
let mut lo = estimate - 20.0;
let mut hi = estimate + 20.0;
for _ in 0..80 {
let mid = (lo + hi) / 2.0;
let lon = solar_longitude(mid);
let diff = signed_angle_diff(lon, target);
if diff < 0.0 {
lo = mid;
} else {
hi = mid;
}
if (hi - lo) < 1e-8 {
break;
}
}
let result = (lo + hi) / 2.0;
if result < jd {
return next_season(season, jd + 300.0);
}
Ok(result)
}
pub fn seasons_of_year(jd: f64) -> Result<Vec<(Season, f64)>> {
let seasons = [
Season::VernalEquinox,
Season::SummerSolstice,
Season::AutumnalEquinox,
Season::WinterSolstice,
];
let start = jd - 180.0;
let mut results = Vec::with_capacity(4);
for season in seasons {
let event_jd = next_season(season, start)?;
if event_jd < jd + 366.0 {
results.push((season, event_jd));
}
}
results.sort_by(|a, b| a.1.total_cmp(&b.1));
Ok(results)
}
pub fn next_conjunction(body1: Planet, body2: Planet, jd: f64) -> Result<SynodicEvent> {
find_synodic_event(body1, body2, jd, false)
}
pub fn next_opposition(body1: Planet, body2: Planet, jd: f64) -> Result<SynodicEvent> {
find_synodic_event(body1, body2, jd, true)
}
fn longitude_of(planet: Planet, jd: f64) -> Result<f64> {
match planet {
Planet::Sun => Ok(solar_longitude(jd)),
Planet::Moon => Ok(crate::moon::lunar_longitude(jd)),
_ => Ok(crate::planetary::compute_position(planet, jd)?.longitude_deg),
}
}
fn find_synodic_event(
body1: Planet,
body2: Planet,
jd: f64,
opposition: bool,
) -> Result<SynodicEvent> {
let target_sep = if opposition { 180.0 } else { 0.0 };
let step = 5.0;
let max_steps = 200;
let mut jd_prev = jd;
let mut diff_prev = synodic_diff(body1, body2, jd_prev, target_sep)?;
for _ in 0..max_steps {
let jd_curr = jd_prev + step;
let diff_curr = synodic_diff(body1, body2, jd_curr, target_sep)?;
if diff_prev * diff_curr < 0.0 {
let mut lo = jd_prev;
let mut hi = jd_curr;
for _ in 0..60 {
let mid = (lo + hi) / 2.0;
let diff_mid = synodic_diff(body1, body2, mid, target_sep)?;
if diff_mid * diff_prev > 0.0 {
lo = mid;
} else {
hi = mid;
}
if (hi - lo) < 1e-8 {
break;
}
}
let event_jd = (lo + hi) / 2.0;
let lon1 = longitude_of(body1, event_jd)?;
let lon2 = longitude_of(body2, event_jd)?;
let sep = crate::aspect::angular_separation(lon1, lon2);
if (sep - target_sep).abs() < 10.0 {
return Ok(SynodicEvent {
body1,
body2,
jd: event_jd,
separation: sep,
is_opposition: opposition,
});
}
}
diff_prev = diff_curr;
jd_prev = jd_curr;
}
Err(JyotishError::InvalidParameter(
"synodic event search did not converge".into(),
))
}
fn synodic_diff(body1: Planet, body2: Planet, jd: f64, target_sep: f64) -> Result<f64> {
let lon1 = longitude_of(body1, jd)?;
let lon2 = longitude_of(body2, jd)?;
let elongation = normalize_degrees(lon1 - lon2);
let mut diff = elongation - target_sep;
if diff > 180.0 {
diff -= 360.0;
} else if diff < -180.0 {
diff += 360.0;
}
Ok(diff)
}
fn signed_angle_diff(lon: f64, target: f64) -> f64 {
let mut diff = lon - target;
if diff > 180.0 {
diff -= 360.0;
} else if diff < -180.0 {
diff += 360.0;
}
diff
}
#[cfg(test)]
mod tests {
use super::*;
const JD_J2000: f64 = 2_451_545.0;
#[test]
fn vernal_equinox_2000() {
let jd = next_season(Season::VernalEquinox, JD_J2000).unwrap();
assert!((jd - 2_451_623.8).abs() < 1.5, "got {jd}");
}
#[test]
fn summer_solstice_2000() {
let jd = next_season(Season::SummerSolstice, JD_J2000).unwrap();
assert!((jd - 2_451_716.0).abs() < 2.0, "got {jd}");
}
#[test]
fn all_four_seasons() {
let seasons = seasons_of_year(JD_J2000).unwrap();
assert_eq!(seasons.len(), 4);
for i in 0..3 {
assert!(seasons[i].1 < seasons[i + 1].1);
}
}
#[test]
fn season_display() {
assert_eq!(Season::VernalEquinox.to_string(), "Vernal Equinox");
assert_eq!(Season::WinterSolstice.to_string(), "Winter Solstice");
}
#[test]
fn season_serde() {
let season = Season::AutumnalEquinox;
let json = serde_json::to_string(&season).unwrap();
let restored: Season = serde_json::from_str(&json).unwrap();
assert_eq!(restored, season);
}
#[test]
fn conjunction_venus_jupiter() {
let event = next_conjunction(Planet::Venus, Planet::Jupiter, JD_J2000).unwrap();
assert!(event.jd > JD_J2000);
assert!(!event.is_opposition);
assert!(event.separation < 5.0, "separation = {}", event.separation);
}
#[test]
fn opposition_jupiter_sun() {
let event = next_opposition(Planet::Jupiter, Planet::Sun, JD_J2000).unwrap();
assert!(event.jd > JD_J2000);
assert!(event.is_opposition);
assert!(
(event.separation - 180.0).abs() < 10.0,
"separation = {}",
event.separation
);
}
}