Blog Archive

Showing posts with label Horizons Web Interface. Show all posts
Showing posts with label Horizons Web Interface. Show all posts

Sunday, February 12, 2017

2016 WF9 - a simulation based on Feb 11th orbital params

This is a simulation of asteroid 2016 WF9 based on orbital parameters published on HORIZONS Web-Interface on Feb 11th, 2017.
 
At this time, the orbital condition code is 4

Orbital Elements at Epoch 2457800.5 (2017-Feb-16.0) TDB
Reference: JPL 13 (heliocentric ecliptic J2000)

 Element Value Uncertainty (1-sigma)   Units 
e .6580360863711519 9.9808e-06
a 2.870895540801432 8.6058e-05 au
q .9817426747520661 8.7514e-07 au
i 14.99562014245548 9.8595e-05 deg
node 125.4263986134548 0.00019086 deg
peri 342.4337260195148 8.1316e-05 deg
M 3.004319337477674 0.00013463 deg
tp 2457785.672494105614
(2017-Feb-01.17249411)
9.1832e-05 JED
period 1776.742590373345
4.86
0.07989
0.0002187
d
yr
n .2026179830159605 9.1105e-06 deg/d
Q 4.760048406850797 0.00014269 au

This asteroid was simulated with the Mercury6 orbit simulator together with 100 virtual clones generated with the package R.

These 100 clones were generated so that their orbital parameters are normally distributed around the nominal value of asteroid 2016 WF9 and their standard deviation is almost equal to the uncertainty shown above.

Virtual Asteroids: summary


mean sd
a 2.8709 8.68E−05
e 0.65804 9.51E−06
i 14.9956 9.48E−05
w 342.434 8.18E−05
om 125.426 0.00019
M 3.00432 0.00014

Simulation parameters
  • period simulated: past 1e8 days 
  • time step 0.1 days
  • ejection distance = 100 au
  • N-body algorithm: Conservative Bulirsch-Stoer

Simulation results
  • 1 out of 101 virtual asteroid was discarded because it would have collided with Jupiter
  • 44 out of 101 virtual asteroids came from the outskirt of the solar system (i.e., there was a time in the past when their distance was more than 100 au - ejection distance from the simulator point of view), so it makes sense to think they were comets.

This is the density distribution of arrival time in the solar system:


 Actually, according to the simulation, the first cometary event was 30084 years ago and the last one occurred 273009 years ago.

Probability of being a comet
I tried to use a R survival package called survminer to display the probability of a virtual asteroid being a comet, a probability that increases as you go more and more in the past.

I built a table with three columns:
  • virtual asteroid id
  • year: time of arrival into solar system, or end-time of the simulation (right censored data)
  • event: in case of arrival from the outskirt of the solar system, this flag is TRUE, otherwise it is FALSE coherently with column year.
The Surv function available in R can be used to convert the table in order to have a formal "survival" object and the survfit function can then be used to display a summary table showing how the various asteroids were progressively lost in favor of a comet - the lost of an asteroid is equivalent to the "death" event in a survival analysis).

These are the first lines of the survfit result:

Year(in the past) n.risk n.event survival std.err lower 95% CI upper 95% CI
30084 101 1 0.990 0.010 0.971 1.000
67281 100 1 0.980 0.014 0.953 1.000
70233 99 1 0.970 0.017 0.938 1.000
73898 98 1 0.960 0.019 0.923 0.999
74074 97 1 0.950 0.022 0.909 0.994
75356 96 1 0.941 0.024 0.896 0.988
80292 95 1 0.931 0.025 0.882 0.982
90790 94 1 0.921 0.027 0.870 0.975

How to read:
we started with 101 virtual asteroids, then at Year 30084 in the past we see one of them coming from the outskirt of the solar system, so we do not count it as an asteroid ...then, at Year 67281 in the past, this event occurs again ... and so on till we arrive at the last virtual asteroid being counted as a comet at Year 273009 in the past (last row not shown).

Based on this, we can draw a plot showing the probability of a virtual asteroid being a comet as a function of time:




Kind Regards,
Alessandro Odasso

