practical_astronomy_rust/
coordinates.rs

1use crate::macros as pa_m;
2use crate::util as pa_u;
3
4/// Convert an Angle (degrees, minutes, and seconds) to Decimal Degrees.
5pub fn angle_to_decimal_degrees(degrees: f64, minutes: f64, seconds: f64) -> f64 {
6    let a = seconds.abs() / 60.0;
7    let b = (minutes.abs() + a) / 60.0;
8    let c = degrees.abs() + b;
9    let d = if degrees < 0.0 || minutes < 0.0 || seconds < 0.0 {
10        -c
11    } else {
12        c
13    };
14
15    return d;
16}
17
18/// Convert Decimal Degrees to an Angle (degrees, minutes, and seconds).
19///
20/// ## Returns
21/// degrees, minutes, seconds
22pub fn decimal_degrees_to_angle(decimal_degrees: f64) -> (f64, f64, f64) {
23    let unsigned_decimal = decimal_degrees.abs();
24    let total_seconds = unsigned_decimal * 3600.0;
25    let seconds_2_dp = pa_u::round_f64(total_seconds % 60.0, 2);
26    let corrected_seconds = if seconds_2_dp == 60.0 {
27        0.0
28    } else {
29        seconds_2_dp
30    };
31    let corrected_remainder = if seconds_2_dp == 60.0 {
32        total_seconds + 60.0
33    } else {
34        total_seconds
35    };
36    let minutes = (corrected_remainder / 60.0).floor() % 60.0;
37    let unsigned_degrees = (corrected_remainder / 3600.0).floor();
38    let signed_degrees = if decimal_degrees < 0.0 {
39        -1.0 * unsigned_degrees
40    } else {
41        unsigned_degrees
42    };
43
44    return (signed_degrees, minutes, corrected_seconds.floor());
45}
46
47/// Convert Right Ascension to Hour Angle.
48pub fn right_ascension_to_hour_angle(
49    ra_hours: f64,
50    ra_minutes: f64,
51    ra_seconds: f64,
52    lct_hours: f64,
53    lct_minutes: f64,
54    lct_seconds: f64,
55    is_daylight_saving: bool,
56    zone_correction: i32,
57    local_day: f64,
58    local_month: u32,
59    local_year: u32,
60    geographical_longitude: f64,
61) -> (f64, f64, f64) {
62    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
63
64    let hour_angle = pa_m::ra_ha(
65        ra_hours,
66        ra_minutes,
67        ra_seconds,
68        lct_hours,
69        lct_minutes,
70        lct_seconds,
71        daylight_saving,
72        zone_correction,
73        local_day,
74        local_month,
75        local_year,
76        geographical_longitude,
77    );
78
79    let hour_angle_hours = pa_m::dh_hour(hour_angle);
80    let hour_angle_minutes = pa_m::dh_min(hour_angle);
81    let hour_angle_seconds = pa_m::dh_sec(hour_angle);
82
83    return (
84        hour_angle_hours as f64,
85        hour_angle_minutes as f64,
86        hour_angle_seconds,
87    );
88}
89
90/// Convert Hour Angle to Right Ascension.
91pub fn hour_angle_to_right_ascension(
92    hour_angle_hours: f64,
93    hour_angle_minutes: f64,
94    hour_angle_seconds: f64,
95    lct_hours: f64,
96    lct_minutes: f64,
97    lct_seconds: f64,
98    is_daylight_saving: bool,
99    zone_correction: i32,
100    local_day: f64,
101    local_month: u32,
102    local_year: u32,
103    geographical_longitude: f64,
104) -> (f64, f64, f64) {
105    let daylight_saving = if is_daylight_saving == true { 1 } else { 0 };
106
107    let right_ascension = pa_m::ha_ra(
108        hour_angle_hours,
109        hour_angle_minutes,
110        hour_angle_seconds,
111        lct_hours,
112        lct_minutes,
113        lct_seconds as f64,
114        daylight_saving,
115        zone_correction,
116        local_day as f64,
117        local_month,
118        local_year,
119        geographical_longitude,
120    );
121
122    let right_ascension_hours = pa_m::dh_hour(right_ascension);
123    let right_ascension_minutes = pa_m::dh_min(right_ascension);
124    let right_ascension_seconds = pa_m::dh_sec(right_ascension);
125
126    return (
127        right_ascension_hours as f64,
128        right_ascension_minutes as f64,
129        right_ascension_seconds,
130    );
131}
132
133/// Convert Equatorial Coordinates to Horizon Coordinates.
134pub fn equatorial_coordinates_to_horizon_coordinates(
135    hour_angle_hours: f64,
136    hour_angle_minutes: f64,
137    hour_angle_seconds: f64,
138    declination_degrees: f64,
139    declination_minutes: f64,
140    declination_seconds: f64,
141    geographical_latitude: f64,
142) -> (f64, f64, f64, f64, f64, f64) {
143    let azimuth_in_decimal_degrees = pa_m::eq_az(
144        hour_angle_hours,
145        hour_angle_minutes,
146        hour_angle_seconds,
147        declination_degrees,
148        declination_minutes,
149        declination_seconds,
150        geographical_latitude,
151    );
152
153    let altitude_in_decimal_degrees = pa_m::eq_alt(
154        hour_angle_hours,
155        hour_angle_minutes,
156        hour_angle_seconds,
157        declination_degrees,
158        declination_minutes,
159        declination_seconds,
160        geographical_latitude,
161    );
162
163    let azimuth_degrees = pa_m::dd_deg(azimuth_in_decimal_degrees);
164    let azimuth_minutes = pa_m::dd_min(azimuth_in_decimal_degrees);
165    let azimuth_seconds = pa_m::dd_sec(azimuth_in_decimal_degrees);
166
167    let altitude_degrees = pa_m::dd_deg(altitude_in_decimal_degrees);
168    let altitude_minutes = pa_m::dd_min(altitude_in_decimal_degrees);
169    let altitude_seconds = pa_m::dd_sec(altitude_in_decimal_degrees);
170
171    return (
172        azimuth_degrees,
173        azimuth_minutes,
174        azimuth_seconds,
175        altitude_degrees,
176        altitude_minutes,
177        altitude_seconds,
178    );
179}
180
181/// Convert Horizon Coordinates to Equatorial Coordinates.
182pub fn horizon_coordinates_to_equatorial_coordinates(
183    azimuth_degrees: f64,
184    azimuth_minutes: f64,
185    azimuth_seconds: f64,
186    altitude_degrees: f64,
187    altitude_minutes: f64,
188    altitude_seconds: f64,
189    geographical_latitude: f64,
190) -> (f64, f64, f64, f64, f64, f64) {
191    let hour_angle_in_decimal_degrees = pa_m::hor_ha(
192        azimuth_degrees,
193        azimuth_minutes,
194        azimuth_seconds,
195        altitude_degrees,
196        altitude_minutes,
197        altitude_seconds,
198        geographical_latitude,
199    );
200
201    let declination_in_decimal_degrees = pa_m::hor_dec(
202        azimuth_degrees,
203        azimuth_minutes,
204        azimuth_seconds,
205        altitude_degrees,
206        altitude_minutes,
207        altitude_seconds,
208        geographical_latitude,
209    );
210
211    let hour_angle_hours = pa_m::dh_hour(hour_angle_in_decimal_degrees);
212    let hour_angle_minutes = pa_m::dh_min(hour_angle_in_decimal_degrees);
213    let hour_angle_seconds = pa_m::dh_sec(hour_angle_in_decimal_degrees);
214
215    let declination_degrees = pa_m::dd_deg(declination_in_decimal_degrees);
216    let declination_minutes = pa_m::dd_min(declination_in_decimal_degrees);
217    let declination_seconds = pa_m::dd_sec(declination_in_decimal_degrees);
218
219    return (
220        hour_angle_hours as f64,
221        hour_angle_minutes as f64,
222        hour_angle_seconds,
223        declination_degrees,
224        declination_minutes,
225        declination_seconds,
226    );
227}
228
229/// Calculate Mean Obliquity of the Ecliptic for a Greenwich Date.
230pub fn mean_obliquity_of_the_ecliptic(
231    greenwich_day: f64,
232    greenwich_month: u32,
233    greenwich_year: u32,
234) -> f64 {
235    let jd = pa_m::cd_jd(greenwich_day, greenwich_month, greenwich_year);
236    let mjd = jd - 2451545.0;
237    let t = mjd / 36525.0;
238    let de1 = t * (46.815 + t * (0.0006 - (t * 0.00181)));
239    let de2 = de1 / 3600.0;
240
241    return 23.439292 - de2;
242}
243
244/// Convert Ecliptic Coordinates to Equatorial Coordinates.
245pub fn ecliptic_coordinate_to_equatorial_coordinate(
246    ecliptic_longitude_degrees: f64,
247    ecliptic_longitude_minutes: f64,
248    ecliptic_longitude_seconds: f64,
249    ecliptic_latitude_degrees: f64,
250    ecliptic_latitude_minutes: f64,
251    ecliptic_latitude_seconds: f64,
252    greenwich_day: f64,
253    greenwich_month: u32,
254    greenwich_year: u32,
255) -> (f64, f64, f64, f64, f64, f64) {
256    let eclon_deg = pa_m::dms_dd(
257        ecliptic_longitude_degrees,
258        ecliptic_longitude_minutes,
259        ecliptic_longitude_seconds,
260    );
261    let eclat_deg = pa_m::dms_dd(
262        ecliptic_latitude_degrees,
263        ecliptic_latitude_minutes,
264        ecliptic_latitude_seconds,
265    );
266    let eclon_rad = eclon_deg.to_radians();
267    let eclat_rad = eclat_deg.to_radians();
268    let obliq_deg = pa_m::obliq(greenwich_day, greenwich_month, greenwich_year);
269    let obliq_rad = obliq_deg.to_radians();
270    let sin_dec =
271        eclat_rad.sin() * obliq_rad.cos() + eclat_rad.cos() * obliq_rad.sin() * eclon_rad.sin();
272    let dec_rad = sin_dec.asin();
273    let dec_deg = pa_m::degrees(dec_rad);
274    let y = eclon_rad.sin() * obliq_rad.cos() - eclat_rad.tan() * obliq_rad.sin();
275    let x = eclon_rad.cos();
276    let ra_rad = y.atan2(x);
277    let ra_deg1 = pa_m::degrees(ra_rad);
278    let ra_deg2 = ra_deg1 - 360.0 * (ra_deg1 / 360.0).floor();
279    let ra_hours = pa_m::dd_dh(ra_deg2);
280
281    let out_ra_hours = pa_m::dh_hour(ra_hours);
282    let out_ra_minutes = pa_m::dh_min(ra_hours);
283    let out_ra_seconds = pa_m::dh_sec(ra_hours);
284    let out_dec_degrees = pa_m::dd_deg(dec_deg);
285    let out_dec_minutes = pa_m::dd_min(dec_deg);
286    let out_dec_seconds = pa_m::dd_sec(dec_deg);
287
288    return (
289        out_ra_hours as f64,
290        out_ra_minutes as f64,
291        out_ra_seconds,
292        out_dec_degrees,
293        out_dec_minutes,
294        out_dec_seconds,
295    );
296}
297
298/// Convert Equatorial Coordinates to Ecliptic Coordinates.
299pub fn equatorial_coordinate_to_ecliptic_coordinate(
300    ra_hours: f64,
301    ra_minutes: f64,
302    ra_seconds: f64,
303    dec_degrees: f64,
304    dec_minutes: f64,
305    dec_seconds: f64,
306    gw_day: f64,
307    gw_month: u32,
308    gw_year: u32,
309) -> (f64, f64, f64, f64, f64, f64) {
310    let ra_deg = pa_m::dh_dd(pa_m::hms_dh(ra_hours, ra_minutes, ra_seconds));
311    let dec_deg = pa_m::dms_dd(dec_degrees, dec_minutes, dec_seconds);
312    let ra_rad = ra_deg.to_radians();
313    let dec_rad = dec_deg.to_radians();
314    let obliq_deg = pa_m::obliq(gw_day, gw_month, gw_year);
315    let obliq_rad = obliq_deg.to_radians();
316    let sin_ecl_lat =
317        dec_rad.sin() * obliq_rad.cos() - dec_rad.cos() * obliq_rad.sin() * ra_rad.sin();
318    let ecl_lat_rad = sin_ecl_lat.asin();
319    let ecl_lat_deg = pa_m::degrees(ecl_lat_rad);
320    let y = ra_rad.sin() * obliq_rad.cos() + dec_rad.tan() * obliq_rad.sin();
321    let x = ra_rad.cos();
322    let ecl_long_rad = y.atan2(x);
323    let ecl_long_deg1 = pa_m::degrees(ecl_long_rad);
324    let ecl_long_deg2 = ecl_long_deg1 - 360.0 * (ecl_long_deg1 / 360.0).floor();
325
326    let out_ecl_long_deg = pa_m::dd_deg(ecl_long_deg2);
327    let out_ecl_long_min = pa_m::dd_min(ecl_long_deg2);
328    let out_ecl_long_sec = pa_m::dd_sec(ecl_long_deg2);
329    let out_ecl_lat_deg = pa_m::dd_deg(ecl_lat_deg);
330    let out_ecl_lat_min = pa_m::dd_min(ecl_lat_deg);
331    let out_ecl_lat_sec = pa_m::dd_sec(ecl_lat_deg);
332
333    return (
334        out_ecl_long_deg,
335        out_ecl_long_min,
336        out_ecl_long_sec,
337        out_ecl_lat_deg,
338        out_ecl_lat_min,
339        out_ecl_lat_sec,
340    );
341}
342
343/// Convert Equatorial Coordinates to Galactic Coordinates.
344pub fn equatorial_coordinate_to_galactic_coordinate(
345    ra_hours: f64,
346    ra_minutes: f64,
347    ra_seconds: f64,
348    dec_degrees: f64,
349    dec_minutes: f64,
350    dec_seconds: f64,
351) -> (f64, f64, f64, f64, f64, f64) {
352    let ra_deg = pa_m::dh_dd(pa_m::hms_dh(ra_hours, ra_minutes, ra_seconds));
353    let dec_deg = pa_m::dms_dd(dec_degrees, dec_minutes, dec_seconds);
354    let ra_rad = ra_deg.to_radians();
355    let dec_rad = dec_deg.to_radians();
356    let sin_b = dec_rad.cos()
357        * (27.4 as f64).to_radians().cos()
358        * (ra_rad - (192.25 as f64).to_radians()).cos()
359        + dec_rad.sin() * (27.4 as f64).to_radians().sin();
360    let b_radians = sin_b.asin();
361    let b_deg = pa_m::degrees(b_radians);
362    let y = dec_rad.sin() - sin_b * (27.4 as f64).to_radians().sin();
363    let x = dec_rad.cos()
364        * (ra_rad - (192.25 as f64).to_radians()).sin()
365        * (27.4 as f64).to_radians().cos();
366    let long_deg1 = pa_m::degrees(y.atan2(x)) + 33.0;
367    let long_deg2 = long_deg1 - 360.0 * (long_deg1 / 360.0).floor();
368
369    let gal_long_deg = pa_m::dd_deg(long_deg2);
370    let gal_long_min = pa_m::dd_min(long_deg2);
371    let gal_long_sec = pa_m::dd_sec(long_deg2);
372    let gal_lat_deg = pa_m::dd_deg(b_deg);
373    let gal_lat_min = pa_m::dd_min(b_deg);
374    let gal_lat_sec = pa_m::dd_sec(b_deg);
375
376    return (
377        gal_long_deg,
378        gal_long_min,
379        gal_long_sec,
380        gal_lat_deg,
381        gal_lat_min,
382        gal_lat_sec,
383    );
384}
385
386/// Convert Galactic Coordinates to Equatorial Coordinates.
387pub fn galactic_coordinate_to_equatorial_coordinate(
388    gal_long_deg: f64,
389    gal_long_min: f64,
390    gal_long_sec: f64,
391    gal_lat_deg: f64,
392    gal_lat_min: f64,
393    gal_lat_sec: f64,
394) -> (f64, f64, f64, f64, f64, f64) {
395    let glong_deg = pa_m::dms_dd(gal_long_deg, gal_long_min, gal_long_sec);
396    let glat_deg = pa_m::dms_dd(gal_lat_deg, gal_lat_min, gal_lat_sec);
397    let glong_rad = glong_deg.to_radians();
398    let glat_rad = glat_deg.to_radians();
399    let sin_dec = glat_rad.cos()
400        * (27.4 as f64).to_radians().cos()
401        * (glong_rad - (33 as f64).to_radians()).sin()
402        + glat_rad.sin() * (27.4 as f64).to_radians().sin();
403    let dec_radians = sin_dec.asin();
404    let dec_deg = pa_m::degrees(dec_radians);
405    let y = glat_rad.cos() * (glong_rad - (33 as f64).to_radians()).cos();
406    let x = glat_rad.sin() * ((27.4 as f64).to_radians()).cos()
407        - (glat_rad).cos()
408            * ((27.4 as f64).to_radians()).sin()
409            * (glong_rad - (33 as f64).to_radians()).sin();
410
411    let ra_deg1 = pa_m::degrees(y.atan2(x)) + 192.25;
412    let ra_deg2 = ra_deg1 - 360.0 * (ra_deg1 / 360.0).floor();
413    let ra_hours1 = pa_m::dd_dh(ra_deg2);
414
415    let ra_hours = pa_m::dh_hour(ra_hours1);
416    let ra_minutes = pa_m::dh_min(ra_hours1);
417    let ra_seconds = pa_m::dh_sec(ra_hours1);
418    let dec_degrees = pa_m::dd_deg(dec_deg);
419    let dec_minutes = pa_m::dd_min(dec_deg);
420    let dec_seconds = pa_m::dd_sec(dec_deg);
421
422    return (
423        ra_hours as f64,
424        ra_minutes as f64,
425        ra_seconds,
426        dec_degrees,
427        dec_minutes,
428        dec_seconds,
429    );
430}
431
432/// Calculate the angle between two celestial objects.
433pub fn angle_between_two_objects(
434    ra_long_1_hour_deg: f64,
435    ra_long_1_min: f64,
436    ra_long_1_sec: f64,
437    dec_lat_1_deg: f64,
438    dec_lat_1_min: f64,
439    dec_lat_1_sec: f64,
440    ra_long_2_hour_deg: f64,
441    ra_long_2_min: f64,
442    ra_long_2_sec: f64,
443    dec_lat_2_deg: f64,
444    dec_lat_2_min: f64,
445    dec_lat_2_sec: f64,
446    hour_or_degree: String,
447) -> (f64, f64, f64) {
448    let ra_long_1_decimal = if hour_or_degree == "H" {
449        pa_m::hms_dh(ra_long_1_hour_deg, ra_long_1_min, ra_long_1_sec)
450    } else {
451        pa_m::dms_dd(ra_long_1_hour_deg, ra_long_1_min, ra_long_1_sec)
452    };
453    let ra_long_1_deg = if hour_or_degree == "H" {
454        pa_m::dh_dd(ra_long_1_decimal)
455    } else {
456        ra_long_1_decimal
457    };
458    let ra_long_1_rad = ra_long_1_deg.to_radians();
459    let dec_lat_1_deg1 = pa_m::dms_dd(dec_lat_1_deg, dec_lat_1_min, dec_lat_1_sec);
460    let dec_lat_1_rad = dec_lat_1_deg1.to_radians();
461    let ra_long_2_decimal = if hour_or_degree == "H" {
462        pa_m::hms_dh(ra_long_2_hour_deg, ra_long_2_min, ra_long_2_sec)
463    } else {
464        pa_m::dms_dd(ra_long_2_hour_deg, ra_long_2_min, ra_long_2_sec)
465    };
466    let ra_long_2_deg = if hour_or_degree == "H" {
467        pa_m::dh_dd(ra_long_2_decimal)
468    } else {
469        ra_long_2_decimal
470    };
471    let ra_long_2_rad = ra_long_2_deg.to_radians();
472    let dec_lat_2_deg1 = pa_m::dms_dd(dec_lat_2_deg, dec_lat_2_min, dec_lat_2_sec);
473    let dec_lat_2_rad = dec_lat_2_deg1.to_radians();
474
475    let cos_d = dec_lat_1_rad.sin() * dec_lat_2_rad.sin()
476        + dec_lat_1_rad.cos() * dec_lat_2_rad.cos() * (ra_long_1_rad - ra_long_2_rad).cos();
477    let d_rad = cos_d.acos();
478    let d_deg = pa_m::degrees(d_rad);
479
480    let angle_deg = pa_m::dd_deg(d_deg);
481    let angle_min = pa_m::dd_min(d_deg);
482    let angle_sec = pa_m::dd_sec(d_deg);
483
484    return (angle_deg, angle_min, angle_sec);
485}
486
487/// Rising and Setting times.
488///
489/// ## Arguments
490/// * `ra_hours` -- Right Ascension, in hours.
491/// * `ra_minutes` -- Right Ascension, in minutes.
492/// * `ra_seconds` -- Right Ascension, in seconds.
493/// * `dec_deg` -- Declination, in degrees.
494/// * `dec_min` -- Declination, in minutes.
495/// * `dec_sec` -- Declination, in seconds.
496/// * `gw_date_day` -- Greenwich Date, day part.
497/// * `gw_date_month` -- Greenwich Date, month part.
498/// * `gw_date_year` -- Greenwich Date, year part.
499/// * `geog_long_deg` -- Geographical Longitude, in degrees.
500/// * `geog_lat_deg` -- Geographical Latitude, in degrees.
501/// * `vert_shift_deg` -- Vertical Shift, in degrees.
502///
503/// ## Returns
504/// * `rise_set_status` -- "Never Rises", "Circumpolar", or "OK".
505/// * `ut_rise_hour` -- Rise time, UT, hour part.
506/// * `ut_rise_min` -- Rise time, UT, minute part.
507/// * `ut_set_hour` -- Set time, UT, hour part.
508/// * `ut_set_min` -- Set time, UT, minute part.
509/// * `az_rise` -- Azimuth angle, at rise.
510/// * `az_set` -- Azimuth angle, at set.
511pub fn rising_and_setting(
512    ra_hours: f64,
513    ra_minutes: f64,
514    ra_seconds: f64,
515    dec_deg: f64,
516    dec_min: f64,
517    dec_sec: f64,
518    gw_date_day: f64,
519    gw_date_month: u32,
520    gw_date_year: u32,
521    geog_long_deg: f64,
522    geog_lat_deg: f64,
523    vert_shift_deg: f64,
524) -> (String, f64, f64, f64, f64, f64, f64) {
525    let ra_hours1 = pa_m::hms_dh(ra_hours, ra_minutes, ra_seconds);
526    let dec_rad = (pa_m::dms_dd(dec_deg, dec_min, dec_sec)).to_radians();
527    let vertical_displ_radians = (vert_shift_deg).to_radians();
528    let geo_lat_radians = (geog_lat_deg).to_radians();
529    let cos_h = -((vertical_displ_radians).sin() + (geo_lat_radians).sin() * (dec_rad).sin())
530        / ((geo_lat_radians).cos() * (dec_rad).cos());
531    let h_hours = pa_m::dd_dh(pa_m::degrees((cos_h).acos()));
532    let lst_rise_hours = (ra_hours1 - h_hours) - 24.0 * ((ra_hours1 - h_hours) / 24.0).floor();
533    let lst_set_hours = (ra_hours1 + h_hours) - 24.0 * ((ra_hours1 + h_hours) / 24.0).floor();
534    let a_deg = pa_m::degrees(
535        (((dec_rad).sin() + (vertical_displ_radians).sin() * (geo_lat_radians).sin())
536            / ((vertical_displ_radians).cos() * (geo_lat_radians).cos()))
537        .acos(),
538    );
539    let az_rise_deg = a_deg - 360.0 * (a_deg / 360.0).floor();
540    let az_set_deg = (360.0 - a_deg) - 360.0 * ((360.0 - a_deg) / 360.0).floor();
541    let ut_rise_hours1 = pa_m::gst_ut(
542        pa_m::lst_gst(lst_rise_hours, 0.0, 0.0, geog_long_deg),
543        0.0,
544        0.0,
545        gw_date_day,
546        gw_date_month,
547        gw_date_year,
548    );
549    let ut_set_hours1 = pa_m::gst_ut(
550        pa_m::lst_gst(lst_set_hours, 0.0, 0.0, geog_long_deg),
551        0.0,
552        0.0,
553        gw_date_day,
554        gw_date_month,
555        gw_date_year,
556    );
557    let ut_rise_adjusted_hours = ut_rise_hours1 + 0.008333;
558    let ut_set_adjusted_hours = ut_set_hours1 + 0.008333;
559
560    let mut rise_set_status = "OK";
561    if cos_h > 1.0 {
562        rise_set_status = "never rises";
563    }
564    if cos_h < -1.0 {
565        rise_set_status = "circumpolar";
566    }
567
568    let ut_rise_hour = if rise_set_status == "OK" {
569        pa_m::dh_hour(ut_rise_adjusted_hours) as f64
570    } else {
571        0.0
572    };
573    let ut_rise_min = if rise_set_status == "OK" {
574        pa_m::dh_min(ut_rise_adjusted_hours) as f64
575    } else {
576        0.0
577    };
578    let ut_set_hour = if rise_set_status == "OK" {
579        pa_m::dh_hour(ut_set_adjusted_hours) as f64
580    } else {
581        0.0
582    };
583    let ut_set_min = if rise_set_status == "OK" {
584        pa_m::dh_min(ut_set_adjusted_hours) as f64
585    } else {
586        0.0
587    };
588    let az_rise = if rise_set_status == "OK" {
589        pa_u::round_f64(az_rise_deg, 2)
590    } else {
591        0.0
592    };
593    let az_set = if rise_set_status == "OK" {
594        pa_u::round_f64(az_set_deg, 2)
595    } else {
596        0.0
597    };
598
599    return (
600        rise_set_status.to_string(),
601        ut_rise_hour,
602        ut_rise_min,
603        ut_set_hour,
604        ut_set_min,
605        az_rise,
606        az_set,
607    );
608}
609
610/// Calculate precession (corrected coordinates between two epochs).
611///
612/// ## Returns
613/// * corrected RA hour
614/// * corrected RA minutes
615/// * corrected RA seconds
616/// * corrected Declination degrees
617/// * corrected Declination minutes
618/// * corrected Declination seconds
619pub fn correct_for_precession(
620    ra_hour: f64,
621    ra_minutes: f64,
622    ra_seconds: f64,
623    dec_deg: f64,
624    dec_minutes: f64,
625    dec_seconds: f64,
626    epoch1_day: f64,
627    epoch1_month: u32,
628    epoch1_year: u32,
629    epoch2_day: f64,
630    epoch2_month: u32,
631    epoch2_year: u32,
632) -> (f64, f64, f64, f64, f64, f64) {
633    let ra_1_rad = (pa_m::dh_dd(pa_m::hms_dh(ra_hour, ra_minutes, ra_seconds))).to_radians();
634    let dec_1_rad = (pa_m::dms_dd(dec_deg, dec_minutes, dec_seconds)).to_radians();
635    let t_centuries = (pa_m::cd_jd(epoch1_day, epoch1_month, epoch1_year) - 2415020.0) / 36525.0;
636    let m_sec = 3.07234 + (0.00186 * t_centuries);
637    let n_arcsec = 20.0468 - (0.0085 * t_centuries);
638    let n_years = (pa_m::cd_jd(epoch2_day, epoch2_month, epoch2_year)
639        - pa_m::cd_jd(epoch1_day, epoch1_month, epoch1_year))
640        / 365.25;
641    let s1_hours =
642        ((m_sec + (n_arcsec * (ra_1_rad).sin() * (dec_1_rad).tan() / 15.0)) * n_years) / 3600.0;
643    let ra_2_hours = pa_m::hms_dh(ra_hour, ra_minutes, ra_seconds) + s1_hours;
644    let s2_deg = (n_arcsec * (ra_1_rad).cos() * n_years) / 3600.0;
645    let dec_2_deg = pa_m::dms_dd(dec_deg, dec_minutes, dec_seconds) + s2_deg;
646
647    let corrected_ra_hour = pa_m::dh_hour(ra_2_hours);
648    let corrected_ra_minutes = pa_m::dh_min(ra_2_hours);
649    let corrected_ra_seconds = pa_m::dh_sec(ra_2_hours);
650    let corrected_dec_deg = pa_m::dd_deg(dec_2_deg);
651    let corrected_dec_minutes = pa_m::dd_min(dec_2_deg);
652    let corrected_dec_seconds = pa_m::dd_sec(dec_2_deg);
653
654    return (
655        corrected_ra_hour as f64,
656        corrected_ra_minutes as f64,
657        corrected_ra_seconds,
658        corrected_dec_deg,
659        corrected_dec_minutes,
660        corrected_dec_seconds,
661    );
662}
663
664/// Calculate nutation for two values: ecliptic longitude and obliquity, for a Greenwich date.
665///
666/// ## Returns
667/// * nutation in ecliptic longitude (degrees)
668/// * nutation in obliquity (degrees)
669pub fn nutation_in_ecliptic_longitude_and_obliquity(
670    greenwich_day: f64,
671    greenwich_month: u32,
672    greenwich_year: u32,
673) -> (f64, f64) {
674    let jd_days = pa_m::cd_jd(greenwich_day, greenwich_month, greenwich_year);
675    let t_centuries = (jd_days - 2415020.0) / 36525.0;
676    let a_deg = 100.0021358 * t_centuries;
677    let l_1_deg = 279.6967 + (0.000303 * t_centuries * t_centuries);
678    let l_deg1 = l_1_deg + 360.0 * (a_deg - (a_deg).floor());
679    let l_deg2 = l_deg1 - 360.0 * (l_deg1 / 360.0).floor();
680    let l_rad = (l_deg2).to_radians();
681    let b_deg = 5.372617 * t_centuries;
682    let n_deg1 = 259.1833 - 360.0 * (b_deg - (b_deg).floor());
683    let n_deg2 = n_deg1 - 360.0 * ((n_deg1 / 360.0).floor());
684    let n_rad = (n_deg2).to_radians();
685    let nut_in_long_arcsec = -17.2 * (n_rad).sin() - 1.3 * (2.0 * l_rad).sin();
686    let nut_in_obl_arcsec = 9.2 * (n_rad).cos() + 0.5 * (2.0 * l_rad).cos();
687
688    let nut_in_long_deg = nut_in_long_arcsec / 3600.0;
689    let nut_in_obl_deg = nut_in_obl_arcsec / 3600.0;
690
691    return (nut_in_long_deg, nut_in_obl_deg);
692}
693
694/// Correct ecliptic coordinates for the effects of aberration.
695///
696/// ## Returns
697/// * apparent ecliptic longitude (degrees, minutes, seconds)
698/// * apparent ecliptic latitude (degrees, minutes, seconds)
699pub fn correct_for_aberration(
700    ut_hour: f64,
701    ut_minutes: f64,
702    ut_seconds: f64,
703    gw_day: f64,
704    gw_month: u32,
705    gw_year: u32,
706    true_ecl_long_deg: f64,
707    true_ecl_long_min: f64,
708    true_ecl_long_sec: f64,
709    true_ecl_lat_deg: f64,
710    true_ecl_lat_min: f64,
711    true_ecl_lat_sec: f64,
712) -> (f64, f64, f64, f64, f64, f64) {
713    let true_long_deg = pa_m::dms_dd(true_ecl_long_deg, true_ecl_long_min, true_ecl_long_sec);
714    let true_lat_deg = pa_m::dms_dd(true_ecl_lat_deg, true_ecl_lat_min, true_ecl_lat_sec);
715    let sun_true_long_deg = pa_m::sun_long(
716        ut_hour, ut_minutes, ut_seconds, 0, 0, gw_day, gw_month, gw_year,
717    );
718    let dlong_arcsec = -20.5 * ((sun_true_long_deg - true_long_deg).to_radians()).cos()
719        / ((true_lat_deg).to_radians()).cos();
720    let dlat_arcsec = -20.5
721        * ((sun_true_long_deg - true_long_deg).to_radians()).sin()
722        * ((true_lat_deg).to_radians()).sin();
723    let apparent_long_deg = true_long_deg + (dlong_arcsec / 3600.0);
724    let apparent_lat_deg = true_lat_deg + (dlat_arcsec / 3600.0);
725
726    let apparent_ecl_long_deg = pa_m::dd_deg(apparent_long_deg);
727    let apparent_ecl_long_min = pa_m::dd_min(apparent_long_deg);
728    let apparent_ecl_long_sec = pa_m::dd_sec(apparent_long_deg);
729    let apparent_ecl_lat_deg = pa_m::dd_deg(apparent_lat_deg);
730    let apparent_ecl_lat_min = pa_m::dd_min(apparent_lat_deg);
731    let apparent_ecl_lat_sec = pa_m::dd_sec(apparent_lat_deg);
732
733    return (
734        apparent_ecl_long_deg,
735        apparent_ecl_long_min,
736        apparent_ecl_long_sec,
737        apparent_ecl_lat_deg,
738        apparent_ecl_lat_min,
739        apparent_ecl_lat_sec,
740    );
741}
742
743/// Calculate corrected RA/Dec, accounting for atmospheric refraction.
744///
745/// NOTE: Valid values for coordinate_type are "TRUE" and "APPARENT".
746///
747/// ## Returns
748/// * corrected RA hours,minutes,seconds
749/// * corrected Declination degrees,minutes,seconds
750pub fn atmospheric_refraction(
751    true_ra_hour: f64,
752    true_ra_min: f64,
753    true_ra_sec: f64,
754    true_dec_deg: f64,
755    true_dec_min: f64,
756    true_dec_sec: f64,
757    coordinate_type: String,
758    geog_long_deg: f64,
759    geog_lat_deg: f64,
760    daylight_saving_hours: i32,
761    timezone_hours: i32,
762    lcd_day: f64,
763    lcd_month: u32,
764    lcd_year: u32,
765    lct_hour: f64,
766    lct_min: f64,
767    lct_sec: f64,
768    atmospheric_pressure_mbar: f64,
769    atmospheric_temperature_celsius: f64,
770) -> (f64, f64, f64, f64, f64, f64) {
771    let ha_hour = pa_m::ra_ha(
772        true_ra_hour,
773        true_ra_min,
774        true_ra_sec,
775        lct_hour,
776        lct_min,
777        lct_sec,
778        daylight_saving_hours,
779        timezone_hours,
780        lcd_day,
781        lcd_month,
782        lcd_year,
783        geog_long_deg,
784    );
785    let azimuth_deg = pa_m::eq_az(
786        ha_hour,
787        0.0,
788        0.0,
789        true_dec_deg,
790        true_dec_min,
791        true_dec_sec,
792        geog_lat_deg,
793    );
794    let altitude_deg = pa_m::eq_alt(
795        ha_hour,
796        0.0,
797        0.0,
798        true_dec_deg,
799        true_dec_min,
800        true_dec_sec,
801        geog_lat_deg,
802    );
803    let corrected_altitude_deg = pa_m::refract(
804        altitude_deg,
805        coordinate_type,
806        atmospheric_pressure_mbar,
807        atmospheric_temperature_celsius,
808    );
809
810    let corrected_ha_hour = pa_m::hor_ha(
811        azimuth_deg,
812        0.0,
813        0.0,
814        corrected_altitude_deg,
815        0.0,
816        0.0,
817        geog_lat_deg,
818    );
819    let corrected_ra_hour1 = pa_m::ha_ra(
820        corrected_ha_hour,
821        0.0,
822        0.0,
823        lct_hour,
824        lct_min,
825        lct_sec,
826        daylight_saving_hours,
827        timezone_hours,
828        lcd_day,
829        lcd_month,
830        lcd_year,
831        geog_long_deg,
832    );
833    let corrected_dec_deg1 = pa_m::hor_dec(
834        azimuth_deg,
835        0.0,
836        0.0,
837        corrected_altitude_deg,
838        0.0,
839        0.0,
840        geog_lat_deg,
841    );
842
843    let corrected_ra_hour = pa_m::dh_hour(corrected_ra_hour1);
844    let corrected_ra_min = pa_m::dh_min(corrected_ra_hour1);
845    let corrected_ra_sec = pa_m::dh_sec(corrected_ra_hour1);
846    let corrected_dec_deg = pa_m::dd_deg(corrected_dec_deg1);
847    let corrected_dec_min = pa_m::dd_min(corrected_dec_deg1);
848    let corrected_dec_sec = pa_m::dd_sec(corrected_dec_deg1);
849
850    return (
851        corrected_ra_hour as f64,
852        corrected_ra_min as f64,
853        corrected_ra_sec,
854        corrected_dec_deg,
855        corrected_dec_min,
856        corrected_dec_sec,
857    );
858}
859
860/// Calculate corrected RA/Dec, accounting for geocentric parallax.
861///
862/// NOTE: Valid values for coordinate_type are "TRUE" and "APPARENT".
863///
864/// ## Returns
865/// * corrected RA hours,minutes,seconds
866/// * corrected Declination degrees,minutes,seconds
867pub fn corrections_for_geocentric_parallax(
868    ra_hour: f64,
869    ra_min: f64,
870    ra_sec: f64,
871    dec_deg: f64,
872    dec_min: f64,
873    dec_sec: f64,
874    coordinate_type: String,
875    equatorial_hor_parallax_deg: f64,
876    geog_long_deg: f64,
877    geog_lat_deg: f64,
878    height_m: f64,
879    daylight_saving: i32,
880    timezone_hours: i32,
881    lcd_day: f64,
882    lcd_month: u32,
883    lcd_year: u32,
884    lct_hour: f64,
885    lct_min: f64,
886    lct_sec: f64,
887) -> (f64, f64, f64, f64, f64, f64) {
888    let ha_hours = pa_m::ra_ha(
889        ra_hour,
890        ra_min,
891        ra_sec,
892        lct_hour,
893        lct_min,
894        lct_sec,
895        daylight_saving,
896        timezone_hours,
897        lcd_day,
898        lcd_month,
899        lcd_year,
900        geog_long_deg,
901    );
902
903    let corrected_ha_hours = pa_m::parallax_ha(
904        ha_hours,
905        0.0,
906        0.0,
907        dec_deg,
908        dec_min,
909        dec_sec,
910        coordinate_type.to_string(),
911        geog_lat_deg,
912        height_m,
913        equatorial_hor_parallax_deg,
914    );
915
916    let corrected_ra_hours = pa_m::ha_ra(
917        corrected_ha_hours,
918        0.0,
919        0.0,
920        lct_hour,
921        lct_min,
922        lct_sec,
923        daylight_saving,
924        timezone_hours,
925        lcd_day,
926        lcd_month,
927        lcd_year,
928        geog_long_deg,
929    );
930
931    let corrected_dec_deg1 = pa_m::parallax_dec(
932        ha_hours,
933        0.0,
934        0.0,
935        dec_deg,
936        dec_min,
937        dec_sec,
938        coordinate_type.to_string(),
939        geog_lat_deg,
940        height_m,
941        equatorial_hor_parallax_deg,
942    );
943
944    let corrected_ra_hour = pa_m::dh_hour(corrected_ra_hours);
945    let corrected_ra_min = pa_m::dh_min(corrected_ra_hours);
946    let corrected_ra_sec = pa_m::dh_sec(corrected_ra_hours);
947    let corrected_dec_deg = pa_m::dd_deg(corrected_dec_deg1);
948    let corrected_dec_min = pa_m::dd_min(corrected_dec_deg1);
949    let corrected_dec_sec = pa_m::dd_sec(corrected_dec_deg1);
950
951    return (
952        corrected_ra_hour as f64,
953        corrected_ra_min as f64,
954        corrected_ra_sec,
955        corrected_dec_deg,
956        corrected_dec_min,
957        corrected_dec_sec,
958    );
959}
960
961/// Calculate heliographic coordinates for a given Greenwich date, with a given heliographic position angle and heliographic displacement in arc minutes.
962///
963/// ## Returns
964/// heliographic longitude and heliographic latitude, in degrees
965pub fn heliographic_coordinates(
966    helio_position_angle_deg: f64,
967    helio_displacement_arcmin: f64,
968    gwdate_day: f64,
969    gwdate_month: u32,
970    gwdate_year: u32,
971) -> (f64, f64) {
972    let julian_date_days = pa_m::cd_jd(gwdate_day, gwdate_month, gwdate_year);
973    let t_centuries = (julian_date_days - 2415020.0) / 36525.0;
974    let long_asc_node_deg = pa_m::dms_dd(74.0, 22.0, 0.0) + (84.0 * t_centuries / 60.0);
975    let sun_long_deg = pa_m::sun_long(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year);
976    let y = ((long_asc_node_deg - sun_long_deg).to_radians()).sin()
977        * ((pa_m::dms_dd(7.0, 15.0, 0.0)).to_radians()).cos();
978    let x = -((long_asc_node_deg - sun_long_deg).to_radians()).cos();
979    let a_deg = pa_m::degrees(y.atan2(x));
980    let m_deg1 = 360.0 - (360.0 * (julian_date_days - 2398220.0) / 25.38);
981    let m_deg2 = m_deg1 - 360.0 * (m_deg1 / 360.0).floor();
982    let l0_deg1 = m_deg2 + a_deg;
983    let _l0_deg2 = l0_deg1 - 360.0 * (l0_deg1 / 360.0).floor();
984    let b0_rad = (((sun_long_deg - long_asc_node_deg).to_radians()).sin()
985        * ((pa_m::dms_dd(7.0, 15.0, 0.0)).to_radians()).sin())
986    .asin();
987    let theta1_rad = (-((sun_long_deg).to_radians()).cos()
988        * ((pa_m::obliq(gwdate_day, gwdate_month, gwdate_year)).to_radians()).tan())
989    .atan();
990    let theta2_rad = (-((long_asc_node_deg - sun_long_deg).to_radians()).cos()
991        * ((pa_m::dms_dd(7.0, 15.0, 0.0)).to_radians()).tan())
992    .atan();
993    let p_deg = pa_m::degrees(theta1_rad + theta2_rad);
994    let rho1_deg = helio_displacement_arcmin / 60.0;
995    let rho_rad = (2.0 * rho1_deg
996        / pa_m::sun_dia(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year))
997    .asin()
998        - (rho1_deg).to_radians();
999    let b_rad = ((b0_rad).sin() * (rho_rad).cos()
1000        + (b0_rad).cos()
1001            * (rho_rad).sin()
1002            * ((p_deg - helio_position_angle_deg).to_radians()).cos())
1003    .asin();
1004    let b_deg = pa_m::degrees(b_rad);
1005    let l_deg1 = pa_m::degrees(
1006        ((rho_rad).sin() * ((p_deg - helio_position_angle_deg).to_radians()).sin() / (b_rad).cos())
1007            .asin(),
1008    ) + l0_deg1;
1009    let l_deg2 = l_deg1 - 360.0 * (l_deg1 / 360.0).floor();
1010
1011    let helio_long_deg = pa_u::round_f64(l_deg2, 2);
1012    let helio_lat_deg = pa_u::round_f64(b_deg, 2);
1013
1014    return (helio_long_deg, helio_lat_deg);
1015}
1016
1017/// Calculate carrington rotation number for a Greenwich date.
1018///
1019/// ## Returns
1020/// carrington rotation number
1021pub fn carrington_rotation_number(gwdate_day: f64, gwdate_month: u32, gwdate_year: u32) -> i32 {
1022    let julian_date_days = pa_m::cd_jd(gwdate_day, gwdate_month, gwdate_year);
1023
1024    let crn = 1690 + pa_u::round_f64((julian_date_days - 2444235.34) / 27.2753, 0) as i32;
1025
1026    return crn;
1027}
1028
1029/// Calculate selenographic (lunar) coordinates (sub-Earth).
1030///
1031/// ## Returns
1032/// * sub-earth longitude
1033/// * sub-earth latitude
1034/// * position angle of pole
1035pub fn selenographic_coordinates_1(
1036    gwdate_day: f64,
1037    gwdate_month: u32,
1038    gwdate_year: u32,
1039) -> (f64, f64, f64) {
1040    let julian_date_days = pa_m::cd_jd(gwdate_day, gwdate_month, gwdate_year);
1041    let t_centuries = (julian_date_days - 2451545.0) / 36525.0;
1042    let long_asc_node_deg = 125.044522 - 1934.136261 * t_centuries;
1043    let f1 = 93.27191 + 483202.0175 * t_centuries;
1044    let f2 = f1 - 360.0 * (f1 / 360.0).floor();
1045    let geocentric_moon_long_deg =
1046        pa_m::moon_long(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year);
1047    let geocentric_moon_lat_rad =
1048        (pa_m::moon_lat(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year)).to_radians();
1049    let inclination_rad = (pa_m::dms_dd(1.0, 32.0, 32.7)).to_radians();
1050    let node_long_rad = (long_asc_node_deg - geocentric_moon_long_deg).to_radians();
1051    let sin_be = -(inclination_rad).cos() * (geocentric_moon_lat_rad).sin()
1052        + (inclination_rad).sin() * (geocentric_moon_lat_rad).cos() * (node_long_rad).sin();
1053    let sub_earth_lat_deg = pa_m::degrees((sin_be).asin());
1054    let a_rad = (-(geocentric_moon_lat_rad).sin() * (inclination_rad).sin()
1055        - (geocentric_moon_lat_rad).cos() * (inclination_rad).cos() * (node_long_rad).sin())
1056    .atan2((geocentric_moon_lat_rad).cos() * (node_long_rad).cos());
1057    let a_deg = pa_m::degrees(a_rad);
1058    let sub_earth_long_deg1 = a_deg - f2;
1059    let sub_earth_long_deg2 = sub_earth_long_deg1 - 360.0 * (sub_earth_long_deg1 / 360.0).floor();
1060    let sub_earth_long_deg3 = if sub_earth_long_deg2 > 180.0 {
1061        sub_earth_long_deg2 - 360.0
1062    } else {
1063        sub_earth_long_deg2
1064    };
1065    let c1_rad = ((node_long_rad).cos() * (inclination_rad).sin()
1066        / ((geocentric_moon_lat_rad).cos() * (inclination_rad).cos()
1067            + (geocentric_moon_lat_rad).sin() * (inclination_rad).sin() * (node_long_rad).sin()))
1068    .atan();
1069    let obliquity_rad = (pa_m::obliq(gwdate_day, gwdate_month, gwdate_year)).to_radians();
1070    let c2_rad = ((obliquity_rad).sin() * ((geocentric_moon_long_deg).to_radians()).cos()
1071        / ((obliquity_rad).sin()
1072            * (geocentric_moon_lat_rad).sin()
1073            * ((geocentric_moon_long_deg).to_radians()).sin()
1074            - (obliquity_rad).cos() * (geocentric_moon_lat_rad).cos()))
1075    .atan();
1076    let c_deg = pa_m::degrees(c1_rad + c2_rad);
1077
1078    let sub_earth_longitude = pa_u::round_f64(sub_earth_long_deg3, 2);
1079    let sub_earth_latitude = pa_u::round_f64(sub_earth_lat_deg, 2);
1080    let position_angle_of_pole = pa_u::round_f64(c_deg, 2);
1081
1082    return (
1083        sub_earth_longitude,
1084        sub_earth_latitude,
1085        position_angle_of_pole,
1086    );
1087}
1088
1089/// Calculate selenographic (lunar) coordinates (sub-Solar).
1090///
1091/// ## Returns
1092/// * sub-solar longitude
1093/// * sub-solar colongitude
1094/// * sub-solar latitude
1095pub fn selenographic_coordinates_2(
1096    gwdate_day: f64,
1097    gwdate_month: u32,
1098    gwdate_year: u32,
1099) -> (f64, f64, f64) {
1100    let julian_date_days = pa_m::cd_jd(gwdate_day, gwdate_month, gwdate_year);
1101    let t_centuries = (julian_date_days - 2451545.0) / 36525.0;
1102    let long_asc_node_deg = 125.044522 - 1934.136261 * t_centuries;
1103    let f1 = 93.27191 + 483202.0175 * t_centuries;
1104    let f2 = f1 - 360.0 * (f1 / 360.0).floor();
1105    let sun_geocentric_long_deg =
1106        pa_m::sun_long(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year);
1107    let moon_equ_hor_parallax_arc_min =
1108        pa_m::moon_hp(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year) * 60.0;
1109    let sun_earth_dist_au =
1110        pa_m::sun_dist(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year);
1111    let geocentric_moon_lat_rad =
1112        (pa_m::moon_lat(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year)).to_radians();
1113    let geocentric_moon_long_deg =
1114        pa_m::moon_long(0.0, 0.0, 0.0, 0, 0, gwdate_day, gwdate_month, gwdate_year);
1115    let adjusted_moon_long_deg = sun_geocentric_long_deg
1116        + 180.0
1117        + (26.4
1118            * (geocentric_moon_lat_rad).cos()
1119            * ((sun_geocentric_long_deg - geocentric_moon_long_deg).to_radians()).sin()
1120            / (moon_equ_hor_parallax_arc_min * sun_earth_dist_au));
1121    let adjusted_moon_lat_rad =
1122        0.14666 * geocentric_moon_lat_rad / (moon_equ_hor_parallax_arc_min * sun_earth_dist_au);
1123    let inclination_rad = (pa_m::dms_dd(1.0, 32.0, 32.7)).to_radians();
1124    let node_long_rad = (long_asc_node_deg - adjusted_moon_long_deg).to_radians();
1125    let sin_bs = -(inclination_rad).cos() * (adjusted_moon_lat_rad).sin()
1126        + (inclination_rad).sin() * (adjusted_moon_lat_rad).cos() * (node_long_rad).sin();
1127    let sub_solar_lat_deg = pa_m::degrees((sin_bs).asin());
1128    let a_rad = (-(adjusted_moon_lat_rad).sin() * (inclination_rad).sin()
1129        - (adjusted_moon_lat_rad).cos() * (inclination_rad).cos() * (node_long_rad).sin())
1130    .atan2((adjusted_moon_lat_rad).cos() * (node_long_rad).cos());
1131    let a_deg = pa_m::degrees(a_rad);
1132    let sub_solar_long_deg1 = a_deg - f2;
1133    let sub_solar_long_deg2 = sub_solar_long_deg1 - 360.0 * (sub_solar_long_deg1 / 360.0).floor();
1134    let sub_solar_long_deg3 = if sub_solar_long_deg2 > 180.0 {
1135        sub_solar_long_deg2 - 360.0
1136    } else {
1137        sub_solar_long_deg2
1138    };
1139    let sub_solar_colong_deg = 90.0 - sub_solar_long_deg3;
1140
1141    let sub_solar_longitude = pa_u::round_f64(sub_solar_long_deg3, 2);
1142    let sub_solar_colongitude = pa_u::round_f64(sub_solar_colong_deg, 2);
1143    let sub_solar_latitude = pa_u::round_f64(sub_solar_lat_deg, 2);
1144
1145    return (
1146        sub_solar_longitude,
1147        sub_solar_colongitude,
1148        sub_solar_latitude,
1149    );
1150}
practical_astronomy_rust/coordinates.rs

practical_astronomy_rust/
coordinates.rs
