LightCurve Class

The LightCurve Class within SORA was created to reduce and analyze stellar occultations’ light curves. Here we present some useful tasks that allows the user: - set the model’s parameters - prepare the light curve to be reduced - fit one or all parameters simultaneously.

The documentation here contains the details about every step.

This Jupyter-Notebook was designed as a tutorial for how to work with the LightCurve Class. Any further question, please contact the core team: Altair Ramos Gomes Júnior, Bruno Eduardo Morgado, Gustavo Benedetti Rossi, and Rodrigo Carlos Boufleur.

The LightCurve Docstring was designed to help the users. Also, each function has its Docstring containing its main purpose and the needed parameters to set (physical description and formats). Please, do not hesitate to use it.

0. Index

1. Instantiating the LightCurve Object

2. Set and/or modify model parameters

3. Preparing the LightCurve

   3.1 Light curve normalization using a polinomial fit

4. Automatic dettection of Occultations

5. Complete occultation light curve model and fitting

   5.1 Automatic mode

   5.2 Fitting only the immersion

   5.3 Fitting the times and the opacity

6. Viewing and saving the results

## SORA package
from sora import LightCurve

## Other main and necessary packages
from astropy.time import Time
import astropy.units as u
import numpy as np
import matplotlib.pylab as pl
import os

SORA version: 0.3

1. Instantiating a LightCurve Object

The LightCurve Class can be instantiated in different ways.

The LightCurve main goal is to reduce and analyse stellar occultations light curves, which is a file composed of two or three columns: time, flux, and flux uncertainty. An ASCII file can be used as input, where the first column is the time (Julian Date or seconds relative to a reference time), the second column is the flux for each time, and the third is the flux uncertainty (optional). Another option is to furnish theses parameters directly as NumPy arrays. If the time provided is in JD the reference time will be set as 00:00:00.000 of the occultation date, otherwise the user should add the user reference time.

Other observational parameters needed are the exposure time for the respective light curve (in seconds) and the filter parameters (central wavelength and bandwidth, in microns). If the filter parameter is not furnished the default values are used: \(\lambda_0\) as 0.7 and \(\Delta \lambda\) as 0.3 microns.

If the parameters of the light curve were already calculated using SORA or other method, the user could instantiate the LightCurve by only providing the instants of immersion and emersion and their uncertainties, with no need to provide a light curve file. This is useful for example in cases where the user does not want to redo the fitting process and go directly to the chords and ellipse fit procedure within the Occultation Class.

The user can also instantiate a LightCurve Object for the negative observations. In this case, it is just needed the initial and the end time of the observational sequence. In this case, an Error (and a warning message) will be raised for all functions in this module that require times and fluxes, but the object will be instantiated and can be used for other purposes.

LightCurve?

Init signature: LightCurve(name='', **kwargs)
Docstring:
Defines a Light Curve.

Parameters
----------
name : `str`
    The name of the LightCurve. Each time an LightCurve object is defined
    the name must be different.

tref : `astropy.time.Time`, `str`, `float`
    Instant of reference.

    Format: `Julian Date`, string in ISO format or Time object.
    Required only if LightCurve have input fluxes and given time is
    not in Julian Date.

central_bandpass : `int`, `float`, otpional, default=0.7
    The center band pass of the detector used in observation. Value in microns.

delta_bandpass : `int`, `float`, optional, default=0.3
    The band pass width of the detector used in observation. Value in microns.

exptime : `int`, `float`
    The exposure time of the observation, in seconds.
    *NOT* required in cases *2*, *3* and *4* below.
    *Required* in case *1* below.

**kwargs: `int`, `float`
    Object velocity, distance, and star diameter.

    Note
    ----
    vel : `int`, `float`
        Velocity in km/s.

    dist : `int`, `float`
        Object distance in AU.

    d_star : `float`
        Star diameter, in km.


Warning
-------
Input data must be one of the 4 options below:

1) Input data from file with time and flux
    `file (str)`: a file with the time and flux. A third column with the error in
    flux can also be given.

    `usecols (int, tuple, array)`: Which columns to read, with the first being the
    time, the seconds the flux and third the flux error (optional).

    **Example:**

    >>> LightCurve(name, file, exptime) # dflux can also be given


