Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement IAU1980 nutation #23

Merged
merged 9 commits into from
Nov 20, 2023
Merged

Implement IAU1980 nutation #23

merged 9 commits into from
Nov 20, 2023

Conversation

AngusGMorrison
Copy link
Contributor

@AngusGMorrison AngusGMorrison commented Nov 16, 2023

Implements nutation calculation for the IAU 1980 model, and provides the skeleton to calculate nutation for any given model, pending their implementation.

I had to extend and refactor a few parts of the time crate to support conversions to the inputs expected by the nutation algorithm. My expectation is that this will make implementing future models much easier.

I additionally extracted a few newly shared constants, types and functions to common modules.

crates/lox_core/src/bodies/nutation/iau1980.rs Outdated Show resolved Hide resolved
crates/lox_core/src/bodies/nutation/iau1980.rs Outdated Show resolved Hide resolved
crates/lox_core/src/bodies/nutation.rs Outdated Show resolved Hide resolved
@@ -35,7 +34,7 @@ pub fn mean_moon_sun_elongation(t: TDBJulianCenturiesSinceJ2000) -> Radians {
-0.00003169,
],
);
arcsec_to_rad_wrapping(arcsec)
arcsec_to_rad_two_pi(arcsec)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Renamed for consistency with other modular conversions I found.

@@ -162,20 +161,6 @@ impl MeanLongitude for Uranus {
}
}

const ARCSECONDS_IN_CIRCLE: f64 = 360.0 * 60.0 * 60.0;
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Extracted to common module.

// match ERFA outputs rather than any target threshold.
const THRESHOLD: f64 = 1e-11;
// match ERFA outputs rather than any target TOLERANCE.
const TOLERANCE: f64 = 1e-11;
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Renamed for consistency with other tests.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A less arbitrary definition based on machine epsilon would be:

const TOLERANCE = f64::EPSILON.sqrt(); # ≈ 1e-8

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The square root is probably a bit aggressive for this particular case. Some of the test outputs are on the order of e-6, and I'd hope to test for more than 3 decimal places of precision.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From reading the Random ASCII post you sent me, it seems like the tolerance does generally have to be tailored to the use case.

#[derive(Debug, Default, Clone, PartialEq)]
pub struct Nutation {
/// δψ
pub longitude: Radians,
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For all names in this PR, I'll let me know whether the plain English names or the Greek shorthand is clearer to you as an astrodynamicist. (The plain English is clearest for when you don't have that experience, let me tell you! But we shouldn't adopt a naming convention that will annoy the target audience)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My rule of thumb: Greek letters and one-letter variable names inside the code to stay closer to the source material (e.g. formula in the paper) but plain names in the API.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like it 👍

@@ -113,10 +114,6 @@ pub fn cartesian_to_keplerian(grav_param: f64, pos: DVec3, vel: DVec3) -> Kepler
)
}

fn normalize_two_pi(a: f64, center: f64) -> f64 {
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Moved to common module.

@AngusGMorrison AngusGMorrison marked this pull request as ready for review November 19, 2023 20:20
Copy link
Member

@helgee helgee left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM. Left some non-blocking comments.

#[derive(Debug, Default, Clone, PartialEq)]
pub struct Nutation {
/// δψ
pub longitude: Radians,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My rule of thumb: Greek letters and one-letter variable names inside the code to stay closer to the source material (e.g. formula in the paper) but plain names in the API.

Comment on lines +69 to +87
pub fn from_datetime(scale: TimeScale, dt: DateTime) -> Self {
Self::from_date_and_time(scale, dt.date(), dt.time())
}

/// Returns the J2000 epoch in the given timescale.
pub fn j2000(scale: TimeScale) -> Self {
Self::from_raw(scale, RawEpoch::default())
}

/// Returns, as an epoch in the given timescale, midday on the first day of the proleptic Julian
/// calendar.
pub fn jd0(scale: TimeScale) -> Self {
// This represents 4713 BC, since there is no year 0 separating BC and AD.
let first_proleptic_day = Date::new_unchecked(ProlepticJulian, -4712, 1, 1);
let midday = Time::new(12, 0, 0).expect("midday should be a valid time");
Self::from_date_and_time(scale, first_proleptic_day, midday)
}

fn from_raw(scale: TimeScale, raw: RawEpoch) -> Self {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not entirely convinced about these API changes but we should revisit this later.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm content for these to be part of the crate API. I can't speak to their usefulness to the outside world, but they make testing more convenient.

@AngusGMorrison AngusGMorrison merged commit db93294 into main Nov 20, 2023
5 checks passed
@AngusGMorrison AngusGMorrison deleted the am/nutation branch November 20, 2023 15:42
This was referenced Jul 19, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants