why do RA and DEC differ from oce values?

[dankelley]

I think oce::moonAngle() and oce::sunAngle() are returning correct values, since they accurately predict eclipses, etc. But a quick test suggests that ephemeris::ephemeris() is producing different results (even after converting from hour-angle to angle).

It seems that my test is wrong, oce is wrong, or ephemeris is wrong.

library(oce)
#> Loading required package: gsw
#> Loading required package: testthat
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.0, PROJ 7.2.0
library(ephemeris)
fix <- function(x)
{
    while (x < 0)
        x <- x + 360
    while (x > 360)
        x <- x - 360
    x
}
t <- as.POSIXct("2020-01-01", tz="UTC")

d <- list()

oce <- moonAngle(t)
eph <- ephemeris("s:Moon", t0=t, nbd=1)
d$oceMoonDEC <- fix(oce$declination)
d$ephMoonDEC <- fix(eph$DECdec * 360 / 24)
d$oceMoonRA <- fix(oce$rightAscension)
d$ephMoonRA <- fix(eph$RAdec * 360 / 24)

oce <- sunAngle(t)
eph <- ephemeris("s:Sun", t0=t, nbd=1)
d$oceSunDEC <- fix(oce$declination)
d$ephSunDEC <- fix(eph$DECdec * 360 / 24)
d$oceSunRA <- fix(oce$rightAscension)
d$ephSunRA <- fix(eph$RAdec * 360 / 24)
str(d)
#> List of 8
#>  $ oceMoonDEC: num 350
#>  $ ephMoonDEC: num 209
#>  $ oceMoonRA : num 349
#>  $ ephMoonRA : num 349
#>  $ oceSunDEC : num 337
#>  $ ephSunDEC : num 350
#>  $ oceSunRA  : num 281
#>  $ ephSunRA  : num 275

^{Created on 2020-12-27 by the reprex package (v0.3.0)}

[dankelley]

Let's check oce::moonAngle() against a website.

oce

> m <- oce::moonAngle("2020-12-28",lon=0,lat=0)
> m$declination
[1] 21.56793
> m$rightAscension
[1] 70.93225

In-the-sky Website

https://in-the-sky.org/ephemeris.php?ird=1&irs=1&ima=1&iob=1&objtype=1&objpl=Moon&objtxt=the+Moon&tz=1&startday=27&startmonth=12&startyear=2020&interval=1&rows=25 gives

Let's convert this to degrees:

> DEC <- 21+31/60+43/3600
> DEC
[1] 21.52861
> RA <- 15*(4+42/60+27/3600) # right ascension is in hours and 15=360/24
> RA
[1] 70.6125

Comparison

> RA-m$rightAscension
[1] -0.3197498
> DEC-m$declination
[1] -0.03931779

Conclusion

I don't think oce::moonAngle() is wildly wrong.

[dankelley]

Oh. I see. I was doing declination as though it was in hour angles, but it is actually in angles.

Now I think the moon calculations are okay.

It makes sense to look into the sun values, though.

library(oce)
#> Loading required package: gsw
#> Loading required package: testthat
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.0, PROJ 7.2.0
library(ephemeris)
fix <- function(x)
{
    while (x < 0)
        x <- x + 360
    while (x > 360)
        x <- x - 360
    x
}
t <- as.POSIXct("2020-01-01", tz="UTC")

d <- list()

oce <- moonAngle(t)
eph <- ephemeris("s:Moon", t0=t, nbd=1)
d$oceMoonDEC <- fix(oce$declination)
d$ephMoonDEC <- fix(eph$DECdec)
d$oceMoonRA <- fix(oce$rightAscension)
d$ephMoonRA <- fix(eph$RAdec * 360 / 24)

oce <- sunAngle(t)
eph <- ephemeris("s:Sun", t0=t, nbd=1)
d$oceSunDEC <- fix(oce$declination)
d$ephSunDEC <- fix(eph$DECdec)
d$oceSunRA <- fix(oce$rightAscension)
d$ephSunRA <- fix(eph$RAdec * 360 / 24)
print(as.data.frame(d))
#>   oceMoonDEC ephMoonDEC oceMoonRA ephMoonRA oceSunDEC ephSunDEC oceSunRA
#> 1   350.0205   349.9189  349.1567  348.9149  336.9421  335.3551 280.9003
#>   ephSunRA
#> 1 274.5263

^{Created on 2020-12-27 by the reprex package (v0.3.0)}

[dankelley]

oce::sunAngle() test. (Note that I figured out how to give the lon and lat, which I had not done in the in-the-sky tests above.)

Conclusion In this test, oce::sunAngle() agreed with the in-the-sky website to within 0.5deg for both RA and DEC.