2) Input data when file is not given:
    `time`: time must be a list of times, in seconds from tref, or Julian Date, or
    a Time object.

    `flux`: flux must be a list of fluxes. It must have the same lenght as time.

    `dflux`: if file not given, dflux must be a list of fluxes errors. It must
    have the same lenght as time. (not required)

    **Example:**

    >>> LightCurve(name, flux, time, exptime) # dflux can also be given


Cases for when `time` and `flux` are not given.


3) Input for a positive occultation:
    `immersion`: The instant of immersion.

    `emersion`: The instant of emersion.

    `immersion_err`: Immersion time uncertainty, in seconds.

    `emersion_err`: Emersion time uncertainty, in seconds.

    **Example:**

    >>> LightCurve(name, immersion, immersion_err, emersion, emersion_err)


4) Input for a negative occultation:
    `initial_time`: The initial time of observation.

    `end_time`: The end time of observation.

    **Example:**

    >>> LightCurve(name, initial_time, end_time)
File:           ~/Documentos/códigos/SORA/sora/lightcurve/core.py
Type:           type
Subclasses:

Example of LightCurve intantiate with an ASCII file

# The file lc_example.dat contains time in JD and flux
lc_1a = LightCurve(name='Example 1a', file='input/lightcurves/lc_example.dat', exptime=0.100)

# The file lc_example_trend.dat contains time in seconds and flux, so it is needed a reference time
lc_1b = LightCurve(name='Example 1b', file='input/lightcurves/lc_example_trend.dat',
                   tref='2017-06-22 00:00', exptime=0.100)

# If the file has more columns, the user can set which columns to use. Always in a sense, time, flux, flux error [optional].
lc_1c = LightCurve(name='Example 1c', file='input/lightcurves/lc_example_colunms.dat',
                   usecols=[2,4], exptime=0.100)

Example of LightCurve intantiate with NumPy arrays with times as Julian Dates

time, flux = np.loadtxt('input/lightcurves/lc_example_2.dat', unpack=True)

print('first time: ', time[0])

lc_2a = LightCurve(name='Example 2a', time=time, flux=flux, exptime=0.075)

first time:  2457927.38872249

Example of LightCurve instantiate with NumPy arrays with times as seconds relative to a reference time

time, flux = np.loadtxt('input/lightcurves/lc_example_2.dat',unpack=True)

tref = Time('2017-06-22 21:00').jd

time_sec = (time - tref)*u.d.to(u.s)

print('first time: ',time_sec[0])

lc_2b = LightCurve(name='Example 2b',time=time_sec,flux=flux,tref='2017-06-22 21:00',exptime=0.075)

first time:  1185.6231406331062

Example of LightCurve intantiate with the required parameters (positive)

lc_3 = LightCurve(name='Example 3',
                  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.010,
                  emersion ='2017-06-22 21:21:19.988',emersion_err=0.040)

Example of LightCurve intantiate with the required parameters (negative)

lc_4 = LightCurve(name='Example 4',
                  initial_time='2017-06-22 21:16:00.094',
                  end_time ='2017-06-22 21:28:00.018')

After the LightCurve was instantiated, some parameters can be easily retrieved

print('Name: {}'.format(lc_1a.name))
print('Initial time: {} UTC'.format(lc_1a.initial_time.iso))
print('Mean time:    {} UTC'.format(lc_1a.time_mean.iso))
print('End time:     {} UTC'.format(lc_1a.end_time.iso))

Name: Example 1a
Initial time: 2017-06-22 21:20:00.056 UTC
Mean time:    2017-06-22 21:21:40.007 UTC
End time:     2017-06-22 21:23:19.958 UTC

For cases where the immersion and emersion time were instanciated (or calculated), the mean time is between these values.

print('Initial time:   {} UTC'.format(lc_3.initial_time.iso))
print('Immersion time: {} UTC'.format(lc_3.immersion.iso))
print('Mean time:      {} UTC'.format(lc_3.time_mean.iso))
print('Emersion time:  {} UTC'.format(lc_3.emersion.iso))
print('End time:       {} UTC'.format(lc_3.end_time.iso))

Initial time:   2017-06-22 21:16:00.094 UTC
Immersion time: 2017-06-22 21:21:15.628 UTC
Mean time:      2017-06-22 21:21:17.808 UTC
Emersion time:  2017-06-22 21:21:19.988 UTC
End time:       2017-06-22 21:28:00.018 UTC

