add AGN_lightcurves_simulation

opt/lsst_stack

opt/vif

opt/AGILE_Mstar_Mbh

opt/AGN_lightcurves_simulations

#opt/AGN_lightcurves_simulations

opt/CEP_lsst_quasar_project

opt/GrowBHs

opt/mock_catalog_SED

#opt/mock_catalog_SED

opt/quasarlf

opt/ivano_feature_extraction

.ipynb_checkpoints

__pycache__

*~

*.swp

.DS_store

*.orig

import numpy as np

from astropy.cosmology import FlatLambdaCDM

# Define the cosmological parameters (you can adjust these as needed)

H0 = 70  # Hubble constant in km/s/Mpc

Om0 = 0.3  # Matter density parameter

cosmo = FlatLambdaCDM(H0=H0, Om0=Om0)

# Authors: Andjelka Kovacevic, Isidora Jankov & Viktor Radovic

# Additional contributions by Vincenzo Petrecca

# Function to convert absolute magnitude to apparent magnitude

def abs_mag_to_apparent(abs_mag, redshift):

    luminosity_distance = cosmo.luminosity_distance(redshift)

    apparent_mag = abs_mag + 5 * np.log10(luminosity_distance.to("pc").value / 10)

    return apparent_mag

def get_logtau_logsf_inf(lambda_rest, mag_i, logMbh):

    # Constants (A,B,C,D from Tab.2 of Suberlak et al. 2020)

    coeff_tau = [2.597, 0.17, 0.035, 0.141]

    coeff_sfinf = [-0.476, -0.479, 0.118, 0.118]

    # Calculate damping time scale (Eq 12, Suberlak et al. 2020)

    logtau = (

        coeff_tau[0]

        + coeff_tau[1] * (np.log10(lambda_rest / 4000))

        + coeff_tau[2] * (mag_i + 23)

        + coeff_tau[3] * (logMbh - 9)

    )

    # Calculate SF at infinity (related to sigma, amplitude of correlation decay)

    logsf_inf = (

        coeff_sfinf[0]

        + coeff_sfinf[1] * (np.log10(lambda_rest / 4000))

        + coeff_sfinf[2] * (mag_i + 23)

        + coeff_sfinf[3] * (logMbh - 9)

    )

    return logtau, logsf_inf