https://in-the-sky.org/ephemeris.php?ird=1&irs=1&ima=1&iob=1&objtype=1&objpl=Sun&objtxt=the+Sun&tz=1&startday=1&startmonth=1&startyear=2020&interval=1&rows=25

DECoce <- oce::sunAngle("2020-01-02",lon=0,lat=0,useRefraction=FALSE)$declination
RAoce <- oce::sunAngle("2020-01-02",lon=0,lat=0,useRefraction=FALSE)$rightAscension
DECits <- -(23+00/60+03/3600)
RAits <- 15*(18+46/60+48/3600)-360
DECoce
#> [1] -22.97744
DECits
#> [1] -23.00083
DECoce - DECits
#> [1] 0.02339639
RAoce
#> [1] -77.99575
RAits
#> [1] -78.3
RAoce - RAits
#> [1] 0.3042516

^{Created on 2020-12-27 by the reprex package (v0.3.0)}

[dankelley]

Since I have verified oce::moonAngle() and oce::sunAngle() to within 0.5deg for both RA and DEC, I think ephemeris() must be setting up its query incorrectly.

[dankelley]

I added teph as an arg but varying it through possible values does not change much

library(ephemeris)
for(teph in 1:4)
    cat("teph=", teph,
        unlist(ephemeris("s:Sun",t0="2020-01-02",nbd=1,teph=teph)[c("RA","DEC")]),"\n")
#> teph= 1  +18 25 03.32428  -24 40 17.1715 
#> teph= 2  +18 26 14.28111  -24 39 32.3612 
#> teph= 3  +18 26 18.96535  -24 39 33.6084 
#> teph= 4  +18 25 05.26670  -24 40 18.4777

^{Created on 2020-12-27 by the reprex package (v0.3.0)}

[dankelley]

Oh, this is interesting. I was ignoring the header, assuming that if the URL responded to the query, I was getting what I wanted. So, I was skipping the header. Bad idea.

As make test5 in sandbox shows, what I thought would give the sun (-name=s:Sun) is actually giving ... Mercury. That's why the RA and DEC are near what I get in oce::sunAngle().