If the LightCurve were instantiated with times and fluxes, the exposure time and the mean cycle can also be obtained

print('Exposure time: {:.4f} s'.format(lc_1a.exptime))
print('Cycle time:    {:.4f} s'.format(lc_1a.cycle))

Exposure time: 0.1000 s
Cycle time:    0.1002 s

2. Set and/or modifing model parameters

In order to create a stellar occultation light curve model, some physical parameters are needed. First, we need to ‘project’ the time axis of the light curve in the sky plane using the event’s velocity (\(vel\) in km/s). Then the Fresnel diffraction should be applied. That diffraction works in a scale that depends on the object distance (\(dist\) in AU) and the observational wavelength (\(\lambda_0\) and \(\Delta \lambda\), in microns). The last parameter needed is the occulted star apparent diameter at the object distance (\(d\_star\) in km). In all cases, the user can provide each parameter with its own unit, otherwise SORA will consider the standard units described above.

Let’s start with the shadow velocity during the occultation, in km/s

lc_1a.set_vel(vel=22.0)
print('Vel: {:3.1f} km/s'.format(lc_1a.vel))

lc_1b.set_vel(vel=22.0*u.km/u.s)
print('Vel: {:3.1f} km/s'.format(lc_1b.vel))

lc_2a.set_vel(vel=22000.0*u.m/u.s)
print('Vel: {:3.1f} km/s'.format(lc_2a.vel))

Vel: 22.0 km/s
Vel: 22.0 km/s
Vel: 22.0 km/s

Now, the object distance at occultation time, in AU

lc_1a.set_dist(dist=15)
print('Distance: {:3.1f} AU'.format(lc_1a.dist))

lc_1b.set_dist(dist=15*u.AU)
print('Distance: {:3.1f} AU'.format(lc_1b.dist))

lc_2a.set_dist(dist=2243968060.5*u.km)
print('Distance: {:3.1f} AU'.format(lc_2a.dist))

Distance: 15.0 AU
Distance: 15.0 AU
Distance: 15.0 AU

And then, the observational wavelength, central value (:math:`lambda_0`) and its bandwidth (:math:`Deltalambda`), in microns

Default value set as \(\lambda_0\) = 0.7 and \(\Delta\lambda\) = 0.3 microns, no filter (Clear).

print('Observational wavelength centred at {:1.3f} with a bandwidth of {:1.3f} microns'
      .format(lc_1a.central_bandpass,lc_1a.delta_bandpass))

lc_1b.set_filter(central_bandpass=0.8, delta_bandpass=0.2)
print('Observational wavelength centred at {:1.3f} with a bandwidth of {:1.3f} microns'
      .format(lc_1b.central_bandpass,lc_1b.delta_bandpass))

lc_2a.set_filter(central_bandpass=0.8*u.micrometer, delta_bandpass=0.2*u.micrometer)
print('Observational wavelength centred at {:1.3f} with a bandwidth of {:1.3f} microns'
      .format(lc_2a.central_bandpass,lc_2a.delta_bandpass))

lc_4.set_filter(central_bandpass=800*u.nm, delta_bandpass=200*u.nm)
print('Observational wavelength centred at {:1.3f} with a bandwidth of {:1.3f} microns'
      .format(lc_4.central_bandpass,lc_4.delta_bandpass))

Observational wavelength centred at 0.700 with a bandwidth of 0.300 microns
Observational wavelength centred at 0.800 with a bandwidth of 0.200 microns
Observational wavelength centred at 0.800 with a bandwidth of 0.200 microns
Observational wavelength centred at 0.800 with a bandwidth of 0.200 microns

With the model parameters, the mean Fresnel Scale (:math:`F_S`) is automatically calculated.

\(F_S = \frac{1}{2}.\left(\sqrt{\frac{(\lambda + 0.5.\Delta\lambda).dist}{2}} + \sqrt{\frac{(\lambda - 0.5.\Delta\lambda).dist}{2}}\right)\)

print('Fresnel Scale: {:1.3f} km'.format(lc_1a.fresnel_scale))

Fresnel Scale: 0.881 km

Last, let’s set the stellar diameter projected at the body distance, in km

lc_1a.set_star_diam(d_star=0.2)
print('Stellar diameter: {:1.2f} km'.format(lc_1a.d_star))

lc_1b.set_star_diam(d_star=0.2*u.km)
print('Stellar diameter: {:1.2f} km'.format(lc_1b.d_star))