Citations:
Survminer
Alboukadel Kassambara and Marcin Kosinski (2016). survminer: Drawing Survival Curves using 'ggplot2'
      Mercury Simulator - Mercury6
      J.E.Chambers (1999) "A Hybrid
      Symplectic Integrator that Permits Close Encounters between
      Massive Bodies''. Monthly Notices of the Royal Astronomical
      Society, vol 304, pp793-799.

Monday, January 9, 2017

2016 WF9 - a simulation based on Jan 5th orbital params

On Jan 5th, I downloaded the orbital parameters of asteroid 2016 WF9 and their uncertainty from the Horizons Web interface:

param     value        1-sigma

a             2.8729806 0.00200340
e             0.6582656 0.00022291
i            14.9994893 0.00329270
w        342.4286190 0.00216310
om      125.4222393 0.00428050
M        353.6898676 0.00394660

In spite of the fact that the uncertainty is still very high, I tried to make a simulation generating 100 virtual asteroids to see how many of them in the past would seem to come from a distance greater than 100 AU (cometary origin).

First, I used the R programming environment to generate the virtual asteroids.
An obvious check is that the mean value of their orbital parameters must be almost equal to the nominal values of asteroid 2016 WF9 and their standard deviation should be also almost equal to the uncertainty shown above.

In fact, looking at the virtual asteroids, I got:

a     mean     2.8729   sd 0.002
e     mean     0.6583   sd 2e-04
i      mean   14.9988   sd 0.0035
w    mean 342.4285   sd 0.0023
om  mean 125.4218   sd 0.0046
M    mean 353.6903   sd 0.004


As a second step I used the Mercury6 simulator by John E. Chambers, to simulate the past 1e8 days, output interval 100 days (hybrid algorithm).
The program simulated the behaviour of 101 asteroids (100 virtual asteroid + the "real" asteroid):

The result was that about 60 virtual asteroids would seem to have a cometary orgin (i.e. there was a time in the past when their distance from the sun was grater than 100 AU).

The earliest date is about 37500 years ago, the oldest date is constrained by the simulation period (about 275000 years ago) and the median value is 155000 years ago.

(for curiosity: the nominal orbital values of 2016 WF9 itself tells us that this can be a comet arrived 76000 years ago).

This is the graph showing the arrival time of the 60 virtual comets:



In the next weeks, when the orbit of asteroid  2016 WF9 will be better defined, I would like to repeat the simulation to see if and how much the simulation results will be different.

Kind Regards,
Alessandro Odasso


Citations

      Mercury Simulator - Mercury6
      J.E.Chambers (1999) ``A Hybrid
      Symplectic Integrator that Permits Close Encounters between
      Massive Bodies''. Monthly Notices of the Royal Astronomical
      Society, vol 304, pp793-799.

Thursday, October 6, 2016

Amor asteroids - H mag vs Tisserand parameter

As seen in the previous blog where I focused on Apollos, I am now going to do the same for Amors: look for any relation between H mag and the Tisserand parameter with respect to a generic body with semi-major axis ap.


T_{P}\ ={\frac  {a_{P}}{a}}+2\cdot {\sqrt  {{\frac  {a}{a_{P}}}(1-e^{2})}}\cos i

The objective is to see what happens as ap changes in a given range.

Two simple ways to do this:
  • method1: calculate the (spearman) correlation between H mag and Tp, plot the result for many values of ap.
  • method2: divide the H mag and Tp range into quartiles; then, calculate the Chi-squared statistics on the resulting contingency table. Plot the result for many values of ap.

These two methods give almost the same result in the case of Apollo asteroids.
What is very strange, at least for me, is that these two methods give a different answer in the case of Amor asteroids.

In the following sections, I will describe in more detail:
  • Apollo - method1
  • Apollo - method2
  • Amor - method1
  • Amor - method2 
  • appendix (more graphs showing H mag vs others orbital parameters)

Apollo - method1
We calculate the spearman correlation and plot it against many values of ap. We search for minimum or maximum.
This method and the graphical result has been described in the previous blog.

Apollo - method2
I did not applied this approach in the previous blog, so I describe it in detail and I show the result.

First, we divide the H mag into quartiles (H-quartiles).
Then we calculate the Tp param for all asteroids and we divide the resulting range into quartiles (Tp-quartiles).
After that, we build a contincency table showing H-quartiles vs Tp-quartiles and, finally, we calculate the Chi-squared statistic.