def LC2(

    T=7300,

    deltatc=1,

    z=0,

    mag_i=-26,

    logMbh=9.0,

    lambda_rest=4000,

    oscillations=False,

    A=0.14,

    noise=0.00005,

    frame="observed",

    type2=False,

    meanmag=None,

):

    """

    Generate one artificial light curve using a stochastic model based on the

    Damped random walk (DRW) proccess and scaling relations from Suberlak+21.

    Parameters describing the model are characteristic amplitude ("logsig2" in

    code) and time scale of the exponentially-decaying variability ("tau" in

    code), both infered from physical quantities such are supermassive black

    hole mass and/or luminosity of the AGN. For further details regarding the

    model see Kovačević et al. (2021) and references therein.

    Parameters:

    -----------

    T: int, default=7300 (20 years)

        Total time span of the light curve in days. It is recommended to generate light curves to be at least

        10 times longer than their characteristic timescale (Kozłowski 2017).

    deltatc: int, default=1

        Cadence (or sampling rate) - time interval between two consecutive samplings of the light

        curve in days.

    z: float, default=0

        Redshift.

    mag_i: float, default=-26

        Mean i-band absolute magnitude

    logMbh: flat, default=9.0

        Log of Black Hole Mass

    lamda_rest: float, default=4000

        Rest frame wavelenght (Angstrom)

    oscillations: bool, default=False

        If True, light curve simulation will take an oscillatory signal into account.

    A: float, default=0.14

        Amplitude of the oscillatory signal in magnitudes (used only if oscillations=True).

    noise: float, default=0.00005

        Amount of noise to include in the light curve simulation.

    frame: {'observed', 'rest'}, default='observed'

        Frame of reference.

    type2: bool, default=False

        If True, we use a Damped-DRW model for type2, reducing both variability

        amplitude and SF slope (through tau)

    Returns:

    --------

    tt: np.array

        Days when the light curve was sampled.

    yy: np.array

        Magnitudes of the simulated light curve.

    References:

    -----------

    Ivezić, Ž., et al. 2019, ApJ, 873, 111 (https://iopscience.iop.org/article/10.3847/1538-4357/ab042c)

    Kelly, B.C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895 (https://iopscience.iop.org/article/10.1088/0004-637X/698/1/895)

    Kovačević, A., et al. 2021, submitted to MNRAS (https://github.com/LSST-sersag/white_paper/blob/main/data/paper.pdf)

    Kozłowski, S. 2017, A&A, 597, A128 (https://www.aanda.org/articles/aa/full_html/2017/01/aa29890-16/aa29890-16.html)

    """

    # print(

    #    "T =", T,

    #    "deltac =", deltatc,

    #    "z =", z,

    #    "mag_i =", mag_i,

    #    "logMbh =", logMbh,

    #    "lambda_rest =", lambda_rest,

    #    "oscillations =", oscillations,

    #    "A =", A,

    #    "noise =", noise,

    #    "frame =", frame,

    #    "type2 =", type2,

    #    "meanmag =", meanmag,

    # )

    # Generating survey days

    tt = np.arange(0, T, int(deltatc))

    times = tt.shape[0]

    # Calculate logtau, logsf_inf

    logtau, logsf_inf = get_logtau_logsf_inf(lambda_rest, mag_i, logMbh)

    # tau \propto lambda_rest ** (-4 / 3)

    if frame == "observed":

        # Convering to observed frame (Eq 17, 18 from Kelly et al. 2009)

        tau = np.power(10, logtau) * (1 + z)

        sf_inf = np.power(10, logsf_inf) / np.sqrt(1 + z)

        # sf_inf = np.power(10,logsf_inf)

    elif frame == "rest":

        tau = np.power(10, logtau)

        sf_inf = np.power(10, logsf_inf)

    # Damping factors on Tau and SF_inf (D-DRW for Type 2 AGNs)

    tau = np.where(type2, tau / 10, tau)

    sf_inf = np.where(type2, sf_inf / 10, sf_inf)

    # if type2 == True:

    #    tau = tau / 10

    #    sf_inf = sf_inf / 10

    # Get the mean magnitude from Mag_i (Cosmological parameters are set above)

    if meanmag is None:

        meanmag = mag_i

        if frame == "observed":

            meanmag = abs_mag_to_apparent(mag_i, z)

    # Calculating light curve magnitudes

    mag_i = np.atleast_1d(mag_i)

    ratio = -deltatc / tau

    sigma = sf_inf / np.sqrt(tau)  # See end of Sec. 2.2 (MacLeod et al. 2010)

    SFCONST2 = np.power(sigma, 2)

    def get_ss1():

        ss = np.zeros((mag_i.size, times))

        ss[:, 0] = meanmag  # light curve is initialized

        for i in range(1, times):

            ss[:, i] = np.random.normal(

                ss[:, i - 1] * np.exp(ratio) + meanmag * (1 - np.exp(ratio)),

                np.sqrt(0.5 * SFCONST2 * tau * (1 - np.exp(2 * ratio))),

            )

        return ss

    ss = get_ss1()

    # OPTIONAL: Calculating error (Ivezic et al. 2019)

    # gamma=0.039

    # m5=24.7

    # x=np.zeros(ss.shape)

    # x=np.power(10, 0.4*(ss-m5))

    # err = (0.005*0.005) + (0.04-gamma)*x + gamma*x*x

    ## Final light curve with oscillations

    # if oscillations:

    #    # Calculate underlying periodicity

    #    conver = 173.145  # convert from LightDays to AU

    #    lightdays = 10

    #    P = np.sqrt(((lightdays * conver) ** 3) / (msmbh))

    #    # Calculating and adding oscillatory signal

    #    sinus = A * np.sin(2 * np.pi * tt / (P * 365))

    #    ss = ss + sinus

    #    for i in range(times):

    #        # Adding error and noise to each magnitude value

    #        yy[:, i] = ss[:, i]  # + np.random.normal(0,((noise*ss[i])),1) + np.sqrt(err[i])

    #    return tt, yy, tau, sf_inf

    # Final light curve without oscillations

    if not oscillations:

        yy = np.zeros((mag_i.size, times))

        yy[:, :] = ss[:, :]

        # AV 20241004

        if False:

            for i in range(times):

                print(i, times)

                # Adding error and noise to each magnitude value

                yy[:, i] = ss[:, i]  # + np.random.normal(0,((noise*ss[i])),1) + np.sqrt(err[i])

        return tt, np.squeeze(yy), tau, sf_inf, meanmag, sigma

############################################################################################################

%% Cell type:markdown id:86d90505 tags:

# Simulation of AGN light curves using Damped Random Walk (DRW)

%% Cell type:code id:11112661 tags:

``` python

# Import standard libraries

# Import simulation code

import DRW_sim

import matplotlib.pyplot as plt

import numpy as np

```

%% Cell type:code id:d500da76 tags:

``` python

# %matplotlib inline

plt.rcParams.update(

    {"xtick.top": True, "ytick.right": True, "xtick.direction": "in", "ytick.direction": "in"}

)

```

