<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from this Jupyter notebook by Heiner Igel (@heinerigel), Lion Krischer (@krischer) and Taufiqurrahman (@git-taufiqurrahman) which is a supplemenatry material to the book Computational Seismology
Step1:  Solving the 1D acoustic wave equation by finite-differences
In the previous lectures we derived acoustic and elastic approximations to the partial differential equations of motion, describing seismic wave propagation in 3D isotropic elastic media. For homogeneous acoustic media we derived analytical solutions in terms of Green's functions. Finally, we dealt with the discretization of material parameters and wavefields in continous media and the finite-difference approximation of partial derivatives. 
In this lesson we wrap up the results of all previous lectures by solving the 1D acoustic wave equation for a homogeneous medium using the finite-difference method.
Finite difference solution
As derived in this and this lecture, the acoustic wave equation in 1D with constant density is
\begin{equation}
\frac{\partial^2 p(x,t)}{\partial t^2} \ = \ vp(x)^2 \frac{\partial^2 p(x,t)}{\partial x^2} + s(x,t) \nonumber
\end{equation}
with pressure $p$, acoustic velocity $vp$, and source term $s$. Both second derivatives can be approximated by a 3-point difference formula. For example for the time derivative, we get
Step2:  Source time function
To excitate wave propagation in our 1D homogenous model, we use the first derivative of the Gaussian
Step3:  Analytical Solution
In this lecture we calculated the Green's functions for the homogenous acoustic wave equation
\begin{equation}
\frac{\partial^2}{\partial t^2} G(x,t;x_s, t_s) \ - \ vp_0^2 \Delta G(x,t;x_s, t_s) \ = \delta (x-x_s) \delta (t-t_s) \nonumber
\end{equation}
where $\Delta$ denotes the Laplace operator and the $\delta-$function is defined as
\begin{equation}
\delta(x) = \left{
\begin{array}{ll}
    \infty &x=0 \
    0 &x\neq 0
\end{array}
\right.\notag
\end{equation}
with the constraint
Step4:  Comparison of numerical with analytical solution
In the code below we solve the homogeneous acoustic wave equation by the 3-point difference operator and compare the numerical results with the analytical solution.
<ASSISTANT_TASK:>
Python Code:
# Execute this cell to load the notebook's style sheet, then ignore it
from IPython.core.display import HTML
css_file = '../../style/custom.css'
HTML(open(css_file, "r").read())
Explanation: Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from this Jupyter notebook by Heiner Igel (@heinerigel), Lion Krischer (@krischer) and Taufiqurrahman (@git-taufiqurrahman) which is a supplemenatry material to the book Computational Seismology: A Practical Introduction,  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi
End of explanation
# Import Libraries (PLEASE RUN THIS CODE FIRST!) 
# ----------------------------------------------
import numpy as np
import matplotlib
# Show Plot in The Notebook
matplotlib.use("nbagg")
import matplotlib.pyplot as plt
# Sub-plot Configuration
# ----------------------
from matplotlib import gridspec 
# Ignore Warning Messages
# -----------------------
import warnings
warnings.filterwarnings("ignore")
# Definition of modelling parameters
# ----------------------------------
         # number of grid points in x-direction
         # grid point distance in x-direction
         # wave speed in medium (m/s)
         # source location in grid in x-direction
         # receiver location in grid in x-direction
         # maximum number of time steps
         # time step
Explanation: Solving the 1D acoustic wave equation by finite-differences
In the previous lectures we derived acoustic and elastic approximations to the partial differential equations of motion, describing seismic wave propagation in 3D isotropic elastic media. For homogeneous acoustic media we derived analytical solutions in terms of Green's functions. Finally, we dealt with the discretization of material parameters and wavefields in continous media and the finite-difference approximation of partial derivatives. 
In this lesson we wrap up the results of all previous lectures by solving the 1D acoustic wave equation for a homogeneous medium using the finite-difference method.
Finite difference solution
As derived in this and this lecture, the acoustic wave equation in 1D with constant density is
\begin{equation}
\frac{\partial^2 p(x,t)}{\partial t^2} \ = \ vp(x)^2 \frac{\partial^2 p(x,t)}{\partial x^2} + s(x,t) \nonumber
\end{equation}
with pressure $p$, acoustic velocity $vp$, and source term $s$. Both second derivatives can be approximated by a 3-point difference formula. For example for the time derivative, we get:
\begin{equation}
\frac{\partial^2 p(x,t)}{\partial t^2} \ \approx \ \frac{p(x,t+dt) - 2 p(x,t) + p(x,t-dt)}{dt^2}, \nonumber
\end{equation}
and equivalently for the spatial derivative: 
\begin{equation}
\frac{\partial^2 p(x,t)}{\partial x^2} \ \approx \ \frac{p(x+dx,t) - 2 p(x,t) + p(x-dx,t)}{dx^2}, \nonumber
\end{equation}
Injecting these approximations into the wave equation allows us to formulate the pressure p(x) for the time step $t+dt$ (the future) as a function of the pressure at time $t$ (now) and $t-dt$ (the past). This is called an explicit scheme allowing the $extrapolation$  of the space-dependent field into the future only looking at the nearest neighbourhood.
After discretization of the P-wave velocity and pressure wavefield at the discrete spatial grid points $i = 0, 1, 2, ..., nx$ and time steps $n = 0, 1, 2, ..., nt$, we can replace the time-dependent part (upper index time, lower index space) by
\begin{equation}
 \frac{p_{i}^{n+1} - 2 p_{i}^n + p_{i}^{n-1}}{\mathrm{d}t^2} \ = \ vp_{i}^2 \biggl( \frac{\partial^2 p}{\partial x^2}\biggr) \ + s_{i}^n \nonumber
\end{equation}
Solving for $p_{i}^{n+1}$ leads to the extrapolation scheme:
\begin{equation}
p_{i}^{n+1} \ = \ vp_i^2 \mathrm{d}t^2 \left( \frac{\partial^2 p}{\partial x^2} \right) + 2p_{i}^n - p_{i}^{n-1} + \mathrm{d}t^2 s_{i}^n
\end{equation}
The spatial derivatives are determined by 
\begin{equation}
\frac{\partial^2 p(x,t)}{\partial x^2} \ \approx \ \frac{p_{i+1}^{n} - 2 p_{i}^n + p_{i-1}^{n}}{\mathrm{d}x^2} \nonumber
\end{equation}
Eq. (1) is the essential core of the FD modelling code. Because we derived analytical solutions for wave propagation in a homogeneous medium, we should test our first code implementation for a similar medium, by setting:
\begin{equation}
vp_i = vp0 \notag
\end{equation}
at each spatial grid point $i = 0, 1, 2, ..., nx$, in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The initial condition is 
\begin{equation}
p_{i}^0 = 0, \nonumber
\end{equation}
so the modelling starts with zero pressure amplitude at each spatial grid point $i = 0, 1, 2, ..., nx$. As boundary conditions, we assume 
\begin{align}
p_{0}^n = 0, \nonumber\
p_{nx}^n = 0, \nonumber\
\end{align}
for all time steps n. This Dirichlet boundary condition, leads to artifical boundary reflections which would obviously not describe a homogeneous medium. For now, we simply extend the model, so that boundary reflections are not recorded at the receiver positions.
Let's implement it ...
End of explanation
# Plot Source Time Function 
# -------------------------
 # dominant frequency of the source (Hz)
 # source time shift
print('Source frequency =', f0, 'Hz')
# Source time function (Gaussian)
# -------------------------------
# 1st derivative of a Gaussian
src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))
# Plot position configuration
# ---------------------------
plt.ion()
fig1 = plt.figure(figsize=(6, 5))
gs1  = gridspec.GridSpec(1, 2, width_ratios=[1, 1], hspace=0.3, wspace=0.3)
# Plot source time function
# -------------------------
ax1  = plt.subplot(gs1[0])
ax1.plot(time, src) # plot source time function
ax1.set_title('Source Time Function')
ax1.set_xlim(time[0], time[-1])
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Amplitude')
# Plot source spectrum
# --------------------
ax2  = plt.subplot(gs1[1])
      # source time function in frequency domain
      # time domain to frequency domain
      # plot frequency and amplitude
      # only display frequency from 0 to 250 Hz
