compara

import numpy as np
import math
import matplotlib.pyplot as plt

c   = 3e8
l0  = 600e-9
t0  = l0/c

k   = 1.0*((2*math.pi)/l0)
w   = c*k

dz  = l0/40
Nz  = 100
zi  = -Nz*dz
zf  =  Nz*dz

z   = np.zeros(2*Nz+1)
Exf = np.zeros(2*Nz+1) # solucion de fourier

for iz in range(0,2*Nz+1):
  z[iz]   = zi+iz*dz
  Exf[iz] = math.sin(k*z[iz]) 

Nt=150
dt=dz/(2*c)
t      = np.zeros(Nt+1)
fuente = np.zeros(Nt+1)

#for it in range(0,Nt+1):
#  t[it]=it*dt
#  fuente[it] = math.sin(w*t[it])
#  for iz in range(0,2*Nz+1):
#    Exa[iz] = math.sin(k*z[iz]-w*t[it])

exn = np.zeros(2*Nz+1)
hyn = np.zeros(2*Nz+1)
Exa = np.zeros(2*Nz+1)
for it in range(0,Nt,2):
  t[it]=it*dt
  fuente[it] = math.sin(-w*t[it])
  for iz in range(0,2*Nz-2,2):
    exn[iz+2] = exn[iz+2]-0.5*(hyn[iz+3]-hyn[iz+1])
  exn[int(Nz)]=fuente[it]
  for iz in range(0,2*Nz,2):
    hyn[iz+1] = hyn[iz+1]-0.5*(exn[iz+2]-exn[iz])
  for iz in range(0,2*Nz+1):
    Exa[iz] = math.sin(k*z[iz]-w*t[it])
    plt.grid()  
  plt.ylim(-1,1)
  plt.plot(z,Exa,'-b',z,exn,'.r')
  if it<  10:            plt.savefig("exa00{}.png".format(it))
  if it>= 10 and it<100: plt.savefig( "exa0{}.png".format(it))
  if it>= 100:           plt.savefig(  "exa{}.png".format(it))
  plt.close()

Nt=80
t     = np.zeros(Nt+1)
fuente= np.zeros(Nt+1)
E_n   = np.zeros(Nz+1)
E_np1 = np.zeros(Nz+1)
E_np2 = np.zeros(Nz+1)
Ea    = np.zeros(Nz+1)
z=np.zeros(Nz+1)

for iz in range(0,Nz+1):
  z[iz]=iz*dz

for it in range(0,Nt+1):
  t[it]  = it*dt
  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 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])
    Ea[iz]=math.sin(k*(z[iz]-(Nz/2)*dz)-w*t[it])
  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,'ob',z,Ea,'-r')
  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()