add ivezic formulae

data/baseline_*.*

# opt software

opt/astro

opt/egg

opt/imsim

opt/lsst_stack

            #"exclude_patterns=notebooks/*,_build", # Exclude notebooks and build dir from pre-commit

            "exclude_patterns=_build", # Exclude notebooks and build dir from pre-commit

            "-j 4", # Use multiple threads

            "-d .sphinx-cache", # Use persistent cache

          ]

    # Run unit tests, verify that they pass. Note that coverage is run against

    # the ./src directory here because that is what will be committed. In the

"""Middleware for reference catalog config."""

import numpy as np

import util

from lsst.meas.algorithms.convertRefcatManager import ConvertGaiaManager

from lsst_inaf_agile import util

class ConvertAgileManager(ConvertGaiaManager):

    """

#!/usr/bin/env python3

"""

Ivezic formulae for calculating LSST errors from Ivano 2025/03/18

"""

import pickle

import numpy as np

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

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

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

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

# Modified script to try to recover the m5 values in the paper.

#

# Original script was in src.bak2

import fitsio

import matplotlib.pyplot as plt

import numpy as np

from lsst_inaf_agile.ivezic import compute_m5, load_Ivezic_19_tab2

from scipy import stats

magmin = 15.0

magmax = 28.0

dbin = 1.00

baseline = fitsio.read("data/baseline/baseline_v4.0_10yrs.fits", ext=1)

def get_baseline_airmass_seeing_sky_exptime(baseline, b, is_cosmos=True, statistic=np.median):

    select_filter = baseline["filter"] == b

    select_cosmos = baseline["target_name"] == "DD:COSMOS"

    select = select_filter * select_cosmos

    airmass = statistic(baseline[select]["airmass"])

    seeing = statistic(baseline[select]["seeingFwhmEff"])

    sky = statistic(baseline[select]["skyBrightness"])

    exptime = baseline["visitExposureTime"][select]

    assert np.all(exptime[0] == exptime[1:])

    exptime = exptime[0]

    assert 1.0 < airmass < 2.0

    assert 0.5 < seeing < 2.0

    assert 15.0 < sky < 25.0

    return airmass, seeing, sky, exptime

def get_delta_mag(flux, dflux):

    return 2.5 * (dflux / flux) / np.log(10)

def main(magmin, magmax, dbin):

    fig, axes = plt.subplots(2, 3, figsize=(3 * 6.4, 2 * 4.8))

    print("# band airmass seeing sky m5")

    for b in "ugrizy":

        for statistic, label in ((np.median, "median"), (lambda a: stats.mode(a).mode, "mode")):

            airmass, seeing, sky, exptime = get_baseline_airmass_seeing_sky_exptime(

                baseline, b, statistic=statistic

            )

            # sigma = (

            #    compute_sigma2(

            #        magvec,

            #        b,

            #        typical_m5=False,

            #        airmass=airmass,

            #        sigma_sys=0.0,

            #        theta_eff=seeing,

            #        m_sky=sky,

            #        exp_time=exptime,

            #    )

            #    ** 0.5

            # )

            m5_ivezic = compute_m5(load_Ivezic_19_tab2()[b], airmass, exptime, seeing, sky)

            print(

                b,

                "%6.2f" % airmass,

                "%6.2f" % seeing,

                "%6.2f" % sky,

                "%6.2f" % m5_ivezic,

                "(%s)" % label,

                flush=True,

            )

        print()

if __name__ == "__main__":

    main(magmin, magmax, dbin)