lc_2a.set_star_diam(d_star=200*u.m)
print('Stellar diameter: {:1.2f} km'.format(lc_2a.d_star))

Stellar diameter: 0.20 km
Stellar diameter: 0.20 km
Stellar diameter: 0.20 km

If the user instantiated the LightCurve object without times and flux and want to do it later (or redo it)

lc_4.set_flux(file='input/lightcurves/lc_example_3.dat',exptime=0.30, tref='2018-05-20 00:00:00.000')

3. Preparing the LightCurve

With the LightCurve object instantiated, the user may want to improve the quality of his curve using some processes.

For SORA v0.2 the only process implemented is the normalization of a light curve using a polinomial fit.

Other processes will be developed for future use, such as the binning of the light curve (sum, mean and median), a Savitzky–Golay filter, among others.

First, let’s check the expected magnitude drop and bottom flux

lc_1b.calc_magnitude_drop(mag_obj=18.6,mag_star=14.3)

print('Expected magnitude drop: {:.3f}'.format(lc_1b.mag_drop))
print('Expected bottom flux:    {:.3f}'.format(lc_1b.bottom_flux))

Expected magnitude drop: 4.320
Expected bottom flux:    0.019

Now, let’s look at this light curve

lc_1b.plot_lc()

There is a linear trend at this light curve.

#### 3.1. Light curve normalization using a polinomial fit

The user can remove this linear trend using the function LightCurve.normalize(). That function can be used automatically or with the desired inputed parameters.

lc_1b.normalize?

Signature:
lc_1b.normalize(
    poly_deg=None,
    mask=None,
    flux_min=0.0,
    flux_max=1.0,
    plot=False,
)
Docstring:
Returns the normalized flux within the flux min and flux max defined scale.

Parameters
----------
poly_deg : `int`
    Degree of the polynomial to be fitted.

mask : `bool` array
    Which values to be fitted.

flux_min : `int`, `float`
    Event flux to be set as 0.

flux_max : `int`, `float`
    Baseline flux to be set as 1.

plot : `bool`
    If True plot the steps for visual aid.
File:      ~/Documentos/códigos/SORA/sora/lightcurve/core.py
Type:      method

Let’s try the automatic mode

lc_1b.normalize()
lc_1b.plot_lc()

Normalization using a 1 degree polynomial
There is no improvement with a 2 degree polynomial

If the user wants to see the steps inside the automatic mode, just set the plot parameter as True

lc_1b.normalize(plot=True)
lc_1b.plot_lc()

Normalization using a 1 degree polynomial
There is no improvement with a 2 degree polynomial

As an standard procedure within the SORA library the user can control all the parameters

This is useful when the automatic mode does not work correctly (particular cases), or for better control of the procedure.

mask = (lc_1b.time < 76875) + (lc_1b.time > 76900)
lc_1b.normalize(poly_deg=5,mask=mask)

lc_1b.plot_lc()
pl.axvline(76875)
pl.axvline(76900)
pl.show()

Also, the user can control the maximum and/or minimum value of the flux

lc_4.plot_lc()
pl.ylim(-0.5,1.5)
pl.show()

In this case the bottom flux goes to 0.5, instead of the expected 0.0. This can be solved setting the flux_min parameter

lc_4.normalize(flux_min=0.50)

lc_4.plot_lc()
pl.ylim(-0.5,1.5)
pl.show()

Normalization using a 0 degree polynomial
There is no improvement with a 1 degree polynomial

If the user is not satisfied with the normalization, it is possible to throw away and redo the normalization process and to return to the initial stage using the LightCurve.reset_flux()

lc_4.reset_flux()

lc_4.plot_lc()
pl.ylim(-0.5,1.5)
pl.show()

4. Automatic detection of Occultations

Before considering any effect (Fresnell scale, star diameter, etc), the user can choose that SORA detect the events in the light curve automatically. This process is done using the BLS algorithm (Kovacs et al., 2002). This function will always find events and give their Signal to Noise Raito (SNR). The user will consider the values for the events found and determine the significance of each event. This is a useful tool for detection of low depth events or to check for secondary events, for example.

The result of the LightCurve.occ_detect() will be a python dictionary containing: a rank for the detection; duration; central, immersion and emersion times; time uncertainty; flux depth and its uncertainty; baseline flux and its uncertainty; SNR; and a boollean mask with the values inside the detection. All the times are in seconds. The automatic mode will determine only one detection (the most significant).

