Occultation Class
The Occultation Class within SORA was created to reduce and analyze stellar occultations, and manage all the other Objects in a SORA occultation. Here we have some useful tasks that allow converting the times for each observatory in positions in the sky plane (\(\xi\), \(\eta\)), fit an ellipse to the points, obtain the astrometrical position resulting, among others.
The documentation here contains the details about every step.
This Jupyter-Notebook was designed as a tutorial for how to work with the Occultation Class. More information about the other classes, please refer to their specif Jupyter-Notebook. Any further question, please contact the core team: Altair Ramos Gomes Júnior, Bruno Eduardo Morgado, Gustavo Benedetti Rossi, and Rodrigo Carlos Boufleur.
The Occultation Docstring was designed to help the users. Also, each function has its Docstring containing its main purpose and the needed parameters (physical description and formats). Please, do not hesitate to use it.
0. Index
[1]:
## SORA package
from sora import Occultation, Body, Star, LightCurve, Observer
## Other main packages
#from astropy.time import Time
import astropy.units as u
## Usual packages
import numpy as np
import matplotlib.pylab as pl
import os
SORA version: 0.3
1. Instantiating an Occultation Object and adding observations
The Occultation Class can be instantiated in only one way. For this it is needed a Star Object, an Body Object with an Ephemeris Object (EphemKernel, EphemPlanet or EphemHorizons), and the occultation time (within 60 minutes of the correct value).
[2]:
Occultation?
Init signature:
Occultation(
star,
body=None,
ephem=None,
time=None,
reference_center='geocenter',
)
Docstring:
Instantiates the Occultation object and performs the reduction of the
occultation.
Attributes
----------
star : `sora.Star`, `str`, required
the coordinate of the star in the same reference frame as the ephemeris.
It must be a Star object or a string with the coordinates of the object
to search on Vizier.
body : `sora.Body`, `str`
Object that will occult the star. It must be a Body object or its name
to search in the Small Body Database.
ephem : `sora.Ephem`, `list`
Object ephemeris. It must be an Ephemeris object or a list.
time : `str`, `astropy.time.Time`, required
Reference time of the occultation. Time does not need to be exact, but
needs to be within approximately 50 minutes of the occultation closest
approach to calculate occultation parameters.
reference_center : `str`, `sora.Observer`, `sora.Spacecraft`
A SORA observer object or a string 'geocenter'.
The occultation parameters will be calculated in respect
to this reference as center of projection.
Important
---------
When instantiating with "body" and "ephem", the user may define the
Occultation in 3 ways:
1. With `body` and `ephem`.
2. With only "body". In this case, the "body" parameter must be a Body
object and have an ephemeris associated (see Body documentation).
3. With only `ephem`. In this case, the `ephem` parameter must be one of the
Ephem Classes and have a name (see Ephem documentation) to search for the
body in the Small Body Database.
File: ~/Documentos/códigos/SORA/sora/occultation/core.py
Type: type
Subclasses:
First let’s instantiate the ``Star`` and the ``Body``
SORA will automatically search for the star information. A warning will raise when any information is missing. In this example, there is no star radius available in Gaia.
[3]:
chariklo = Body(name='Chariklo',
ephem=['input/bsp/Chariklo.bsp', 'input/bsp/de438_small.bsp'])
star_occ = Star(coord='18 55 15.65250 -31 31 21.67051') # Occ Chariklo 22-06-2017
Obtaining data for Chariklo from SBDB
1 GaiaDR3 star found band={'G': 14.223702}
star coordinate at J2016.0: RA=18h55m15.65210s +/- 0.0197 mas, DEC=-31d31m21.6676s +/- 0.018 mas
Downloading star parameters from I/297/out
Now, we can instantiate the Occultation
[4]:
occ = Occultation(star=star_occ, body=chariklo, time='2017-06-22 21:18')
print(occ)
Stellar occultation of star GaiaDR3 6760223758801661440 by 10199 Chariklo (1997 CU26).
Geocentric Closest Approach: 0.049 arcsec
Instant of CA: 2017-06-22 21:18:48.200
Position Angle: 359.72 deg
Geocentric shadow velocity: -22.00 km / s
Sun-Geocenter-Target angle: 166.42 deg
Moon-Geocenter-Target angle: 149.11 deg
No observations reported
###############################################################################
STAR
###############################################################################
GaiaDR3 star Source ID: 6760223758801661440
ICRS star coordinate at J2016.0:
RA=18h55m15.65210s +/- 0.0197 mas, DEC=-31d31m21.6676s +/- 0.0180 mas
pmRA=3.556 +/- 0.025 mas/yr, pmDEC=-2.050 +/- 0.020 mas/yr
GaiaDR3 Proper motion corrected as suggested by Cantat-Gaudin & Brandt (2021)
Plx=0.2121 +/- 0.0228 mas, Rad. Vel.=-40.49 +/- 3.73 km/s
Magnitudes: G: 14.224, B: 14.320, V: 13.530, R: 14.180, J: 12.395, H: 11.781,
K: 11.627
Apparent diameter from Kervella et. al (2004):
V: 0.0216 mas, B: 0.0216 mas
Apparent diameter from van Belle (1999):
sg: B: 0.0238 mas, V: 0.0244 mas
ms: B: 0.0261 mas, V: 0.0198 mas
vs: B: 0.0350 mas, V: 0.0315 mas
Geocentric star coordinate at occultation Epoch (2017-06-22 21:18:48.200):
RA=18h55m15.65251s +/- 0.0327 mas, DEC=-31d31m21.6706s +/- 0.0341 mas
###############################################################################
10199 Chariklo (1997 CU26)
###############################################################################
Object Orbital Class: Centaur
Spectral Type:
SMASS: D [Reference: EAR-A-5-DDR-TAXONOMY-V4.0]
Relatively featureless spectrum with very steep red slope.
Discovered 1997-Feb-15 by Spacewatch at Kitt Peak
Physical parameters:
Diameter:
302 +/- 30 km
Reference: Earth, Moon, and Planets, v. 89, Issue 1, p. 117-134 (2002),
Rotation:
7.004 +/- 0 h
Reference: LCDB (Rev. 2021-June); Warner et al., 2009, [Result based on less than full coverage, so that the period may be wrong by 30 percent or so.] REFERENCE LIST:[Fornasier, S.; Lazzaro, D.; Alvarez-Candal, A.; Snodgrass, C.; et al. (2014) Astron. Astrophys. 568, L11.], [Leiva, R.; Sicardy, B.; Camargo, J.I.B.; Desmars, J.; et al. (2017) Astron. J. 154, A159.]
Absolute Magnitude:
6.54 +/- 0 mag
Reference: MPO691682,
Albedo:
0.045 +/- 0.01
Reference: Earth, Moon, and Planets, v. 89, Issue 1, p. 117-134 (2002),
Ellipsoid: 151.0 x 151.0 x 151.0
----------- Ephemeris -----------
EphemKernel: CHARIKLO/DE438_SMALL (SPKID=2010199)
Ephem Error: RA*cosDEC: 0.000 arcsec; DEC: 0.000 arcsec
Offset applied: RA*cosDEC: 0.0000 arcsec; DEC: 0.0000 arcsec
If the occultation was observed by a spacecraft, the Occultation parameters should be refered to this observer. For that, it is necessary to pass the observer as the center of reference for the Occultation
[5]:
# from sora import Spacecraft
# nh = Spacecraft(name='New Horizons', spkid='-98', ephem='horizons')
# nh_occ = Occultation(star=star, body=body, time=time, reference_center=nh)
Note that at the print(Occultation) there are no observations added to this Occultation yet.
It is needed one Observer and one LightCurve to define an Chord.
Let’s instanciate the Observers
Here we will give 5 locations
[6]:
out = Observer(name='Outeniqua' ,lon='+16 49 17.710', lat='-21 17 58.170', height =1416)
ond = Observer(name='Onduruquea' ,lon='+15 59 33.750', lat='-21 36 26.040', height =1220)
tiv = Observer(name='Tivoli' ,lon='+18 01 01.240', lat='-23 27 40.190', height =1344)
whc = Observer(name='Windhoek' ,lon='+17 06 31.900', lat='-22 41 55.160', height =1902)
hak = Observer(name='Hakos' ,lon='+16 21 41.320', lat='-23 14 11.040', height =1843)
If a spacecraft observed the occultation, a Spacecraft object could be defined and added to the Occultation normally.
Let’s instanciate the LightCurves
Here we give 6 light curves. Note that 2 are for the same site observed with different telescopes.
[7]:
#THE ERRORS ARE INCREASED FOR BETTER VISUALIZATION
out_lc = LightCurve(name='Outeniqua lc',
initial_time='2017-06-22 21:20:00.056',
end_time ='2017-06-22 21:29:59.963',
immersion='2017-06-22 21:21:20.329',immersion_err=0.320,
emersion ='2017-06-22 21:21:30.343',emersion_err=0.340)
ond_lc = LightCurve(name='Onduruquea lc',
initial_time='2017-06-22 21:11:52.175',
end_time ='2017-06-22 21:25:13.389',
immersion='2017-06-22 21:21:22.213',immersion_err=0.100,
emersion ='2017-06-22 21:21:33.824',emersion_err=0.110)
tiv_lc = LightCurve(name='Tivoli lc',
initial_time='2017-06-22 21:16:00.094',
end_time ='2017-06-22 21:28:00.018',
immersion='2017-06-22 21:21:15.628',immersion_err=0.700,
emersion ='2017-06-22 21:21:19.988',emersion_err=0.700)
whc_c14_lc = LightCurve(name='Windhoek C14 lc',
initial_time='2017-06-22 21:12:48.250',
end_time ='2017-06-22 21:32:47.963',
immersion='2017-06-22 21:21:17.609',immersion_err=0.240,
emersion ='2017-06-22 21:21:27.564',emersion_err=0.260)
whc_d16_lc = LightCurve(name='Windhoek D16 lc',
initial_time='2017-06-22 21:20:01.884',
end_time ='2017-06-22 21:22:21.894',
immersion='2017-06-22 21:21:17.288',immersion_err=0.280,
emersion ='2017-06-22 21:21:27.228',emersion_err=0.340)
hak_lc = LightCurve(name='Hakos lc',
initial_time='2017-06-22 21:10:19.461',
end_time ='2017-06-22 21:30:19.345')
To add the observation to the Occultation just use the Occultation.chords.add_chord()
During this step the LightCurve model parameters are automaticaly updated using the Occultation parameters: the Shadow’s velocity, the star diameter and the distance to the occulting body. This means that is no need for the user to do the LightCurve.set_dist(), LightCurve.set_vel() and LightCurve.set_star_diam().
[8]:
occ.chords.add_chord(observer=out, lightcurve=out_lc)
occ.chords.add_chord(observer=ond, lightcurve=ond_lc)
occ.chords.add_chord(observer=tiv, lightcurve=tiv_lc)
/home/altair/Documentos/códigos/SORA/sora/body/core.py:332: UserWarning: H and/or G is not defined for 10199 Chariklo. Searching into JPL Horizons service
warnings.warn('H and/or G is not defined for {}. Searching into JPL Horizons service'.format(self.shortname))
[8]:
<Chord: Tivoli>
Note that a warning comes up if the flux drop is not calculated automatically.
The same goes for two observation at with the same Observer , however is important to define their names as different values
[9]:
occ.chords.add_chord(name='Windhoek C14', observer=whc, lightcurve=whc_c14_lc)
occ.chords.add_chord(name='Windhoek D16', observer=whc, lightcurve=whc_d16_lc)
[9]:
<Chord: Windhoek D16>
Also, the same goes for negative observations
[10]:
occ.chords.add_chord(observer=hak, lightcurve=hak_lc)
[10]:
<Chord: Hakos>
If an spacecraft observed the occultation, just define the Spacecraft object and add it to the Occutation as usual.
[11]:
# from sora import Spacecraft
# nh = Spacecraft(name='New Horizons', spkid='-98', ephem='horizons')
# occ.chords.add_chord(observer=nh, lightcurve=lightcurve)
To check which observers were added to the Occultation just use Occultation.chords
[12]:
occ.chords
[12]:
<ChordList:
0: Chord(Outeniqua)
1: Chord(Onduruquea)
2: Chord(Tivoli)
3: Chord(Windhoek C14)
4: Chord(Windhoek D16)
5: Chord(Hakos)>
If any of them was wrongly added just remove it using Occultation.chords.remove_chord()
[13]:
occ.chords.remove_chord(name='Outeniqua')
[14]:
occ.chords
[14]:
<ChordList:
0: Chord(Onduruquea)
1: Chord(Tivoli)
2: Chord(Windhoek C14)
3: Chord(Windhoek D16)
4: Chord(Hakos)>
If all the chords were wrongly added the user can clear the Occultation using Occultation.chords.clear()
[15]:
occ.chords.clear()
[16]:
occ.chords
[16]:
<ChordList:>
[17]:
occ.chords.add_chord(observer=out, lightcurve=out_lc)
occ.chords.add_chord(observer=ond, lightcurve=ond_lc)
occ.chords.add_chord(observer=tiv, lightcurve=tiv_lc)
occ.chords.add_chord(name='Windhoek C14', observer=whc, lightcurve=whc_c14_lc)
occ.chords.add_chord(name='Windhoek D16', observer=whc, lightcurve=whc_d16_lc)
occ.chords.add_chord(observer=hak, lightcurve=hak_lc)
[17]:
<Chord: Hakos>
2. Projecting the times in the sky plane and the Chords
This step is done automatically without the user having to ask for it. The user can see it automatically using Occultation.chords.summary() that contains all positions for all observations added.
[18]:
occ.chords.summary()
Name Longitude Latitude status time f g
------------ ------------ ------------- ------------ ----------------------- --------- -------
Outeniqua 16d49m17.71s -21d17m58.17s Initial Time 2017-06-22 21:20:00.056 -1912.82 45.76
Immersion 2017-06-22 21:21:20.329 -118.28 65.39
Emersion 2017-06-22 21:21:30.343 105.60 67.83
End Time 2017-06-22 21:29:59.963 11501.40 190.23
Onduruquea 15d59m33.75s -21d36m26.04s Initial Time 2017-06-22 21:11:52.175 -12876.38 -135.47
Immersion 2017-06-22 21:21:22.213 -138.25 7.19
Emersion 2017-06-22 21:21:33.824 121.28 10.05
End Time 2017-06-22 21:25:13.389 5029.49 63.76
Tivoli 18d01m01.24s -23d27m40.19s Initial Time 2017-06-22 21:16:00.094 -7123.66 -202.63
Immersion 2017-06-22 21:21:15.628 -70.56 -126.77
Emersion 2017-06-22 21:21:19.988 26.91 -125.73
End Time 2017-06-22 21:28:00.018 8971.45 -31.61
Windhoek C14 17d06m31.9s -22d41m55.16s Initial Time 2017-06-22 21:12:48.250 -11505.42 -198.50
Immersion 2017-06-22 21:21:17.609 -121.72 -73.52
Emersion 2017-06-22 21:21:27.564 100.82 -71.11
End Time 2017-06-22 21:32:47.963 15314.98 90.02
Windhoek D16 17d06m31.9s -22d41m55.16s Initial Time 2017-06-22 21:20:01.884 -1814.42 -91.87
Immersion 2017-06-22 21:21:17.288 -128.89 -73.59
Emersion 2017-06-22 21:21:27.228 93.31 -71.19
End Time 2017-06-22 21:22:21.894 1315.35 -58.00
Hakos 16d21m41.32s -23d14m11.04s Initial Time 2017-06-22 21:10:19.461 -14875.34 -316.83
End Time 2017-06-22 21:30:19.345 11939.87 -23.52
Let’s give a look to an specific Chord
[19]:
chord = occ.chords['Outeniqua']
The user can easily acess the f and g for the immersion and emersion times using the Chord.get_fg()
[20]:
chord.get_fg()
[20]:
(array([-118.27875072, 105.59793973]), array([65.3915666 , 67.83395663]))
To see the f and g projection for the immersion and emersion times used in the fitting process, the user can use the following function
[21]:
names, fg, error = occ.chords.get_limb_points()
print(names)
print(fg)
print(error)
['Outeniqua_immersion' 'Outeniqua_emersion' 'Onduruquea_immersion'
'Onduruquea_emersion' 'Tivoli_immersion' 'Tivoli_emersion'
'Windhoek C14_immersion' 'Windhoek C14_emersion' 'Windhoek D16_immersion'
'Windhoek D16_emersion']
[[-118.27875072 65.3915666 ]
[ 105.59793973 67.83395663]
[-138.25381491 7.1933912 ]
[ 121.275169 10.05237144]
[ -70.56164003 -126.77190052]
[ 26.91022935 -125.73348183]
[-121.71711407 -73.51645459]
[ 100.81847461 -71.10987823]
[-128.89276574 -73.59407832]
[ 93.30744752 -71.19108221]]
[[ 7.15400526 0.07806865]
[ 7.60120413 0.08290053]
[ 2.23518846 0.02463189]
[ 2.45872898 0.0270767 ]
[15.64912805 0.16673645]
[15.64918323 0.16669294]
[ 5.36497525 0.05803441]
[ 5.81210573 0.06283513]
[ 6.25913618 0.06770901]
[ 7.60044325 0.08217062]]
The user can also set a specific time
[22]:
chord.get_fg(time='2017-06-22 21:00:00.000')
[22]:
(-28724.98453460884, -258.35753813600354)
Finally, the user can visually see the chord in the Sky-plane using the Chord.plot_chord().
[23]:
chord.plot_chord(segment='positive', color='blue')
chord.plot_chord(segment='error', color='red')
pl.xlim(+250,-250)
pl.ylim(-250,+250)
pl.show()
Let’s define a LightCurve with times and fluxes
[24]:
chord2 = occ.chords['Hakos']
chord2.lightcurve.set_flux(file='input/lightcurves/lc_example_5.dat', exptime=1.000, usecols=[0,1])
The Chord.plot_chord() can also be used with an linestyle = "exposure" to see the times where the the chord in fact holds information and the readout time (without information).
[25]:
chord2.plot_chord(segment='negative', linestyle='exposure', color='green')
pl.xlim(+250,-250)
pl.ylim(-250,+250)
pl.show()
Now let’s consider all the chords
[26]:
occ.chords.plot_chords()
occ.chords.plot_chords(segment='error', color='red')
pl.legend(loc=1)
pl.xlim(+170,-330)
pl.ylim(-250,+250)
pl.show()
If there are known time shifts, this can be easily solved using LightCurve.dt
[27]:
out_lc.dt = -0.150
ond_lc.dt = -0.190
tiv_lc.dt = -0.150
whc_c14_lc.dt = -0.375
whc_d16_lc.dt = +0.000
[28]:
occ.chords.plot_chords()
occ.chords.plot_chords(segment='error', color='red')
pl.legend(loc=1)
pl.xlim(+170,-330)
pl.ylim(-250,+250)
pl.show()
The user can save these positions
[29]:
occ.to_file()
3. Ellipse fit
The next step is the fitting of an ellipse. The five parameters that describe the ellipse are:
1. and 2. The centre position (\(f_0\), \(g_0\))
3. The apparent equatorial radius, semi-major axix (\(a'\))
4. The oblatness (\(\epsilon' = \frac{a' - b'}{a'}\))
5. The position angle of the pole, semi-minor axis (\(P\))
The result of the fit is a ChiSquare Object, and its functions can be found at its specific Jupyter-Notebook.
Here, there is only the manual method, the user should provide the parameters to the fit and also the region for searching each parameters.
The equation to be minimize is:
\(\chi^2 = \sum_{i}^{N} \frac{(r_{i} - r'_{i})^2}{\sigma_i^2 + \sigma_{model}^2}\)
where: - \(r_i\) is the radial distance between the \(i^{th}\) observed point and the ellipse centre; - \(r'_i\) is the radial distance between the modelled ellipse’s \(i^{th}\) point and the ellipse centre; - \(\sigma_i\) is the unceartainty of the \(i^{th}\) observed point - \(\sigma_{model}\) is the model uncertainty, that is releated to the difference between the real apparent shape of the occultating object and the ellipse model.
[30]:
occ.fit_ellipse?
Signature: occ.fit_ellipse(**kwargs)
Docstring:
Fits an ellipse to given occultation using given parameters.
Parameters
----------
center_f : `int`, `float`, default=0
The coordinate in f of the ellipse center.
center_g : `int`, `float`, default=0
The coordinate in g of the ellipse center.
equatorial_radius : `int`, `float`
The Equatorial radius (semi-major axis) of the ellipse.
oblateness : `int`, `float`, default=0
The oblateness of the ellipse.
position_angle : `int`, `float`, default=0
The pole position angle of the ellipse in degrees.
Zero is in the North direction ('g-positive'). Positive clockwise.
dcenter_f : `int`, `float`
Interval for coordinate f of the ellipse center.
dcenter_g : `int`, `float`
Interval for coordinate g of the ellipse center.
dequatorial_radius `int`, `float`
Interval for the Equatorial radius (semi-major axis) of the ellipse.
doblateness : `int`, `float`
Interval for the oblateness of the ellipse
dposition_angle : `int`, `float`
Interval for the pole position angle of the ellipse in degrees.
loop : `int`, default=10000000
The number of ellipses to attempt fitting.
dchi_min : `int`, `float`
If given, it will only save ellipsis which chi square are smaller than
chi_min + dchi_min. By default `None` when used with `method='chisqr`, and
`3` for other methods.
number_chi : `int`, default=10000
In the `chisqr` method if `dchi_min` is given, the procedure is repeated until
`number_chi` is reached.
In other methods it is the number of values (simulations) that should lie within
the provided `sigma_result`.
verbose : `bool`, default=False
If True, it prints information while fitting.
ellipse_error : `int`, `float`
Model uncertainty to be considered in the fit, in km.
sigma_result : `int`, `float`
Sigma value to be considered as result.
method : `str`, default=`least_squares`
Method used to perform the fit. Available methods are:
`chisqr` : monte carlo computation method used in versions of SORA <= 0.2.1.
`fastchi` : monte carlo computation method, allows multithreading .
`least_squares` or `ls`: best fit done used levenberg marquardt convergence algorithm.
`differential_evolution` or ``de`: best fit done using genetic algorithms.
All methods return a Chisquare object.
threads : `int`
Number of threads/workers used to perform parallel computations of the chi square
object. It works with all methods except `chisqr`, by default 1.
Returns
-------
chisquare : `sora.ChiSquare`
A ChiSquare object with all parameters.
Important
---------
Each occultation is added as the first argument(s) directly.
Mandatory input parameters: 'center_f', 'center_g', 'equatorial_radius',
'oblateness', and 'position_angle'.
Parameters fitting interval: 'dcenter_f', 'dcenter_g', 'dequatorial_radius',
'doblateness', and 'dposition_angle'. Default values are set to zero.
Search done between (value - dvalue) and (value + dvalue).
Examples
--------
To fit the ellipse to the chords of occ1 Occultation object:
>>> fit_ellipse(occ1, **kwargs)
To fit the ellipse to the chords of occ1 and occ2 Occultation objects together:
>>> fit_ellipse(occ1, occ2, **kwargs)
File: ~/Documentos/códigos/SORA/sora/occultation/core.py
Type: method
[31]:
ellipse_chi2 = occ.fit_ellipse(center_f=-15.046, center_g=-2.495, dcenter_f=10, dcenter_g=20,
equatorial_radius=138, dequatorial_radius=50, oblateness=0.093,
doblateness=0.20, position_angle=90, dposition_angle=90 ,loop=100000,
dchi_min=10,number_chi=20000)
Note that here we only used 100000 loops, but it is recommended at least 10 millions for a good sample for each parameter.
[32]:
print(ellipse_chi2)
Minimum chi-square: 0.062
Number of fitted points: 10
Number of fitted parameters: 5
Minimum chi-square per degree of freedom: 0.012
center_f:
1-sigma: -13.429 +/- 0.789
3-sigma: -13.488 +/- 4.379
center_g:
1-sigma: -1.969 +/- 3.795
3-sigma: 0.912 +/- 16.823
equatorial_radius:
1-sigma: 139.351 +/- 5.457
3-sigma: 152.721 +/- 25.282
oblateness:
1-sigma: 0.085 +/- 0.037
3-sigma: 0.140 +/- 0.140
position_angle:
1-sigma: 126.475 +/- 18.116
3-sigma: 125.088 +/- 89.978
[33]:
ellipse_chi2.plot_chi2()
Now, besides the chords we can plot the fitted ellipse
[34]:
from sora.extra import draw_ellipse
occ.chords.plot_chords()
occ.chords.plot_chords(segment='error', color='red')
#plotting the best fitted ellipse, in black
draw_ellipse(**ellipse_chi2.get_values())
pl.legend(loc=1)
pl.xlim(+170,-330)
pl.ylim(-250,+250)
pl.show()
Not just the fitted ellipse, but the user can also plot all the ellipses inside an sigma region, for instance within :math:`3sigma`
This step will, possible, plot large number of ellipses, this can spend some CPU time.
[35]:
occ.chords.plot_chords()
occ.chords.plot_chords(segment='error', color='red')
#plotting the best fitted ellipse, in black
draw_ellipse(**ellipse_chi2.get_values())
# ploting all the ellipses within 3-sigma, in gray
draw_ellipse(**ellipse_chi2.get_values(sigma=3),alpha=1.0,lw=1)
pl.legend(loc=1)
pl.xlim(+170,-330)
pl.ylim(-250,+250)
pl.show()
As can be seen, there are some ellipses that cross the Hakos negative chord. They are impossible solutions as the Hakos lightcurve show no evidence of an occultation. The User can filter this solutions using the filter_negative_chord() function.
[36]:
for from sora.occultation import filter_negative_chord
[37]:
filter_negative_chord?
Signature: filter_negative_chord(chord, chisquare, step=1, sigma=0)
Docstring:
Get points for the ellipse with the given input parameters.
Parameters
----------
chord : `sora.observer.Chord`
Chord object, must be associated to an Occultation to work.
chisquare : `sora.extra.ChiSquare`
Resulted ChiSquare object of fit_ellipse.
sigma : `int`, `float`
Uncertainty of the expected ellipse, in km.
step : `int`, `float`, `str`
If a number, it corresponds to the step, in seconds, for each point of
the chord path. The step can also be equal to ``'exposure'``. In this
case, the chord path will consider the lightcurve individual times and
exptime.
File: ~/Documentos/códigos/SORA/sora/occultation/utils.py
Type: function
[38]:
filter_chi2 = filter_negative_chord(chord=occ.chords['Hakos'], chisquare=ellipse_chi2, step=0.5)
print(filter_chi2)
Filter chord: Hakos |█████...................................| - 14%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |██████████████..........................| - 37%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |████████████████████████................| - 60%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |█████████████████████████████████.......| - 84%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |████████████████████████████████████████| - 100%
Minimum chi-square: 0.062
Number of fitted points: 10
Number of fitted parameters: 5
Minimum chi-square per degree of freedom: 0.012
center_f:
1-sigma: -13.429 +/- 0.789
3-sigma: -13.488 +/- 4.379
center_g:
1-sigma: -1.969 +/- 3.795
3-sigma: 0.912 +/- 16.823
equatorial_radius:
1-sigma: 139.351 +/- 5.457
3-sigma: 149.280 +/- 21.842
oblateness:
1-sigma: 0.085 +/- 0.037
3-sigma: 0.131 +/- 0.131
position_angle:
1-sigma: 126.475 +/- 18.116
3-sigma: 125.088 +/- 89.978
[39]:
occ.chords.plot_chords()
occ.chords.plot_chords(segment='error', color='red')
#plotting the best fitted ellipse, in black
draw_ellipse(**filter_chi2.get_values())
# ploting all the ellipses within 3-sigma, in gray
draw_ellipse(**filter_chi2.get_values(sigma=3),alpha=1.0,lw=1)
pl.legend(loc=4,ncol=2)
pl.xlim(+90,-90)
pl.ylim(-180,-100)
pl.show()
The user can set step = "exposure" and consider only the region where data was acquired, allowing ellipses that crosses the negative chord during its readout time.
[40]:
filter_2_chi2 = filter_negative_chord(chord = occ.chords['Hakos'], chisquare = ellipse_chi2, step='exposure')
print(filter_2_chi2)
Filter chord: Hakos |█████...................................| - 14%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |███████████████.........................| - 38%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |█████████████████████████...............| - 62%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |███████████████████████████████████.....| - 88%
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Filter chord: Hakos |████████████████████████████████████████| - 100%
Minimum chi-square: 0.062
Number of fitted points: 10
Number of fitted parameters: 5
Minimum chi-square per degree of freedom: 0.012
center_f:
1-sigma: -13.429 +/- 0.789
3-sigma: -13.488 +/- 4.379
center_g:
1-sigma: -1.969 +/- 3.795
3-sigma: 0.912 +/- 16.823
equatorial_radius:
1-sigma: 139.351 +/- 5.457
3-sigma: 149.280 +/- 21.842
oblateness:
1-sigma: 0.085 +/- 0.037
3-sigma: 0.131 +/- 0.131
position_angle:
1-sigma: 126.475 +/- 18.116
3-sigma: 125.088 +/- 89.978
[41]:
occ.chords.plot_chords(segment='positive')
occ.chords.plot_chords(segment='error', color='red')
chord2.plot_chord(segment='negative', linestyle='exposure', color='m', label='Hakos')
#plotting the best fitted ellipse, in black
draw_ellipse(**filter_2_chi2.get_values())
# ploting all the ellipses within 3-sigma, in gray
draw_ellipse(**filter_2_chi2.get_values(sigma=3),alpha=1.0,lw=1)
pl.legend(loc=4,ncol=2)
pl.xlim(+90,-90)
pl.ylim(-180,-100)
pl.show()
The radial velocity issue
During the LightCurve fitting, the user has to add the velocity for that light curve. In the initial process, only the geocentric velocity was determined, and only after the fit, a radial velocity can be correctly calculated.
Usually, the light curves’ features are dominated by the exposure time. However, in the cases that the fresnel diffraction plays a significant role in it, we suggest that the LightCurve fitting procedure should be redone with the correct velocities.
[42]:
occ.check_velocities()
Outeniqua - Velocity used: 22.004
Immersion Radial Velocity: 17.800
Emersion Radial Velocity: 18.284
Onduruquea - Velocity used: 22.004
Immersion Radial Velocity: 22.252
Emersion Radial Velocity: 22.246
Tivoli - Velocity used: 22.004
Immersion Radial Velocity: 8.626
Emersion Radial Velocity: 5.195
Windhoek C14 - Velocity used: 22.004
Immersion Radial Velocity: 18.170
Emersion Radial Velocity: 17.470
Windhoek D16 - Velocity used: 22.004
Immersion Radial Velocity: 18.107
Emersion Radial Velocity: 17.526
4. Viewing and saving the results
In the end, we can see the results using Python dictionaries. The Occultation.fitted_params will have fitted parameters and their \(1\sigma\) uncertainties. The Occultation.chi2_params will have some information about the fit and its quality.
[43]:
occ.fitted_params
[43]:
{'equatorial_radius': [139.3506109843088, 5.4567114516244],
'center_f': [-13.428622523637122, 0.7891892760699086],
'center_g': [-1.9686948831613296, 3.795102365340013],
'oblateness': [0.08488597520674623, 0.03706668452200804],
'position_angle': [126.47517825871333, 18.116387367826775]}
[44]:
occ.chi2_params
[44]:
{'chord_name': ['Outeniqua_immersion',
'Outeniqua_emersion',
'Onduruquea_immersion',
'Onduruquea_emersion',
'Tivoli_immersion',
'Tivoli_emersion',
'Windhoek C14_immersion',
'Windhoek C14_emersion',
'Windhoek D16_immersion',
'Windhoek D16_emersion'],
'radial_dispersion': array([-0.14982613, -2.3388093 , -1.21884535, -1.3706693 , -0.17265982,
-2.22234731, -0.58491025, -1.06365652, -1.66000562, -0.34293474]),
'position_angle': array([301.88964659, 58.90450147, 274.03961814, 84.75566445,
205.85090919, 163.36646323, 238.44925955, 123.18339268,
238.18760663, 122.96493764]),
'radial_error': array([ 7.15442539, 7.60165495, 2.23532561, 2.45888165, 15.65002209,
15.65007494, 5.36529313, 5.81244286, 6.25950441, 7.60088706]),
'chi2_min': 0.06220013599688508,
'nparam': 5,
'npts': 10}
The astrometrical positions obtained can be seen using Occultation.new_astrometric_position()
[45]:
occ.new_astrometric_position()
Ephemeris offset (km): X = -13.4 km +/- 0.8 km; Y = -2.0 km +/- 3.8 km
Ephemeris offset (mas): da_cos_dec = -1.263 +/- 0.074; d_dec = -0.185 +/- 0.357
Astrometric object position at time 2017-06-22 21:18:48.200 for reference 'geocenter'
RA = 18 55 15.6523925 +/- 0.081 mas; DEC = -31 31 21.622083 +/- 0.359 mas
If the user wants the position at a specific time
[46]:
occ.new_astrometric_position('2017-06-22 21:21:00.000')
Ephemeris offset (km): X = -13.4 km +/- 0.8 km; Y = -2.0 km +/- 3.8 km
Ephemeris offset (mas): da_cos_dec = -1.263 +/- 0.074; d_dec = -0.185 +/- 0.357
Astrometric object position at time 2017-06-22 21:21:00.000 for reference 'geocenter'
RA = 18 55 15.6310596 +/- 0.081 mas; DEC = -31 31 21.623477 +/- 0.359 mas
Also a ‘post-fit’ occultation map can be created using Occultation.plot_occ_map()
The function that plots the map contains a large number of kwargs. They can be used to completely control the map and its tutorial can be found at SORA documentation here.
[47]:
occ.plot_occ_map(centermap_delta=[-3500,+400],zoom=20,nameimg='output/map_posfit')
Projected shadow radius = 135.2 km
output/map_posfit.png generated
All this information can also be seen using print(Occultation)
[48]:
print(occ)
Stellar occultation of star GaiaDR3 6760223758801661440 by 10199 Chariklo (1997 CU26).
Geocentric Closest Approach: 0.049 arcsec
Instant of CA: 2017-06-22 21:18:48.200
Position Angle: 359.72 deg
Geocentric shadow velocity: -22.00 km / s
Sun-Geocenter-Target angle: 166.42 deg
Moon-Geocenter-Target angle: 149.11 deg
5 positive observations
1 negative observations
###############################################################################
STAR
###############################################################################
GaiaDR3 star Source ID: 6760223758801661440
ICRS star coordinate at J2016.0:
RA=18h55m15.65210s +/- 0.0197 mas, DEC=-31d31m21.6676s +/- 0.0180 mas
pmRA=3.556 +/- 0.025 mas/yr, pmDEC=-2.050 +/- 0.020 mas/yr
GaiaDR3 Proper motion corrected as suggested by Cantat-Gaudin & Brandt (2021)
Plx=0.2121 +/- 0.0228 mas, Rad. Vel.=-40.49 +/- 3.73 km/s
Magnitudes: G: 14.224, B: 14.320, V: 13.530, R: 14.180, J: 12.395, H: 11.781,
K: 11.627
Apparent diameter from Kervella et. al (2004):
V: 0.0216 mas, B: 0.0216 mas
Apparent diameter from van Belle (1999):
sg: B: 0.0238 mas, V: 0.0244 mas
ms: B: 0.0261 mas, V: 0.0198 mas
vs: B: 0.0350 mas, V: 0.0315 mas
Geocentric star coordinate at occultation Epoch (2017-06-22 21:18:48.200):
RA=18h55m15.65251s +/- 0.0327 mas, DEC=-31d31m21.6706s +/- 0.0341 mas
###############################################################################
10199 Chariklo (1997 CU26)
###############################################################################
Object Orbital Class: Centaur
Spectral Type:
SMASS: D [Reference: EAR-A-5-DDR-TAXONOMY-V4.0]
Relatively featureless spectrum with very steep red slope.
Discovered 1997-Feb-15 by Spacewatch at Kitt Peak
Physical parameters:
Diameter:
302 +/- 30 km
Reference: Earth, Moon, and Planets, v. 89, Issue 1, p. 117-134 (2002),
Rotation:
7.004 +/- 0 h
Reference: LCDB (Rev. 2021-June); Warner et al., 2009, [Result based on less than full coverage, so that the period may be wrong by 30 percent or so.] REFERENCE LIST:[Fornasier, S.; Lazzaro, D.; Alvarez-Candal, A.; Snodgrass, C.; et al. (2014) Astron. Astrophys. 568, L11.], [Leiva, R.; Sicardy, B.; Camargo, J.I.B.; Desmars, J.; et al. (2017) Astron. J. 154, A159.]
Absolute Magnitude:
6.54 +/- 0 mag
Reference: JPL Horizons,
Phase Slope:
0.15 +/- 0
Reference: JPL Horizons,
Albedo:
0.045 +/- 0.01
Reference: Earth, Moon, and Planets, v. 89, Issue 1, p. 117-134 (2002),
Ellipsoid: 151.0 x 151.0 x 151.0
----------- Ephemeris -----------
EphemKernel: CHARIKLO/DE438_SMALL (SPKID=2010199)
Ephem Error: RA*cosDEC: 0.000 arcsec; DEC: 0.000 arcsec
Offset applied: RA*cosDEC: 0.0000 arcsec; DEC: 0.0000 arcsec
###############################################################################
POSITIVE OBSERVATIONS
###############################################################################
-------------------------------------------------------------------------------
Site: Outeniqua
Geodetic coordinates: Lon: 16d49m17.71s, Lat: -21d17m58.17s, height: 1.416 km
Target altitude: 56.6 deg
Target azimuth: 115.3 deg
Light curve name: Outeniqua lc
Initial time: 2017-06-22 21:20:00.056 UTC
End time: 2017-06-22 21:29:59.963 UTC
Duration: 9.998 minutes
Time offset: -0.150 seconds
Object LightCurve was not instantiated with time and flux.
Bandpass: 0.700 +/- 0.300 microns
Object Distance: 14.66 AU
Used shadow velocity: 22.004 km/s
Fresnel scale: 0.040 seconds or 0.87 km
Stellar size effect: 0.010 seconds or 0.23 km
Object LightCurve model was not fitted.
Immersion time: 2017-06-22 21:21:20.179 UTC +/- 0.320 seconds
Emersion time: 2017-06-22 21:21:30.193 UTC +/- 0.340 seconds
-------------------------------------------------------------------------------
Site: Onduruquea
Geodetic coordinates: Lon: 15d59m33.75s, Lat: -21d36m26.04s, height: 1.220 km
Target altitude: 56.1 deg
Target azimuth: 114.7 deg
Light curve name: Onduruquea lc
Initial time: 2017-06-22 21:11:52.175 UTC
End time: 2017-06-22 21:25:13.389 UTC
Duration: 13.354 minutes
Time offset: -0.190 seconds
Object LightCurve was not instantiated with time and flux.
Bandpass: 0.700 +/- 0.300 microns
Object Distance: 14.66 AU
Used shadow velocity: 22.004 km/s
Fresnel scale: 0.040 seconds or 0.87 km
Stellar size effect: 0.010 seconds or 0.23 km
Object LightCurve model was not fitted.
Immersion time: 2017-06-22 21:21:22.023 UTC +/- 0.100 seconds
Emersion time: 2017-06-22 21:21:33.634 UTC +/- 0.110 seconds
-------------------------------------------------------------------------------
Site: Tivoli
Geodetic coordinates: Lon: 18d01m01.24s, Lat: -23d27m40.19s, height: 1.344 km
Target altitude: 58.5 deg
Target azimuth: 112.4 deg
Light curve name: Tivoli lc
Initial time: 2017-06-22 21:16:00.094 UTC
End time: 2017-06-22 21:28:00.018 UTC
Duration: 11.999 minutes
Time offset: -0.150 seconds
Object LightCurve was not instantiated with time and flux.
Bandpass: 0.700 +/- 0.300 microns
Object Distance: 14.66 AU
Used shadow velocity: 22.004 km/s
Fresnel scale: 0.040 seconds or 0.87 km
Stellar size effect: 0.010 seconds or 0.23 km
Object LightCurve model was not fitted.
Immersion time: 2017-06-22 21:21:15.478 UTC +/- 0.700 seconds
Emersion time: 2017-06-22 21:21:19.838 UTC +/- 0.700 seconds
-------------------------------------------------------------------------------
Site: Windhoek
Geodetic coordinates: Lon: 17d06m31.9s, Lat: -22d41m55.16s, height: 1.902 km
Target altitude: 57.4 deg
Target azimuth: 113.4 deg
Light curve name: Windhoek C14 lc
Initial time: 2017-06-22 21:12:48.250 UTC
End time: 2017-06-22 21:32:47.963 UTC
Duration: 19.995 minutes
Time offset: -0.375 seconds
Object LightCurve was not instantiated with time and flux.
Bandpass: 0.700 +/- 0.300 microns
Object Distance: 14.66 AU
Used shadow velocity: 22.004 km/s
Fresnel scale: 0.040 seconds or 0.87 km
Stellar size effect: 0.010 seconds or 0.23 km
Object LightCurve model was not fitted.
Immersion time: 2017-06-22 21:21:17.234 UTC +/- 0.240 seconds
Emersion time: 2017-06-22 21:21:27.189 UTC +/- 0.260 seconds
-------------------------------------------------------------------------------
Site: Windhoek
Geodetic coordinates: Lon: 17d06m31.9s, Lat: -22d41m55.16s, height: 1.902 km
Target altitude: 57.4 deg
Target azimuth: 113.4 deg
Light curve name: Windhoek D16 lc
Initial time: 2017-06-22 21:20:01.884 UTC
End time: 2017-06-22 21:22:21.894 UTC
Duration: 2.333 minutes
Time offset: 0.000 seconds
Object LightCurve was not instantiated with time and flux.
Bandpass: 0.700 +/- 0.300 microns
Object Distance: 14.66 AU
Used shadow velocity: 22.004 km/s
Fresnel scale: 0.040 seconds or 0.87 km
Stellar size effect: 0.010 seconds or 0.23 km
Object LightCurve model was not fitted.
Immersion time: 2017-06-22 21:21:17.288 UTC +/- 0.280 seconds
Emersion time: 2017-06-22 21:21:27.228 UTC +/- 0.340 seconds
###############################################################################
NEGATIVE OBSERVATIONS
###############################################################################
-------------------------------------------------------------------------------
Site: Hakos
Geodetic coordinates: Lon: 16d21m41.32s, Lat: -23d14m11.04s, height: 1.843 km
Target altitude: 56.8 deg
Target azimuth: 112.5 deg
Light curve name: Hakos lc
Initial time: 2017-06-22 21:10:19.961 UTC
End time: 2017-06-22 21:30:16.955 UTC
Duration: 19.950 minutes
Time offset: 0.000 seconds
Exposure time: 1.0000 seconds
Cycle time: 2.9998 seconds
Num. data points: 400
Bandpass: 0.700 +/- 0.300 microns
Object Distance: 14.66 AU
Used shadow velocity: 22.004 km/s
Fresnel scale: 0.040 seconds or 0.87 km
Stellar size effect: 0.010 seconds or 0.23 km
Object LightCurve model was not fitted.
Immersion and emersion times were not fitted or instantiated.
###############################################################################
RESULTS
###############################################################################
Fitted Ellipse:
equatorial_radius: 139.351 +/- 5.457
center_f: -13.429 +/- 0.789
center_g: -1.969 +/- 3.795
oblateness: 0.085 +/- 0.037
position_angle: 126.475 +/- 18.116
polar_radius: 127.522 km
equivalent_radius: 133.305 km
geometric albedo (V): 0.061 (6.1%)
Minimum chi-square: 0.062
Number of fitted points: 10
Number of fitted parameters: 5
Minimum chi-square per degree of freedom: 0.012
Radial dispersion: -1.112 +/- 0.803 km
Radial error: 7.579 +/- 4.655 km
Ephemeris offset (km): X = -13.4 km +/- 0.8 km; Y = -2.0 km +/- 3.8 km
Ephemeris offset (mas): da_cos_dec = -1.263 +/- 0.074; d_dec = -0.185 +/- 0.357
Astrometric object position at time 2017-06-22 21:18:48.200 for reference 'geocenter'
RA = 18 55 15.6523925 +/- 0.081 mas; DEC = -31 31 21.622083 +/- 0.359 mas
And this can be saved to an ASCII file using Occultation.to_log
[49]:
occ.to_log('output/Test_occ.log')
!ls output/*.log
output/LC_test_1a.log output/Test_occ.log
This Jupyter-Notebook was designed as a tutorial for how to work with the Occultation Class. More information about the other classes, please refer to their specif Jupyter-Notebook. Any further question, please contact the core team: Altair Ramos Gomes Júnior, Bruno Eduardo Morgado, Gustavo Benedetti Rossi, and Rodrigo Carlos Boufleur.
The End