Rafael-Ayala / asteRisk

astrodynamics in R
5 stars 0 forks source link

targetTime in sgdp4 #3

Closed helgasoft closed 2 years ago

helgasoft commented 2 years ago

From the documentation, about targetTime in sgdp4: It can be provided in two different ways: either as a number corresponding to the time in minutes counting from epoch at which the orbit should be propagated, or as a date-time string in UTC, in which case the date-time string for epoch must be provided through initialDateTime. I see different orbits for target1 and target2 in code below

mins <- 15  # minutes into the future
target1 <- mins
target2 <- as.character(as.POSIXlt(Sys.time(), tz='UTC') + 60*mins)
sgdp4(..., targetTime=target1 )

Could you please double check whether they are indeed equivalent ? UPDATE: actually here is some comparison code:

library(asteRisk)
    # The following orbital parameters correspond to an object with NORAD catalogue
    # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC.
n0 <- 1.007781*((2*pi)/(1440))  # Multiplication by 2pi/1440 to convert to radians/min
e0 <- 0.002664 # mean eccentricity at epoch
i0 <- 3.8536*pi/180 # mean inclination at epoch in radians
M0 <- 48.3*pi/180 # mean anomaly at epoch in radians
omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians
OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians
Bstar <- 1e-04 # drag coefficient
epochDateTime <- "2006-06-26 00:58:29.34"

# Let´s calculate the position and velocity of the satellite 1 day later

state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0,
            Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440)
target2 <- as.character(as.POSIXlt(epochDateTime, tz='UTC') + 60*1440)
state_1day_TEME2 <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0,
            Bstar=Bstar, initialDateTime=epochDateTime, targetTime=target2)
cat('1) number:',state_1day_TEME$position, '\n2) date-time:',state_1day_TEME2$position)

# 1) number:    5500.959 41590.28 138.3334 
# 2) date-time: 5501.997 41590.14 138.2627

SOLVED

as.POSIXlt cuts off the seconds and 2006-06-26 00:58:29.34 becomes 2006-06-26 00:58:29, hence the slight difference. An extra line options("digits.secs"=2) needs to be added on top of the code forcing as.POSIXlt to include seconds and then results become identical.

Rafael-Ayala commented 2 years ago

Very interesting observation and glad to know you managed to solve it! I am a bit surprised that that was required thought. The following code for example, uses as.POSIXlt to convert a date-time to a POSIXlt object, which are then converted to Julian dates. The difference between the 2 is calculated and then converted back to seconds, resulting in the expected 0.2 seconds difference:

date1=as.numeric(julian(as.POSIXlt("2020-04-19 14:00:00.5", tz="UTC")))
date2=as.numeric(julian(as.POSIXlt("2020-04-19 14:00:00.3", tz="UTC")))
abs(date1 - date2)*86400

I did not have to set up any options to get the correct result. It might be due to some differences in the system local environment? Would you mind trying the above and letting me know the results? I use as.POSIXlt/as.POSIXct on some other places in the package, and it would be great to know if indeed I should ensure expected behavior in all locales by setting that option explicitly.

Thanks!

helgasoft commented 2 years ago

Here is what I get:

> date1=as.numeric(julian(as.POSIXlt("2020-04-19 14:00:00.5", tz="UTC")))
> date2=as.numeric(julian(as.POSIXlt("2020-04-19 14:00:00.3", tz="UTC")))
> abs(date1 - date2)*86400
[1] 0.2000002

Your example above produces same results regardless of the value of options("digits.secs"). Maybe because it deals with numbers only. But command as.POSIXlt returns a character vector without sub-seconds if options("digits.secs") is not set:

> epochDateTime <- "2006-06-26 00:58:29.34"

> options("digits.secs"=NULL)
> as.POSIXlt(epochDateTime, tz='UTC')
[1] "2006-06-26 00:58:29 UTC"

> options("digits.secs"=2)
> as.POSIXlt(epochDateTime, tz='UTC')
[1] "2006-06-26 00:58:29.34 UTC

From the ?POSIXlt help docs: Sub-second Accuracy Classes "POSIXct" and "POSIXlt" are able to express fractions of a second. (Conversion of fractions between the two forms may not be exact, but will have better than microsecond accuracy.) Fractional seconds are printed only if options("digits.secs") is set: see strftime.

Rafael-Ayala commented 2 years ago

I see, that makes sense! I have done some more tests, and indeed the loss of second fractions happens only when converting to character string (either explicitly or implicitly). For example

library(nanotime)
as.nanotime(as.POSIXct(900 * 604800 + 678.5, origin="1980-01-06 00:00:00", tz="UTC"))
## Correctly keeps the fractional seconds
as.character(as.POSIXct(900 * 604800 + 678.5, origin="1980-01-06 00:00:00", tz="UTC"))
## Fractional seconds are lost due to conversion to character type

Will certainly keep this in mind!