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

Elliptical orbit can't compute true anomaly consistently #340

Closed
translunar opened this issue Oct 15, 2024 · 2 comments · Fixed by #342
Closed

Elliptical orbit can't compute true anomaly consistently #340

translunar opened this issue Oct 15, 2024 · 2 comments · Fixed by #342
Labels
good first issue Good for newcomers Kind: improvement New feature or request

Comments

@translunar
Copy link

translunar commented Oct 15, 2024

Bug report

Describe the bug

Can't compute true anomaly for an elliptical orbit consistently (or, consequently, get Cartesian coordinates at various times)

To Reproduce

def make_elliptical_orbit(
        direction: str = 'north',
        epoch: anise.time.Epoch = anise.time.Epoch.system_now(),
        true_anomaly: float = 0.0) -> anise.astro.Orbit:
    """Generate a northern or southern elliptical orbit (12 hour period).

    From Table 9.
    """
    almanac = spice_transforms.almanac
    mj2k = Frames.MOON_J2000
    loaded_mj2k = almanac.frame_info(mj2k)

    if direction == 'north':
        asc_node_d = 270.0
        aop_d = 270.0
    elif direction == 'south':
        asc_node_d = 0.0
        aop_d = 90.0
    else:
        raise ValueError(f"expected north or south, got {direction}")
    return anise.astro.Orbit.from_keplerian(
        sma=6142.400,  # sma
        ecc=0.6,  # ecc
        inc=57.7,  # inc
        aop=aop_d,  # argument of periselene
        raan=asc_node_d,  # ascending node
        ta=true_anomaly,
        epoch=epoch,
        frame=loaded_mj2k)

orbit = make_elliptical_orbit('north')
rs = []

for t in range(0, 12 * 60 * 60, dt_s):
    # Get PV throughout the orbit
    delta = Unit.Second * t
    print(t)
    x = orbit.at_epoch(orbit.epoch + delta)
    rs.append(np.array([x.x_km, x.y_km, x.z_km]))

I get the following output:

(Pdb) c
0
60
120
180
240
300
360
420
480
540
600
660
720
780
840
900
960
1020
1080
1140
1200
1260
1320
1380
1440
1500
1560
Traceback (most recent call last):
    x = orbit.at_epoch(orbit.epoch + delta)
Exception: max iterations reached (1001) when computing the true anomaly from the mean anomaly

Expected behavior

It should be really easy to output the Cartesian coordinates of this elliptical orbit for various epochs.

Code to reproduce the issue

See above.

Platform

Linux, Anise 0.4.4

Additional context

Add any other context about the problem here.

@ChristopherRabotin ChristopherRabotin added Kind: improvement New feature or request good first issue Good for newcomers labels Oct 16, 2024
@ChristopherRabotin
Copy link
Member

Thanks for the report. The function that computes the true anomaly from the mean anomaly loops until the difference between two consecutive true anomaly computations is less than 1e-16 radians:

. This check is too strict. An MA_EPSILON of 1e-10 is probably sufficient.

@ChristopherRabotin
Copy link
Member

I was able to reproduce your bug, and fix it by changing that precision to 1e-12 degrees (from 1e-16). I was just about to release 0.4.5, so the fix will be in that version tomorrow morning.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
good first issue Good for newcomers Kind: improvement New feature or request
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants