Create a bending power law light curve

import numpy as np
import matplotlib.pyplot as plt
from slsim.Util import astro_util
from slsim.Sources.SourceVariability import accretion_disk_reprocessing
from slsim.Sources.SourceVariability import variability
import warnings

warnings.filterwarnings(action="ignore")

# define light curve parameters
length_of_light_curve = 1000  # Days
dt = 1  # Days
mean_magnitude = 10  # magnitude
standard_deviation = 0.2
# I am defining a random seed so the light curve
# is consistent through all realizations
random_seed = 14

# define bending power law parameters
low_frequency_index = 1
high_frequency_index = 3
log_breakpoint_frequency = -1

bpl_kwarg_dict = {
    "length_of_light_curve": length_of_light_curve,
    "time_resolution": dt,
    "log_breakpoint_frequency": log_breakpoint_frequency,
    "low_frequency_slope": low_frequency_index,
    "high_frequency_slope": high_frequency_index,
    "mean_magnitude": mean_magnitude,
    "standard_deviation": standard_deviation,
    "seed": random_seed,
}

time_stamps, light_curve = astro_util.generate_signal_from_bending_power_law(
    **bpl_kwarg_dict
)

print("mean: " + str(light_curve.mean()))
print("standard deviation: " + str(light_curve.std()))

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_stamps, light_curve)
ax.set_xlabel("time [days]")
ax.set_ylabel("flux [mag.]")
plt.show()

mean: 10.0
standard deviation: 0.19999999116369654

Create identical light curve using variability module

variable_light_curve = variability.Variability(
    variability_model="bending_power_law", **bpl_kwarg_dict
)

some_time_stamps = np.linspace(-50, 1100, 20000)
interpolated_light_curve = variable_light_curve.variability_at_time(some_time_stamps)
# interpolation to timestamps outside defined light curve default to 0

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(some_time_stamps, interpolated_light_curve)
ax.set_xlabel("time [days]")
ax.set_ylabel("flux [mag.]")
plt.show()

Create a light curve using a defined power spectrum density

frequencies = astro_util.define_frequencies(length_of_light_curve, dt)
power_spectrum_1 = frequencies ** (
    -3
)  # Note the slope of (-3) suppresses high frequency signals.
power_spectrum_2 = frequencies ** (-1)  # A slope of (-1) becomes very noisy.

other_time_stamps = np.linspace(200, 800, 1000)

psd_kwarg_dict_1 = {
    "length_of_light_curve": length_of_light_curve,
    "time_resolution": dt,
    "input_frequencies": frequencies,
    "input_psd": power_spectrum_1,
    "mean_magnitude": mean_magnitude,
    "standard_deviation": standard_deviation,
    "seed": random_seed,
}
psd_kwarg_dict_2 = {
    "length_of_light_curve": length_of_light_curve,
    "time_resolution": dt,
    "input_frequencies": frequencies,
    "input_psd": power_spectrum_2,
    "mean_magnitude": mean_magnitude,
    "standard_deviation": standard_deviation,
    "seed": random_seed,
}

psd_variable_signal = variability.Variability(
    variability_model="user_defined_psd", **psd_kwarg_dict_1
)
psd_signal_1 = psd_variable_signal.variability_at_time(other_time_stamps)
psd_variable_signal = variability.Variability(
    variability_model="user_defined_psd", **psd_kwarg_dict_2
)
psd_signal_2 = psd_variable_signal.variability_at_time(other_time_stamps)

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(other_time_stamps, psd_signal_1)
ax.plot(other_time_stamps, psd_signal_2, alpha=0.5)
ax.set_xlabel("time [days]")
ax.set_ylabel("amplitude [arb.]")
plt.show()

Use these as an intrinsic signal to be reprocessed

outer_radius_of_accretion_disk = 1000  # R_g = G M_bh / c^2
corona_height = 30  # R_g
r_resolution = 1000  # pixels per outer radius
inclination_angle = 45  # degrees
black_hole_mass_exponent = 9.0
black_hole_spin = 0.3
eddington_ratio = 0.15

agn_kwarg_dict = {
    "r_out": outer_radius_of_accretion_disk,
    "corona_height": corona_height,
    "r_resolution": r_resolution,
    "inclination_angle": inclination_angle,
    "black_hole_mass_exponent": black_hole_mass_exponent,
    "black_hole_spin": black_hole_spin,
    "eddington_ratio": eddington_ratio,
}

agn_reprocessor = accretion_disk_reprocessing.AccretionDiskReprocessing(
    "lamppost", **agn_kwarg_dict
)

agn_reprocessor.define_intrinsic_signal(
    time_array=time_stamps, magnitude_array=light_curve
)
reprocessed_200nm = agn_reprocessor.reprocess_signal(
    rest_frame_wavelength_in_nanometers=200
)
reprocessed_500nm = agn_reprocessor.reprocess_signal(
    rest_frame_wavelength_in_nanometers=500
)

# Note the earliest times do not have enough signal history built up to generate a reprocessed signal

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_stamps, light_curve, color="black", linewidth=0.5, label="driving signal")
ax.plot(time_stamps, reprocessed_200nm, alpha=0.7, label="200 nm")
ax.plot(time_stamps, reprocessed_500nm, alpha=0.7, label="500 nm")
ax.set_xlabel("time [days]")
ax.set_ylabel("amplitude [arb.]")
ax.set_xlim(50, 200)  # Zoom in on a small region to see delays
ax.set_ylim(9, 11)
ax.legend()
plt.show()

