Calculate Fourier parameters from shape data

In this tutorial we take a look at how to construct Harmonica’s Fourier parameterisation of a transmission string from shape data. This is useful for forward modelling transit light curves of various objects with known shapes. You only need to determine some high-order Fourier representation of the shape, and then you may easily simulate high-precision light curves.

Let us start by loading some shape data. Here is some data made earlier, based on the transit of Phobos as viewed from the Perseverance rover on Mars.

import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt


r_data = np.array([ 0.25299594, 0.25049962, 0.2555872 , 0.24131962, 0.23499399, 0.2281387 , 0.21770261, 0.21355986, 0.21756298, 0.22216417, 0.215748  , 0.21737008, 0.22181932, 0.22220379, 0.22215942, 0.2191899 , 0.22060524, 0.22316191, 0.22795516, 0.23203883, 0.23715028, 0.24817015, 0.26129245, 0.2701211 , 0.27631858, 0.279554  , 0.27550068, 0.26970037, 0.25844681, 0.2551263 , 0.2564253 , 0.25649795, 0.2619272 , 0.26137834, 0.24959362, 0.23925815, 0.2218478 , 0.21541933, 0.21006674, 0.1952882 , 0.18672611, 0.19014673, 0.20044651, 0.2134876 , 0.21243847, 0.21826262, 0.23046709, 0.24290151, 0.24455098, 0.24425059, 0.25006916])
theta_data = np.array([-3.05626382, -2.92458806, -2.81512117, -2.71398938, -2.56361121, -2.37500762, -2.19475882, -2.05340439, -1.89915229, -1.78311113, -1.68792003, -1.54576688, -1.45830389, -1.36707739, -1.26004262, -1.12498761, -1.03931859, -0.91754355, -0.8460734 , -0.72241586, -0.64825102, -0.51378753, -0.40700302, -0.30187693, -0.20512046, -0.1247507 , -0.0469304 , 0.04035871, 0.12437005, 0.17993957, 0.24397848, 0.35182409, 0.48800481, 0.57248726, 0.72834884, 0.84510295, 0.98657143, 1.12488908, 1.2518673 , 1.36829014, 1.56557413, 1.70233616, 1.84736078, 1.96211635, 2.0600249 , 2.31436316, 2.4535093 , 2.60749337, 2.78292412, 2.95552791, 3.08844006])

plt.figure(figsize=(10, 7))
plt.gca().set_aspect("equal", "datalim")
plt.scatter(r_data * np.cos(theta_data), r_data * np.sin(theta_data),
            color=cm.inferno(0.1), label="Shape data")
plt.plot(0.2324 * np.cos(theta_data), 0.2324 * np.sin(theta_data),
         color="#b9b9b9", ls="--", label="Reference circle")
plt.xlabel("x / stellar radii", fontsize=13)
plt.ylabel("y / stellar radii", fontsize=13)
plt.legend(loc="lower left", fontsize=12)
plt.show()

In polar coordinates:

plt.figure(figsize=(10, 7))
plt.scatter(theta_data * 180. / np.pi, r_data,
            color=cm.inferno(0.1), label="Shape data")
plt.plot(theta_data * 180. / np.pi, 0.2324 * np.ones(theta_data.shape[0]),
         color="#b9b9b9", ls="--", label="Reference circle")
plt.xlabel("$\\theta$ / degrees", fontsize=13)
plt.ylabel("$r_{\\rm{p}}$ / stellar radii", fontsize=13)
plt.legend(loc="lower left", fontsize=12)
plt.show()

The Fourier parameters of Harmonica’s transmission strings can be calculated from these data. This involves solving:

\[a_0 = \frac{1}{2 \pi} \int_{-\pi}^{\pi} r(\theta) \,d \theta,\]

\[a_n = \frac{1}{\pi} \int_{-\pi}^{\pi} r(\theta) \cos{(n \theta)} \,d \theta,\]

\[b_n = \frac{1}{\pi} \int_{-\pi}^{\pi} r(\theta) \sin{(n \theta)} \,d \theta,\]

where each integrand is a generated interpolation function based on the data. Here let us compute up to \(n = 15\).

from scipy import interpolate, integrate


