practical_astronomy_rust/
comet.rs

1use crate::cometdata as pa_c;
2use crate::macros as pa_m;
3use crate::util as pa_u;
4
5/// Calculate position of an elliptical comet.
6///
7/// ## Arguments
8/// * `lct_hour` -- Local civil time, hour part.
9/// * `lct_min` -- Local civil time, minutes part.
10/// * `lct_sec` -- Local civil time, seconds part.
11/// * `is_daylight_saving` -- Is daylight savings in effect?
12/// * `zone_correction_hours` -- Time zone correction, in hours.
13/// * `local_date_day` -- Local date, day part.
14/// * `local_date_month` -- Local date, month part.
15/// * `local_date_year` -- Local date, year part.
16/// * `comet_name` -- Name of comet, e.g., "Halley".
17///
18/// ## Returns
19/// * `comet_ra_hour` -- Right ascension of comet (hour part)
20/// * `comet_ra_min` -- Right ascension of comet (minutes part)
21/// * `comet_dec_deg` -- Declination of comet (degrees part)
22/// * `comet_dec_min` -- Declination of comet (minutes part)
23/// * `comet_dist_earth` -- Comet's distance from Earth (AU)
24pub fn position_of_elliptical_comet(
25    lct_hour: f64,
26    lct_min: f64,
27    lct_sec: f64,
28    is_daylight_saving: bool,
29    zone_correction_hours: i32,
30    local_date_day: f64,
31    local_date_month: u32,
32    local_date_year: u32,
33    comet_name: String,
34) -> (f64, f64, f64, f64, f64) {
35    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
36
37    let greenwich_date_day = pa_m::lct_gday(
38        lct_hour,
39        lct_min,
40        lct_sec,
41        daylight_saving,
42        zone_correction_hours,
43        local_date_day,
44        local_date_month,
45        local_date_year,
46    );
47    let greenwich_date_month = pa_m::lct_gmonth(
48        lct_hour,
49        lct_min,
50        lct_sec,
51        daylight_saving,
52        zone_correction_hours,
53        local_date_day,
54        local_date_month,
55        local_date_year,
56    );
57    let greenwich_date_year = pa_m::lct_gyear(
58        lct_hour,
59        lct_min,
60        lct_sec,
61        daylight_saving,
62        zone_correction_hours,
63        local_date_day,
64        local_date_month,
65        local_date_year,
66    );
67
68    let (comet_info, _comet_info_status) = pa_c::get_comet_info_elliptical_vector(comet_name);
69
70    let time_since_epoch_years = (pa_m::cd_jd(
71        greenwich_date_day,
72        greenwich_date_month,
73        greenwich_date_year,
74    ) - pa_m::cd_jd(0.0, 1, greenwich_date_year))
75        / 365.242191
76        + greenwich_date_year as f64
77        - comet_info.epoch;
78    let mc_deg = 360.0 * time_since_epoch_years / comet_info.period;
79    let mc_rad = (mc_deg - 360.0 * (mc_deg / 360.0).floor()).to_radians();
80    let eccentricity = comet_info.ecc;
81    let true_anomaly_deg = pa_m::degrees(pa_m::true_anomaly(mc_rad, eccentricity));
82    let lc_deg = true_anomaly_deg + comet_info.peri;
83    let r_au = comet_info.axis * (1.0 - eccentricity * eccentricity)
84        / (1.0 + eccentricity * ((true_anomaly_deg).to_radians()).cos());
85    let lc_node_rad = (lc_deg - comet_info.node).to_radians();
86    let psi_rad = ((lc_node_rad).sin() * ((comet_info.incl).to_radians()).sin()).asin();
87
88    let y = (lc_node_rad).sin() * ((comet_info.incl).to_radians()).cos();
89    let x = (lc_node_rad).cos();
90
91    let ld_deg = pa_m::degrees(y.atan2(x)) + comet_info.node;
92    let rd_au = r_au * (psi_rad).cos();
93
94    let earth_longitude_le_deg = pa_m::sun_long(
95        lct_hour,
96        lct_min,
97        lct_sec,
98        daylight_saving,
99        zone_correction_hours,
100        local_date_day,
101        local_date_month,
102        local_date_year,
103    ) + 180.0;
104    let earth_radius_vector_au = pa_m::sun_dist(
105        lct_hour,
106        lct_min,
107        lct_sec,
108        daylight_saving,
109        zone_correction_hours,
110        local_date_day,
111        local_date_month,
112        local_date_year,
113    );
114
115    let le_ld_rad = (earth_longitude_le_deg - ld_deg).to_radians();
116    let a_rad = if rd_au < earth_radius_vector_au {
117        (rd_au * (le_ld_rad).sin()).atan2(earth_radius_vector_au - rd_au * (le_ld_rad).cos())
118    } else {
119        (earth_radius_vector_au * (-le_ld_rad).sin())
120            .atan2(rd_au - earth_radius_vector_au * (le_ld_rad).cos())
121    };
122
123    let comet_long_deg1 = if rd_au < earth_radius_vector_au {
124        180.0 + earth_longitude_le_deg + pa_m::degrees(a_rad)
125    } else {
126        pa_m::degrees(a_rad) + ld_deg
127    };
128    let comet_long_deg = comet_long_deg1 - 360.0 * (comet_long_deg1 / 360.0).floor();
129    let comet_lat_deg = pa_m::degrees(
130        (rd_au * (psi_rad).tan() * ((comet_long_deg1 - ld_deg).to_radians()).sin()
131            / (earth_radius_vector_au * (-le_ld_rad).sin()))
132        .atan(),
133    );
134    let comet_ra_hours1 = pa_m::dd_dh(pa_m::ec_ra(
135        comet_long_deg,
136        0.0,
137        0.0,
138        comet_lat_deg,
139        0.0,
140        0.0,
141        greenwich_date_day,
142        greenwich_date_month,
143        greenwich_date_year,
144    ));
145    let comet_dec_deg1 = pa_m::ec_dec(
146        comet_long_deg,
147        0.0,
148        0.0,
149        comet_lat_deg,
150        0.0,
151        0.0,
152        greenwich_date_day,
153        greenwich_date_month,
154        greenwich_date_year,
155    );
156    let comet_distance_au = (num::pow(earth_radius_vector_au, 2) + num::pow(r_au, 2)
157        - 2.0
158            * earth_radius_vector_au
159            * r_au
160            * ((lc_deg - earth_longitude_le_deg).to_radians()).cos()
161            * (psi_rad).cos())
162    .sqrt();
163
164    let comet_ra_hour = pa_m::dh_hour(comet_ra_hours1 + 0.008333);
165    let comet_ra_min = pa_m::dh_min(comet_ra_hours1 + 0.008333);
166    let comet_dec_deg = pa_m::dd_deg(comet_dec_deg1 + 0.008333);
167    let comet_dec_min = pa_m::dd_min(comet_dec_deg1 + 0.008333);
168    let comet_dist_earth = pa_u::round_f64(comet_distance_au, 2);
169
170    return (
171        comet_ra_hour as f64,
172        comet_ra_min as f64,
173        comet_dec_deg,
174        comet_dec_min,
175        comet_dist_earth,
176    );
177}
178
179/// Calculate position of a parabolic comet.
180///
181/// ## Arguments
182/// * `lct_hour` -- Local civil time, hour part.
183/// * `lct_min` -- Local civil time, minutes part.
184/// * `lct_sec` -- Local civil time, seconds part.
185/// * `is_daylight_saving` -- Is daylight savings in effect?
186/// * `zone_correction_hours` -- Time zone correction, in hours.
187/// * `local_date_day` -- Local date, day part.
188/// * `local_date_month` -- Local date, month part.
189/// * `local_date_year` -- Local date, year part.
190/// * `comet_name` -- Name of comet, e.g., "Kohler".
191///
192/// ## Returns
193/// * `comet_ra_hour` -- Right ascension of comet (hour part)
194/// * `comet_ra_min` -- Right ascension of comet (minutes part)
195/// * `comet_ra_sec` -- Right ascension of comet (seconds part)
196/// * `comet_dec_deg` -- Declination of comet (degrees part)
197/// * `comet_dec_min` -- Declination of comet (minutes part)
198/// * `comet_dec_sec` -- Declination of comet (seconds part)
199/// * `comet_dist_earth` -- Comet's distance from Earth (AU)
200pub fn position_of_parabolic_comet(
201    lct_hour: f64,
202    lct_min: f64,
203    lct_sec: f64,
204    is_daylight_saving: bool,
205    zone_correction_hours: i32,
206    local_date_day: f64,
207    local_date_month: u32,
208    local_date_year: u32,
209    comet_name: String,
210) -> (f64, f64, f64, f64, f64, f64, f64) {
211    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
212
213    let greenwich_date_day = pa_m::lct_gday(
214        lct_hour,
215        lct_min,
216        lct_sec,
217        daylight_saving,
218        zone_correction_hours,
219        local_date_day,
220        local_date_month,
221        local_date_year,
222    );
223    let greenwich_date_month = pa_m::lct_gmonth(
224        lct_hour,
225        lct_min,
226        lct_sec,
227        daylight_saving,
228        zone_correction_hours,
229        local_date_day,
230        local_date_month,
231        local_date_year,
232    );
233    let greenwich_date_year = pa_m::lct_gyear(
234        lct_hour,
235        lct_min,
236        lct_sec,
237        daylight_saving,
238        zone_correction_hours,
239        local_date_day,
240        local_date_month,
241        local_date_year,
242    );
243
244    let _ut_hours = pa_m::lct_ut(
245        lct_hour,
246        lct_min,
247        lct_sec,
248        daylight_saving,
249        zone_correction_hours,
250        local_date_day,
251        local_date_month,
252        local_date_year,
253    );
254
255    let (comet_info, _comet_info_status) = pa_c::get_comet_info_parabolic_vector(comet_name);
256
257    let perihelion_epoch_day = comet_info.epoch_peri_day;
258    let perihelion_epoch_month = comet_info.epoch_peri_month;
259    let perihelion_epoch_year = comet_info.epoch_peri_year;
260    let q_au = comet_info.peri_dist;
261    let inclination_deg = comet_info.incl;
262    let perihelion_deg = comet_info.arg_peri;
263    let node_deg = comet_info.node;
264
265    let (comet_long_deg, comet_lat_deg, comet_dist_au) = pa_m::p_comet_long_lat_dist(
266        lct_hour,
267        lct_min,
268        lct_sec,
269        daylight_saving,
270        zone_correction_hours,
271        local_date_day,
272        local_date_month,
273        local_date_year,
274        perihelion_epoch_day,
275        perihelion_epoch_month,
276        perihelion_epoch_year,
277        q_au,
278        inclination_deg,
279        perihelion_deg,
280        node_deg,
281    );
282
283    let comet_ra_hours = pa_m::dd_dh(pa_m::ec_ra(
284        comet_long_deg,
285        0.0,
286        0.0,
287        comet_lat_deg,
288        0.0,
289        0.0,
290        greenwich_date_day,
291        greenwich_date_month,
292        greenwich_date_year,
293    ));
294    let comet_dec_deg1 = pa_m::ec_dec(
295        comet_long_deg,
296        0.0,
297        0.0,
298        comet_lat_deg,
299        0.0,
300        0.0,
301        greenwich_date_day,
302        greenwich_date_month,
303        greenwich_date_year,
304    );
305
306    let comet_ra_hour = pa_m::dh_hour(comet_ra_hours);
307    let comet_ra_min = pa_m::dh_min(comet_ra_hours);
308    let comet_ra_sec = pa_m::dh_sec(comet_ra_hours);
309    let comet_dec_deg = pa_m::dd_deg(comet_dec_deg1);
310    let comet_dec_min = pa_m::dd_min(comet_dec_deg1);
311    let comet_dec_sec = pa_m::dd_sec(comet_dec_deg1);
312    let comet_dist_earth = pa_u::round_f64(comet_dist_au, 2);
313
314    return (
315        comet_ra_hour as f64,
316        comet_ra_min as f64,
317        comet_ra_sec,
318        comet_dec_deg,
319        comet_dec_min,
320        comet_dec_sec,
321        comet_dist_earth,
322    );
323}
practical_astronomy_rust/comet.rs

practical_astronomy_rust/
comet.rs
