thomasorb / orcs

ORCS is an analysis engine for SITELLE spectral cubes.
GNU General Public License v3.0
9 stars 6 forks source link

RVCorrect.py #18

Closed Klaus130161 closed 6 years ago

Klaus130161 commented 6 years ago

Dear Thomas,

I found your Python code RVCorrect.py on the Internet and have a question about the accuracy of the result and another one about the sign of the geographical length of the observation site.

First of all, I searched the IRAF. net forum for an example that was calculated with IRAF: http://iraf.net/forum/viewtopic.php?showtopic=138738

The example was:

The object: (year, month, day, UT, RA, DEC) 2007 4 30 08:41:59 15:46:43 37:49:45:45

latitude = 19:49.7 longitude = 155:28.7 !!!! sign wrong ??? altitude = 4160.

The result of IRAF (Frank Valdes): HJD VOBS VHELIO VLSR VDIURNAL VLUNAR VANNUAL VS. SOLAR 2454220.86577 0.00 0.54 17.96 0.239 -0.002 0.299 17.428

Then I calculated the same example with your code and I received the following result:

VOBS: 0.000000 km. s-1 HJD: 2454220.865771 VHELIO: 0.566202 km. s-1 VLSR: 18.003718 km. s-1 VDIURNAL: 0.239563 km. s-1 VLUNAR: -0.002047 km. s-1 VANNUAL: 0.328686 km. s-1 VSOLAR: 17,437515 km. s-1

The difference between the two calculations is 0.566-0.54 = 0,026km/sec.

Do you know the reason for this discrepancy?

I have also done a comparison calculation with barycorr http://astroutils.astronomy.ohio-state.edu/exofast/barycorr.html Here is the result: 567.467035794 m/sec. The difference in this case is only 1.265 m/sec !!

The problem, however, is that the sign of longitude has to be accepted positively for this result. Is that correct?

Thank you in advance for your effort Klaus

thomasorb commented 6 years ago

Dear @Klaus130161,

I have counter checked the results returned by my module with the results returned by the module astropy (http://docs.astropy.org/en/stable/coordinates/velocities.html) with the following code

import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation, FK5
keck = EarthLocation.from_geodetic(lat=19.828333*u.deg, lon=-155.47833333*u.deg, height=4160*u.m) # those are exactly the coordinates of the observatory returned by IRAF (observ.lat, observ.lon)
sc = SkyCoord('15:46:43 +37:49:45', FK5, unit=(u.hourangle, u.deg))
obstime = Time(2454220.865771, format='jd') # exactly the number of JD returned by IRAF (and my module)
barycorr = sc.radial_velocity_correction(obstime=obstime, location=keck)
print barycorr.to(u.km/u.s)
heliocorr = sc.radial_velocity_correction('heliocentric', obstime=obstime, location=keck)
print heliocorr.to(u.km/u.s)

which returns

0.565774937224 km / s # barycentric  (barrycorr.html: 0.567)
0.563097179679 km / s # heliocentric (rvcorrect.py: 0.566 / IRAF: 0.54)

It means that, in terms of barycentric velocity, astropy and the barycorr function (http://astroutils.astronomy.ohio-state.edu/exofast/barycorr.html) return comparable values within the 3 m/s precision error reported in both implementations. And, in terms of heliocentric velocity, both astropy and my implementation are giving similar results. I thus believe that there might be a precision problem in IRAF's calculation or some other subtlety I am not aware of.

For the sign problem, in IRAF's function, the longitude is positive towards west which is the opposite for other conventions.

You pointed out a good problem which is the use of this module rvcorrect.py in ORCS which might have not been completely checked and compared for a large range of input. In fact I will remove it and consider using the astropy implementation which is certainly far more robust as it is checked by numerous people ;)

Thank you !

@Klaus130161 Please tell me if it answers your questions so that I can close this issue

thomasorb commented 6 years ago

I have changed the heliocentric radial velocity correction method. It has been replaced with the much more reliable astropy.coordinates.SkyCoord.radial_velocity_correction method (https://github.com/thomasorb/orcs/commit/f89c9b99892e20148d31480bc803dd316b3c1b0f). rvcorrect.py module has been deleted (https://github.com/thomasorb/orcs/commit/993a0be388de10113b03b20161be0c3e4c01e708).

If you want to work with it feel free to make a branch from https://github.com/thomasorb/orcs/commit/f89c9b99892e20148d31480bc803dd316b3c1b0f (it is under GPL3 licence)

Klaus130161 commented 6 years ago

Dear Thomas,

thank you for your answer and your effort. You can close the issue.

Best regards Klaus

—————————————————————————— Dr. rer. nat. Klaus Vollmann

Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR)

Raumfahrtmanagement / Space Agency Technik für Raumfahrtsysteme und Robotik / General Technologies and Robotics Königswinterer Str. 522-524 53227 Bonn

Telefon +49 (0)228 447-212 Telefax +49 (0)228 447-718 Internet www.dlr.dehttp://www.dlr.de/

Von: Thomas Martin [mailto:notifications@github.com] Gesendet: Dienstag, 9. Januar 2018 23:31 An: thomasorb/orcs Cc: Vollmann, Klaus; Mention Betreff: Re: [thomasorb/orcs] RVCorrect.py (#18)

Dear @Klaus130161https://github.com/klaus130161,

I have counter checked the results returned by my module with the results returned by the module astropy (http://docs.astropy.org/en/stable/coordinates/velocities.html) with the following code

import astropy.units as u

from astropy.time import Time

from astropy.coordinates import SkyCoord, EarthLocation, FK5

keck = EarthLocation.from_geodetic(lat=19.828333u.deg, lon=-155.47833333u.deg, height=4160*u.m) # those are exactly the coordinates of the observatory returned by IRAF (observ.lat, observ.lon)

sc = SkyCoord('15:46:43 +37:49:45', FK5, unit=(u.hourangle, u.deg))

obstime = Time(2454220.865771, format='jd') # exactly the number of JD returned by IRAF (and my module)

barycorr = sc.radial_velocity_correction(obstime=obstime, location=keck)

print barycorr.to(u.km/u.s)

heliocorr = sc.radial_velocity_correction('heliocentric', obstime=obstime, location=keck)

print heliocorr.to(u.km/u.s)

which returns

0.565774937224 km / s # barycentric (barrycorr.html: 0.567)

0.563097179679 km / s # heliocentric (rvcorrect.py: 0.566 / IRAF: 0.54)

It means that, in terms of barycentric velocity, astropy and the barycorr function (http://astroutils.astronomy.ohio-state.edu/exofast/barycorr.html) return comparable values within the 3 m/s precision error reported in both implementations. And, in terms of heliocentric velocity, both astropy and my implementation are giving similar results. I thus believe that there might be a precision problem in IRAF's calculation or some other subtlety I am not aware of.

For the sign problem, in IRAF's function, the longitude is positive towards west which is the opposite for other conventions.

You pointed out a good problem which is the use of this module rvcorrect.py in ORCS which might have not been completely checked and compared for a large range of input. In fact I will remove it and consider using the astropy implementation which is certainly far more robust as it is checked by numerous people ;)

Thank you !

@Klaus130161https://github.com/klaus130161 Please tell me if it answers your questions so that I can close this issue

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/thomasorb/orcs/issues/18#issuecomment-356435727, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AhDf6XW-Do6eTZAZ-oM5ps6svX870QeLks5tI-iVgaJpZM4RF8WP.