n_harmonics = 15
transmission_string_params = []
for n in range(n_harmonics + 1):

    if n == 0:
        # a_0 term.
        integrand_cos0_func = interpolate.interp1d(
            theta_data, r_data,
            kind="cubic", bounds_error=False, fill_value="extrapolate")
        cn = integrate.quad(
            integrand_cos0_func, -np.pi, np.pi,
            epsabs=1.e-7, epsrel=1.e-7, limit=500)[0] / (2. * np.pi)
        transmission_string_params.append(cn)
    else:
        # a_n terms.
        integrand_cosn_func = interpolate.interp1d(
            theta_data, r_data * np.cos(n * theta_data),
            kind="cubic", bounds_error=False, fill_value="extrapolate")
        cn = integrate.quad(
            integrand_cosn_func, -np.pi, np.pi,
            epsabs=1.e-7, epsrel=1.e-7, limit=500)[0] / (1. * np.pi)
        transmission_string_params.append(cn)

        # b_n terms.
        integrand_sinn_func = interpolate.interp1d(
            theta_data, r_data * np.sin(n * theta_data),
            kind="cubic", bounds_error=False, fill_value="extrapolate")
        sn = integrate.quad(
            integrand_sinn_func, -np.pi, np.pi,
            epsabs=1.e-7, epsrel=1.e-7, limit=500)[0] / (1. * np.pi)
        transmission_string_params.append(sn)

transmission_string_params = np.array(transmission_string_params)

In practice the above method tends to produce numerical ringing when a large number of harmonics are computed. As an alternative, you may find using numpy’s discrete Fourier transform is a more stable approach.

theta_data_even_spacing = np.linspace(-np.pi, np.pi, 101, endpoint=True)
r_data_even_spacing = np.interp(
    theta_data_even_spacing, theta_data, r_data)  # Interpolate onto evenly-spaced grid.
r_data_even_spacing = np.concatenate(
    [r_data_even_spacing[theta_data_even_spacing >= 0.],
     r_data_even_spacing[theta_data_even_spacing < 0.]])  # Reorder for the interval 0 - 2*pi.

n_terms = 2 * n_harmonics + 1
n_data = len(r_data_even_spacing)

c = np.fft.fft(r_data_even_spacing) / n_data
r_coeffs = np.empty(n_data)
r_coeffs[0] = c[0].real
r_coeffs[1::2] = 2*c[1:n_data//2+1].real
r_coeffs[2::2] = -2*c[1:n_data//2+1].imag
transmission_string_params = r_coeffs[:n_terms]

Now plot the resulting transmission string and check you have a good approximation of your input shape:

from harmonica import HarmonicaTransit


ht = HarmonicaTransit()
ht.set_planet_transmission_string(transmission_string_params)
theta_model = np.linspace(-np.pi, np.pi, 1000, endpoint=False)

plt.figure(figsize=(10, 7))
plt.gca().set_aspect("equal", "datalim")
plt.scatter(r_data * np.cos(theta_data), r_data * np.sin(theta_data),
            color=cm.inferno(0.1), label="Shape data")
plt.plot(ht.get_planet_transmission_string(theta_model) * np.cos(theta_model),
         ht.get_planet_transmission_string(theta_model) * np.sin(theta_model),
         c=cm.Purples(0.8), lw=2.5, label="Transmission string")
plt.plot(0.2324 * np.cos(theta_data), 0.2324 * np.sin(theta_data),
         color="#b9b9b9", ls="--", label="Reference circle")
plt.xlabel("x / stellar radii", fontsize=13)
plt.ylabel("y / stellar radii", fontsize=13)
plt.legend(loc="lower left", fontsize=12)
plt.show()

In polar coordinates:

plt.figure(figsize=(10, 7))
plt.scatter(theta_data * 180. / np.pi, r_data,
            color=cm.inferno(0.1), label="Shape data")
plt.plot(theta_model * 180. / np.pi, ht.get_planet_transmission_string(theta_model),
         c=cm.Purples(0.8), lw=2.5, label="Transmission string")
plt.plot(theta_data * 180. / np.pi, 0.2324 * np.ones(theta_data.shape[0]),
         color="#b9b9b9", ls="--", label="Reference circle")
plt.xlabel("$\\theta$ / degrees", fontsize=13)
plt.ylabel("$r_{\\rm{p}}$ / stellar radii", fontsize=13)
plt.legend(loc="lower left", fontsize=12)
plt.show()

These parameters can now be used to generate transit light curves of your object. The only thing left to do is specify the orbit and stellar limb darkening.

Note

When using many Fourier coefficients you may need to increase the number of terms used to compute the light curve solutions, in order to generate a sufficiently precise result for your use case. This can easily be accomplished by manually setting pnl_c and pnl_e when you instantiate the HarmonicaTransit class. You can also check the estimated precision with the get_precision_estimate() method.

times = np.linspace(-0.2, 0.2, 500)
ht = HarmonicaTransit(times, pnl_c=100, pnl_e=100)
ht.set_orbit(t0=0., period=4., a=7., inc=88. * np.pi / 180.)
ht.set_stellar_limb_darkening(u=np.array([0.074, 0.193]), limb_dark_law="quadratic")
ht.set_planet_transmission_string(transmission_string_params)
light_curve = ht.get_transit_light_curve()

plt.figure(figsize=(10, 7))
plt.plot(times, light_curve, c="#000000")
plt.xlabel('Time / days')
plt.ylabel('Relative flux')
plt.show()