Skip to content

Fix proper motions not propagated correctly#14

Merged
agabrown merged 1 commit intoagabrown:masterfrom
henrysky:pm_propagation
Dec 14, 2024
Merged

Fix proper motions not propagated correctly#14
agabrown merged 1 commit intoagabrown:masterfrom
henrysky:pm_propagation

Conversation

@henrysky
Copy link
Copy Markdown
Contributor

@henrysky henrysky commented Dec 14, 2024

Thanks for this amazing package to show how astrometric data are propagated.

But I've found that the speed of a stars is not conserved when propagating to a different epoch. For example Acturus went from ~121.7 km/s at J2000 to ~95.3 km/s at J100000. Since the propagation assumes straight-line motion at constant speed, something is wrong. I've found a typo in the code leading to speed not conserved.

Here is a code snippet that I've used to sanity check

from pygaia.astrometry.coordinates import EpochPropagation
from astroquery.simbad import Simbad
import numpy as np

MAS2RAD = 4.84813681109536e-9
PC_TO_KM = 3.086e13  # parsec to km
AU = 149597870.691  # astronomical unit in km
JYEAR_SECONDS = 31557600.0

customSimbad = Simbad()
customSimbad.add_votable_fields(
    "ra(d2000)", "dec(d2000)", "pmra", "pmdec", "otype", "plx", "plx_error", "rv_value"
)
result = customSimbad.query_object("Arcturus")

ra = result["RA_d2000"][0]
dec = result["DEC_d2000"][0]
pmra = result["PMRA"][0]
pmdec = result["PMDEC"][0]
plx = result["PLX_VALUE"][0]
rv = result["RV_VALUE"][0]

totvel = (
    np.sqrt(
        (
            np.sqrt(pmra**2 + pmdec**2)
            * (MAS2RAD / JYEAR_SECONDS)
            * PC_TO_KM
            * (1000 / plx)
        )
        ** 2
        + rv**2
    ),
)

ra_new, dec_new, parallax_new, pmra_new, pmdec_new, pmr_new = (
    EpochPropagation().propagate_astrometry(ra, dec, plx, pmra, pmdec, rv, 2000, 100000)
)

totvel_new = (
    np.sqrt(
        (
            np.sqrt(pmra_new**2 + pmdec_new**2)
            * (MAS2RAD / JYEAR_SECONDS)
            * PC_TO_KM
            * (1000 / parallax_new)
        )
        ** 2
        + ((pmr_new / parallax_new) * (AU / JYEAR_SECONDS)) ** 2
    ),
)

# these two values should be similar
print(totvel, totvel_new)

which will print (np.float64(121.7675972466994),) (np.float64(95.28250254167489),) and with this PR will print (np.float64(121.7675972466994),) (np.float64(121.7609982767927),)

Copy link
Copy Markdown
Owner

@agabrown agabrown left a comment

Choose a reason for hiding this comment

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

Thanks for spotting! I never included this obvious test case.

@agabrown agabrown merged commit 4a68e00 into agabrown:master Dec 14, 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