Konrad Komorowski - Orbital Determination Of Asteroid 1998 Yp11

ORBITAL DETERMINATION OF ASTEROID 1998 YP11 Konrad Komorowski, 7/23/08

ABSTRACT This document describes the process of determining the orbital elements of asteroid 1998 YP11. The observations were carried out by “Cretaceous Crustaceans” team (Juhee Bae, Marie Fulmer, Konrad Komorowski) at Etscorn Observatory (New Mexico Tech, 34° 4' 21" N; 106° 54' 52" W, IAU code: 719) between 6/17/08 and 7/12/08 as a part of the Summer Science Program 2008 curriculum. Four observations were performed with the following results: # UT time Right ascension Declination s h m s s 1 2008-06-18 05:44:18.359 ± 60 19 23 27.082 ± 0.051 + 36° 14' 57.29" ± 0.74" 2 2008-06-22 05:47:58.234 ± 60 s 19h 22m 14.749s ± 0.052s + 33° 5' 42.47" ± 0.41" s h m s s 3 2008-07-01 06:35:44.546 ± 60 19 18 48.490 ± 0.045 + 25° 10' 47.96" ± 0.59" s h m s s 4 2008-07-12 09:01:12.171 ± 45 19 14 4.994 ± 0.073 + 14° 26' 2.99" ± 0.44" Fig. 1 – equatorial coordinates of 1998 YP11 asteroid as observed at four different nights from Frank T. Etscorn Campus Observatory in Socorro, NM (IAU code: 719) Orbital elements calculated using Gaussian orbit determination method on data from Fig. 1: Element Value Uncertainty a [AU] 1.686156 0.078135 e 0.378604 0.025090 i [°] 14.567170 1.051213 Ω [°] 145.393673 0.986805 73.650606 1.604025 ω [°] M [°] 0 0 d Perihelion 2008-04-27 5 19h passage time 07:54:01 28m 35.23s Fig. 2 – calculated orbital elements of the asteroid 1998 YP11

-1-

INTRODUCTION The goal of the project is to determine the orbital elements of a near-Earth asteroid 1998 YP11. It will be achieved by performing observations of the asteroids at 4 different nights. Various Python programming language scripts (the core one performing Gaussian orbit determination) will be then employed to determine a set of 6 standard orbital elements of the asteroid.

MATERIALS AND METHODS Materials All the data was collected using CCD imaging with a Takahashi 6” refractor. The properties of the CCD chip are summarized in Fig. 3. Resolution: 2184 px × 1472 px ē max absorption: 77 000 ē/px Read noise: 8.8 ē Pixel size: 6.8 μm Gain: 1.3 ē/ADU Plate scale: 169.2 “/mm Fig. 3 – properties of the CCD chip used Observations were normally carried out with a time span of 8 ± 4 days. Each of the analyzed images was taken with a 2×2 binning in order to reduce noise and simplify analysis, which reduced the images resolution to 1092 px × 736 px. Usual exposure time was 120 sec, with the exception of the last observation when a shorter exposure of 90 sec was needed to avoid recording asteroid’s trail – an effect of its faster motion against the background of stars.

Method Data reduction and observations preparation were performed using a collection of 4 scripts written in the Python programming language. These included: 1. Asteroid’s position interpolation script (own algorithm). 2. Plate reduction script: proper motion data application (own algorithm), image rotation and asteroid’s position calculation (algorithm by Dr. Kevin Krisciunas) with a full-blown uncertainty estimation (own algorithm). 3. Orbital elements determination script (Gaussian method) with a Monte Carlo error analysis and output normalization (own algorithm). 4. Ephemeris generator script (algorithm by Dr. Agnes Kim) with a Monte Carlo error analysis.

-2-

Positions of the asteroid on consecutive observation nights were predicted using a Python interpolation script (linear interpolation between the previous and next ephemeris data). Ephemeris data was supplied by the Summer Science Program 2008 at Socorro faculty. Accurate positions of the asteroid on the obtained CCD images were then determined using a plate reduction script written in the Python programming language. The script’s algorithm was prepared by Dr. Kevin Krisciunas. The script itself included modules reading data from an external CSV1 file, weighting the objects’ positions based on the brightness of adjacent pixels and accounting for the proper motion of the stars. Gaussian method was the core of the orbital determination script, also written in Python. Its computer implementation had been prepared by Dr. Agnes Kim. Monte Carlo method was adopted for error analysis – the probability distribution function used for randomizing input variables was a standard Gaussian curve; and the same p.d.f. was assumed for the results interpretation. Gaussian method was refined by topocentric to geocentric coordinates correction and account for the travel time of light. The orbital elements were then inserted into the ephemeris generator script. The ephemeris was then generated for the time of the other observation, to make sure that the calculated elements are correct.

                                                             1

Comma Separated Value

-3-

DATA AND ANALYSIS Data collection The positions of the asteroid (see Appendix 5) were determined accurately from CDD images taken at four different nights (see Appendices 1 – 4) using the plate reduction script:

#1

40.00 35.00 30.00 25.00 20.00 #4 15.00 10.00 5.00 0.00

#2 #3

19.400

19.375

19.350

19.325

19.300

19.275

19.250

Dec (°)

Dec vs RA diagram for 1999 YP11

19.225

RA (h)

Fig. 4 – declination vs. right ascension graph for 1998 YP11 observations

Perihelion passage time Classical orbital elements are given with mean anomaly, M, at the time of perihelion. However, the crude mean anomaly received from the Gaussian method corresponds to the time of the middle observation of the object and does not have to be equal to 0° (perihelion). The OD founds the time (t2) of previous passage though the perihelion (M2 = 0) using the following relation: ‫ܯ‬ଵ ൌ ‫ܯ‬ଶ ൅ ට

ߤ ሺ‫ ݐ‬െ ‫ݐ‬ଶ ሻ ܽଷ ଵ

‫ݐ‬ଶ ൌ ‫ݐ‬ଵ െ

‫ܯ‬ଵ ට

ߤ ܽଷ

Data analysis Four different orbit determinations were performed – each for a different set of observations. The results of each of them are summarized in Appendices 6 – 9.

-4-

In each of the 4 cases Monte Carlo analysis was performed. It turned out that in some iterations of the Monte Carlo analysis the script would encounter errors – non-converging series, significantly different from the one of others or values of cosine/sine greater than one. In such cases the particular partial results were ignored. Unavoidably, such omissions did contribute to the skewness of the general result.2 To counteract this unwanted result, the script mirrored each data (here: orbital elements) against the normal value after obtaining the primary Monte Carlo set of results. Obtained in such a way new set of partial results was then added to the old one. As the skewness of the new set was opposite to the skewness of the old one, their sum added to zero. Such a set of results was called normalized. It was then noticed that for calculations performed in a shorter time span or with irregular dates rate of such error omissions is high. What directly corresponds with the skewness: Skew/s.d.m. ratio for partial orbital data for 4 obs. sets 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

a e i Ω

ω M #1 : #2 : #3 (5.8%)

#1 : #2 : #4 (38.6%)

#1 : #3 : #4 (0.2%)

#2 : #3 : #4 (10.5%)

Fig. 5 – skew as a fraction of uncertainty (absolute value) for primary partial orbital data for 4 observations sets; in brackets – percentage of failed iterations Sets for which the results were significantly distorted (as #1:#2:#4 and #2:#3:#4) are not scientifically valuable. Therefore, they should not be included in the final result. The uncertainty itself is also a very important factor. Standard deviation from the mean of the partial data was assumed to be also the result’s uncertainty:

                                                             2

Non-skewed result was considered to be the one obtained with zero error analysis, i.e. after performing the OD for input variables as they were – with no randomization.

-5-

S.d.m. of partial orbital data for 4 observations sets 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

a [AU] e i [°] Ω [°]

ω [°] M [°] #1 : #2 : #3

#1 : #2 : #4

#1 : #3 : #4

#2 : #3 : #4

Fig. 6 – Standard deviation from the mean of partial orbital data for 4 observations sets Skewness check excluded sets #1:#2:#4 and #2:#3:#4 (to a lesser extend). Uncertainty check excluded sets #1:#2:4 and #2:#3:#4. The only relatively precise and unskewed orbital elements calculation had been performed based on observations set #1:#3:#4. This was to be expected, as observations #1, #3 and #4 were performed against the longest time span. Also, they were the most regular ones. These two factors minimized statistical errors resulting from distribution of observations. Consequently, the orbital elements will be calculated basing on observations #1, #3 and #4. Observation #2 will be used for sanity check of the elements. Element Value Uncertainty a [AU] 1.686156 0.078135 e 0.378604 0.025090 i [°] 14.567170 1.051213 Ω [°] 145.393673 0.986805 73.650606 1.604025 ω [°] M(t) [°] 29.146860 2.226832 Perihelion 2008-04-27 5d 19h 28m passage time 07:54:00.851276 35.232159s Fig. 7 – orbital elements calculated from observations #1, #3, #4, 5 000 iterations, 99.8% success rate, middle observation on 2008-07-01 at 06:35:44.546 ± 60 s

-6-

Data consistency check The orbital elements obtained were then inserted into the ephemeris generator script, to determine the position of the asteroid at the time of observation #2. The following result was obtained: Coordinate RA Dec

Ephemeris Ephemeris Observed Observed value value uncertainty uncertainty 19h 22m 15.747s 1h 43m 7.223s 19h 22m 14.749s 0.052s + 33° 7’ 1.91” 10° 31’ 24.49” + 33° 5' 42.47" 0.41" Fig. 8 – coordinates of the asteroid for the time of observation #2

Δ 0.988s 1’ 19.44”

Notice the very high ephemeris uncertainty as compared to the actual discrepancy between the ephemeris and observation. It is because uncertainty in the ephemeris generator is calculated regardless of the source of orbital elements data. It is only a coincidence that the ephemeris point was very close in time to the data points used to calculate orbital elements. This radical closeness in time minimized the error. However, if the ephemeris point was a few months apart from the time of measurements taken at the Etscorn observatory, then the actual discrepancy between the prediction and observation could have been as much as 1h 43m in RA and 10° 31’ in Dec.

-7-

CONCLUSION Data accuracy check The orbital elements are certainly correct within the specified error bonds according to every possible internal measure of verification (i.e. consistency of orbital elements received from different combinations of observations, ephemeris check for the other observation). The most reliable external method of checking the calculated orbital elements is Jet Propulsion Laboratory’s Horizons ephemeris computation system.3 Query for 1998 YP11’s orbital elements on 2008-04-27 07:54 (calculated perihelion passage revealed): JPL Horizons Calculated Calculated Δ value value uncertainty a [AU] 1.720593 1.686156 -0.034437 0.078135 e 0.388646 0.378604 -0.010042 0.025090 i [°] 15.030003 14.567170 -0.462833 1.051213 Ω [°] 144.928843 145.393673 0.464830 0.986805 74.493758 73.650606 -0.843152 1.604025 ω [°] M [°] 359.835702 0 0.164298 Fig. 9 – calculated orbital elements of the asteroid 1998 YP11 compared

Element

It can be seen that the difference between the JPL Horizons system values (assumed to be very precise) and the ones calculated in this paper is always smaller than the estimated uncertainty. It means that measurements, calculations were accurate and error estimation was honest. This means that the results are scientifically significant and worth of submission to the Minor Planet Center. Note on uncertainty of M: in my orbital elements the uncertainty was calculated for the time of perihelion passage, i.e. when M = 0° (without any uncertainty). However, there is no possibility of querying the JPL Horizons system by the value of M to return a specific date, so this uncertainty could not be directly tested. However in the primary calculation of orbital elements based on observations #1, #3 and #4 the uncertainty in M was 2.226832°. Assuming a similar value for 2008-04-27 07:54, this orbital element is also fully confirmed by JPL Horizons system.

Source of errors Lack of precision of the observational techniques was an evident source of error. It was revealed in the discrepancies between calculated orbital elements and the ones taken from JPL Horizons system. However, it was accurately predicted in the Monte Carlo error analysis and                                                              3

http://ssd.jpl.nasa.gov/?horizons

-8-

confirmed in its magnitude by the Horizons system. The solution for a smaller error would be to use 1x1 image binning, do shorter exposures and more advanced pixels brightness weighting to obtain more accurate positions of the asteroid and reference stars on the CCD images. Also, the fact that RA lines converge at the celestial poles, what distorts the interpolated system of coordinates, should be taken into account to achieve higher measurements precision. Another source of error is the very short time of observations – less than 6 weeks in total. It has been shown by skewness/s.d.m. partial data analysis how profoundly unevenly spaced/shorttermed observations affect the uncertainty in the result. The only solution to that problem would be to observe the asteroid in a greater time span. Orbital elements have been determined using only one set of three observations, as the rest of possible combinations have been found to distort the results/lower their accuracy in an unacceptable way. More observations would make it possible to perform combination of multiple results to achieve higher accuracy.

-9-

APPENDICES # 2UCAC RA J2000 Dec J2000 μRA μDec X Y 1 44518844 19 23 10.784 +36 14 22.32 -6.0 3.8 227.01 140.01 2 44518824 19 23 07.809 +36 14 07.08 -3.1 -1.9 243.00 146.00 3 44518837 19 23 09.358 +36 13 07.63 -3.4 -2.8 235.01 171.97 4 44518849 19 23 12.253 +36 13 25.07 0.2 -7.8 220.02 164.00 5 44518914 19 23 23.614 +36 13 03.11 9.7 10.7 160.00 174.00 6 44518897 19 23 20.383 +36 14 10.28 0.8 -10.9 177.01 144.99 7 44518881 19 23 17.832 +36 15 03.62 1.0 -9.0 191.01 121.00 8 44518903 19 23 21.223 +36 15 33.66 10.6 2.6 173.00 108.00 9 44518884 19 23 18.175 +36 16 18.54 6.9 -18.4 189.00 89.00 10 44518833 19 23 08.945 +36 15 38.86 8.0 -4.1 237.00 105.98 11 44518826 19 23 07.904 +36 16 03.20 4.5 -1.1 243.05 96.06 12 44518778 19 23 00.498 +36 15 56.23 3.8 -12.1 282.01 99.00 As 1998 YP11 221.00 139.00 Appendix 1 – contents of the CCD image from observation #1 on 2008-06-18 at 05:44:18.359 ±60 s (proper motion is measured in mas/yr) # 2UCAC RA J2000 Dec J2000 μRA μDec X Y 1 43479541 19 22 26.571 +33 05 55.48 -10.7 -2.7 77.97 118.00 2 43479545 19 22 26.989 +33 07 09.51 7.0 2.0 75.98 86.00 3 43479478 19 22 15.779 +33 07 41.27 0.4 -5.1 118.00 62.99 4 43479434 19 22 07.820 +33 07 50.00 3.0 4.7 148.03 61.97 5 43479464 19 22 13.645 +33 08 06.22 -1.6 -3.3 180.00 68.99 6 43479441 19 22 09.212 +33 05 40.49 -3.9 -2.7 172.03 124.99 7 43479427 19 22 06.911 +33 05 29.21 0.0 -12.1 184.83 129.97 8 43479406 19 22 03.782 +33 05 05.90 -0.4 0.8 201.97 140.01 9 43652885 19 22 06.080 +33 04 35.32 13.5 3.9 189.01 153.04 10 43479499 19 22 19.234 +33 04 57.60 1.6 0.2 117.04 143.03 11 43479536 19 22 26.027 +33 04 49.36 6.5 1.0 80.05 146.96 12 43479515 19 22 22.224 +33 03 31.15 -5.7 -7.6 101.00 180.97 As 1998 YP11 141.00 124.00 Appendix 2 – contents of the CCD image from observation #2 on 2008-06-22 at 05:47:58.234 ±60 s (proper motion is measured in mas/yr) # 1 2 3 4 5

2UCAC 40732161 40732143 40732134 40732120 40732097

RA J2000 19 18 50.895 19 18 49.370 19 18 48.786 19 18 46.016 19 18 43.675

Dec J2000 +25 11 01.82 +25 11 55.42 +25 11 26.91 +25 12 13.83 +25 11 51.89

- 10 -

μRA -6.0 3.8 -0.5 11.4 10.9

μDec -6.9 -1.1 -4.2 -2.9 -10.2

X 63.99 72.99 76.00 92.99 106.99

Y 57.00 33.99 46.01 25.99 35.01

# 2UCAC RA J2000 Dec J2000 μRA μDec X Y 6 40732020 19 18 5.580 +25 12 21.43 0.7 -9.0 153.99 22.98 7 40732125 19 18 47.330 +25 10 31.25 -3.3 -9.8 85.00 70.00 8 40732083 19 18 42.479 +25 10 02.94 -5.9 3.1 113.00 82.96 9 40732099 19 18 44.011 +25 09 58.51 1.1 1.8 104.00 84.99 10 40732167 19 18 51.887 +25 09 06.37 2.4 -1.2 57.98 107.00 11 40732194 19 18 54.721 +25 09 33.12 -0.5 -6.4 41.00 95.00 12 40732197 19 18 55.097 +25 10 38.32 4.6 6.9 38.98 66.99 As 1998 YP11 77.00 63.00 Appendix 3 – contents of the CCD image from observation #3 on 2008-06-18 at 05:44:18.359 ±60 s (proper motion is measured in mas/yr) # 2UCAC RA J2000 Dec J2000 μRA μDec X Y 1 36877103 19 14 10.694 +14 25 01.74 9.6 0.4 135.00 78.00 2 36877065 19 14 04.850 +14 25 34.85 -2.5 -0.4 98.01 91.99 3 36877044 19 14 02.788 +14 25 19.35 4.6 6.2 86.00 84.99 4 36877039 19 14 01.680 +14 25 56.42 4.9 -1.8 78.01 100.98 5 36877038 19 14 01.672 +14 27 01.48 -4.6 -5.1 78.00 129.00 6 36877030 19 14 00.893 +14 27 32.11 6.7 -4.1 72.01 142.00 7 36877067 19 14 05.025 +14 26 54.16 13.8 0.0 99.00 126.00 8 36877079 19 14 06.168 +14 26 49.09 5.9 -16.4 106.00 124.00 9 36877089 19 14 07.810 +14 26 44.91 10.2 -1.9 116.00 122.00 10 36877147 19 14 17.605 +14 25 52.51 4.8 -0.1 179.00 100.99 As 1998 YP11 99.00 104.00 Appendix 4 – contents of the CCD image from observation #4 on 2008-06-18 at 05:44:18.359 ±60 s (proper motion is measured in mas/yr) # UT time Right ascension Declination s h m s s 1 2008-06-18 05:44:18.359 ± 60 19 23 27.082 ± 0.051 + 36° 14' 57.29" ± 0.74" 2 2008-06-22 05:47:58.234 ± 60 s 19h 22m 14.749s ± 0.052s + 33° 5' 42.47" ± 0.41" s h m s s 3 2008-07-01 06:35:44.546 ± 60 19 18 48.490 ± 0.045 + 25° 10' 47.96" ± 0.59" 4 2008-07-12 09:01:12.171 ± 45 s 19h 14m 4.994s ± 0.073s + 14° 26' 2.99" ± 0.44" Appendix 5 – equatorial coordinates of 1998 YP11 asteroid as observed at four different nights from Frank T. Etscorn Campus Observatory in Socorro, NM (IAU code: 719) Element

Initial value

Uncertainty

Skew

a [AU] e i [°] Ω [°] ω [°]

1.655172 0.361543 14.065205 145.649424 73.191064

0.167479 0.058366 2.306594 1.977235 3.395554

-0.001292 -0.006200 -0.238857 +0.272037 -0.441414

- 11 -

Normalized value 1.656465 0.367742 14.304062 145.377387 73.632477

Normalized uncertainty 0.167475 0.058691 2.318807 1.995760 3.423949

Normalized Normalized value uncertainty M(t) [°] 26.469491 4.676070 0.679451 25.790040 4.724935 Perihelion 2008-04-22 13d 18h 5m -5d 1h 32m 2008-04-21 16d 3h 8m passage time 12:30:11.081461 34.585347s 47.068368s 03:07:24.730182 46.839596s Appendix 6 – orbital elements calculated from observations #1, #2, #3, 5 000 iterations, 94.2% success rate, middle observation on 2008-06-22 at 05:48:58.000234 ± 60 s Element

Initial value

Uncertainty

Skew

Normalized Normalized value uncertainty a [AU] 1.563310 0.076771 -0.071829 1.635139 0.105133 e 0.333427 0.030312 -0.026868 0.360295 0.040505 i [°] 12.980434 1.221062 -1.087472 14.067906 1.635096 Ω [°] 146.443520 1.140067 0.983712 145.459808 1.505785 71.796584 1.886706 -1.626096 73.422680 2.490722 ω [°] M(t) [°] 28.515636 2.534486 2.178086 26.337550 3.341766 d h m Perihelion 2008-04-27 7d 6h 35m -1 16 42 2008-04-27 9d 1h 44m passage time 01:34:22.388478 7.322609s 29.897610s 16:32:57.769458 04.944411s Appendix 7 – orbital elements calculated from observations #1, #2, #4, 5 000 iterations, 61.4% success rate, middle observation on 2008-06-22 at 05:48:58.000234 ± 60 s Element

Initial value

Uncertainty

Skew

Normalized Normalized value uncertainty a [AU] 1.685026 0.078131 -0.001129 1.686156 0.078135 e 0.377115 0.025047 -0.001488 0.378604 0.025090 i [°] 14.506337 1.049504 -0.060833 14.567170 1.051213 Ω [°] 145.476670 0.983357 +0.082997 145.393673 0.986805 73.533011 1.599789 -0.117595 73.650606 1.604025 ω [°] M(t) [°] 29.315447 2.220552 0.168586 29.146860 2.226832 d h m h m Perihelion 2008-04-28 4 14 14 22 58 2008-04-27 5d 19h 28m s s passage time 11:36:59.627291 00.851720 40.543812 07:54:00.851276 35.232159s Appendix 8 – orbital elements calculated from observations #1, #3, #4, 5 000 iterations, 99.8% success rate, middle observation on 2008-07-01 at 06:35:44.546 ± 60 s Element

Initial value

Uncertainty

Skew

Element

Initial value

Uncertainty

Skew

a [AU] e i [°] Ω [°] ω [°] M(t) [°]

1.585063 0.340312 13.031847 146.916947 71.282602 32.676157

0.128419 0.048136 2.000277 2.092853 3.286904 4.549375

-0.015778 -0.009486 -0.382916 +0.492963 -0.730204 +1.010619

- 12 -

Normalized value 1.600841 0.349798 13.414763 146.423985 72.012806 31.665538

Normalized uncertainty 0.129377 0.049059 2.036493 2.150019 3.366866 4.660039

Normalized Normalized value uncertainty d h m d h m Perihelion 2008-04-26 14 2 09 -2 20 46 2008-04-26 14d 5h 31m s s passage time 01:37:09.664462 11.877144 36.182565 18:08:46.086318 12.435960s Appendix 9 – orbital elements calculated from observations #2, #3, #4, 5 000 iterations, 89.5% success rate, middle observation on 2008-07-01 at 06:35:44.546 ± 60 s Element

Initial value

Uncertainty

Skew

Appendices 10 - 25 included electronically: 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.

1_2008_06_18_i1_cropped_prfv160.csv (application/vnd.ms-excel) 4K 1_2008_06_22_i2_cropped_prfv160.csv (application/vnd.ms-excel) 2K 1_2008_07_01_i2_cropped_prfv160.csv (application/vnd.ms-excel) 2K 1_2008_07_12_i6_cropped_prfv160.csv (application/vnd.ms-excel) 2K 1_2008_06_18_i1_cropped_prv160_output.txt (text/plain) 4K 1_2008_06_22_i2_cropped_prv160_output.txt (text/plain) 4K 1_2008_07_01_i2_cropped_prv160_output.txt (text/plain) 4K 1_2008_07_12_i6_cropped_prv160_output.txt (text/plain) 3K 1998_YP11_OD_123.csv (application/vnd.ms-excel) 1K 1998_YP11_OD_124.csv (application/vnd.ms-excel) 1K 1998_YP11_OD_134.csv (application/vnd.ms-excel) 1K 1998_YP11_OD_234.csv (application/vnd.ms-excel) 1K KKAsteroidsPositionInterpolation_1.0.0.py (text/plain) 8K KKEphemerisGenerator_1.1.0.py (text/plain) 20K PlateReduction_1.6.0_Team1.py (text/plain) 31K KKOrbitalDetermination_1.4.1.py (text/plain) 52K

- 13 -