#!/usr/bin/env python3
"""
Ivezic formulae for calculating LSST errors from Ivano 2025/03/18
"""
import pickle
import numpy as np
[docs]
def compute_sigma2(
magnitude,
band,
sigma_sys=0.005,
typical_m5=True,
airmass=1,
exp_time=30,
theta_eff=None,
m_sky=None,
*args,
**kwargs,
):
assert band in ["u", "g", "r", "i", "z", "y"]
dizionario = load_Ivezic_19_tab2(*args, **kwargs)
sigma2_rand = compute_sigma2_rand(
magnitude,
dizionario[band],
typical_m5=typical_m5,
airmass=airmass,
exp_time=exp_time,
theta_eff=theta_eff,
m_sky=m_sky,
)
sigma2 = sigma_sys**2 + sigma2_rand
return sigma2
[docs]
def load_Ivezic_19_tab2(dirname="./"):
with open(f"{dirname}data/Ivezic_19_tab_2.pickle", "rb") as f:
dizionario = pickle.load(f)
return dizionario
[docs]
def compute_sigma2_rand(
magnitude, dizionario, typical_m5=True, airmass=1, exp_time=30, theta_eff=None, m_sky=None
):
if typical_m5:
m5 = dizionario["m_5"]
else:
m5 = compute_m5(dizionario, airmass=airmass, exp_time=exp_time, theta_eff=theta_eff, m_sky=m_sky)
x = 10 ** (0.4 * (magnitude - m5))
gamma = dizionario["gamma"]
sigma2_rand = (0.04 - gamma) * x + gamma * x * x # Eq. 5 Ivezic+19
return sigma2_rand
[docs]
def compute_m5(dizionario, airmass=1, exp_time=30, theta_eff=None, m_sky=None):
## Computes the 5sigma depth magnitude, Eq 6 in Ivezic+19
theta_eff = dizionario["theta_eff"] if theta_eff is None else theta_eff
m_sky = dizionario["m_sky"] if m_sky is None else m_sky
m_5 = dizionario["C_m"] + 0.5 * (m_sky - 21)
m_5 += 2.50 * np.log10(0.7 / theta_eff)
m_5 += 1.25 * np.log10(exp_time / 30)
m_5 -= dizionario["k_m"] * (airmass - 1)
return m_5