lc_2a.occ_detect?

Signature:
lc_2a.occ_detect(
    maximum_duration=None,
    dur_step=None,
    snr_limit=None,
    n_detections=None,
    plot=False,
)
Docstring:
Detects automatically the occultation event in the light curve.

Detects a 'square well' shaped transit. All parameters are optional.

Parameters
----------
maximum_duration : `float`, default: light curve time span
    Maximum duration of the occultation event.

dur_step : `float`, default: 1/2 of sampling rate
    Step size to sweep occultation duration event.

snr_limit : `float`, default=None
    Minimum occultation SNR.

n_detections : `int`, default=1
    Number of detections regardless the SNR. `n_detections` is
    superseded by `snr_limit`.

plot : `bool`
    True if output plots are desired.


Returns
-------
OrderedDict : `dict`
    An ordered dictionary of :attr:`name`::attr:`value` pairs for each
    parameter.

Examples
--------
>>> lc = LightCurve(time=time, flux=flux, exptime=0.0, name='lc_example')
>>> params = lc.occ_detect()
>>> params
{'rank': 1,
'occultation_duration': 40.1384063065052,
'central_time': 7916.773870512843,
'immersion_time': 7896.7046673595905,
'emersion_time': 7936.843073666096,
'time_err': 0.05011036992073059,
'depth': 0.8663887801707082,
'depth_err': 0.10986223384336465,
'baseline': 0.9110181732552853,
'baseline_err': 0.19045768512595365,
'snr': 7.886138392251848,
'occ_mask': array([False, False, False, ..., False, False, False])}
File:      ~/Documentos/códigos/SORA/sora/lightcurve/core.py
Type:      method

For visual aid the user can set the plot parameter as True

preliminary = lc_2a.occ_detect(plot=True)
preliminary

{'rank': 1,
 'occultation_duration': 11.52958944439888,
 'central_time': 76888.023570925,
 'immersion_time': 76882.2587762028,
 'emersion_time': 76893.7883656472,
 'time_err': 0.04017278552055359,
 'depth': 0.9797601789722584,
 'depth_err': 0.06536566413075075,
 'baseline': 1.0271920160757066,
 'baseline_err': 0.11532581974159437,
 'snr': 14.988911869883966,
 'occ_mask': array([False, False, False, ..., False, False, False])}

The user can choose to detect more than one event, two parameters can be used to control this n_detections and snr_limit.

With the n_detections, the user will set the number of detections they want. On the other hand snr_limit will obtain all the detection within the set limit. In all cases is up to the user to confirm if these detections have physical meaning.

preliminary = lc_2a.occ_detect(n_detections=5, plot=True)

preliminary = lc_2a.occ_detect(snr_limit=5, plot=True)

Other two parameters can be used for otimization purpose

The first is the maximum_duration, and this will set an upper limit to the size of the occultation (in seconds). The user can also set the dur_step, which is the step size in seconds of the occultation search. Those parameters will decrease the CPU time of these searches.

preliminary = lc_2a.occ_detect(maximum_duration=12,dur_step=1, plot=True)
preliminary

{'rank': 1,
 'occultation_duration': 11.0,
 'central_time': 76888.3231406331,
 'immersion_time': 76882.8231406331,
 'emersion_time': 76893.8231406331,
 'time_err': 0.04017278552055359,
 'depth': 0.9802503938156636,
 'depth_err': 0.062074094327636595,
 'baseline': 1.0241496004823303,
 'baseline_err': 0.1266541830883394,
 'snr': 15.791618136895426,
 'occ_mask': array([False, False, False, ..., False, False, False])}

5. Complete occultation light curve model and fitting

The model considers a sharp-edge occultation model (geometric) convolved with Fresnel diffraction, stellar diameter (projected at the body distance), CCD bandwidth, and finite integration time. It is created using a smaller time resolution than the observational light curve. By default, the model will have a resolution that is 1/10 of the mean Fresnel Scale (in km) divided by the event’s velocity (in km/s), if the exposure time is shorter than this value the software will use 1/10 of the exposure time. The parameters of interest are, mostly, the immersion and emersion times. Moreover, the opacity of the object can also be of interest for particular cases (occultations by rings, for instance), where 1.0 means an opaque body and 0.0 a transparent one.

