Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat(hfbcat): Move get_doppler_shifted_freq from ephemeris #73

Merged
merged 1 commit into from
Jul 21, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
90 changes: 9 additions & 81 deletions ch_util/ephemeris.py
Original file line number Diff line number Diff line change
Expand Up @@ -820,90 +820,18 @@ def peak_RA(body, date=None, deg=False):
return ra


def get_doppler_shifted_freq(
source: Union[skyfield.starlib.Star, str],
date: Union[float, list],
freq_rest: Union[float, list] = None,
obs: Observer = chime,
) -> np.array:
"""Calculate Doppler shifted frequency of spectral feature with rest
frequency `freq_rest`, seen towards source `source` at time `date`, due to
Earth's motion and rotation, following the relativistic Doppler effect.

Parameters
----------
source
Position(s) on the sky. If the input is a `str`, attempt to resolve this
from `ch_util.hfbcat.HFBCatalog`.
date
Unix time(s) for which to calculate Doppler shift.
freq_rest
Rest frequency(ies) in MHz. If None, attempt to obtain rest frequency
of absorption feature from `ch_util.hfbcat.HFBCatalog.freq_abs`.
obs
An Observer instance to use. If not supplied use `chime`. For many
calculations changing from this default will make little difference.

Returns
-------
freq_obs
Doppler shifted frequencies in MHz. Array where rows correspond to the
different input rest frequencies and columns correspond either to input
times or to input sky positions (whichever contains multiple entries).

Notes
-----
Only one of `source` and `date` can contain multiple entries.

Example
-------
To get the Doppler shifted frequencies of a feature with a rest frequency
of 600 MHz for two positions on the sky at a single point in time (Unix
time 1717809534 = 2024-06-07T21:18:54+00:00), run:

>>> from skyfield.starlib import Star
>>> from skyfield.positionlib import Angle
>>> from ch_util.tools import get_doppler_shifted_freq
>>> coord = Star(ra=Angle(degrees=[100, 110]), dec=Angle(degrees=[45, 50]))
>>> get_doppler_shifted_freq(coord, 1717809534, 600)
"""

from scipy.constants import c as speed_of_light

from ch_util.fluxcat import _ensure_list
from ch_util.hfbcat import HFBCatalog

# For source string inputs, get skyfield Star object from HFB catalog
if isinstance(source, str):
try:
source = HFBCatalog[source].skyfield
except KeyError:
raise ValueError(f"Could not find source {source} in HFB catalog.")

# Get rest frequency from HFB catalog
if freq_rest is None:
if not source.names or source.names not in HFBCatalog:
raise ValueError(
"Rest frequencies must be supplied unless source can be found "
"in ch_util.hfbcat.HFBCatalog. "
f"Got source {source} with names {source.names}"
)
else:
freq_rest = HFBCatalog[source.names].freq_abs

# Prepare rest frequencies for broadcasting
freq_rest = np.asarray(_ensure_list(freq_rest))[:, np.newaxis]
# For backwards compatibility
def get_doppler_shifted_freq(*args, **kwargs):
"""Deprecated. Use `ch_util.hfbcat.get_doppler_shifted_freq`."""
from . import hfbcat

# Get rate at which the distance between the observer and source changes
# (positive for observer and source moving appart)
range_rate = get_range_rate(source, date, obs)
import warnings

# Compute observed frequencies from rest frequencies
# using relativistic Doppler effect
beta = range_rate / speed_of_light
freq_obs = freq_rest * np.sqrt((1.0 - beta) / (1.0 + beta))
warnings.warn(
"Use `ch_util.hfbcat.get_doppler_shifted_freq` instead.", DeprecationWarning
)

return freq_obs
return hfbcat.get_doppler_shifted_freq(*args, **kwargs)


def get_range_rate(
Expand Down
92 changes: 91 additions & 1 deletion ch_util/hfbcat.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,16 @@
Catalog of HFB test targets
"""

from __future__ import annotations
import os
import numpy as np
from typing import TYPE_CHECKING, Union

from .fluxcat import FluxCatalog
from . import ephemeris
from .fluxcat import FluxCatalog, _ensure_list

if TYPE_CHECKING:
import skyfield.starlib.Star

# Define the source collection that should be loaded when this module is imported.
HFB_COLLECTION = os.path.join(
Expand Down Expand Up @@ -79,5 +86,88 @@ def __init__(
self.freq_abs = freq_abs


def get_doppler_shifted_freq(
source: Union[skyfield.starlib.Star, str],
date: Union[float, list],
freq_rest: Union[float, list] = None,
obs: ephemeris.Observer = ephemeris.chime,
) -> np.array:
"""Calculate Doppler shifted frequency of spectral feature with rest
frequency `freq_rest`, seen towards source `source` at time `date`, due to
Earth's motion and rotation, following the relativistic Doppler effect.

Parameters
----------
source
Position(s) on the sky. If the input is a `str`, attempt to resolve this
from `ch_util.hfbcat.HFBCatalog`.
date
Unix time(s) for which to calculate Doppler shift.
freq_rest
Rest frequency(ies) in MHz. If None, attempt to obtain rest frequency
of absorption feature from `ch_util.hfbcat.HFBCatalog.freq_abs`.
obs
An Observer instance to use. If not supplied use `chime`. For many
calculations changing from this default will make little difference.

Returns
-------
freq_obs
Doppler shifted frequencies in MHz. Array where rows correspond to the
different input rest frequencies and columns correspond either to input
times or to input sky positions (whichever contains multiple entries).

Notes
-----
Only one of `source` and `date` can contain multiple entries.

Example
-------
To get the Doppler shifted frequencies of a feature with a rest frequency
of 600 MHz for two positions on the sky at a single point in time (Unix
time 1717809534 = 2024-06-07T21:18:54+00:00), run:

>>> from skyfield.starlib import Star
>>> from skyfield.units import Angle
>>> from ch_util.hfbcat import get_doppler_shifted_freq
>>> coord = Star(ra=Angle(degrees=[100, 110]), dec=Angle(degrees=[45, 50]))
>>> get_doppler_shifted_freq(coord, 1717809534, 600)
"""

from scipy.constants import c as speed_of_light

# For source string inputs, get skyfield Star object from HFB catalog
if isinstance(source, str):
try:
source = HFBCatalog[source].skyfield
except KeyError:
raise ValueError(f"Could not find source {source} in HFB catalog.")

# Get rest frequency from HFB catalog
if freq_rest is None:
if not source.names or source.names not in HFBCatalog:
raise ValueError(
"Rest frequencies must be supplied unless source can be found "
"in ch_util.hfbcat.HFBCatalog. "
f"Got source {source} with names {source.names}"
)
else:
freq_rest = HFBCatalog[source.names].freq_abs

# Prepare rest frequencies for broadcasting
freq_rest = np.asarray(_ensure_list(freq_rest))[:, np.newaxis]

# Get rate at which the distance between the observer and source changes
# (positive for observer and source moving appart)
range_rate = ephemeris.get_range_rate(source, date, obs)

# Compute observed frequencies from rest frequencies
# using relativistic Doppler effect
beta = range_rate / speed_of_light
freq_obs = freq_rest * np.sqrt((1.0 - beta) / (1.0 + beta))

return freq_obs


# Load the HFB target list
HFBCatalog.load(HFB_COLLECTION)
Loading