%% Cell type:markdown id:3d0d70d1 tags:

#### Simulate a DRW light curve

%% Cell type:code id:d7661124 tags:

``` python

np.random.seed(9)

# np.random.seed(434)

```

%% Cell type:raw id:2f89c0aa tags:

from astropy.cosmology import FlatLambdaCDM

# Define the cosmological parameters (you can adjust these as needed)

H0 = 70  # Hubble constant in km/s/Mpc

Om0 = 0.3  # Matter density parameter

cosmo = FlatLambdaCDM(H0=H0, Om0=Om0)

# Function to convert absolute magnitude to apparent magnitude

def abs_mag_to_apparent(abs_mag, redshift):

    luminosity_distance = cosmo.luminosity_distance(redshift)

    apparent_mag = abs_mag + 5 * np.log10(luminosity_distance.to('pc').value / 10)

    return apparent_mag

%% Cell type:code id:5b510c34 tags:

``` python

redshift = 1.0

mass = 9.0

wv_rest = 4000

abs_mag_i = -25

```

%% Cell type:code id:461ddb92 tags:

``` python

# New scaling relations (including flag for AGN type)

tt, yy, tau, sf_inf, meanmag, sigma = DRW_sim.LC2(

    T=7300,

    deltatc=1,

    z=redshift,

    mag_i=abs_mag_i,

    logMbh=mass,

    lambda_rest=wv_rest,

    oscillations=False,

    frame="observed",

    type2=False,

)

# Take only 10 years (full LSST lenght)

tt = tt[:3650]

yy = yy[:3650]

```

%% Cell type:code id:17666179 tags:

``` python

print("Damping time:", tau)

print("Amplitude:", sf_inf)

print("Sigma:", sigma)

print("\nMean magnitude (input):", meanmag)

print("Mean magnitude (output):", np.mean(yy))

```

%% Cell type:code id:29797ac8 tags:

``` python

# Plot the obtained light curve

fig = plt.figure(figsize=(6, 4), dpi=120, tight_layout=True)

plt.scatter(tt, yy, marker=".", s=2, label="Continuos")

plt.minorticks_on()

plt.gca().invert_yaxis()

plt.xlabel("t [days]")

plt.ylabel("magnitude")

# plt.legend()

# plt.savefig('sim_drw.png')

plt.show()

```

%% Cell type:markdown id:791caef3 tags:

### Define cadence

%% Cell type:code id:9de5a802 tags:

``` python

# Add a typical survey cadence (ugrizy == 0,1,2,3,4,5)

# cad_tt, cad_yy = DRW_sim.var_cad(tt, yy, m3=mar[2],m4=apr[2],m5=may[2],m6=jun[2])

# cad_tt, cad_yy = DRW_sim.var_cad(tt, yy, m3=5,m4=5)

cad_tt, cad_yy = DRW_sim.var_cad(tt, yy, m3=4, m4=5, m5=5, m6=3)

```

%% Cell type:code id:2a3bd507 tags:

``` python

# Plot the obtained light curve

fig = plt.figure(figsize=(6, 4), dpi=120, tight_layout=True)

plt.scatter(tt, yy, marker=".", s=2, label="Continuos")

plt.scatter(cad_tt, cad_yy, marker="o", s=4, label="Typical survey")

plt.minorticks_on()

plt.gca().invert_yaxis()

plt.xlabel("t [days]")

plt.ylabel("magnitude")

plt.legend(markerscale=2, frameon=False)

# plt.savefig('sim_drw.png')

plt.show()

```

%% Cell type:markdown id:77e90bbc tags:

### Simulate different filters

%% Cell type:code id:70074d52 tags:

``` python

# LSST central wavelenght for each filter (u,g,r,i,z,y in Angstrom)

lambda_lsst = [3560, 4767, 6216, 7545, 8708, 10042]

```

%% Cell type:code id:8835a62a tags:

``` python

labels = ["u", "g", "r", "i", "z", "y"]

# Plot multi-wavelenght light curves

fig = plt.figure(figsize=(7, 5), dpi=120, tight_layout=True)

for i in range(6):

    np.random.seed(9)

    # Get rest-frame wavelenght for each filter

    flt = lambda_lsst[i] / (1 + redshift)

    # Simulate DRW light curves

    tt, yy, tau, sf_inf, meanmag, sigma = DRW_sim.LC2(

        T=7300,

        deltatc=1,

        z=redshift,

        mag_i=abs_mag_i,

        logMbh=mass,

        lambda_rest=flt,

        oscillations=False,

        frame="observed",

        type2=False,

    )

    # Take only 10 years (full LSST lenght)

    tt = tt[:3650]

    yy = yy[:3650]

    print(tau, sf_inf)

    # Plot

    plt.scatter(tt, yy, marker="o", s=2, alpha=0.7, label=labels[i])

    # plt.plot(tt, yy, linewidth=0.6, alpha=0.9, label=labels[i])

plt.minorticks_on()

plt.gca().invert_yaxis()

plt.xlabel("t [days]")

plt.ylabel("magnitude")

plt.legend(markerscale=5, frameon=False)

# plt.legend(bbox_to_anchor=(1, 1))

plt.show()

```

