Getting apparent magnitude of a star for a specific date

[gloglo17]

Hi. Can this library be used to find out what the apparent magnitude of a particular variable star was/will be at a particular time? If so, would you please post an example source code for Algol? Thank you very much!

[waqasbhatti]

Hello!

This can be done in a few steps. The prediction of the apparent magnitudes on future dates won't be as accurate as those from a dedicated eclipsing binary modeling code, but probably enough to plan observations.

First, you need a light curve for Algol. I used STScI's MAST service to get one from the TESS archive [a full-frame-image (FFI) high level science product (HLSP)]: hlsp_qlp_tess_ffi_s0018-0000000346783960_tess_v01_llc.txt

This is a CSV file that you can read in:

import pandas as pd
import numpy as np

lc = pd.read_csv('hlsp_qlp_tess_ffi_s0018-0000000346783960_tess_v01_llc.txt', usecols=[0, 1], names=['BTJD', 'flux'], comment='#')

You can then use some astrobase functions to find the epoch of the primary eclipse phase (time at eclipse minimum), given the period of the eclipsing binary (2.867328 days from Wikipedia):

from astrobase.lcfit import nonphysical
from astrobase import plotbase

# use a Savitsky-Golay smoothing 'fit' to find the epoch
fit_results = nonphysical.savgol_fit_magseries(np.array(lc['BTJD']), np.array(lc['flux']), np.full_like(np.array(lc['BTJD']), 0.001), 2.867328, plotfit='savgol-fit.png', magsarefluxes=True, sigclip=50.0)

epoch = fit_results['fitinfo']['fitepoch']

Then, you can use this to plot a phased flux time-series to estimate some parameters for an eclipsing binary model fit:

times, fluxes, errs = np.array(lc['BTJD']), np.array(lc['flux']), np.full_like(np.array(lc['BTJD']), 0.001)
plotbase.plot_phased_magseries(times, fluxes, 2.867328, magsarefluxes=True, sigclip=50.0, epoch=epoch, outfile='phased-lc.png')

That gives you a plot like the one below:

Next, you can fit a simple eclipsing binary model:

from astrobase.lcfit import eclipses

ebfit = eclipses.gaussianeb_fit_magseries(times, fluxes, errs, [2.867328, 1796.84939402, 0.4, 0.1, 0.2, 0.5], plotfit='ebfit.png', magsarefluxes=True, sigclip=50.0, param_bounds={'period': 'fixed'})

This makes a fit plot:

Using this model, you can predict the magnitudes of the star into the future. First, convert the units of time in the light curve from TESS to BJD (Baryocentric Julian Date):

ref_bjd = 2457000.0
times_bjd = ref_bjd + times

# get the eclipse fit params
fit_params = ebfit['fitinfo']['finalparams']

# convert the fit epoch into BJD as well
fit_params[1] = ref_bjd + fit_params[1]

Generate an array of BJDs into the future:

from astrobase import timeutils
jd_now = timeutils.jd_now()
jd_to_oneyear = np.linspace(jd_now - 3.0, jd_now + 365.0, num=365*24)

Use the eclipse model to predict fluxes on these dates:

from astrobase.lcmodels import eclipses as ebmodels

model_fluxes, model_phases, model_times, _, _ = ebmodels.invgauss_eclipses_func(fit_params, jd_to_oneyear, np.full_like(jd_to_oneyear, 1.0), np.full_like(jd_to_oneyear, 0.0001))

Convert model fluxes to TESS apparent magnitudes. The median TESS apparent magnitude of Algol is 3.3684 and the relation between the normalized flux in the TESS light curve and apparent magnitude is:

apparent_mag - object_tess_mag = -2.5 log(normalized_flux)

So we can do the following to get an array of apparent magnitudes:

model_mags = 3.3684 - 2.5*np.log10(model_fluxes)

To figure out the apparent brightness of Algol on specific dates, you can convert these to Julian dates and then use the model output to see the magnitudes on those dates. For example, let's get the apparent magnitude light curve for the next two weeks.

from datetime import datetime, timezone, timedelta

# get today's date and two weeks in the future
date_now = datetime.now(tz=timezone.utc)
date_twoweeks = date_now + timedelta(days=14)