The fit consists into create models and compare them with the observational light curve, through a standard \(\chi^2\) test. The parameters that minimizes the \(\chi^2\) are the fitted parameters and their uncertainties are obtained from the values within \(1\sigma = \chi^2_{min} + 1\) and \(3\sigma = \chi^2_{min} + 9\). The result of the fit is a ChiSquare object, and its functions can be found at its specific Jupyter-Notebook.

lc_1a.occ_model?

Signature:
lc_1a.occ_model(
    immersion_time,
    emersion_time,
    opacity,
    mask,
    npt_star=12,
    time_resolution_factor=10,
    flux_min=0,
    flux_max=1,
)
Docstring:
Returns the modelled light curve.

The modelled light curve takes into account the fresnel diffraction, the
star diameter and the instrumental response.

Parameters
----------
immersion_time : `int`, `float`
    Immersion time, in seconds.

emersion_time : `int`, `float`
    Emersion time, in seconds.

opacity : `int`, `float`
    Opacity. Opaque = 1.0, transparent = 0.0,

mask : `bool` array
    Mask with True values to be computed.

npt_star : `int`, default=12
    Number of subdivisions for computing the star size effects.

time_resolution_factor : `int`, `float`, default: 10*fresnel scale
    Steps for fresnel scale used for modelling the light curve.

flux_min : `int`, `float`, default=0
    Bottom flux (only object).

flux_max : `int`, `float`, default=1
    Base flux (object plus star).
File:      ~/Documentos/códigos/SORA/sora/lightcurve/core.py
Type:      method

Let’s find some parameters to create the model

immersion_time = 76880.000  # seconds relative to tref
emersion_time  = 76891.000  # seconds relative to tref
opacity = 1.0            # 1.0 == opaque; 0.0 == transparent

tmin = immersion_time  - 20 # seconds relative to tref
tmax = emersion_time   + 20 # seconds relative to tref

mask = (lc_1a.time > tmin) & (lc_1a.time < tmax)

Let’s give a look at the parameters we created (visual aid)

pl.plot(lc_1a.time,lc_1a.flux,'k.-',zorder=1,label='Observation')
pl.plot(lc_1a.time[mask],lc_1a.flux[mask],'bo',zorder=0,label='Points inside the mask')
pl.axvline(immersion_time,color='r',linestyle='-',label='Immersion')
pl.axvline(emersion_time,color='r',linestyle='--',label='Emersion')
pl.legend()
pl.show()

Now, let’s create the model

lc_1a.occ_model(immersion_time, emersion_time, opacity, mask)

lc_1a.plot_lc()
pl.xlim(76870,76900)
pl.show()

Let’s start with the fitting procedure

lc_1a.occ_lcfit?

Signature: lc_1a.occ_lcfit(**kwargs)
Docstring:
Monte Carlo chi square fit for occultations lightcurve.

Parameters
----------
tmin : `int`, `float`
    Minimum time to consider in the fit procedure, in seconds.

tmax : `int`, `float`
    Maximum time to consider in the fit procedure, in seconds.

flux_min : `int`, `float`, default=0
    Bottom flux (only object).

flux_max :`int`, `float`, default=1
    Base flux (object plus star).

immersion_time : `int`, `float`
    Initial guess for immersion time, in seconds.

emersion_time : `int`, `float`
    Initial guess for emersion time, in seconds.

opacity : `int`, `float`, default=1
    Initial guess for opacity. Opaque = 1, Transparent = 0.

delta_t : `int`, `float`
    Interval to fit immersion or emersion time.

dopacity : `int`, `float`, default=0
    Interval to fit opacity.

sigma : `int`, `float`, `array`, 'auto'
    Fluxes errors. If None it will use the `self.dflux`. If 'auto' it
    will calculate using the region outside the event.

loop : `int`, default=10000
    Number of tests to be done.

sigma_result : `int`, `float`
    Sigma value to be considered as result.

method : `str`, default=`chisqr`
    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
-------
chi2 : `sora.extra.ChiSquare`
    ChiSquare object.
File:      ~/Documentos/códigos/SORA/sora/lightcurve/core.py
Type:      method

5.1. Automatic mode

This step is one of the most time consuming processes within the SORA package. It is expected that it would take between 1-10 minutes to complete 10000 loops. To spend less computational time with this tutorial, we will use only 1000 loops. However it is suggested that the user always set at least 10000 loops.

