import numpy as np
import matplotlib.pyplot as plt
import math
c=3e8
l0 = 1e-9
t0 = l0/c
tc = 4*t0
sigma = l0
w = 1.0*((2*math.pi)/t0)
dz = l0/10
dt = dz/(2*c)
Nt=100
fuente=np.zeros(Nt+1)
t =np.zeros(Nt+1)
Nz=100
z=np.zeros(Nz+1)
E_n=np.zeros(Nz+1)
E_np1=np.zeros(Nz+1)
E_np2=np.zeros(Nz+1)
for it in range(0,Nt+1):
t[it] = it*dt
# fuente[it] = math.exp(-(0.5/(sigma*sigma))*(c*(t[it]-tc))*(c*(t[it]-tc)))
fuente[it] = math.sin(w*t[it])
it = 0
E_n[int(Nz/2)]=fuente[it]
it = 1
E_np1[int(Nz/2)]=fuente[it]
for iz in range(0,Nz+1):
z[iz]=iz*dz
plt.ylim(-1,1)
plt.plot(z,E_n,'o')
plt.savefig("E_00{}.png".format(0))
plt.close()
plt.ylim(-1,1)
plt.plot(z,E_np1,'o')
plt.savefig("E_00{}.png".format(1))
plt.close()
for it in range(2,Nt):
for iz in range(0,Nz-1):
E_np2[iz+1]=2*E_np1[iz+1]-E_n[iz+1]+0.25*(E_np1[iz+2]-2*E_np1[iz+1]+E_np1[iz])
E_np2[int(Nz/2)]= fuente[it]
for iz in range(0,Nz-1):
E_n[iz]=E_np1[iz]
E_np1[iz] = E_np2[iz]
plt.ylim(-1,1)
plt.plot(z,E_np2,'-b')
if it < 10: plt.savefig("E_00{}.png".format(it))
if it >= 10 and it < 100: plt.savefig( "E_0{}.png".format(it))
if it >= 100 and it < 1000: plt.savefig( "E_{}.png".format(it))
plt.close()