# convert dates to Julian dates
jd_now = timeutils.datetime_to_jd(date_now)
jd_twoweeks = timeutils.datetime_to_jd(date_twoweeks)

# sort the model output by time first
sorted_time_idx = np.argsort(model_times)
sorted_model_times = model_times[sorted_time_idx]
sorted_model_mags = model_mags[sorted_time_idx]

# then apply the date filter
date_filtered_idx = np.where((sorted_model_times >= jd_now) & (sorted_model_times <= jd_twoweeks))
date_filtered_times = sorted_model_times[date_filtered_idx]
date_filtered_model_mags = sorted_model_mags[date_filtered_idx]

Now, you can plot the apparent magnitudes for the next two weeks:

import matplotlib.pyplot as plt
plt.figure()
plt.scatter(date_filtered_times, date_filtered_model_mags, marker='.')
plt.ylim(4, 3)
plt.savefig('algol-next-twoweeks.png')

[gloglo17]

Thank you so much for your wonderfully elaborated response. Unfortunately, I'm pretty much a layman when it comes to astronomy and don't really understand it. Actually I need a function to calculate the apparent magnitude of Algol for a given time because I want to include it as a parameter in my AI project. (Later I will need the same for other variable stars, but if I understand correctly, the principle is exactly the same, just a different dataset will be needed.)

I tried to make a function based on what you wrote, but I'm missing a step, so it doesn't work, could you please take a look?

import numpy as np
import juliandate

def get_algol_magnitude(date): # python datetime
    # read data from text file
    data = np.genfromtxt('hlsp_qlp_tess_ffi_s0018-0000000346783960_tess_v01_llc.txt', delimiter=',')

    # calculate Julian Date and BTJD offset
    jd = juliandate.from_gregorian(date.year,date.month,date.day,date.hour,date.minute,date.second)
    btjd_offset = 2454833.0 - jd

    # find index of closest time to datetime_obj
    index = np.absolute(data[0] - btjd_offset).argmin()

    # calculate Algol magnitude
    magnitude = -2.5 * math.log10(data[1][index] / data[2][index]) 

    return magnitude

[waqasbhatti]

Ah okay, I think the problems I can see are:

• btjd_offset is using the JD of the start date from the Kepler telescope instead of TESS. Since the light curve above is from TESS, you probably want to use the TESS start date of 2457000.0 in that calculation: btjd_offset = jd - 2457000.0
• For magnitude, you only need column index 1 from the light curve, because that is already normalized flux (flux / flux_median). You also need the apparent magnitude corresponding to the flux_median, which is the apparent magnitude of the star in the TESS band.
• The calculation of the magnitude would be: magnitude = tess_algol_mag - 2.5 * math.log10(data[1][index]) where tess_algol_mag is 3.3684. You can get TESS magnitudes for stars using the TESS input catalog.
• The method here will only work if your intended query date is in the range of the observation dates of TESS for this object. If not, you'll need a eclipsing binary model to predict magnitudes on arbitary dates.

[gloglo17]

Someone advised me how to calculate the phase of Algol for a given date. This should give a number between 0.0 and 1.0 where 0.0 means Algol is at a minimum, 0.5 means Algol is half way between minima and 1.0 means Algol is at a minimum again.

import juliandate
import datetime

def get_algol_phase(time):
     time_jd = juliandate.from_gregorian(date.year,date.month,date.day,date.hour,date.minute,date.second)

     epoch_minimum_jd=2460072.826389
     algol_period=2.867328

     days_from_epoch = abs(epoch minimum_jd - time_jd)
     algol_phase = (days_from_epoch / algol_period) % 1.0   

     return algol_phase

time = datetime.datetime.now()
algol_phase = get_algol_phase(time)

Is there a way how to use this together with your code so that I can get the magnitude for the given phase? :)

[waqasbhatti]

If you have series of times and associated magnitudes, you can get the phases and associated magnitudes at each phase using astrobase.lcmath.phase_magseries() function. You'll need the period and epoch.

from astrobase import lcmath

# times_jd and mags are from an existing light curve
phased_lc = lcmath.phase_magseries(times_jd, mags, period, epoch, wrap=False, sort=True)

# returned phased mags are sorted in order of phase
phases = phased_lc["phase"]
phased_mags = phased_lc["mags"]

Then, you can use numpy's array indexing or np.where() to find the mags at any phase value.