practical_astronomy_rust/
planet.rs

1use crate::macros as pa_m;
2use crate::planetdata as pa_pd;
3use crate::util as pa_u;
4
5/// Calculate approximate position of a planet.
6///
7/// ## Arguments
8/// * `lct_hour` -- Local civil time, in hours.
9/// * `lct_min` -- Local civil time, in minutes.
10/// * `lct_sec` -- Local civil time, in seconds.
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/// * `planet_name` -- Name of planet, e.g., "Jupiter".
17///
18/// ## Returns
19/// * `planet_ra_hour` -- Right ascension of planet (hour part)
20/// * `planet_ra_min` -- Right ascension of planet (minutes part)
21/// * `planet_ra_sec` -- Right ascension of planet (seconds part)
22/// * `planet_dec_deg` -- Declination of planet (degrees part)
23/// * `planet_dec_min` -- Declination of planet (minutes part)
24/// * `planet_dec_sec` -- Declination of planet (seconds part)
25pub fn approximate_position_of_planet(
26    lct_hour: f64,
27    lct_min: f64,
28    lct_sec: f64,
29    is_daylight_saving: bool,
30    zone_correction_hours: i32,
31    local_date_day: f64,
32    local_date_month: u32,
33    local_date_year: u32,
34    planet_name: String,
35) -> (f64, f64, f64, f64, f64, f64) {
36    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
37
38    let (planet_info, _planet_info_status) = pa_pd::get_planet_info_vector(planet_name);
39
40    let planet_tp_from_table = planet_info.tp;
41    let planet_long_from_table = planet_info.long;
42    let planet_peri_from_table = planet_info.peri;
43    let planet_ecc_from_table = planet_info.ecc;
44    let planet_axis_from_table = planet_info.axis;
45    let planet_incl_from_table = planet_info.incl;
46    let planet_node_from_table = planet_info.node;
47
48    let gdate_day = pa_m::lct_gday(
49        lct_hour,
50        lct_min,
51        lct_sec,
52        daylight_saving,
53        zone_correction_hours,
54        local_date_day,
55        local_date_month,
56        local_date_year,
57    );
58    let gdate_month = pa_m::lct_gmonth(
59        lct_hour,
60        lct_min,
61        lct_sec,
62        daylight_saving,
63        zone_correction_hours,
64        local_date_day,
65        local_date_month,
66        local_date_year,
67    );
68    let gdate_year = pa_m::lct_gyear(
69        lct_hour,
70        lct_min,
71        lct_sec,
72        daylight_saving,
73        zone_correction_hours,
74        local_date_day,
75        local_date_month,
76        local_date_year,
77    );
78
79    let ut_hours = pa_m::lct_ut(
80        lct_hour,
81        lct_min,
82        lct_sec,
83        daylight_saving,
84        zone_correction_hours,
85        local_date_day,
86        local_date_month,
87        local_date_year,
88    );
89    let d_days = pa_m::cd_jd(gdate_day + (ut_hours / 24.0), gdate_month, gdate_year)
90        - pa_m::cd_jd(0.0, 1, 2010);
91    let np_deg1 = 360.0 * d_days / (365.242191 * planet_tp_from_table);
92    let np_deg2 = np_deg1 - 360.0 * (np_deg1 / 360.0).floor();
93    let mp_deg = np_deg2 + planet_long_from_table - planet_peri_from_table;
94    let lp_deg1 = np_deg2
95        + (360.0 * planet_ecc_from_table * mp_deg.to_radians().sin() / std::f64::consts::PI)
96        + planet_long_from_table;
97    let lp_deg2 = lp_deg1 - 360.0 * (lp_deg1 / 360.0).floor();
98    let planet_true_anomaly_deg = lp_deg2 - planet_peri_from_table;
99    let r_au = planet_axis_from_table * (1.0 - num::pow(planet_ecc_from_table, 2))
100        / (1.0 + planet_ecc_from_table * planet_true_anomaly_deg.to_radians().cos());
101
102    let (earth_info, _earth_info_status) = pa_pd::get_planet_info_vector("Earth".to_string());
103
104    let earth_tp_from_table = earth_info.tp;
105    let earth_long_from_table = earth_info.long;
106    let earth_peri_from_table = earth_info.peri;
107    let earth_ecc_from_table = earth_info.ecc;
108    let earth_axis_from_table = earth_info.axis;
109
110    let ne_deg1 = 360.0 * d_days / (365.242191 * earth_tp_from_table);
111    let ne_deg2 = ne_deg1 - 360.0 * (ne_deg1 / 360.0).floor();
112    let me_deg = ne_deg2 + earth_long_from_table - earth_peri_from_table;
113    let le_deg1 = ne_deg2
114        + earth_long_from_table
115        + 360.0 * earth_ecc_from_table * me_deg.to_radians().sin() / std::f64::consts::PI;
116    let le_deg2 = le_deg1 - 360.0 * (le_deg1 / 360.0).floor();
117    let earth_true_anomaly_deg = le_deg2 - earth_peri_from_table;
118    let r_au2 = earth_axis_from_table * (1.0 - num::pow(earth_ecc_from_table, 2))
119        / (1.0 + earth_ecc_from_table * earth_true_anomaly_deg.to_radians().cos());
120    let lp_node_rad = (lp_deg2 - planet_node_from_table).to_radians();
121    let psi_rad = ((lp_node_rad).sin() * planet_incl_from_table.to_radians().sin()).asin();
122    let y = lp_node_rad.sin() * planet_incl_from_table.to_radians().cos();
123    let x = lp_node_rad.cos();
124    let ld_deg = pa_m::degrees(y.atan2(x)) + planet_node_from_table;
125    let rd_au = r_au * psi_rad.cos();
126    let le_ld_rad = (le_deg2 - ld_deg).to_radians();
127    let atan2_type_1 = (rd_au * le_ld_rad.sin()).atan2(r_au2 - rd_au * le_ld_rad.cos());
128    let atan2_type_2 = (r_au2 * (-le_ld_rad).sin()).atan2(rd_au - r_au2 * le_ld_rad.cos());
129    let a_rad = if rd_au < 1.0 {
130        atan2_type_1
131    } else {
132        atan2_type_2
133    };
134    let lamda_deg1 = if rd_au < 1.0 {
135        180.0 + le_deg2 + pa_m::degrees(a_rad)
136    } else {
137        pa_m::degrees(a_rad) + ld_deg
138    };
139    let lamda_deg2 = lamda_deg1 - 360.0 * (lamda_deg1 / 360.0).floor();
140    let beta_deg = pa_m::degrees(
141        (rd_au * psi_rad.tan() * ((lamda_deg2 - ld_deg).to_radians()).sin()
142            / (r_au2 * (-le_ld_rad).sin()))
143        .atan(),
144    );
145    let ra_hours = pa_m::dd_dh(pa_m::ec_ra(
146        lamda_deg2,
147        0.0,
148        0.0,
149        beta_deg,
150        0.0,
151        0.0,
152        gdate_day,
153        gdate_month,
154        gdate_year,
155    ));
156    let dec_deg = pa_m::ec_dec(
157        lamda_deg2,
158        0.0,
159        0.0,
160        beta_deg,
161        0.0,
162        0.0,
163        gdate_day,
164        gdate_month,
165        gdate_year,
166    );
167
168    let planet_ra_hour = pa_m::dh_hour(ra_hours);
169    let planet_ra_min = pa_m::dh_min(ra_hours);
170    let planet_ra_sec = pa_m::dh_sec(ra_hours);
171    let planet_dec_deg = pa_m::dd_deg(dec_deg);
172    let planet_dec_min = pa_m::dd_min(dec_deg);
173    let planet_dec_sec = pa_m::dd_sec(dec_deg);
174
175    return (
176        planet_ra_hour as f64,
177        planet_ra_min as f64,
178        planet_ra_sec,
179        planet_dec_deg,
180        planet_dec_min,
181        planet_dec_sec,
182    );
183}
184
185/// Calculate precise position of a planet.
186///
187/// ## Arguments
188/// * `lct_hour` -- Local civil time, hour part.
189/// * `lct_min` -- Local civil time, minutes part.
190/// * `lct_sec` -- Local civil time, seconds part.
191/// * `is_daylight_saving` -- Is daylight savings in effect?
192/// * `zone_correction_hours` -- Time zone correction, in hours.
193/// * `local_date_day` -- Local date, day part.
194/// * `local_date_month` -- Local date, month part.
195/// * `local_date_year` -- Local date, year part.
196/// * `planet_name` -- Name of planet, e.g., "Jupiter".
197///
198/// ## Returns
199/// * `planet_ra_hour` -- Right ascension of planet (hour part)
200/// * `planet_ra_min` -- Right ascension of planet (minutes part)
201/// * `planet_ra_sec` -- Right ascension of planet (seconds part)
202/// * `planet_dec_deg` -- Declination of planet (degrees part)
203/// * `planet_dec_min` -- Declination of planet (minutes part)
204/// * `planet_dec_sec` -- Declination of planet (seconds part)
205pub fn precise_position_of_planet(
206    lct_hour: f64,
207    lct_min: f64,
208    lct_sec: f64,
209    is_daylight_saving: bool,
210    zone_correction_hours: i32,
211    local_date_day: f64,
212    local_date_month: u32,
213    local_date_year: u32,
214    planet_name: String,
215) -> (f64, f64, f64, f64, f64, f64) {
216    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
217
218    let _gdate_day = pa_m::lct_gday(
219        lct_hour,
220        lct_min,
221        lct_sec,
222        daylight_saving,
223        zone_correction_hours,
224        local_date_day,
225        local_date_month,
226        local_date_year,
227    );
228    let _gdate_month = pa_m::lct_gmonth(
229        lct_hour,
230        lct_min,
231        lct_sec,
232        daylight_saving,
233        zone_correction_hours,
234        local_date_day,
235        local_date_month,
236        local_date_year,
237    );
238    let _gdate_year = pa_m::lct_gyear(
239        lct_hour,
240        lct_min,
241        lct_sec,
242        daylight_saving,
243        zone_correction_hours,
244        local_date_day,
245        local_date_month,
246        local_date_year,
247    );
248
249    let (
250        planet_ecl_long_deg,
251        planet_ecl_lat_deg,
252        _planet_distance_au,
253        _planet_h_long1,
254        _planet_h_long2,
255        _planet_h_lat,
256        _planet_r_vect,
257    ) = pa_m::planet_coordinates(
258        lct_hour,
259        lct_min,
260        lct_sec,
261        daylight_saving,
262        zone_correction_hours,
263        local_date_day,
264        local_date_month,
265        local_date_year,
266        planet_name,
267    );
268
269    let planet_ra_hours = pa_m::dd_dh(pa_m::ec_ra(
270        planet_ecl_long_deg,
271        0.0,
272        0.0,
273        planet_ecl_lat_deg,
274        0.0,
275        0.0,
276        local_date_day,
277        local_date_month,
278        local_date_year,
279    ));
280    let planet_dec_deg1 = pa_m::ec_dec(
281        planet_ecl_long_deg,
282        0.0,
283        0.0,
284        planet_ecl_lat_deg,
285        0.0,
286        0.0,
287        local_date_day,
288        local_date_month,
289        local_date_year,
290    );
291
292    let planet_ra_hour = pa_m::dh_hour(planet_ra_hours);
293    let planet_ra_min = pa_m::dh_min(planet_ra_hours);
294    let planet_ra_sec = pa_m::dh_sec(planet_ra_hours);
295    let planet_dec_deg = pa_m::dd_deg(planet_dec_deg1);
296    let planet_dec_min = pa_m::dd_min(planet_dec_deg1);
297    let planet_dec_sec = pa_m::dd_sec(planet_dec_deg1);
298
299    return (
300        planet_ra_hour as f64,
301        planet_ra_min as f64,
302        planet_ra_sec,
303        planet_dec_deg,
304        planet_dec_min,
305        planet_dec_sec,
306    );
307}
308
309/// Calculate several visual aspects of a planet.
310///
311/// ## Arguments
312/// * `lct_hour` -- Local civil time, hour part.
313/// * `lct_min` -- Local civil time, minutes part.
314/// * `lct_sec` -- Local civil time, seconds part.
315/// * `is_daylight_saving` -- Is daylight savings in effect?
316/// * `zone_correction_hours` -- Time zone correction, in hours.
317/// * `local_date_day` -- Local date, day part.
318/// * `local_date_month` -- Local date, month part.
319/// * `local_date_year` -- Local date, year part.
320/// * `planet_name` -- Name of planet, e.g., "Jupiter".
321///
322/// ## Returns
323/// * `distance_au` -- Planet's distance from Earth, in AU.
324/// * `ang_dia_arcsec` -- Angular diameter of the planet, in arcseconds.
325/// * `phase` -- Illuminated fraction of the planet.
326/// * `light_time_hour` -- Light travel time from planet to Earth, hour part.
327/// * `light_time_minutes` -- Light travel time from planet to Earth, minutes part.
328/// * `light_time_seconds` -- Light travel time from planet to Earth, seconds part.
329/// * `pos_angle_bright_limb_deg` -- Position-angle of the bright limb.
330/// * `approximate_magnitude` -- Apparent brightness of the planet.
331pub fn visual_aspects_of_a_planet(
332    lct_hour: f64,
333    lct_min: f64,
334    lct_sec: f64,
335    is_daylight_saving: bool,
336    zone_correction_hours: i32,
337    local_date_day: f64,
338    local_date_month: u32,
339    local_date_year: u32,
340    planet_name: String,
341) -> (f64, f64, f64, f64, f64, f64, f64, f64) {
342    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
343
344    let greenwich_date_day = pa_m::lct_gday(
345        lct_hour,
346        lct_min,
347        lct_sec,
348        daylight_saving,
349        zone_correction_hours,
350        local_date_day,
351        local_date_month,
352        local_date_year,
353    );
354    let greenwich_date_month = pa_m::lct_gmonth(
355        lct_hour,
356        lct_min,
357        lct_sec,
358        daylight_saving,
359        zone_correction_hours,
360        local_date_day,
361        local_date_month,
362        local_date_year,
363    );
364    let greenwich_date_year = pa_m::lct_gyear(
365        lct_hour,
366        lct_min,
367        lct_sec,
368        daylight_saving,
369        zone_correction_hours,
370        local_date_day,
371        local_date_month,
372        local_date_year,
373    );
374
375    let (
376        planet_ecl_long_deg,
377        planet_ecl_lat_deg,
378        planet_dist_au,
379        planet_h_long1,
380        _temp3,
381        _temp4,
382        planet_r_vect,
383    ) = pa_m::planet_coordinates(
384        lct_hour,
385        lct_min,
386        lct_sec,
387        daylight_saving,
388        zone_correction_hours,
389        local_date_day,
390        local_date_month,
391        local_date_year,
392        planet_name.to_string(),
393    );
394
395    let planet_ra_rad = (pa_m::ec_ra(
396        planet_ecl_long_deg,
397        0.0,
398        0.0,
399        planet_ecl_lat_deg,
400        0.0,
401        0.0,
402        local_date_day,
403        local_date_month,
404        local_date_year,
405    ))
406    .to_radians();
407    let planet_dec_rad = (pa_m::ec_dec(
408        planet_ecl_long_deg,
409        0.0,
410        0.0,
411        planet_ecl_lat_deg,
412        0.0,
413        0.0,
414        local_date_day,
415        local_date_month,
416        local_date_year,
417    ))
418    .to_radians();
419
420    let light_travel_time_hours = planet_dist_au * 0.1386;
421    let (planet_info, _planet_info_status) = pa_pd::get_planet_info_vector(planet_name.to_string());
422    let angular_diameter_arcsec = planet_info.theta0 / planet_dist_au;
423    let phase1 = 0.5 * (1.0 + ((planet_ecl_long_deg - planet_h_long1).to_radians()).cos());
424
425    let sun_ecl_long_deg = pa_m::sun_long(
426        lct_hour,
427        lct_min,
428        lct_sec,
429        daylight_saving,
430        zone_correction_hours,
431        local_date_day,
432        local_date_month,
433        local_date_year,
434    );
435    let sun_ra_rad = (pa_m::ec_ra(
436        sun_ecl_long_deg,
437        0.0,
438        0.0,
439        0.0,
440        0.0,
441        0.0,
442        greenwich_date_day,
443        greenwich_date_month,
444        greenwich_date_year,
445    ))
446    .to_radians();
447    let sun_dec_rad = (pa_m::ec_dec(
448        sun_ecl_long_deg,
449        0.0,
450        0.0,
451        0.0,
452        0.0,
453        0.0,
454        greenwich_date_day,
455        greenwich_date_month,
456        greenwich_date_year,
457    ))
458    .to_radians();
459
460    let y = (sun_dec_rad).cos() * (sun_ra_rad - planet_ra_rad).sin();
461    let x = (planet_dec_rad).cos() * (sun_dec_rad).sin()
462        - (planet_dec_rad).sin() * (sun_dec_rad).cos() * (sun_ra_rad - planet_ra_rad).cos();
463
464    let chi_deg = pa_m::degrees(y.atan2(x));
465    let radius_vector_au = planet_r_vect;
466    let approximate_magnitude1 =
467        5.0 * (radius_vector_au * planet_dist_au / (phase1).sqrt()).log10() + planet_info.v0;
468
469    let distance_au = pa_u::round_f64(planet_dist_au, 5);
470    let ang_dia_arcsec = pa_u::round_f64(angular_diameter_arcsec, 1);
471    let phase = pa_u::round_f64(phase1, 2);
472    let light_time_hour = pa_m::dh_hour(light_travel_time_hours);
473    let light_time_minutes = pa_m::dh_min(light_travel_time_hours);
474    let light_time_seconds = pa_m::dh_sec(light_travel_time_hours);
475    let pos_angle_bright_limb_deg = pa_u::round_f64(chi_deg, 1);
476    let approximate_magnitude = pa_u::round_f64(approximate_magnitude1, 1);
477
478    return (
479        distance_au,
480        ang_dia_arcsec,
481        phase,
482        light_time_hour as f64,
483        light_time_minutes as f64,
484        light_time_seconds,
485        pos_angle_bright_limb_deg,
486        approximate_magnitude,
487    );
488}
practical_astronomy_rust/planet.rs

practical_astronomy_rust/
planet.rs
