import sfs
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import HTML

"""
Explanation: Animations of a Pulsating Sphere
End of explanation
"""


import animations_pulsating_sphere as animation

# Pulsating sphere
center = [0, 0, 0]
radius = 0.25
amplitude = 0.05
f = 1000  # frequency
omega = 2 * np.pi * f  # angular frequency

# Axis limits
figsize = (6, 6)
xmin, xmax = -1, 1
ymin, ymax = -1, 1

# Animations
frames = 20  # frames per period

"""
Explanation: In this example, the sound field of a pulsating sphere is visualized.
Different acoustic variables, such as sound pressure,
particle velocity, and particle displacement, are simulated.
The first two quantities are computed with

sfs.fd.source.pulsating_sphere() and 
sfs.fd.source.pulsating_sphere_velocity()

while the last one can be obtained by using

sfs.fd.displacement()

which converts the particle velocity into displacement.
A couple of additional functions are implemented in

animations_pulsating_sphere.py

in order to help creating animating pictures, which is fun!
End of explanation
"""


grid = sfs.util.xyz_grid([xmin, xmax], [ymin, ymax], 0, spacing=0.025)
ani = animation.particle_displacement(
        omega, center, radius, amplitude, grid, frames, figsize, c='Gray')
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Particle Displacement
End of explanation
"""


def hex_grid(xlim, ylim, hex_edge, align='horizontal'):
    if align is 'vertical':
        umin, umax = ylim
        vmin, vmax = xlim
    else:
        umin, umax = xlim
        vmin, vmax = ylim
    du = np.sqrt(3) * hex_edge
    dv = 1.5 * hex_edge
    num_u = int((umax - umin) / du)
    num_v = int((vmax - vmin) / dv)
    u, v = np.meshgrid(np.linspace(umin, umax, num_u),
                       np.linspace(vmin, vmax, num_v))
    u[::2] += 0.5 * du

    if align is 'vertical':
        grid = v, u, 0
    elif align is 'horizontal':
        grid = u, v, 0
    return  grid

grid = hex_grid([xmin, xmax], [ymin, ymax], 0.0125, 'vertical')
ani = animation.particle_displacement(
        omega, center, radius, amplitude, grid, frames, figsize, c='Gray')
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Click the arrow button to start the animation.
to_jshtml() allows you to play with the animation,
e.g. speed up/down the animation (+/- button).
Try to reverse the playback by clicking the left arrow.
You'll see a sound sink.
You can also show the animation by using to_html5_video().
See the documentation for more detail.
Of course, different types of grid can be chosen.
Below is the particle animation using the same parameters
but with a hexagonal grid.
End of explanation
"""


grid = [np.random.uniform(xmin, xmax, 4000),
        np.random.uniform(ymin, ymax, 4000), 0]
ani = animation.particle_displacement(
        omega, center, radius, amplitude, grid, frames, figsize, c='Gray')
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Another one using a random grid.
End of explanation
"""


amplitude = 1e-3
grid = sfs.util.xyz_grid([xmin, xmax], [ymin, ymax], 0, spacing=0.04)
ani = animation.particle_velocity(
        omega, center, radius, amplitude, grid, frames, figsize)
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Each grid has its strengths and weaknesses. Please refer to the
on-line discussion.
Particle Velocity
End of explanation
"""


amplitude = 0.05
impedance_pw = sfs.default.rho0 * sfs.default.c
max_pressure = omega * impedance_pw * amplitude

grid = sfs.util.xyz_grid([xmin, xmax], [ymin, ymax], 0, spacing=0.005)
ani = animation.sound_pressure(
        omega, center, radius, amplitude, grid, frames, pulsate=True,
        figsize=figsize, vmin=-max_pressure, vmax=max_pressure)
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Please notice that the amplitude of the pulsating motion is adjusted
so that the arrows are neither too short nor too long.
This kind of compromise is inevitable since
$$
\text{(particle velocity)} = \text{i} \omega \times (\text{amplitude}),
$$
thus the absolute value of particle velocity is usually
much larger than that of amplitude.
It should be also kept in mind that the hole in the middle
does not visualizes the exact motion of the pulsating sphere.
According to the above equation, the actual amplitude should be
much smaller than the arrow lengths.
The changing rate of its size is also two times higher than the original frequency.
Sound Pressure
End of explanation
"""


amplitude = 1e-6
impedance_pw = sfs.default.rho0 * sfs.default.c
max_pressure = omega * impedance_pw * amplitude

grid = sfs.util.xyz_grid([xmin, xmax], [ymin, ymax], 0, spacing=0.005)
ani = animation.sound_pressure(
        omega, center, radius, amplitude, grid, frames, pulsate=True,
        figsize=figsize, vmin=-max_pressure, vmax=max_pressure)
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Notice that the sound pressure exceeds
the atmospheric pressure ($\approx 10^5$ Pa), which of course makes no sense.
This is due to the large amplitude (50 mm) of the pulsating motion.
It was chosen to better visualize the particle movements
in the earlier animations.
For 1 kHz, the amplitude corresponding to a moderate sound pressure,
let say 1 Pa, is in the order of micrometer.
As it is very small compared to the corresponding wavelength (0.343 m),
the movement of the particles and the spatial structure of the sound field
cannot be observed simultaneously.
Furthermore, at high frequencies, the sound pressure
for a given particle displacement scales with the frequency.
The smaller wavelength (higher frequency) we choose,
it is more likely to end up with a prohibitively high sound pressure.
In the following examples, the amplitude is set to a realistic value 1 $\mu$m.
Notice that the pulsating motion of the sphere is no more visible.
End of explanation
"""


L = 10 * amplitude
xmin_zoom, xmax_zoom = radius - L, radius + L
ymin_zoom, ymax_zoom = -L, L

grid = sfs.util.xyz_grid([xmin_zoom, xmax_zoom], [ymin_zoom, ymax_zoom], 0, spacing=L / 100)
ani = animation.sound_pressure(
        omega, center, radius, amplitude, grid, frames, pulsate=True,
        figsize=figsize, vmin=-max_pressure, vmax=max_pressure)
plt.close()
HTML(ani.to_jshtml())

"""
Explanation: Let's zoom in closer to the boundary of the sphere.
End of explanation
"""


