From 376fc8f888783c170b3780f11e7d46d2dfc85ab4 Mon Sep 17 00:00:00 2001 From: Noah Biederbeck Date: Mon, 8 May 2023 07:15:46 +0000 Subject: [PATCH 1/2] Add first working version with sin(delta) and cos(zenith). --- link-paths.py | 58 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 47 insertions(+), 11 deletions(-) diff --git a/link-paths.py b/link-paths.py index 0c462a7..20af4b7 100644 --- a/link-paths.py +++ b/link-paths.py @@ -5,7 +5,7 @@ import astropy.units as u import numpy as np -from astropy.coordinates import AltAz, EarthLocation +from astropy.coordinates import AltAz, EarthLocation, angular_separation from astropy.table import Table from astropy.time import Time from tqdm import tqdm @@ -45,6 +45,24 @@ template_linkname_model = outdir_model.resolve().as_posix() +def sin_delta(altaz: AltAz): + """Delta is the angle between pointing and magnetic field.""" + b_inc = u.Quantity(-37.36, u.deg) + b_dec = u.Quantity(-4.84, u.deg) + + delta = angular_separation( + lon1=altaz.az, + lat1=altaz.alt, + lon2=b_dec, + lat2=b_inc, + ) + return np.sin(delta.to_value(u.rad)) + + +def cos_zenith(altaz: AltAz): + return np.cos(altaz.zen.to_value(u.rad)) + + @u.quantity_input def build_altaz(*, alt: u.deg = None, zd: u.deg = None, az: u.deg = None) -> AltAz: """Build AltAz from zenith distance and azimuth. @@ -95,6 +113,10 @@ def get_pointings_of_irfs(filelist) -> AltAz: return build_altaz(zd=theta * u.deg, az=az * u.deg) +def euclidean_distance(x1, y1, x2, y2): + return np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) + + def main() -> None: prod = args.prod dec = args.dec @@ -105,6 +127,9 @@ def main() -> None: filelist = [p.name for p in path.parent.iterdir()] irf_pointings: AltAz = get_pointings_of_irfs(filelist) + irf_sin_delta = sin_delta(irf_pointings) + irf_cos_zenith = cos_zenith(irf_pointings) + progress = tqdm(total=n_runs) for night, run_ids in runs.items(): @@ -117,7 +142,15 @@ def main() -> None: alt=run["mean_altitude"] * u.rad, az=run["mean_azimuth"] * u.rad, ) - nearest_irf = pointing.separation(irf_pointings).argmin() + sindelta = sin_delta(pointing) + coszenith = cos_zenith(pointing) + + nearest_irf = euclidean_distance( + x1=sindelta, + y1=coszenith, + x2=irf_sin_delta, + y2=irf_cos_zenith, + ).argmin() node = filelist[nearest_irf] linkname_dl1 = Path( @@ -133,17 +166,20 @@ def main() -> None: target_dl2 = Path(template_target_dl2.format(prod=prod, dec=dec, node=node)) linkname_dl2 = Path(template_linkname_dl2.format(run_id=run_id)) - if not linkname_dl1.exists(): - linkname_dl1.parent.mkdir(exist_ok=True, parents=True) - linkname_dl1.symlink_to(target_dl1) + linkname_dl1.parent.mkdir(exist_ok=True, parents=True) + if linkname_dl1.exists() and linkname_dl1.is_symlink(): + linkname_dl1.unlink() + linkname_dl1.symlink_to(target_dl1) - if not linkname_dl2.exists(): - linkname_dl2.parent.mkdir(exist_ok=True, parents=True) - linkname_dl2.symlink_to(target_dl2) + linkname_dl2.parent.mkdir(exist_ok=True, parents=True) + if linkname_dl2.exists() and linkname_dl2.is_symlink(): + linkname_dl2.unlink() + linkname_dl2.symlink_to(target_dl2) - if not linkname_model.exists(): - linkname_model.parent.mkdir(exist_ok=True, parents=True) - linkname_model.symlink_to(target_model) + linkname_model.parent.mkdir(exist_ok=True, parents=True) + if linkname_model.exists() and linkname_model.is_symlink(): + linkname_model.unlink() + linkname_model.symlink_to(target_model) progress.update() From 0eba74ec7ee8d8e251285404c14fdbfbf90754a7 Mon Sep 17 00:00:00 2001 From: Noah Biederbeck Date: Tue, 9 May 2023 15:31:31 +0200 Subject: [PATCH 2/2] address review comments --- link-paths.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/link-paths.py b/link-paths.py index 20af4b7..1483ace 100644 --- a/link-paths.py +++ b/link-paths.py @@ -47,6 +47,12 @@ def sin_delta(altaz: AltAz): """Delta is the angle between pointing and magnetic field.""" + # Values from + # https://geomag.bgs.ac.uk/data_service/models_compass/igrf_calc.html + # for La Palma coordinates and date=2021-12-01. + # + # Also used here: + # https://github.com/cta-observatory/lst-sim-config/issues/2 b_inc = u.Quantity(-37.36, u.deg) b_dec = u.Quantity(-4.84, u.deg) @@ -127,8 +133,8 @@ def main() -> None: filelist = [p.name for p in path.parent.iterdir()] irf_pointings: AltAz = get_pointings_of_irfs(filelist) - irf_sin_delta = sin_delta(irf_pointings) - irf_cos_zenith = cos_zenith(irf_pointings) + irf_sindelta = sin_delta(irf_pointings) + irf_coszenith = cos_zenith(irf_pointings) progress = tqdm(total=n_runs) @@ -148,8 +154,8 @@ def main() -> None: nearest_irf = euclidean_distance( x1=sindelta, y1=coszenith, - x2=irf_sin_delta, - y2=irf_cos_zenith, + x2=irf_sindelta, + y2=irf_coszenith, ).argmin() node = filelist[nearest_irf]