Data Analysis Functions

import numpy as np

# Load data from space/tab-separated file with header
datos = np.loadtxt("posicion_vs_tiempo.txt", skiprows=1)

t = datos[:, 0]  # Time column
x = datos[:, 1]  # Position X
y = datos[:, 2]  # Position Y

times = []
with open("mediciones_80cm.txt", "r") as f:
    for line in f:
        line = line.strip()
        if not line:
            continue  # Skip empty lines
        try:
            # Split by comma and take second column
            t_str = line.split(",")[1].strip()
            times.append(float(t_str))
        except IndexError:
            print("Línea malformada:", line)
        except ValueError:
            print("Valor no numérico:", line)

times = np.array(times)

import numpy as np

# Dataset from experimental measurements
N = len(times)

# Mean value
t_mean = np.mean(times)

# Standard deviation (sample)
t_std = np.std(times, ddof=1)

# Standard error of the mean
t_err = t_std / np.sqrt(N)

print(f"Tiempo medio: {t_mean:.5f} s")
print(f"Desviación estándar: {t_std:.5f} s")
print(f"Error de la media: {t_err:.5f} s")

# Calculate g from free fall: g = 2s / t²
s = 0.80  # distance in meters
g_values = 2 * s / times**2
g_mean = np.mean(g_values)

# Propagate uncertainty: σ_g = g · 2 · (σ_t / t)
g_err = g_mean * 2 * (t_err / t_mean)

print(f"g promedio: {g_mean:.3f} m/s^2")
print(f"Incertidumbre en g: ±{g_err:.3f} m/s^2")

def smooth(arr, w=5):
    """Apply moving average with window size w."""
    kernel = np.ones(w) / w
    return np.convolve(arr, kernel, mode='same')

# Smooth position data
y_smoothed = smooth(y_data, w=5)
x_smoothed = smooth(x_data, w=5)

import numpy as np

# First derivative: velocity from position
vx = np.gradient(x_smoothed, t)
vy = np.gradient(y_smoothed, t)

# Apply additional smoothing to velocity
vx = smooth(vx, w=5)
vy = smooth(vy, w=5)

# Second derivative: acceleration from velocity
ax = np.gradient(vx, t)
ay = np.gradient(vy, t)

# Smooth acceleration
ax = smooth(ax, w=5)
ay = smooth(ay, w=5)

# Total velocity magnitude
v_total = np.sqrt(vx**2 + vy**2)

import numpy as np

# Quadratic fit for uniformly accelerated motion
# Position: s = s₀ + v₀t + ½at²
coef = np.polyfit(t, x, 2)

# Extract acceleration from coefficient
a = 2 * coef[0]  # Coefficient of t²

# Generate fitted curve
x_fit = np.polyval(coef, t)

print(f"Aceleración estimada: {a:.4f} m/s^2")

from scipy.optimize import curve_fit

def sinusoide_amortiguada(t, A, omega, phi, gamma, offset):
    """Damped sinusoid: A·exp(-γt)·cos(ωt + φ) + offset"""
    return A * np.exp(-gamma * t) * np.cos(omega * t + phi) + offset

# Initial parameter estimates
A0     = (np.max(y_data) - np.min(y_data)) / 2
omega0 = 2 * np.pi / T_estimated
p0     = [A0, omega0, 0.0, 0.05, 0.0]

# Parameter bounds
bounds = ([0, 0.5, -np.pi, 0, -0.5],
          [2,  50, np.pi,  5,  0.5])

# Perform fit
popt, pcov = curve_fit(sinusoide_amortiguada, t_raw, y_data,
                       p0=p0, bounds=bounds, maxfev=20000)

A_fit, omega_fit, phi_fit, gamma_fit, off_fit = popt
perr = np.sqrt(np.diag(pcov))  # Parameter uncertainties

# Extract physical parameters
T_fit = 2 * np.pi / omega_fit  # Period
f_fit = 1.0 / T_fit            # Frequency
omega0 = np.sqrt(omega_fit**2 + gamma_fit**2)  # Natural frequency
k_fit = masa_kg * omega0**2    # Spring constant

from scipy.optimize import curve_fit

def tray_x(t, x0, vx0):
    """Horizontal motion: x = x₀ + vx₀·t"""
    return x0 + vx0 * t

