Unverified Commit 60f07d8e authored by Akke Esa Tapio Viitanen's avatar Akke Esa Tapio Viitanen
Browse files

refactor tau, sf_inf logic

parent 27c75fac

opt/AGN_lightcurves_simulations/AGN_sim/DRW_sim.py

+27 −26
Original line number Diff line number Diff line
@@ -18,7 +18,8 @@ def abs_mag_to_apparent(abs_mag, redshift): 
    return apparent_mag





def get_logtau_logsf_inf(lambda_rest, mag_i, logMbh):

def get_tau_sf_inf(lambda_rest, mag_i, logMbh, z, type2, frame, divide_sf_inf_by_z1):



    # 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]
@@ -39,7 +40,29 @@ def get_logtau_logsf_inf(lambda_rest, mag_i, logMbh): 
        + coeff_sfinf[3] * (logMbh - 9)

    )



    return logtau, logsf_inf

    # 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)



        # NOTE: AGILE DR1 used the 1 / np.sqrt(1 + z)

        sf_inf = np.power(10, logsf_inf)

        if divide_sf_inf_by_z1:

            sf_inf /= np.sqrt(1 + z)



    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



    return tau, sf_inf





def LC2(
@@ -128,30 +151,8 @@ def LC2( 
    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)



        # NOTE: AGILE DR1 used the 1 / np.sqrt(1 + z)

        sf_inf = np.power(10, logsf_inf)

        if divide_sf_inf_by_z1:

            sf_inf /= np.sqrt(1 + z)



    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

    # Calculate tau, sf_inf

    tau, sf_inf = get_tau_sf_inf(lambda_rest, mag_i, logMbh, z, type2, frame, divide_sf_inf_by_z1)



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

    if meanmag is None: