prayers/
astronomy.rs

1use crate::dmath;
2use chrono::{Date, Datelike, Utc};
3
4/// Latitude, Longitude, Altitude (default to 0, in meters)
5///
6/// # Example
7/// ~~~~
8/// use prayers::*;
9///
10/// Coordinates(46.0, 69.0, 0.0);
11/// Coordinates(46.0, 69.0, 25.0);
12/// ~~~~
13#[derive(PartialEq, Debug, Copy, Clone)]
14pub struct Coordinates(pub f64, pub f64, pub f64);
15
16pub fn get_julian_day(date: &Date<Utc>) -> f64 {
17	let mut year = date.year() as f64;
18	let mut month = date.month() as f64;
19	let day = date.day() as f64;
20
21	if month < 3.0 {
22		year -= 1.0;
23		month += 12.0;
24	}
25
26	(365.2425 * year + 30.6001 * month).floor() + day + 1721027.5
27}
28
29pub fn mid_day(julian_date: f64, time: f64) -> f64 {
30	let eqt = sun_position(julian_date + time).1;
31	dmath::fix_hour(12.0 - eqt)
32}
33
34pub fn sun_angle_time(julian_date: f64, latitude: f64, angle: f64, time: f64, ccw: bool) -> f64 {
35	let (decl, eqt) = sun_position(julian_date + time);
36	let t = 1.0 / 15.0
37		* dmath::arccos(
38			&((-dmath::sin(&angle) - dmath::sin(&decl) * dmath::sin(&latitude))
39				/ (dmath::cos(&decl) * dmath::cos(&latitude))),
40		);
41
42	dmath::fix_hour(12.0 - eqt) + if ccw { -t } else { t }
43}
44
45/// (decl, eqt)
46#[allow(non_snake_case)]
47pub fn sun_position(jd: f64) -> (f64, f64) {
48	let D = jd - 2451545.0;
49
50	let q = dmath::fix_angle(280.46061837 + 0.98564736 * D);
51	let g = dmath::fix_angle(357.528 + 0.98560028 * D);
52	let L = dmath::fix_angle(q + 1.915 * dmath::sin(&g) + 0.020 * dmath::sin(&(2.0 * g)));
53	let e = 23.439 - 0.00000036 * D;
54
55	let decl = dmath::arcsin(&(dmath::sin(&e) * dmath::sin(&L)));
56	let eqt = q / 15.0
57		- dmath::fix_hour(dmath::arctan2(&(dmath::cos(&e) * dmath::sin(&L)), &dmath::cos(&L)) / 15.0);
58
59	return (decl, eqt);
60}
61
62pub fn rise_set_angle(elevation: f64) -> f64 {
63	let earth_radius = 6371008.7714; // in meters
64	let angle = dmath::arccos(&(earth_radius / (earth_radius + elevation)));
65	return 0.833 + angle;
66}