Merge pull request #64 from nbiederbeck/allsky-mc-sindelta-coszenith

Map IRFs with sin(delta) and cos(zenith)

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,30 @@
 template_linkname_model = outdir_model.resolve().as_posix()
 
 
+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)
+
+ 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 +119,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 +133,9 @@ def main() -> None:
  filelist = [p.name for p in path.parent.iterdir()]
  irf_pointings: AltAz = get_pointings_of_irfs(filelist)
 
+ irf_sindelta = sin_delta(irf_pointings)
+ irf_coszenith = cos_zenith(irf_pointings)
+
  progress = tqdm(total=n_runs)
 
  for night, run_ids in runs.items():
@@ -117,7 +148,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_sindelta,
+ y2=irf_coszenith,
+ ).argmin()
  node = filelist[nearest_irf]
 
  linkname_dl1 = Path(
@@ -133,17 +172,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()