def tray_y(t, y0, vy0, g):
    """Vertical motion: y = y₀ + vy₀·t - ½g·t²"""
    return y0 + vy0 * t - 0.5 * g * t**2

# Fit X component (uniform motion)
popt_x, pcov_x = curve_fit(tray_x, t_raw, x_raw, p0=[x_raw[0], vx0])
x0_fit, vx0_fit = popt_x

# Fit Y component (uniformly accelerated motion)
popt_y, pcov_y = curve_fit(tray_y, t_raw, y_raw,
                            p0=[y_raw[0], vy0, 9.8],
                            bounds=([-np.inf, -50, 0], [np.inf, 50, 20]))
y0_fit, vy0_fit, g_fit = popt_y

# Calculate initial velocity and angle
v0_fit = np.sqrt(vx0_fit**2 + vy0_fit**2)
angulo_fit = np.degrees(np.arctan2(vy0_fit, vx0_fit))

# Center signal around zero
y_centrado = y_data - np.mean(y_data)

# Find zero crossings
cruces = np.where(np.diff(np.sign(y_centrado)))[0]

if len(cruces) >= 2:
    # Average time between crossings
    T_est = 2 * np.mean(np.diff(t_raw[cruces]))
    omega_est = 2 * np.pi / T_est if T_est > 0 else 5.0
else:
    omega_est = 5.0

from scipy.signal import welch, find_peaks