%% Cell type:markdown id:24f79802 tags:

### Type1 and Typ2 togheter

%% Cell type:code id:9268cd72 tags:

``` python

np.random.seed(9)

```

%% Cell type:code id:b3dc88be tags:

``` python

tt1, yy1, tau1, sf_inf1, meanmag1, sigma1 = DRW_sim.LC2(

    T=7300,

    deltatc=1,

    z=redshift,

    mag_i=abs_mag_i,

    logMbh=mass,

    lambda_rest=wv_rest,

    oscillations=False,

    frame="observed",

    type2=False,

)

tt2, yy2, tau2, sf_inf2, meanmag2, sigma2 = DRW_sim.LC2(

    T=7300,

    deltatc=1,

    z=redshift,

    mag_i=abs_mag_i,

    logMbh=mass,

    lambda_rest=wv_rest,

    oscillations=False,

    frame="observed",

    type2=True,

)

# Take only 10 years (full LSST lenght)

tt1 = tt1[:3650]

yy1 = yy1[:3650]

tt2 = tt2[:3650]

yy2 = yy2[:3650]

```

%% Cell type:code id:7ce8b3e4 tags:

``` python

# Plot the obtained light curve

fig = plt.figure(figsize=(6, 4), dpi=120, tight_layout=True)

plt.scatter(tt1, yy1, marker=".", s=2, color="blue", label="Type1")

plt.scatter(tt2, yy2, marker=".", s=2, color="red", label="Type2 (D-DRW)")

plt.minorticks_on()

plt.gca().invert_yaxis()

plt.xlabel("t [days]")

plt.ylabel("magnitude")

plt.legend(markerscale=5, frameon=False)

plt.show()

```

%% Cell type:code id:ff9485a9 tags:

``` python

print("----- TYPE 1 AGN ------\n")

print("Damping time:", tau1)

print("Amplitude:", sf_inf1)

print("Sigma:", sigma1)

print("Mean magnitude (input):", meanmag1)

print("Mean magnitude (output):", np.mean(yy1))

print("\n----- TYPE 2 AGN (D-DRW) ------\n")

print("Damping time:", tau2)

print("Amplitude:", sf_inf2)

print("Sigma:", sigma2)

print("Mean magnitude (input):", meanmag2)

print("Mean magnitude (output):", np.mean(yy2))

```

%% Cell type:markdown id:e0e957ee tags:

### Try a loop (Optional)

%% Cell type:code id:50095df6 tags:

``` python

for i in range(100):

    tt, yy, tau, sf_inf, meanmag, sigma = DRW_sim.LC2(

        T=7300,

        deltatc=1,

        z=redshift,

        mag_i=abs_mag_i,

        logMbh=mass,

        lambda_rest=wv_rest,

        oscillations=False,

        frame="observed",

        type2=False,

    )

    tt2, yy2, tau2, sf_inf2, meanmag2, sigma2 = DRW_sim.LC2(

        T=7300,

        deltatc=1,

        z=redshift,

        mag_i=abs_mag_i,

        logMbh=mass,

        lambda_rest=wv_rest,

        oscillations=False,

        frame="observed",

        type2=True,

    )

    tt = tt[:3650]

    yy = yy[:3650]

    tt2 = tt2[:3650]

    yy2 = yy2[:3650]

    print("Obj:", i)

    print("Damping time:", tau, "----", tau2)

    print("SF_infty:", sf_inf, "----", sf_inf2)

    # print("Damping time (obserbed):",tau*(1+redshift))

    print("Mean magnitude:", meanmag, "----", np.mean(yy), "\n\n")

    # Plot the obtained light curve

    fig = plt.figure(figsize=(6, 4), dpi=120, tight_layout=True)

    plt.scatter(tt, yy, marker=".", s=2, label="Type 1", color="blue")

    plt.scatter(tt2, yy2, marker=".", s=2, label="Type 2", color="red")

    plt.minorticks_on()

    plt.gca().invert_yaxis()

    plt.xlabel("t [days]")

    plt.ylabel("magnitude")

    # plt.title('Object ',+str(i))

    plt.legend(markerscale=5, frameon=False)

    plt.show()

```

%% Cell type:code id:6173b0a9 tags:

``` python

```