ax2.set_title('Source Spectrum')
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Amplitude')
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
plt.show()
Explanation: Source time function
To excitate wave propagation in our 1D homogenous model, we use the first derivative of the Gaussian:
\begin{equation}
s(t) = -2 (t-t_0) f_0^2 exp(-f_0^2 (t-t_0)^2) \nonumber
\end{equation}
as source time function at the discrete source position $isrc$. Where $t_0$ denotes a time shift and $f_0$ the dominant frequency of the source.
End of explanation
# Analytical solution
# -------------------
G    = time * 0.
# Initialize coordinates
# ----------------------
x    = np.arange(nx)
x    = x * dx       # coordinate in x-direction
for it in range(nt): # Calculate Green's function (Heaviside function)
    if (time[it] - np.abs(x[ir] - x[isrc]) / vp0) >= 0:
        G[it] = 1. / (2 * vp0)
        
    # get limit value from the maximum amplitude
# Plotting convolution of Green's function with source wavelet
plt.plot(time, Gc)
plt.title("Analytical solution" )
plt.xlabel("Time, s")
plt.ylabel("Amplitude")
plt.grid()
plt.show()
Explanation: Analytical Solution
In this lecture we calculated the Green's functions for the homogenous acoustic wave equation
\begin{equation}
\frac{\partial^2}{\partial t^2} G(x,t;x_s, t_s) \ - \ vp_0^2 \Delta G(x,t;x_s, t_s) \ = \delta (x-x_s) \delta (t-t_s) \nonumber
\end{equation}
where $\Delta$ denotes the Laplace operator and the $\delta-$function is defined as
\begin{equation}
\delta(x) = \left{
\begin{array}{ll}
    \infty &x=0 \
    0 &x\neq 0
\end{array}
\right.\notag
\end{equation}
with the constraint:
\begin{equation}
\int_{-\infty}^{\infty} \delta(x)\; dx = 1.\notag
\end{equation}
When comparing numerical with analytical solutions the functions that - in the limit - lead to the $\delta-$function will become very important. An example is the boxcar function
$$
\delta_{bc}(x) = \left{
\begin{array}{ll}
    1/dx &|x|\leq dx/2 \
    0 &\text{elsewhere}   
\end{array}
\right.
$$
fulfilling these properties as $dx\rightarrow0$. These functions are used to properly scale the source terms to obtain correct absolute amplitudes.
To describe analytical solutions for the acoustic wave equation we also make use of the unit step function, also known as the Heaviside function, defined as 
$$
H(x) = \left{
\begin{array}{ll}
    0 &x<0 \
    1 &x \geq 0   
\end{array}
\right.
$$
The Heaviside function is the integral of the $\delta-$function (and vice-versa the $\delta$-function is defined as the derivative of the Heaviside function). In 1D case, the Greens function is proportional to a Heaviside function. 
$$
G=\frac{1}{2vp_0}H\left(t-\frac{|x|}{vp_0}\right)
$$
As the response to an arbitrary source time function can be obtained by convolution 
$$
G_{seis} = G(x,t;x_s,t_s) * s(t),
$$
this implies that the propagating waveform is the integral of the source time function. The response is shown for a source time function with a 1st derivative of a Gaussian.
End of explanation
# Plot Snapshot & Seismogram (PLEASE RERUN THIS CODE AGAIN AFTER SIMULATION!) 
# ---------------------------------------------------------------------------
# Initialize empty pressure arrays
# --------------------------------
 # p at time n (now)
 # p at time n-1 (past)
 # p at time n+1 (present)
 # 2nd spatial derivative of p
# Initialize model (assume homogeneous model)
# -------------------------------------------
      # initialize wave velocity in model
# Initialize empty seismogram
# ---------------------------
 
# Plot position configuration
# ---------------------------
plt.ion()
fig2  = plt.figure(figsize=(6, 4))
gs2   = gridspec.GridSpec(1, 2, width_ratios=[1, 1], hspace=0.3, wspace=0.3)
# Plot 1D wave propagation
# ------------------------
# Note: comma is needed to update the variable
ax3  = plt.subplot(gs2[0])
leg1,= ax3.plot(isrc, 0, 'r*', markersize=11) # plot position of the source in snapshot
leg2,= ax3.plot(ir, 0, 'k^', markersize=8) # plot position of the receiver in snapshot
up31,= ax3.plot(p) # plot pressure update each time step
ax3.set_xlim(0, nx)
ax3.set_ylim(-lim, lim)
ax3.set_title('Time Step (nt) = 0')
ax3.set_xlabel('nx')
ax3.set_ylabel('Amplitude')
ax3.legend((leg1, leg2), ('Source', 'Receiver'), loc='upper right', fontsize=10, numpoints=1)
# Plot seismogram 
# ---------------
# Note: comma is needed to update the variable
ax4  = plt.subplot(gs2[1])
leg3,= ax4.plot(0,0,'r--',markersize=1) # plot analytical solution marker
leg4,= ax4.plot(0,0,'b-',markersize=1) # plot numerical solution marker
up41,= ax4.plot(time, seis) # update recorded seismogram each time step
up42,= ax4.plot([0], [0], 'r|', markersize=15) # update time step position
ax4.yaxis.tick_right()
ax4.yaxis.set_label_position("right")
ax4.set_xlim(time[0], time[-1])
ax4.set_title('Seismogram')
ax4.set_xlabel('Time (s)')
ax4.set_ylabel('Amplitude')
ax4.legend((leg3, leg4), ('Analytical', 'FD'), loc='upper right', fontsize=10, numpoints=1)
plt.plot(time,Gc,'r--') # plot analytical solution
plt.show()
# 1D Wave Propagation (Finite Difference Solution) 
# ------------------------------------------------
# Calculate Partial Derivatives
# -----------------------------
for it in range(nt):
    
    # FD approximation of spatial derivative by 3 point operator
    # ----------------------------------------------------------
    # Explicit Time Extrapolation by 3 point operator
    # -----------------------------------------------
    # Add Source Term at isrc
    # -----------------------
    # Absolute pressure w.r.t analytical solution
                
    # Remap Time Levels
    # -----------------
    
    # Output Seismogram
    # -----------------
    
    # Update Data for Wave Propagation Plot
    # -------------------------------------
    idisp = 2 # display frequency
    if (it % idisp) == 0:
        ax3.set_title('Time Step (nt) = %d' % it)
        up31.set_ydata(p)
        up41.set_ydata(seis)
        up42.set_data(time[it], seis[it])
        plt.gcf().canvas.draw()
Explanation: Comparison of numerical with analytical solution
In the code below we solve the homogeneous acoustic wave equation by the 3-point difference operator and compare the numerical results with the analytical solution.
End of explanation
<END_TASK>