Inclusion of obliquity in the disk-integrated spectra analysis

[frescigno]

While an obliquity can be assigned to the star in the code, no obliquity is ever included in the computation of the stellar disk-integrated spectra. While this would not affect the disk-integration in most cases, it starts to become relevant if any non-axisymmetric mask is applied to the star.

The full rv formula for a position on a rotating sphere as seen by a static observer is defined as:

RV[m/s] = radius[m] 2pi/Period[s] * sin(inc) * sin(colatitude) * sin(longitude – obliquity – phase)

This is similar to the formula used in the code, although not fully.
Starting from within the function spectrum it is defined:

phi, theta = hp.pix2ang(star.sides, range(hp.nside2npix(star.sides)))

Following healpix definitions, this would make phi=colatitude, theta=longitude. This is in opposition to what the strings actually say, but it is not a problem because it is treated consistently within the computation.

Now moving deeper into the spectrum function, it called integrated_spectrum which itself calls doppler_shift, which calls radial_velocity. In this last function the rv is defined as:

rv = radius_m * omega * jnp.sin(theta - phase) * sin_phi
if inc is not None:
rv = rv * jnp.sin(inc)

This formula (assuming the theta=long, phi=colat) is similar to the full RV formula, but misses the obliquity component (obliquity itself is not even an allowed variable and it is not propagated through).
This could be an easy fix as:

if obl is not None:
phase += obl
rv = radius_m * omega * jnp.sin(theta - phase) * sin_phi
if inc is not None:
rv = rv * jnp.sin(inc)

This should not affect the final disk-integrated profile. However, when to test it I implement the change, it does create a shift in the overall profile. This shift is not dependent on phase, but is dependent on the obliquity value.

Why is it so? Could this be caused by how masked_projected_limb is structured, which defines the visible surface based on the distance from the equator?

I have done the changes required to include obliquity in the disk-integrated spectra (the small change in the radial_velocity function, as well changes to allow the propagation down the tree of obl as a variable).
Some results for different obliquities shown below following the example in the tutorials:

star = Star.from_sides(16, u=(0.5, 0.2), period=.02, inc=np.deg2rad(90))
spot = core.soft_spot(star.sides, 0.0, 0.0, 0.15)
wv = jnp.linspace(400, 410, 500) * 1e-9
star_spectrum = 1 - 0.5 * gaussian(wv, wv.mean(), 0.5e-9)
spectra = star_spectrum[:, None] * (star - spot).y
star = star.set(y=spectra, wv=wv)

Obliquity = 0 deg

Obliquity = 20 deg

Obliquity = 50 deg

Obliquity = -20 deg

Obliquity = -50 deg