Now I just have to figure out how to specify the sun in a query. I've looked a bit in the docs but cannot find it so far. Maybe there's no such thing as ephemeris data for the sun? (I'm out of my depth here.)

Details

make test5 does

curl "https://ssp.imcce.fr/webservices/miriade/api/ephemcc.php?-name=s:Sun&-type=&-ep='2020-01-01 12h'&-nbd=5&-step=1.5d&-tscale=UTC&-observer=@500&-theory=INPOP&-teph=1&-tcoor=1&-oscelem=astorb&-mime=text&-from=MiriadeDoc" > test5.txt

and the result is

# Flag: 1
# Ticket: 168731610823280439
#################################################################################
# Solar system object ephemeris by IMCCE/OBSPM/CNRS
# Satellite : 1 Mercury
# Planetary theory: INPOP19A
# Frame: J2000 astrometric
# Coordinates: equatorial
# Frame center: geocenter
# Precession/nutation model: IAU2000
# Relativistic perturbations: yes
# PPN coordinate system: beta=gamma=1, alpha=0
#################################################################################
#         Date UTC        |       RA      |      DEC      |     Dobs     |  VMag   |  Phase  |  Elong.  |    dRAcosDEC    |      dDEC      |     RV
#              h  m  s    |  h  m  s      |   o  '  "     |      au      |   mag   |   deg   |   deg    |   arcsec/min    |   arcsec/min   |    km/s
  2020-01-01T00:00:00.000  +18 18  6.32143  -24 38 41.7913    1.433992355     -0.93     12.24       5.77        0.3941E+01      -0.9462E-01      3.99999
  2020-01-02T12:00:00.000  +18 28 32.36802  -24 40 33.9972    1.436839687     -1.00     10.59       4.99        0.3961E+01      -0.8985E-02      2.57291
  2020-01-04T00:00:00.000  +18 39  1.51678  -24 39 19.3861    1.438448215     -1.07      8.98       4.22        0.3980E+01       0.7834E-01      1.13889
  2020-01-05T12:00:00.000  +18 49 33.41978  -24 34 54.4062    1.438809012     -1.15      7.43       3.48        0.3998E+01       0.1673E+00     -0.30880
  2020-01-07T00:00:00.000  +19 00  7.71665  -24 27 15.7109    1.437907187     -1.22      6.00       2.79        0.4016E+01       0.2577E+00     -1.77722

[dankelley]

Oh, hang on. The notation is that the sun is ... a planet. So -name=p:Sun works. I guess my grade-school teachers had it wrong!

[dankelley]

Here's what I get now. Note that oce and ephemeris now agree to under 0.5deg on both RA and DEC, for both sun and moon. Since my tests with websites suggest that this disagreement is close to that of the disagreement between oce and websites, it is clear that I am on to something.

I will look into this some more later today. I commented out the test, which was just a consistency test anyway. I'll replace that with a test against the website.

In the meantime, the commit is a5ab7c85fd44cabcd9a0392478b248c7e6e03837 (clicking that gives changes from previous commit).

library(oce)
#> Loading required package: gsw
#> Loading required package: testthat
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.0, PROJ 7.2.0
library(ephemeris)
fix <- function(x)
{
    while (x < 0)
        x <- x + 360
    while (x > 360)
        x <- x - 360
    x
}
t <- as.POSIXct("2020-01-01", tz="UTC")

d <- list()

oce <- moonAngle(t)
eph <- ephemeris("s:Moon", t0=t, nbd=1)
d$oceMoonDEC <- fix(oce$declination)
d$ephMoonDEC <- fix(eph$DECdec)
d$oceMoonRA <- fix(oce$rightAscension)
d$ephMoonRA <- fix(eph$RAdec * 360 / 24)

oce <- sunAngle(t)
eph <- ephemeris("p:Sun", t0=t, nbd=1)
d$oceSunDEC <- fix(oce$declination)
d$ephSunDEC <- fix(eph$DECdec)
d$oceSunRA <- fix(oce$rightAscension)
d$ephSunRA <- fix(eph$RAdec * 360 / 24)
print(as.data.frame(d))
#>   oceMoonDEC ephMoonDEC oceMoonRA ephMoonRA oceSunDEC ephSunDEC oceSunRA
#> 1   350.0205   349.9189  349.1567  348.9149  336.9421  336.9207 280.9003
#>   ephSunRA
#> 1 280.5973

cat(sprintf("Moon  RA: %10.3f (oce) %10.3f (ephemeris); %10.3f (diff)\n",
            d$oceMoonRA, d$ephMoonRA, d$oceMoonRA - d$ephMoonRA))
#> Moon  RA:    349.157 (oce)    348.915 (ephemeris);      0.242 (diff)
cat(sprintf("Moon DEC: %10.3f (oce) %10.3f (ephemeris); %10.3f (diff)\n",
            d$oceMoonDEC, d$ephMoonDEC, d$oceMoonDEC - d$ephMoonDEC))
#> Moon DEC:    350.021 (oce)    349.919 (ephemeris);      0.102 (diff)
cat(sprintf("Sun   RA: %10.3f (oce) %10.3f (ephemeris); %10.3f (diff)\n",
            d$oceSunRA, d$ephSunRA, d$oceSunRA - d$ephSunRA))
#> Sun   RA:    280.900 (oce)    280.597 (ephemeris);      0.303 (diff)
cat(sprintf("Sun  DEC: %10.3f (oce) %10.3f (ephemeris); %10.3f (diff)\n",
            d$oceSunDEC, d$ephSunDEC, d$oceSunDEC - d$ephSunDEC))
#> Sun  DEC:    336.942 (oce)    336.921 (ephemeris);      0.021 (diff)

^{Created on 2020-12-28 by the reprex package (v0.3.0)}

[dankelley]

Sun tests: RA and DEC agree to all digits shown on the in-the-sky.org site.

ephemeris::ephemeris(name="p:Sun",t0="2020-01-02",nbd=1)
#>         time               RA             DEC      Dobs   VMag Phase Elong.
#> 1 2020-01-02  +18 46 48.43576  -23 00 03.1820 0.9832732 -26.78  2.54 0.2056
#>   dRAcosDEC dDEC RV    RAdec    DECdec
#> 1  -0.03015   NA NA 18.78012 -23.00088

https://in-the-sky.org/ephemeris.php?ird=1&irs=1&ima=1&iob=1&objtype=1&objpl=Sun&objtxt=the+Sun&tz=1&startday=1&startmonth=1&startyear=2020&interval=1&rows=25

^{Created on 2020-12-28 by the reprex package (v0.3.0)}

Moon

[dankelley]

This is working now in commit 1621aac2144ec3f8aedc5cb3de909273345309a5, so I am closing the issue.

I added a test suite that shows that the sun values agree with the in-the-sky website to within 1 arc-second (i.e. to the precision of in-the-sky), whereas the moon values agree to within 3 arc-seconds for RA and 15 arc-seconds for DEC.

I am not particularly worried about these differences (note that DEC disagreement of 15 arc-seconds is 15/3600 = 0.004 degrees, not a large value for e.g. a tidal computation). If I had the interest and time, I suppose it might be worth exploring how the results vary with e.g. the teph value (see ?ephemeris).