import math
import random
import wave
import numpy as np

sample_rate = 48000
duration = 20.0
total_samples = int(duration * sample_rate)
t = np.linspace(0, duration, total_samples, endpoint=False)

np.random.seed(42)
random.seed(42)

print('1. Synthesizing Flowing Water Stream...')
# A. Base water turbulence (Brownian / Pink noise)
white = np.random.normal(0, 1, total_samples)
# Filter white noise into water turbulence: integrate for Brownian
brownian = np.cumsum(white)
brownian = brownian - np.mean(brownian)
brownian = brownian / np.max(np.abs(brownian))

# Simple multi-band resonance for water body
# Fast convolution or IIR via difference equation
b_water_l = np.zeros(total_samples)
b_water_r = np.zeros(total_samples)

# Multi-rate filter for water body
lp1 = 0.0
lp2 = 0.0
lp_slow_l = np.zeros(total_samples)
lp_slow_r = np.zeros(total_samples)

white_l = np.random.normal(0, 1, total_samples)
white_r = np.random.normal(0, 1, total_samples)

# B. Realistic bubbling stream (thousands of micro bubble grains)
num_bubbles = 3500
bubble_audio_l = np.zeros(total_samples)
bubble_audio_r = np.zeros(total_samples)

bubble_times = np.random.uniform(0.1, duration - 0.2, num_bubbles)
bubble_freqs = np.random.uniform(500, 2400, num_bubbles)
bubble_durs = np.random.uniform(0.012, 0.038, num_bubbles)
bubble_pans = np.random.uniform(0.15, 0.85, num_bubbles)
bubble_amps = np.random.uniform(0.015, 0.065, num_bubbles)

for i in range(num_bubbles):
    t_start = bubble_times[i]
    dur = bubble_durs[i]
    f0 = bubble_freqs[i]
    amp = bubble_amps[i]
    pan = bubble_pans[i]
    
    idx_start = int(t_start * sample_rate)
    n_samples = int(dur * sample_rate)
    if idx_start + n_samples >= total_samples:
        continue
        
    t_b = np.linspace(0, dur, n_samples, endpoint=False)
    # Frequency rises slightly during bubble life (Minnaert pitch rise)
    freq_curve = f0 * (1.0 + 0.15 * (t_b / dur))
    # Damped exponential decay
    env = np.exp(-t_b * (50.0 / dur + 20.0))
    wave_b = amp * np.sin(2.0 * np.pi * freq_curve * t_b) * env
    
    bubble_audio_l[idx_start:idx_start+n_samples] += wave_b * (1.0 - pan)
    bubble_audio_r[idx_start:idx_start+n_samples] += wave_b * pan

# C. Water lapping / wave modulations (swells of water flow)
swell_mod = 0.70 + 0.30 * np.sin(2.0 * np.pi * 0.22 * t + 0.4)
stream_l = (bubble_audio_l * 1.5 + brownian * 0.08) * swell_mod
stream_r = (bubble_audio_r * 1.5 + brownian * 0.08) * swell_mod

# D. Soft droplets (cánh hoa rơi chạm mặt nước)
for t_drop in [2.4, 5.6, 9.1, 13.8, 17.5]:
    idx = int(t_drop * sample_rate)
    dur_d = 0.08
    n_d = int(dur_d * sample_rate)
    t_d = np.linspace(0, dur_d, n_d, endpoint=False)
    f_drop = 1650.0
    wave_drop = 0.08 * np.sin(2.0 * np.pi * f_drop * (1.0 + 0.3 * (t_d/dur_d)) * t_d) * np.exp(-t_d * 45.0)
    stream_l[idx:idx+n_d] += wave_drop * 0.6
    stream_r[idx:idx+n_d] += wave_drop * 0.4

# Normalize stream
max_val = max(np.max(np.abs(stream_l)), np.max(np.abs(stream_r)))
if max_val > 0:
    stream_l = stream_l / max_val * 0.75
    stream_r = stream_r / max_val * 0.75

print('Flowing water synthesized successfully!')