Skip to content
Snippets Groups Projects
Commit 6701ee6c authored by Maaijke Mevius's avatar Maaijke Mevius
Browse files

recover accidentally removed file

parent f43e38c7
No related branches found
No related tags found
3 merge requests!12Rmcalc,!8Draft: Magnetic,!6simple ipp model
Pipeline #101454 passed
# Copyright (C) 2024 ASTRON (Netherlands Institute for Radio Astronomy)
# SPDX-License-Identifier: Apache-2.0
"""Module for getting the Ionospheric Piercepoints"""
from __future__ import annotations
from typing import NamedTuple
import astropy.units as u
import numpy as np
from astropy.coordinates import ITRS, AltAz, EarthLocation, SkyCoord
from astropy.time import Time
r_earth = 6364.62 * u.km
class IPP(NamedTuple):
"""ionospheric piercepoint"""
loc: EarthLocation
"""location of the piercepoints, dimension: times x altitudes. All altitudes are assumed to be equal"""
times: Time
"""array of times"""
los: SkyCoord
"""Line of sight direction in ITRS coordinates"""
airmass: np.ndarray
"""airmass factor to convert to slant values"""
def get_ipp_from_skycoord(
loc: EarthLocation, times: Time, source: SkyCoord, height_array: u.Quantity
) -> IPP:
"""Get ionospheric piercepoints
Parameters
----------
loc : EarthLocation
observer location
times : Time
observation times
source : SkyCoord
source location
height_array : u.Quantity
array of altitudes
Returns
-------
IPP
Ionospheric piercepoints
"""
# Note: at the moment we calculate ipp per station. I think this is ok,
# but maybe we need to include a many stations option
aa = AltAz(location=loc, obstime=times)
altaz = source.transform_to(aa)
return get_ipp_from_altaz(loc, altaz, height_array)
def get_ipp_from_altaz(
loc: EarthLocation, altaz: AltAz, height_array: u.Quantity
) -> IPP:
"""get ionospheric piercepoints from azimuth elevations
Parameters
----------
loc : EarthLocation
observer location
altaz : AltAz
azimuth and elevations for all times
height_array : u.Quantity
array of altitudes
Returns
-------
IPP
_description_
"""
los_dir = altaz.transform_to(ITRS)
ipp, airmass = _get_ipp_simple(height_array=height_array, loc=loc, los_dir=los_dir)
return IPP(
loc=EarthLocation.from_geocentric(*ipp),
times=altaz.obstime,
los=los_dir,
airmass=airmass,
)
def _get_ipp_simple(
height_array: u.Quantity, loc: EarthLocation, los_dir: SkyCoord
) -> list[u.Quantity]:
"""helper function to calculate ionospheric piercepoints using a simple spherical earth model
|loc + alphas * los_dir| = r_earth + height_array, solve for alphas using abc formula
Parameters
----------
height_array : u.Quantity
array of altitudes
loc : EarthLocation
observer location
los_dir : SkyCoord
line of sight (ITRS)
Returns
-------
list[u.Quantity]
ipp.x, ipp.y, ipp.z positions
"""
c_value = r_earth**2 - (r_earth + height_array) ** 2
ipp_vector = los_dir.cartesian[np.newaxis] # unit vector
b_value = (
loc.x * ipp_vector.x.value
+ loc.y * ipp_vector.y.value
+ loc.z * ipp_vector.z.value
)[np.newaxis] # inproduct, add axis for altitudes
alphas = -b_value + np.sqrt(b_value**2 - c_value[:, np.newaxis])
ipp = [
loc.x + alphas * los_dir.x,
loc.y + alphas * los_dir.y,
loc.z + alphas * los_dir.z,
]
inv_airmass = (
ipp[0] * ipp_vector.x.value
+ ipp[1] * ipp_vector.y.value
+ ipp[2] * ipp_vector.z.value
) # inproduct
inv_airmass /= r_earth + height_array[:, np.newaxis] # normalized
airmass = 1.0 / inv_airmass.value
return ipp, airmass
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment