Use get_rm_from_skycoord and compare with values from RMextract
This example script shows how to get ionospheric rotation measures for a give source, observer location and array of times. This script uses the default ionospheric settings, which are good for most purposes.
[1]:
from __future__ import annotations
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import EarthLocation, SkyCoord
from astropy.time import Time
from astropy.visualization import quantity_support, time_support
from spinifex import get_rm
[2]:
### Required to load local data for example - not needed for normal use
from importlib import resources
with resources.as_file(resources.files("spinifex.data.tests")) as datapath:
spinifex_data = datapath
###
[3]:
times = Time("2020-01-08T01:13:00") + np.arange(10) * 1 * u.hr
# create source object
source = SkyCoord(ra=350.85 * u.deg, dec=58.815 * u.deg)
# create Earth Location
lon = 6.367 * u.deg
lat = 52.833 * u.deg
dwingeloo = EarthLocation(lon=lon, lat=lat, height=0 * u.km)
# get rotation measures
rm = get_rm.get_rm_from_skycoord(
loc=dwingeloo,
times=times,
source=source,
# We set these options to use the data packaged with spinifex
# Unsetting them will cause the function to download the data from the internet
output_directory=spinifex_data,
)
INFO models - parse_iono_kwargs: Using ionospheric model <function get_density_ionex_single_layer at 0x7c8704932340> with options server=<Servers.CHAPMAN: 'chapman'> prefix='uqr' url_stem=None time_resolution=None solution='final' output_directory=PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/spinifex/checkouts/latest/spinifex/data/tests') correct_uqrg_rms=True height=<Quantity 350. km> remove_midnight_jumps=True
INFO get_ipp - _get_ipp_simple: Calculating ionospheric piercepoints
INFO get_rm - _get_rm: Calculating rotation measure
INFO download - download_file_http: Downloading from http://chapman.upc.es/tomion/rapid/2020/008_200108.15min/uqrg0080.20i.Z
INFO download - download_file_http: Downloading from http://chapman.upc.es/tomion/rapid/2020/009_200109.15min/uqrg0090.20i.Z
To get the values from the same server, source and time from RMextract, install RMextract and run the following lines:
[4]:
have_rmextract = False
if have_rmextract:
import RMextract
starttime = times[0].mjd * 24 * 3600 # MJD in seconds (casacore definition)
endtime = times[-1].mjd * 24 * 3600
statpos = [dwingeloo.get_itrs().cartesian.xyz.to(u.m).value]
pointing = [source.ra.rad, source.dec.rad] # 3C196 Ra, Dec in radians
rm_dict = RMextract.getRM.getRM(
ionexPath="./ionex_data/",
radec=pointing,
timestep=3600,
timerange=[starttime, endtime],
stat_positions=statpos,
earth_rot=1,
prefix="uqrg",
server="http://chapman.upc.es",
)
Here we directly take the values we got from the above code, to avoid installing RMextract
[5]:
rm_dict = {}
rm_dict["RM"] = {}
rm_dict["RM"]["st1"] = np.array(
[
[
*[0.10242756, 0.08804738, 0.07332737, 0.08251355],
*[0.07816817, 0.08137425, 0.11095578, 0.15794444],
*[0.26828872, 0.36951645, 0.44763038, 0.49045776],
]
]
).T
rm_dict["times"] = np.array(
[
*[5.08515918e09, 5.08516278e09, 5.08516638e09, 5.08516998e09],
*[5.08517358e09, 5.08517718e09, 5.08518078e09, 5.08518438e09],
*[5.08518798e09, 5.08519158e09, 5.08519518e09, 5.08519878e09],
]
)
rm_dict["azimuth"] = {}
rm_dict["azimuth"]["st1"] = np.array(
[
*[-0.5065002, -0.37364739, -0.23410709, -0.09012949],
*[0.05561559, 0.20026593, 0.34108016, 0.475697],
*[0.60222255, 0.7190818, 0.82458167, 0.9159993],
]
)
rm_dict["elev"] = {}
rm_dict["elev"]["st1"] = np.array(
[
*[0.52492472, 0.45731976, 0.40983145, 0.38421338],
*[0.38147449, 0.40172543, 0.44416086, 0.50718446],
*[0.58862976, 0.68600301, 0.79668206, 0.91802653],
]
)
[6]:
rm_rme = rm_dict["RM"]["st1"]
times_rme = Time(rm_dict["times"] / (3600 * 24), format="mjd")
az_rme = np.degrees(rm_dict["azimuth"]["st1"])
elev_rme = np.degrees(rm_dict["elev"]["st1"])
[7]:
# print to screen
rotation_measures = rm.rm
rm_times = rm.times
elevations = rm.elevation
azimuths = rm.azimuth
print("time RM (rad/lambda^2) azimuth (deg) elevation (deg)")
for myrm, tm, az, el in zip(rotation_measures, rm_times, azimuths, elevations):
print(f"{tm.isot} {myrm} {az:3.2f} {el:2.2f}")
print("Values from RMextract:")
for myrm, tm, az, el in zip(rm_rme, times_rme, az_rme, elev_rme):
print(f"{tm.isot} {myrm} {az:3.2f} {el:2.2f}")
time RM (rad/lambda^2) azimuth (deg) elevation (deg)
2020-01-08T01:13:00.000 0.08888034508198857 338.59 26.20
2020-01-08T02:13:00.000 0.0746200275021149 346.59 23.48
2020-01-08T03:13:00.000 0.08350599600850579 354.84 22.01
2020-01-08T04:13:00.000 0.07929786212900167 3.19 21.86
2020-01-08T05:13:00.000 0.08239537924587174 11.47 23.02
2020-01-08T06:13:00.000 0.11311698047238501 19.54 25.45
2020-01-08T07:13:00.000 0.15941962756095648 27.26 29.06
2020-01-08T08:13:00.000 0.26922353394519827 34.50 33.73
2020-01-08T09:13:00.000 0.3721337116368212 41.20 39.30
2020-01-08T10:13:00.000 0.45037082306460624 47.24 45.65
Values from RMextract:
2020-01-08T00:13:00.000 [0.10242756] -29.02 30.08
2020-01-08T01:13:00.000 [0.08804738] -21.41 26.20
2020-01-08T02:13:00.000 [0.07332737] -13.41 23.48
2020-01-08T03:13:00.000 [0.08251355] -5.16 22.01
2020-01-08T04:13:00.000 [0.07816817] 3.19 21.86
2020-01-08T05:13:00.000 [0.08137425] 11.47 23.02
2020-01-08T06:13:00.000 [0.11095578] 19.54 25.45
2020-01-08T07:13:00.000 [0.15794444] 27.26 29.06
2020-01-08T08:13:00.000 [0.26828872] 34.50 33.73
2020-01-08T09:13:00.000 [0.36951645] 41.20 39.31
2020-01-08T10:13:00.000 [0.44763038] 47.25 45.65
2020-01-08T11:13:00.000 [0.49045776] 52.48 52.60
[8]:
with time_support(), quantity_support():
fig, ax = plt.subplots()
im = ax.errorbar(rm.times.datetime, rm.rm, rm.rm_error, fmt=".", label="spinifex")
im2 = ax.plot(times_rme.datetime, rm_rme, "o-", label="RMextract")
ax.set_ylabel("RM / rad/m$^2$")
ax.legend()
plt.xticks(rotation=20)