# Calculate PSD using Welch's method
f_psd, Pxx = welch(y_data, fs=fps, nperseg=min(256, len(y_data)//2))

# Find dominant frequency
peaks_idx, _ = find_peaks(Pxx, height=np.max(Pxx)*0.1)

if len(peaks_idx) > 0:
    f_peak = f_psd[peaks_idx[np.argmax(Pxx[peaks_idx])]]
    print(f"Frecuencia dominante: {f_peak:.3f} Hz")

def _estimate_freq(self, data: np.ndarray):
    """Simple frequency detection via ascending zero crossings."""
    if len(data) < 10:
        return
    
    # Calculate signal midpoint
    mid = (data.max() + data.min()) / 2
    above = data > mid
    
    # Find ascending crossings
    crossings = np.where(np.diff(above.astype(np.int8)) == 1)[0]
    
    if len(crossings) >= 2:
        periods = np.diff(crossings)
        avg_period_samples = periods.mean()
        freq = SAMPLE_RATE / avg_period_samples
        self.lbl_freq.setText(f"{freq:.1f} Hz")
    else:
        self.lbl_freq.setText("< detectable")

import matplotlib.pyplot as plt

plt.style.use('dark_background')

fig, axs = plt.subplots(3, 1, figsize=(13, 10), sharex=True)
fig.suptitle("Sistema Masa-Resorte — Cinemática", 
             fontsize=14, fontweight='bold')

# Position subplot
axs[0].scatter(t_raw, y_data*100, s=4, alpha=0.4, 
               color='#00BFFF', label='Datos')
axs[0].plot(t_raw, y_s*100, lw=1.5, color='#00BFFF', 
            label='Suavizado')
axs[0].set_ylabel("y (cm)")
axs[0].legend(fontsize=8)
axs[0].grid(alpha=0.25)

# Velocity subplot
axs[1].plot(t_raw, vy*100, lw=1.8, color='#7FFF00')
axs[1].axhline(0, color='gray', lw=0.6, linestyle=':')
axs[1].set_ylabel("vy (cm/s)")
axs[1].grid(alpha=0.25)

# Acceleration subplot
axs[2].plot(t_raw, ay*100, lw=1.8, color='#FF6347')
axs[2].set_ylabel("ay (cm/s²)")
axs[2].set_xlabel("t (s)")
axs[2].grid(alpha=0.25)

fig.tight_layout()
fig.savefig("cinematica.png", dpi=150, bbox_inches='tight')

fig, ax = plt.subplots(figsize=(7, 7))
fig.suptitle("Espacio de Fases (y, vy)", fontsize=14, fontweight='bold')

# Scatter plot colored by time
sc = ax.scatter(y_data*100, vy*100, c=t_raw, cmap='plasma', 
                s=10, alpha=0.8)
plt.colorbar(sc, ax=ax, label='t (s)')

ax.set_xlabel("Desplazamiento y (cm)")
ax.set_ylabel("Velocidad vy (cm/s)")
ax.axhline(0, color='gray', lw=0.5, linestyle=':')
ax.axvline(0, color='gray', lw=0.5, linestyle=':')
ax.grid(alpha=0.25)

fig.savefig("espacio_fases.png", dpi=150, bbox_inches='tight')

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
fig.suptitle("Diagrama de Bode  H(s)=1/(ms²+cs+k)", 
             fontsize=13, fontweight='bold')

# Calculate transfer function
omega_range = np.logspace(-1, 2, 1000)  # rad/s
s_vals = 1j * omega_range
H_vals = 1.0 / (masa_kg * s_vals**2 + c_fit * s_vals + k_fit)

# Magnitude plot
H_mag_dB = 20 * np.log10(np.abs(H_vals))
ax1.semilogx(omega_range, H_mag_dB, color='#00BFFF', lw=2)
ax1.axvline(omega0, color='white', lw=1.2, linestyle='--', 
            label=f'ω₀={omega0:.3f} rad/s')
ax1.set_ylabel("|H(jω)| (dB)")
ax1.legend(fontsize=9)
ax1.grid(alpha=0.25, which='both')

# Phase plot
H_fase_deg = np.degrees(np.angle(H_vals))
ax2.semilogx(omega_range, H_fase_deg, color='#FF6347', lw=2)
ax2.set_ylabel("∠H(jω) (°)")
ax2.set_xlabel("ω (rad/s)")
ax2.grid(alpha=0.25, which='both')

fig.savefig("bode.png", dpi=150, bbox_inches='tight')

plt.figure(figsize=(8, 8))

# Time distribution
plt.subplot(2, 1, 1)
plt.hist(times, bins=10)
plt.xlabel("Tiempo de caída (s)")
plt.ylabel("Frecuencia")
plt.title("Distribución de tiempos de caída")

# g value distribution
plt.subplot(2, 1, 2)
plt.hist(g_values, bins=10)
plt.xlabel("Aceleración g (m/s²)")
plt.ylabel("Frecuencia")
plt.title("Distribución de valores de g")

plt.tight_layout()
plt.show()

# Save data with descriptive header
header = (
    "t(s)  x(m)  y(m)  y_desp(m)  vx(m/s)  vy(m/s)  ax(m/s2)  ay(m/s2)\n"
    f"# k={k_fit:.4f} N/m  c={c_fit:.6f} N·s/m  "
    f"omega0={omega0:.4f} rad/s  zeta={zeta:.4f}"
)

np.savetxt("datos_masa_resorte.txt",
           np.column_stack((t_raw, x_data, y_raw, y_data, 
                           vx, vy, ax, ay)),
           header=header)

np.savetxt("posicion_vs_tiempo.txt", datos,
           header="t(s) x(m) y(m)")

# Spring-mass-damper system parameters
masa_kg = 0.150  # kg

# From fitted sinusoid
omega_fit = popt[1]    # Damped angular frequency
gamma_fit = popt[3]    # Damping coefficient

# Natural frequency: ω₀ = sqrt(ω² + γ²)
omega0 = np.sqrt(omega_fit**2 + gamma_fit**2)

# Spring constant: k = m·ω₀²
k_fit = masa_kg * omega0**2

# Damping coefficient: c = 2·m·γ
c_fit = 2 * masa_kg * gamma_fit

# Damping ratio: ζ = γ/ω₀
zeta = gamma_fit / omega0

# Period and frequency
T_fit = 2 * np.pi / omega_fit
f_fit = 1.0 / T_fit

print(f"Constante resorte k: {k_fit:.4f} N/m")
print(f"Coef. amort. c: {c_fit:.6f} N·s/m")
print(f"Razón amort. ζ: {zeta:.4f}")

if len(datos) == 0:
    print("No se detectaron datos. Revisa los filtros.")
    exit()

import warnings
warnings.filterwarnings('ignore')

# Good: constrained parameters
bounds = ([0, 0.5, -np.pi, 0, -0.5],  # lower bounds
          [2,  50,  np.pi, 5,  0.5])  # upper bounds

popt, pcov = curve_fit(model, x, y, p0=p0, bounds=bounds)