Script 05

↓ output ×
 ############################################################################
 #
 #   Name:           plane_earth.py
 #   Python:         3.3.2 on win32
 #   Autor:          Adrian Haas
 #   
 ############################################################################
 #plane earth model
 ############################################################################
 #
 # Numerical animation for the plane earth model with typical rays of
 # interference. Theory i.e:
 # http://www.waves.utoronto.ca/prof/svhum/ece422/notes/18a-planeearth.pdf.
 # The script is a modification of Grundlagen_2D_3.py.
 #
 ############################################################################

 import numpy
 import matplotlib.pyplot as plt
 import matplotlib.animation as animation
 from matplotlib import cm

 ############################################################################
 #define parameters and grid
 ############################################################################

 #constants
 eps0=8.854e-12
 mu0=1.256e-6
 c0=1/numpy.sqrt(eps0*mu0)
 c=c0

 #define grid
 delta_x=5e-2
 delta_y=delta_x
 S=1/2**0.5
 delta_t=delta_x/c*S
 Nx=200
 Ny=500
 Nt=1000

 x_0 = numpy.arange(0,Nx*delta_x,delta_x)
 y_0 = numpy.arange(0,Ny*delta_y,delta_y)
 t_0 = numpy.arange(0,Nt*delta_t,delta_t)

 Ez=numpy.zeros(Nx*Ny*Nt).reshape(Nx,Ny,Nt)

 #source parameter
 freq=850e6
 wave_length=c/freq

 #layer
 cz=numpy.ones(Nx*Ny).reshape(Nx,Ny)*c0

 #ABC constants
 k1=(c*delta_t-delta_x)/(c*delta_t+delta_x)
 k2=(2*delta_x)/(c*delta_t+delta_x)
 k3=(c*delta_t)**2/2/delta_x/(c*delta_t+delta_x)

 ############################################################################
 #Iteration 
 ############################################################################

 for n in range(1,Nt):

      #wave equation

      Ez[1:Nx-1,1:Ny-1,n]=(cz[1:Nx-1,1:Ny-1]**2*delta_t**2/delta_x**2
                           *(Ez[2:Nx,1:Ny-1,n-1]-2*Ez[1:Nx-1,1:Ny-1,n-1]
                             +Ez[0:Nx-2,1:Ny-1,n-1])
                          +cz[1:Nx-1,1:Ny-1]**2*delta_t**2/delta_y**2
                           *(Ez[1:Nx-1,2:Ny,n-1]-2*Ez[1:Nx-1,1:Ny-1,n-1]
                             +Ez[1:Nx-1,0:Ny-2,n-1])
                          +2*Ez[1:Nx-1,1:Ny-1,n-1]-Ez[1:Nx-1,1:Ny-1,n-2])

      #none Mur Abc x=Nx-1

      #full reflection for x=Nx-1, Ez[Nx-1,1:Ny-1,n]=0 by default

      #Mur Abc x=0

      Ez[0,1:Ny-1,n]=(-Ez[1,1:Ny-1,n-2]
                         +k1*(Ez[1,1:Ny-1,n]+Ez[0,1:Ny-1,n-2])
                         +k2*(Ez[1,1:Ny-1,n-1]+Ez[0,1:Ny-1,n-1])
                         +k3*(Ez[0,0:Ny-2,n-1]-2*Ez[0,1:Ny-1,n-1]
                              +Ez[0,2:Ny,n-1]
                              +Ez[1,0:Ny-2,n-1]-2*Ez[1,1:Ny-1,n-1]
                              +Ez[1,2:Ny,n-1]))

      #Mur Abc y=Ny-1

      Ez[1:Nx-1,Ny-1,n]=(-Ez[1:Nx-1,Ny-2,n-2]
                          +k1*(Ez[1:Nx-1,Ny-2,n]+Ez[1:Nx-1,Ny-1,n-2])
                          +k2*(Ez[1:Nx-1,Ny-1,n-1]+Ez[1:Nx-1,Ny-2,n-1])
                          +k3*(Ez[0:Nx-2,Ny-2,n-1]-2*Ez[1:Nx-1,Ny-2,n-1]
                               +Ez[2:Nx,Ny-2,n-1]
                               +Ez[0:Nx-2,Ny-1,n-1]-2*Ez[1:Nx-1,Ny-1,n-1]
                               +Ez[2:Nx,Ny-1,n-1]))

      #Mur Abc y=0

      Ez[1:Nx-1,0,n]=(-Ez[1:Nx-1,1,n-2]
                     +k1*(Ez[1:Nx-1,1,n]+Ez[1:Nx-1,0,n-2])
                     +k2*(Ez[1:Nx-1,1,n-1]+Ez[1:Nx-1,0,n-1])
                     +k3*(Ez[0:Nx-2,0,n-1]-2*Ez[1:Nx-1,0,n-1]
                          +Ez[2:Nx,0,n-1]
                          +Ez[0:Nx-2,1,n-1]-2*Ez[1:Nx-1,1,n-1]
                          +Ez[2:Nx,1,n-1]))

      #corner condition

      Ez[0,0,n]=Ez[1,1,n-2]
      Ez[0,Ny-1,n]=Ez[1,Ny-2,n-2]
      Ez[Nx-1,0,n]=Ez[Nx-2,1,n-2]
      Ez[Nx-1,Ny-1,n]=Ez[Nx-2,Ny-2,n-2]


      #cont. source
      Ez[159:160,59:60,n]=numpy.cos(numpy.pi*freq*n*delta_t)

 ############################################################################
 #Animation
 ############################################################################

 #animation function

 def ani(i):
   wave=Ez[:,:,i+1]
   time="Time: "+str(round(i*numpy.max(t_0)*1e9/len(t_0),1))+" ns"
   pic.set_data(abs(wave))
   textz.set_text(time)
   return pic,textz,

 #init parameters

 zeit="Time: 0 ns"
 wave=Ez[:,:,0]

 #create plot

 fig=plt.figure("Plane earth model")
 plt.title("$u_{tt}-c^2 (u_{xx}+ u_{yy})=0$")
 ax=fig.add_subplot(111)

 pic=ax.imshow(abs(wave),extent=(0, 25, 0,10),interpolation='nearest',
               cmap=cm.Blues,  vmin=0, vmax=0.1)

 textf=fig.text(0.1, 0.8, "Frequency: "+str(int(freq*1e-6))+" MHz")
 textwf=fig.text(0.7, 0.8,
                 "Wavelength: "+str(numpy.round(wave_length,2))+" m")
 texti=fig.text(0.1, 0.15, "Grid size: "+str(Nx)+"*"+str(Ny)+" cells")
 textc=fig.text(0.4, 0.15,
                "Cell size: "+str(1000*delta_x)+"*"+str(1000*delta_y)+" mm")
 textz=fig.text(0.75, 0.15, zeit)

 ax.set_xlabel('[m]')
 ax.set_ylabel('[m]',rotation=0)


 #Animate
 wave_ani = animation.FuncAnimation(fig, ani, frames=Nt-1,interval=1,
                                    blit=False)

 plt.show()
↑ script ×