Script 02

↓ output ×
 ############################################################################
 #
 #   Name:           Grundlagen_2D_3.py
 #   Python:         3.3.2 on win32
 #   Autor:          Adrian Haas
 #   
 ############################################################################
 #2D wave equation
 ############################################################################
 #
 # The script animates a dirac puls hitting a triangle and a circle in two
 # dimension, using the explicit FDM (Finite Difference Method) for the wave
 # equation. It shows the splitting and focusing effect. A short mathematical
 # description of the wave model for electromagnetic waves is in the file
 # http://www.engineering-tools.ch/wellengleichung.pdf
 #
 ############################################################################

 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=1e-3
 delta_y=delta_x
 S=1/2**0.5
 delta_t=delta_x/c*S
 Nx=200
 Ny=200
 Nt=400

 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
 func0=(lambda tx,ty,a,meanx,varx,meany,vary:
      a*numpy.exp(-((tx-meanx)**2/2/varx+(ty-meany)**2/2/vary)))
 puls=numpy.fromfunction(lambda i,j:
                         func0(i,j,1,10,5,10,5),(20,20),dtype=int)

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

 #add objects
 eps_obj=3
 c_obj=1/numpy.sqrt(eps_obj)*c0

 # a) circle
 circle_obj=cz[50:90,100:140]
 circle=numpy.fromfunction(lambda i,j:
                           (i-20)**2+(j-20)**2<20**2,(40,40),dtype=int)
 circle_obj[circle]=c_obj

 # b) triangle
 triangle_obj=cz[120:150,30:60]
 triangle=numpy.fromfunction(lambda i,j: (i-j)<2,(30,30),dtype=int)
 triangle_obj[triangle]=c_obj

 #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])

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

      #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]


      #source
      if n==3:
        Ez[40:60,40:60,3]=puls

 ############################################################################
 #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(wave)
   textz.set_text(time)
   return pic,textz,

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

 #create plot
 fig=plt.figure("Split and Focus")
 plt.title("$u_{tt}-c^2 (u_{xx}+ u_{yy})=0$")
 ax=fig.add_subplot(111)

 pic=ax.imshow(wave,interpolation='nearest', cmap=cm.jet,
               vmin=-0.5, vmax=0.5)
 device=ax.imshow(cz,interpolation='nearest', cmap=cm.gray,alpha=0.2)

 texti=fig.text(0.5, 0.2,
                "Grid size: "+str(Nx)+"*"+str(Ny)+" cells", color='green')
 textc=fig.text(0.5, 0.17,
                "Cell size: "+str(1000*delta_x)+"*"+str(1000*delta_y)+" mm",
                color='green')
 textz=fig.text(0.5, 0.14, zeit, color='green')

 ax.axes.get_xaxis().set_visible(False)
 ax.axes.get_yaxis().set_visible(False)

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

 plt.show()
↑ script ×