# -*- coding: utf-8 -*-
"""
Created on Tue Oct 15 16:46:07 2024

@author: ernesto311
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
from matplotlib.animation import FuncAnimation, PillowWriter
import useful_functions as f

# Generate distances and corresponding frequencies
distances = np.arange(-400, 2000, 1)
f1 = 800
f2 = f1 + distances


# Time parameters
dt = 0.00001
t = np.arange(0, 2, dt)

# Set up the figure and axis for the animation
fig, axs = plt.subplots(1, 2, layout='constrained')
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Amplitude')
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power [1/Hz]')
axs[1].set_xlim(0, 1500)



 # Generate the signal as a sum of two sine waves
def plot_powerspectrum_2(i):
    # Clear the previous plots
    axs[0].cla()
    axs[1].cla()
    
    # Generate the signal as a sum of two sine waves
    x = np.sin(2 * np.pi * f1 * t) + 0.8 * np.sin(2 * np.pi * f2[i] * t)  # Second wave is 20% as strong
    
    
    # Plot the signal (first 20 ms for clarity)
    axs[0].plot(t[t < 0.02], x[t < 0.02])
    axs[0].set_title(f"Signal (f2={f2[i]} Hz)")
    axs[0].set_xlabel('Time [s]')
    axs[0].set_ylabel('Amplitude')
    axs[0].set_ylim(-2, 2)

    x[x < 0] = 0  # Apply half-wave rectification (optional)
    # Compute power spectrum
    freq, power = welch(x, fs=1/dt, nperseg=2**16)
    pref = np.max(power)
    decibel_power = 10 * np.log10(power/pref)
    AM = f.find_AM(f1, 0.5 * f1, f2[i])
    # Plot the power spectrum
    axs[1].plot(freq, power)
    axs[1].set_xlim(0, 3000)
    axs[1].set_title(f'Power Spectrum (f2={f2[i]} Hz)')
    axs[1].set_xlabel('Frequency [Hz]')
    axs[1].set_ylabel('Power [1/Hz]')
    #axs[1].set_ylim(0, 0.00007)
    axs[1].plot(f1, power[np.argmin(np.abs(freq-f1))], 'o')
    axs[1].plot(f2[i], power[np.argmin(np.abs(freq-f2[i]))], 'd')
    axs[1].plot(AM, power[np.argmin(np.abs(freq-AM))], '*')
    axs[1].axvline(AM, alpha = 0.5, color = 'r')

# Create the animation
ani = FuncAnimation(fig, plot_powerspectrum_2, frames=len(distances), interval=500)

# Save the animation as a GIF file (optional)
ani.save("sum_of_sinewaves.gif", writer=PillowWriter(fps=30))