import numpy as np

def newmark_sdof(m, c, k, p, dt, u0=0, v0=0, gamma=0.5, beta=0.25):
    """Newmark-beta (varsayilan: ortalama ivme, kosulsuz kararli)"""
    n = len(p)
    u = np.zeros(n); v = np.zeros(n); a = np.zeros(n)
    u[0], v[0] = u0, v0
    a[0] = (p[0] - c*v0 - k*u0) / m

    a1 = m/(beta*dt**2) + gamma*c/(beta*dt)
    a2 = m/(beta*dt) + (gamma/beta - 1)*c
    a3 = (1/(2*beta) - 1)*m + dt*(gamma/(2*beta) - 1)*c
    k_hat = k + a1

    for i in range(n-1):
        p_hat = p[i+1] + a1*u[i] + a2*v[i] + a3*a[i]
        u[i+1] = p_hat / k_hat
        v[i+1] = (gamma/(beta*dt))*(u[i+1]-u[i]) + (1-gamma/beta)*v[i] \
                 + dt*(1-gamma/(2*beta))*a[i]
        a[i+1] = (u[i+1]-u[i])/(beta*dt**2) - v[i]/(beta*dt) \
                 - (1/(2*beta)-1)*a[i]
    return u, v, a

# Bilinen sistem: analitik cozumu TAM olarak bilinen serbest titresim
m, wn, zeta = 1000.0, 10.0, 0.05
k = m*wn**2
c = 2*zeta*wn*m
u0, v0 = 0.01, 0.0
wd = wn*np.sqrt(1 - zeta**2)

print("dt(s)    T/dt    maks hata(%)   hata orani")
prev_err = None
for dt in [0.05, 0.02, 0.01, 0.005, 0.002]:
    T_end = 3.0
    n = int(T_end/dt) + 1
    t = np.linspace(0, T_end, n)
    p = np.zeros(n)                       # serbest titresim: dis kuvvet yok

    u, _, _ = newmark_sdof(m, c, k, p, dt, u0, v0)

    # Analitik cozum
    ua = u0*np.exp(-zeta*wn*t)*(np.cos(wd*t) + (zeta*wn/wd)*np.sin(wd*t))
    err = np.max(np.abs(u - ua))/u0*100

    ratio = f"{prev_err/err:.2f}x" if prev_err else "-"
    print(f"{dt:.3f}   {2*np.pi/wn/dt:6.1f}   {err:9.4f}     {ratio}")
    prev_err = err