To understand what these reprocessing functions look like, we may plot them

response_function_200nm = agn_reprocessor.define_new_response_function(
    rest_frame_wavelength_in_nanometers=200
)
response_function_500nm = agn_reprocessor.define_new_response_function(
    rest_frame_wavelength_in_nanometers=500
)
# Note the sum of these response function is defined to be 1
# so the convolution doesn't change the amplitudes
print(np.sum(response_function_200nm))

# similar plotting, but changing time axis to units of [days]
grav_rad = astro_util.calculate_gravitational_radius(
    agn_kwarg_dict["black_hole_mass_exponent"]
)
maximum_lag = len(response_function_200nm) * grav_rad / (3e8 * 60 * 60 * 24)
time_lag_axis = np.linspace(0, maximum_lag, len(response_function_200nm))

fig, ax = plt.subplots(2, figsize=(10, 6))
ax[0].plot(response_function_200nm, alpha=0.7, label="200 nm")
ax[0].plot(response_function_500nm, alpha=0.7, label="500 nm")
ax[0].set_xlabel("time delay [R_g]")

ax[1].plot(time_lag_axis, response_function_200nm, alpha=0.7, label="200 nm")
ax[1].plot(time_lag_axis, response_function_500nm, alpha=0.7, label="500 nm")
ax[1].set_xlabel("time delay [days]")

ax[0].legend()
ax[1].legend()
fig.supylabel("response [arb.]")
plt.subplots_adjust(hspace=0.3)
plt.show()

1.0

We can define our own reprocessing kernels using this agn reprocessor as well

time_lags = np.linspace(0, 25, 25)  # days if defined like this

predefined_kernel = np.zeros(np.shape(time_lags))
predefined_kernel[20] = 1  # define a 20 day delay

reprocessed_signal_with_delay = agn_reprocessor.reprocess_signal(
    response_function_time_lags=time_lags,
    response_function_amplitudes=predefined_kernel,
)

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_stamps, light_curve, color="black", linewidth=0.5, label="driving signal")
ax.plot(time_stamps, reprocessed_signal_with_delay, label="delayed 'reprocessing'")
ax.set_xlabel("time [days]")
ax.set_ylabel("amplitude [arb.]")
ax.set_xlim(0, 200)  # Note we do not have a 'reprocessed signal' until time = 20 days.
ax.set_ylim(9, 11)
ax.legend()
plt.show()

Comparing the difference between normal variance in magnitude and amplitude

mean_magnitude = 10
magnitude_variance = -0.4
bpl_kwarg_dict_mag = {
    "length_of_light_curve": length_of_light_curve,
    "time_resolution": dt,
    "log_breakpoint_frequency": log_breakpoint_frequency,
    "low_frequency_slope": low_frequency_index,
    "high_frequency_slope": high_frequency_index,
    "mean_magnitude": mean_magnitude,
    "standard_deviation": magnitude_variance,
    "normal_magnitude_variance": True,
    "seed": random_seed,
}

variable_light_curve = variability.Variability(
    variability_model="bending_power_law", **bpl_kwarg_dict_mag
)

time_axis = np.linspace(0, 999, 1000)
magnitude_light_curve = variable_light_curve.variability_at_time(time_axis)
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_axis, magnitude_light_curve)
ax.plot(
    [time_axis[0], time_axis[-1]], [mean_magnitude, mean_magnitude], "--", color="black"
)
ax.invert_yaxis()
ax.set_xlabel("time [days]")
ax.set_ylabel("magnitude")
plt.show()

Creating a signal with normal flux variance

mean_magnitude = 10
magnitude_variance = 0.4
bpl_kwarg_dict_mag_lognormal = {
    "length_of_light_curve": length_of_light_curve,
    "time_resolution": dt,
    "log_breakpoint_frequency": log_breakpoint_frequency,
    "low_frequency_slope": low_frequency_index,
    "high_frequency_slope": high_frequency_index,
    "mean_magnitude": mean_magnitude,
    "standard_deviation": magnitude_variance,
    "normal_magnitude_variance": False,
    "seed": random_seed,
}

variable_light_curve_lognormal = variability.Variability(
    variability_model="bending_power_law", **bpl_kwarg_dict_mag_lognormal
)

time_axis = np.linspace(0, 999, 1000)
magnitude_light_curve = variable_light_curve.variability_at_time(time_axis)
magnitude_light_curve_lognormal = variable_light_curve_lognormal.variability_at_time(
    time_axis
)
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_axis, magnitude_light_curve, label="normal magnitude variability")
ax.plot(
    time_axis,
    magnitude_light_curve_lognormal,
    label="normal amplitude variability",
    alpha=0.5,
)
ax.plot(
    [time_axis[0], time_axis[-1]],
    [mean_magnitude, mean_magnitude],
    "--",
    color="black",
    label="mean",
)
ax.invert_yaxis()
ax.set_xlabel("time [days]")
ax.set_ylabel("flux [mag.]")
ax.legend()
plt.show()

Note that by converting fluctuations into amplitudes suppresses the brighter fluctuations and accentuates dimmer fluctuations.