lc_chi2 = lc_1a.occ_lcfit(loop=1000)

LightCurve fit: |████████████████████████████████████████|  - 100%

The result of the fitting is a ``ChiSquare`` Object

That is an auxiliary class, and its main information can be found at its specific Jupyter-Notebook, for the purpose of this tutorial only two main functions will be cover: print(ChiSquare) and ChiSquare.plot_chi2().

print(lc_chi2)

Minimum chi-square: 473.974
Number of fitted points: 496
Number of fitted parameters: 2
Minimum chi-square per degree of freedom: 0.959

immersion:
    1-sigma: 76880.326 +/- 0.029
    3-sigma: 76880.353 +/- 0.136

emersion:
    1-sigma: 76890.347 +/- 0.028
    3-sigma: 76890.354 +/- 0.102

lc_chi2.plot_chi2()

These gaps within the chisquare curve are expected as we did a small number of loops (1000)

5.2. Fitting only the immersion

The user can set the parameters of choice to the fitting. That will allow using a more complex set of parameters for the model or do specific fits. One typical example is fitting one side of the chord (only the immersion time, for instance).

immersion_time = 76880.338
tmin = 76873.000
tmax = 76887.000

mask = (lc_1b.time > tmin) & (lc_1b.time < tmax)

Let’s give a look at the parameters we created (visual aid)

pl.plot(lc_1b.time,lc_1b.flux,'k.-',zorder=1,label='Observation')
pl.plot(lc_1b.time[mask],lc_1b.flux[mask],'bo',zorder=0,label='Points inside the mask')
pl.axvline(immersion_time,color='r',linestyle='-',label='Immersion')
pl.legend()
pl.show()

lc_chi2_1b_imm = lc_1b.occ_lcfit(tmin=tmin, tmax=tmax, immersion_time=immersion_time,
                              delta_t=0.5, loop=1000)

LightCurve fit: |████████████████████████████████████████|  - 100%

print(lc_chi2_1b_imm)

Minimum chi-square: 136.864
Number of fitted points: 140
Number of fitted parameters: 1
Minimum chi-square per degree of freedom: 0.985

immersion:
    1-sigma: 76880.323 +/- 0.031
    3-sigma: 76880.347 +/- 0.124

lc_chi2_1b_imm.plot_chi2()

5.3. Fitting the times and the opacity

Let’s give a look in another LightCurve.

lc_2a.plot_lc()

In this light curve, aside the occultation by the main body there are two secondary events (rings), now let’s determine the times and the opacity for the first detection.

lc_chi2_2a_ring = lc_2a.occ_lcfit(tmin=76872.676, tmax=76874.892, immersion_time=76873.676,
                                  emersion_time=76873.892, opacity=0.50, delta_t=0.2,
                                  dopacity=0.5, loop=10000)

LightCurve fit: |████████████████████████████████████████|  - 100%

print(lc_chi2_2a_ring)

Minimum chi-square: 25.746
Number of fitted points: 27
Number of fitted parameters: 3
Minimum chi-square per degree of freedom: 1.073

immersion:
    1-sigma: 76873.676 +/- 0.004
    3-sigma: 76873.670 +/- 0.029

emersion:
    1-sigma: 76873.896 +/- 0.009
    3-sigma: 76873.898 +/- 0.032

opacity:
    1-sigma: 0.488 +/- 0.052
    3-sigma: 0.544 +/- 0.225

lc_chi2_2a_ring.plot_chi2()

lc_2a.plot_lc()
pl.xlim(76869,76878)
pl.legend()
pl.show()

6. Viewing and saving the results

At the end, the immersion and emersion times can be easily accessed using LightCurve.immersion and LightCurve.emersion. Their uncertainties can be found at LightCurve.immersion_err and LightCurve.emersion_err. Plot functions were created to facilitate the user to created the desired plots, an automatic plot can be created using Light_curve.plot_lc() and LightCurve.plot_model(), and the user can control the plot parameter using matplotlib functions. These plots can be saved using matplotlib.pyplot.savefig().

print('Imersion time {} UTC'.format(lc_1a.immersion))
print('Imersion time err: {:.3f} seconds'.format(lc_1a.immersion_err))
print('\n')
print('Emersion time {} UTC'.format(lc_1a.emersion))
print('Emersion time err: {:.3f} seconds'.format(lc_1a.emersion_err))

