Commit ed0cf542 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

Refactor AGN catalog into a separate python module

parent ddc07f7e

.gitignore

+2 −1
Original line number Diff line number Diff line 


*.png



._*
@@ -12,3 +11,5 @@ opt/egg 
opt/imsim

opt/lsst_stack

opt/vif



data/baseline_*.*

src/python/catalog_agn.py

0 → 100644
+790 −0

File added.

Preview size limit exceeded, changes collapsed.

src/python/util.py

+15 −11
Original line number Diff line number Diff line
@@ -11,6 +11,7 @@ General utilities 
import logging

logger = logging.getLogger()



from functools import cache

import os



import astropy.units as u
@@ -18,8 +19,20 @@ import matplotlib.pyplot as plt 
import numpy as np

from scipy.stats import binned_statistic



from astropy.cosmology import FlatLambdaCDM

from astropy import constants

from astropy.coordinates import SkyCoord

from astropy.cosmology import FlatLambdaCDM

from astropy.time import Time

import zou2024

from igm_absorption import my_get_IGM_absorption

from scipy.interpolate import interp1d

from scipy.stats import gaussian_kde



try:

    from quasarlf.pubtools.load_observations import my_get_data

    from quasarlf.pubtools.load_observations import get_data

except (ImportError, OSError):

    logger.warning("Could not load quasarlf")





ROOT = "/staff2/viitanen/agile"
@@ -100,8 +113,6 @@ def egg_band_to_index(egg, band): 
    return bands.index(band)





from astropy.time import Time

from astropy.coordinates import SkyCoord

def get_ra_dec(ra0, dec0, pm_ra_cosdec, pm_dec, mjd, mjd0=51544.5):



    """Get ra, dec for current epoch"""
@@ -250,7 +261,6 @@ def get_E_BV( 
    return np.squeeze(ebv)





from igm_absorption import my_get_IGM_absorption

def luminosity_to_flux(wavlen, luminosity, redshift, distance_in_cm, use_igm=True):



    """
@@ -418,7 +428,6 @@ def get_logn_logs2( 


    # NOTE: optimization on the calculation of the volume/distance

    print("Estimating dV and dl...")

    from scipy.interpolate import interp1d

    zbin = np.arange(z.min() - dz, z.max() + dz + 1e-6, dz)

    fun_dV = interp1d(zbin, cosmo.comoving_volume(zbin).value)

    fun_dl = interp1d(zbin, cosmo.luminosity_distance(zbin).cgs.value)
@@ -521,7 +530,6 @@ def get_star_binary_fbin(star, binary, fbin=0.40, nrepeat=4, seed=1206): 
    return star[is_star], binary2





from scipy.stats import gaussian_kde

def plot_contour2d(fig, ax, x, y, ps=[0.99, 0.95, 0.68], *args, **kwargs):



    print("NOTE, selecting 10,000 points at random...")
@@ -625,7 +633,6 @@ def get_fraction_catastrophic_error(value_estimated, value_true, limit=0.15): 
@np.vectorize

def compare_catalog_to_shen2020(zbin, catalog, dz, dlog_lx, debug=False, *args, **kwargs):



    from quasarlf.pubtools.load_observations import my_get_data

    x1, y1, dy1, dx1 = my_get_data(-4, zbin)



    bins_log_lx = np.arange(32, 52, dlog_lx)
@@ -719,14 +726,12 @@ def get_logN_B( 
    return np.sum(to_sum) * dB * dz





from functools import cache

def get_phi_shen2020(z):



    @cache

    def get(dataid, lam, z):



        # Return the Shen luminosity function for a single dataid

        from quasarlf.pubtools.load_observations import get_data

        x, y, dy, dx = get_data(dataid, z)



        if dataid == -4:
@@ -753,10 +758,9 @@ def get_phi_shen2020(z): 
    return [get(dataid, lam, z) for dataid, lam in dataset]





import zou2024

def get_log_lambda_SAR(i, N, m, z, t):

    ret = zou2024.get_log_lambda_SAR(

        m, z, t, add_ctk=True, log_mstar_lim=(9.50, 12.00), z_lim=(0.00, 4.00)

    )

    print("%8d" % i, "%8d" % N, "%8.2f" % (i / N * 100), "%8.2f" % m, "%8.2f" % z, "%20s" % t, "%8.2f" % ret, end='\r')

    #print("%8d" % i, "%8d" % N, "%8.2f" % (i / N * 100), "%8.2f" % m, "%8.2f" % z, "%20s" % t, "%8.2f" % ret, end='\r')

    return ret

src/tests/test_catalog_agn.py

0 → 100644
+30 −0
Original line number Diff line number Diff line 
#!/usr/bin/env python3

# encoding: utf-8

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-06-21 21:57:44



"""

Test AGN catalog

"""



import logging

logger = logging.getLogger()





from egg import Egg

from catalog_agn import CatalogAGN



from astropy.io import fits

egg = fits.open("data/catalog/test/egg.fits")[1].data

c = CatalogAGN(

    "data/test/test_agn",

    egg,

    "zou+2024",

    1,

    20250621,

    1,

    0,

    0.05,

    0.95

)