Skip to content
Snippets Groups Projects

simple ipp model

Closed Maaijke Mevius requested to merge ipp_calculation into main
All threads resolved!
Files
7
+ 93
13
@@ -11,37 +11,117 @@ import astropy.units as u
import numpy as np
from astropy.coordinates import ITRS, AltAz, EarthLocation, SkyCoord
from astropy.time import Time
from numpy.typing import ArrayLike
r_earth = 6364.62 * u.km
class IPP(NamedTuple):
"""Ionospheric Piercepoints"""
"""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: ArrayLike
loc: EarthLocation, times: Time, source: SkyCoord, height_array: u.Quantity
) -> IPP:
"""Get the ionospheric piercepoints as EarhtLocations for a given EarthLocation, time, SkyCoord"""
# height_array = np.arange(100, 1000, 10) * u.km
"""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, times)
return get_ipp_from_altaz(loc, altaz, height_array)
def get_ipp_from_altaz(
loc: EarthLocation, altaz: AltAz, height_array: ArrayLike, times: Time
loc: EarthLocation, altaz: AltAz, height_array: u.Quantity
) -> IPP:
altaz_dir = altaz.transform_to(ITRS)
ipp_vector = altaz_dir.cartesian[np.newaxis] * height_array[:, np.newaxis]
"""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.to(u.m).x + ipp_vector.x,
loc.to(u.m).y + ipp_vector.y,
loc.to(u.m).z + ipp_vector.z,
loc.x + alphas * los_dir.x,
loc.y + alphas * los_dir.y,
loc.z + alphas * los_dir.z,
]
return IPP(loc=EarthLocation.from_geocentric(*ipp), times=times, los=altaz_dir)
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
Loading