import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

t_start, t_end, num_points = 0.0, 40.0, 4001
t1, t_w = 3.0, 10.0
G_h, G_c = 1.0, 0.06
lam_s = 0.3
w1, lam1 = 0.65, 1.2
w2, lam2 = 0.35, 0.06
D0_target = 0.03

t = np.linspace(t_start, t_end, num_points)
dt = t[1]-t[0]
idx_w = np.argmin(np.abs(t - t_w))

G = np.zeros_like(t)
G[(t>=0)&(t<t1)] = G_h
G[t>=t_w] += G_c

# Single
y = np.zeros_like(t)
for k in range(1, len(t)):
    y[k] = y[k-1] + dt * (-lam_s * y[k-1] + 1.0 * G[k-1])
N_before = y[idx_w]
N_eq_s = G_c / lam_s
A_s = (N_eq_s + D0_target) - N_before
y[idx_w] += A_s
for k in range(idx_w+1, len(t)):
    y[k] = y[k-1] + dt * (-lam_s * y[k-1] + 1.0 * G[k-1])
N_s = y.copy()

# Double
y1 = np.zeros_like(t); y2 = np.zeros_like(t)
for k in range(1, len(t)):
    y1[k] = y1[k-1] + dt * (-lam1 * y1[k-1] + w1 * G[k-1])
    y2[k] = y2[k-1] + dt * (-lam2 * y2[k-1] + w2 * G[k-1])
N_before_d = y1[idx_w]+y2[idx_w]
N_eq_d = w1*G_c/lam1 + w2*G_c/lam2
A_d = (N_eq_d + D0_target) - N_before_d
y1[idx_w] += w1*A_d; y2[idx_w] += w2*A_d
for k in range(idx_w+1, len(t)):
    y1[k] = y1[k-1] + dt * (-lam1*y1[k-1] + w1*G[k-1])
    y2[k] = y2[k-1] + dt * (-lam2*y2[k-1] + w2*G[k-1])
N_d = y1+y2

t_rel = t[idx_w:] - t_w
D_s = N_s[idx_w:] - N_eq_s
D_d = N_d[idx_w:] - N_eq_d

plt.figure(figsize=(9,5))
plt.plot(t_rel, D_s, label="Single exponential (monotonic)")
plt.plot(t_rel, D_d, label="Double exponential (hump)")
plt.axhline(0.0, linestyle="--", linewidth=1, label="Zero deviation")
plt.xlabel(r"$t - t_w$")
plt.ylabel(r"$D(t) = N(t) - N_{\mathrm{eq}}$")
plt.title("Closer curves with same initial deviation; hump only in double")
plt.legend(); plt.grid(True); plt.tight_layout()
plt.savefig("kovacs_close_same_D0_fixed_resaved.png", dpi=200)
plt.show()

df = pd.DataFrame({"t_minus_tw": t_rel, "D_single": D_s, "D_double": D_d})
df.to_csv("kovacs_close_same_D0_fixed_resaved.csv", index=False)