Imersion time 2017-06-22 21:21:20.326 UTC
Imersion time err: 0.029 seconds


Emersion time 2017-06-22 21:21:30.347 UTC
Emersion time err: 0.028 seconds

A complete log can be created by printing the ``LightCurve``

print(lc_1a)

Light curve name: Example 1a
Initial time: 2017-06-22 21:20:00.056 UTC
End time:     2017-06-22 21:23:19.958 UTC
Duration:     3.332 minutes
Time offset:  0.000 seconds

Exposure time:    0.1000 seconds
Cycle time:       0.1002 seconds
Num. data points: 2000

Bandpass:             0.700 +/- 0.300 microns
Object Distance:      15.00 AU
Used shadow velocity: 22.000 km/s
Fresnel scale:        0.040 seconds or 0.88 km
Stellar size effect:  0.009 seconds or 0.20 km
Inst. response:       0.100 seconds or 2.20 km
Dead time effect:     0.000 seconds or 0.00 km
Model resolution:     0.004 seconds or 0.09 km
Modelled baseflux:    1.029
Modelled bottomflux:  0.109
Light curve sigma:    0.336

Immersion time: 2017-06-22 21:21:20.326 UTC +/- 0.029 seconds
Emersion time:  2017-06-22 21:21:30.347 UTC +/- 0.028 seconds

Monte Carlo chi square fit.

Minimum chi-square: 473.974
Number of fitted points: 496
Number of fitted parameters: 2
Minimum chi-square per degree of freedom: 0.959

immersion:
    1-sigma: 76880.326 +/- 0.029
    3-sigma: 76880.353 +/- 0.136

emersion:
    1-sigma: 76890.347 +/- 0.028
    3-sigma: 76890.354 +/- 0.102

And this log can be saved as a file using LightCurve.to_log()

lc_1a.to_log('output/LC_test_1a.log')

!ls output/*.log

output/LC_test_1a.log

The post fitted light curve can be plotted using the LightCurve.plot_lc().

lc_1a.plot_lc()

The post fitted light curve complete model can be plotted using the LightCurve.plot_model()

lc_1a.plot_lc()
lc_1a.plot_model()

pl.xlim(76879.8,76880.8)
pl.ylim(-0.5,1.5)

pl.show()

The LightCurve data values can be saved using LightCurve.to_file().

If the model was correctly fitted with SORA, 4 files should be created: (i) One containing the observational light curve and its convoluted model; (ii) Containing the complete light curve model (geometric, affected by fresnel diffraction and by the star apparent diameter); and (iii and iv) with the labels for each colunm for other two files. Otherwise only 2 files will be created: the observational one and its labels.

lc_1a.to_file?

Signature: lc_1a.to_file(namefile=None)
Docstring:
Saves the light curve to a file.

Parameters
----------
namefile : `str`
    Filename to save the data.
File:      ~/Documentos/códigos/SORA/sora/lightcurve/core.py
Type:      method

lc_1a.to_file(namefile='output/LC_test_1a.dat')

!ls output/*.dat
print('\n')
!ls output/*.dat.label

output/example_chi2.dat  output/LC_test_1a.dat  output/model_LC_test_1a.dat


output/example_chi2.dat.label  output/model_LC_test_1a.dat.label
output/LC_test_1a.dat.label

One last thing, a common issue within the stellar occultation light curves is a time offset for a specific light curve

To fit the body size and shape the user should consider the times (immersion and emersion) for all the observers. It is not unusual to have the different chords projected at the sky plane not be appropriately ‘aligned’. Those shifts are often caused by instrumental issues (lack of GPS, software’s features, internet protocols, among many others). To deal with this issue, the user can add a time offset (in seconds) manually at the LightCurve Object using the LightCurve.dt. Then the SORA package will automatically consider this offset to the fitted (or instantiated) parameters.

lc_1a.dt = 0 #seconds

print('Imersion time {} UTC'.format(lc_1a.immersion))

lc_1a.dt = 30 #seconds

print('Imersion time {} UTC'.format(lc_1a.immersion))

Imersion time 2017-06-22 21:21:20.326 UTC
Imersion time 2017-06-22 21:21:50.326 UTC

This Jupyter-Notebook was designed as a tutorial for how to work with the LightCurve 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