gistools/space/util/time.rs
1use crate::util::Date;
2use core::f64::consts::TAU;
3use libm::floor;
4
5/// Time object to track month-day-hour-minute-second
6#[derive(Debug, Default, Clone, Copy, PartialEq)]
7pub struct TimeStamp {
8 /// month
9 pub mon: f64,
10 /// day
11 pub day: f64,
12 /// hour
13 pub hr: f64,
14 /// minute
15 pub min: f64,
16 /// second
17 pub sec: f64,
18}
19
20/// procedure days2mdhms
21///
22/// this procedure converts the day of the year, days, to the equivalent month
23/// day, hour, minute and second.
24///
25/// algorithm set up array for the number of days per month
26/// find leap year - use 1900 because 2000 is a leap year
27/// loop through a temp value while the value is < the days
28/// perform int conversions to the correct day and month
29/// convert remainder into h m s using type conversions
30///
31/// author david vallado 719-573-2600 1 mar 2001
32///
33/// ## Parameters
34/// - `year`: year to date
35/// - `days`: day of year
36///
37/// ## Returns
38/// Decomposed information into year, month, day, hour, minute and second
39pub fn days2mdhms(year: u16, days: f64) -> TimeStamp {
40 let lmonth = [
41 31.,
42 if year.is_multiple_of(4) { 29. } else { 28. },
43 31.,
44 30.,
45 31.,
46 30.,
47 31.,
48 31.,
49 30.,
50 31.,
51 30.,
52 31.,
53 ];
54 let dayofyr = floor(days);
55
56 // ----------------- find month and day of month ----------------
57 let mut i = 1;
58 let mut inttemp = 0.;
59 while dayofyr > (inttemp + lmonth[i - 1]) && i < 12 {
60 inttemp += lmonth[i - 1];
61 i += 1;
62 }
63
64 let mon = i;
65 let day = dayofyr - inttemp;
66
67 // ----------------- find hours minutes and seconds -------------
68 let mut temp = (days - dayofyr) * 24.0;
69 let hr = floor(temp);
70 temp = (temp - hr) * 60.0;
71 let min = floor(temp);
72 let sec = (temp - min) * 60.0;
73
74 TimeStamp { mon: mon as f64, day, hr, min, sec }
75}
76
77/// procedure jday
78///
79/// this procedure finds the julian date given the year, month, day, and time.
80/// the julian date is defined by each elapsed day since noon, jan 1, 4713 bc.
81///
82/// algorithm calculate the answer in one step for efficiency
83///
84/// author david vallado 719-573-2600 1 mar 2001
85///
86/// ## Parameters
87/// - `year`: year
88/// - `mon`: month
89/// - `day`: day
90/// - `hr`: hour
91/// - `min`: minute
92/// - `sec`: seconds
93/// - `msec`: milliseconds
94///
95/// ## Returns
96/// Julian date
97pub fn jday_internal(
98 year: f64,
99 mon: f64,
100 day: f64,
101 hr: f64,
102 min: f64,
103 sec: f64,
104 msec: Option<f64>,
105) -> f64 {
106 let msec = msec.unwrap_or(0.);
107 367.0 * year - floor(7. * (year + floor((mon + 9.) / 12.0)) * 0.25)
108 + floor((275. * mon) / 9.0)
109 + day
110 + 1721013.5
111 + ((msec / 60000. + sec / 60.0 + min) / 60.0 + hr) / 24.0
112}
113
114/// Builds a julian date from the given parameters
115///
116/// ## Parameters
117/// - `date`: date
118///
119/// ## Returns
120/// Julian date
121pub fn jday(date: &Date) -> f64 {
122 jday_internal(
123 date.year as f64,
124 date.month as f64,
125 date.day as f64,
126 date.hour as f64,
127 date.minute as f64,
128 date.second as f64,
129 None,
130 )
131}
132
133/// function gstime
134///
135/// this function finds the greenwich sidereal time.
136///
137/// author david vallado 719-573-2600 1 mar 2001
138///
139/// ## Parameters
140/// - `jdut1`: julian date
141///
142/// ## Returns
143/// Greenwich sidereal time
144fn gstime_internal(jdut1: f64) -> f64 {
145 let tut1 = (jdut1 - 2451545.0) / 36525.0;
146
147 let mut temp = -6.2e-6 * tut1 * tut1 * tut1
148 + 0.093104 * tut1 * tut1
149 + (876600.0 * 3600.0 + 8640184.812866) * tut1
150 + 67310.54841; // # sec
151 temp = ((temp.to_radians()) / 240.0) % TAU; // 360/86400 = 1/240, to deg, to rad
152
153 // ------------------------ check quadrants ---------------------
154 if temp < 0.0 {
155 temp += TAU;
156 }
157
158 temp
159}
160
161/// find greenwich sidereal time
162///
163/// ## Parameters
164/// - `time`: julian date
165///
166/// ## Returns
167/// Greenwich sidereal time
168pub fn gstime(time: &Date) -> f64 {
169 gstime_internal(jday(time))
170}