The idea is just this: let's imagine for a moment that there is absolutely no relation between H mag and Tp. In such a case the Chi-squared statistic would be zero because there would be no difference between the observed values and the expected values, the latter being, by definition, the values that would be displayed when no relation exists. On the contrary, the greater the Chi-squared statistic and the greater the difference of H mag distributions between the Tp quartiles.

In order to implement all the above in R, besides reading the orbital parameters values downloaded from Horizons Web-Interface in a dataframe, we need a couple of user defined functions:

func_t
function(ap,a,e,i){ap/a+2*(a/ap*(1-e^2))^(1/2)*cos(i/180*pi)}

The above function just implements the Tisserand formula.

func_t_h_chisq
function(ap){p$ts<-func_t(ap,p$a,p$e,p$i);p$tsquartile<-factor(cut(p$ts, fivenum(p$ts), include.lowest = T));return(chisq.test(with(p,table(tsquartile,hquartile)))$statistic)}

All asteroids parameters are in a dataframe called p.
The above function calculates the Tp values (calling the first function) for all asteroids in the dataframe p: then it calculates the quartiles (cut and fivenum instructions), the contincency table (table instruction) and the Chi-squared statistic.

In order to see the result, all we have to do is define a range for ap (lets say 0 <= ap<=30 AU), loop calling the function and plot it:

therange<-seq(0.01,30,0.01)
plot(therange,sapply(therange,function(x) {func_t_h_chisq(x)}),type='l',xlab="ap (tisserand param)",ylab="H-quartiles vs Tp-quartiles - Chi-squared statistic",main="Apollo asteroids")



The maximum correlation is found for ap=1.03 AU
I think that this result is almost the same we achieved with method1 where we got ap = 1.01 AU.
More important: the plot seems smooth, there is a clear maximum.

Amor - method1
See previous blog for R instructions.
This is the plot showing how the correlation varies depending on the choiche of ap:
The maximum correlation in the case of Amors is found for ap = 2.07 AU
However, the overall behaviour is "flat" compared to the bigger variations that we saw in the case of Apollos.
In fact, the boxplot shows some differences (but not so many!) in the distribution of H depending on Tp quartiles:



Amor - method2
Following the same approach already described for Apollo asteroids, in order to see the result, all we have to do is define a range for ap (lets say 0 <= ap<=30 AU), loop calling the function and plot it:

therange<-seq(0.01,30,0.01)

plot(therange,sapply(therange,function(x) {func_t_h_chisq(x)}),type='l',xlab="ap (tisserand param)",ylab="H-quartiles vs Tp-quartiles - Chi-squared statistic",main="Amor asteroids")


This is surprising!
There is a relative minimum for ap = 1.29 AU
There are two relative maxima for ap = 1.0 AU and ap = 1.94 AU

We do not get once again the value ap=2.07 AU that was found before using the correlation method.

Let's see the corresponding boxplots to see if we have smaller differences for ap = 1.29 AU and much greater differences for ap = 1.0 AU and ap = 1.94 AU

boxplot ap=1.29 AU (minimum)
This boxplot below is done for ap=1.29 AU and, in fact, the differences of H distributions depending on different quartiles are small (with a notable exception for the 1st Tp quartile, where Amor asteroids are brighter.

boxplot ap=1.00 AU (maximum)
This boxplot below shows what happens for ap =1 .00 AU where we can see more differences compared to the previous one (although it seems that Q1 and Q4 are more similar):


boxplot ap=1.94 AU (maximum)
This boxplot below shows what happens for ap =1 .94 AU

In this other case, there is still some difference although this time Q2 and Q3 are more similar.
The trend seems consistent with what we have see before for ap=2.07 that was found using the correlation method.

appendix
As already done for Apollo in the previous blog, let's see other graphs for Amors:
  • H mag vs inclination
  • H mag vs perihelium

H mag vs inclination
Higher inclined Amors are brighter


H mag vs perihelium
Amors with greater perihelium (Q4 quartile) are brighter.
This is the opposite compared with Apollo belonging to their Q4 quartile that, on the contrary, were typically darker.


Cheers,
Alessandro Odasso