############################################################################ # # Name: Grundlagen_3D_2.py # Python: 3.3.2 on win32 # Autor: Adrian Haas # ############################################################################ #3D wave equation ############################################################################ # # The script animates a dirac puls hitting a slab in three dimension, # using the explicit FDM (Finite Difference Method) for the wave equation. # A short mathematical description of the wave model for electromagnetic # waves is in the file http://www.engineering-tools.ch/wellengleichung.pdf # (if memory problem -> wrote direct to view/video) # ############################################################################ 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 delta_z=delta_x S=1/3**0.5 delta_t=delta_x/c*S Nx=100 Ny=100 Nz=100 Nt=200 x_0 = numpy.arange(0,Nx*delta_x,delta_x) y_0 = numpy.arange(0,Ny*delta_y,delta_y) z_0 = numpy.arange(0,Ny*delta_z,delta_z) t_0 = numpy.arange(0,Nt*delta_t,delta_t) Ez=numpy.zeros(Nx*Ny*Nz*Nt).reshape(Nx,Ny,Nz,Nt) #source puls=20 #dirac #layers cz=numpy.ones(len(x_0)*len(y_0)*len(z_0)).reshape(len(x_0),len(y_0) ,len(z_0))*c0 sigma_eps=numpy.zeros(len(x_0)*len(y_0)*len(z_0)).reshape(len(x_0),len(y_0) ,len(z_0)) #add objects (with loss) sigma_obj=50 sigma_eps[:,:,60:80]=sigma_obj/eps0 sigma_eps[45:55,45:55,60:80]=0 #ABC constants coeff1=(c*delta_t-delta_x)/(c*delta_t+delta_x) coeff2=(2*delta_x)/(c*delta_t+delta_x) coeff3=(c*delta_t)**2/(2*delta_x*(c*delta_t+delta_x)) ############################################################################ #Iteration ############################################################################ print("Iterate...") for n in range(1,Nt): #wave equation Ez[1:Nx-1,1:Ny-1,1:Nz-1,n]=( (1/(1/delta_t**2+sigma_eps[1:Nx-1,1:Ny-1,1:Nz-1]/delta_t)) *(cz[1:Nx-1,1:Ny-1,1:Nz-1]**2/delta_x**2* (Ez[2:Nx,1:Ny-1,1:Nz-1,n-1]-2*Ez[1:Nx-1,1:Ny-1,1:Nz-1,n-1] +Ez[0:Nx-2,1:Ny-1,1:Nz-1,n-1]) +cz[1:Nx-1,1:Ny-1,1:Nz-1]**2/delta_y**2 *(Ez[1:Nx-1,2:Ny,1:Nz-1,n-1]-2*Ez[1:Nx-1,1:Ny-1,1:Nz-1,n-1] +Ez[1:Nx-1,0:Ny-2,1:Nz-1,n-1]) +cz[1:Nx-1,1:Ny-1,1:Nz-1]**2/delta_z**2 *(Ez[1:Nx-1,1:Ny-1,2:Nz,n-1]-2*Ez[1:Nx-1,1:Ny-1,1:Nz-1,n-1] +Ez[1:Nx-1,1:Ny-1,0:Nz-2,n-1]) +(2/delta_t**2+sigma_eps[1:Nx-1,1:Ny-1,1:Nz-1]/delta_t) *Ez[1:Nx-1,1:Ny-1,1:Nz-1,n-1] -1/delta_t**2*Ez[1:Nx-1,1:Ny-1,1:Nz-1,n-2])) #mur abc for x=0 Ez[0,1:Ny-1,1:Nz-1,n]=( -Ez[1,1:Ny-1,1:Nz-1,n-2] +coeff1*(Ez[1,1:Ny-1,1:Nz-1,n]+Ez[0,1:Ny-1,1:Nz-1,n-2]) +coeff2*(Ez[1,1:Ny-1,1:Nz-1,n-1]+Ez[0,1:Ny-1,1:Nz-1,n-1]) +coeff3*(Ez[0,0:Ny-2,1:Nz-1,n-1]-4*Ez[0,1:Ny-1,1:Nz-1,n-1] +Ez[0,2:Ny,1:Nz-1,n-1] +Ez[1,0:Ny-2,1:Nz-1,n-1]-4*Ez[1,1:Ny-1,1:Nz-1,n-1] +Ez[1,2:Ny,1:Nz-1,n-1] +Ez[0,1:Ny-1,2:Nz,n-1]+Ez[0,1:Ny-1,0:Nz-2,n-1] +Ez[1,1:Ny-1,2:Nz,n-1]+Ez[1,1:Ny-1,0:Nz-2,n-1])) #mur abc for x=Nx-1 Ez[Nx-1,1:Ny-1,1:Nz-1,n]=( -Ez[Nx-2,1:Ny-1,1:Nz-1,n-2] +coeff1*(Ez[Nx-2,1:Ny-1,1:Nz-1,n]+Ez[Nx-1,1:Ny-1,1:Nz-1,n-2]) +coeff2*(Ez[Nx-2,1:Ny-1,1:Nz-1,n-1]+Ez[Nx-1,1:Ny-1,1:Nz-1,n-1]) +coeff3*(Ez[Nx-1,0:Ny-2,1:Nz-1,n-1]-4*Ez[Nx-1,1:Ny-1,1:Nz-1,n-1] +Ez[Nx-1,2:Ny,1:Nz-1,n-1] +Ez[Nx-2,0:Ny-2,1:Nz-1,n-1]-4*Ez[Nx-2,1:Ny-1,1:Nz-1,n-1] +Ez[Nx-2,2:Ny,1:Nz-1,n-1] +Ez[Nx-1,1:Ny-1,2:Nz,n-1]+Ez[Nx-1,1:Ny-1,0:Nz-2,n-1] +Ez[Nx-2,1:Ny-1,2:Nz,n-1]+Ez[Nx-2,1:Ny-1,0:Nz-2,n-1])) #mur abc for y=0 Ez[1:Nx-1,0,1:Nz-1,n]=( -Ez[1:Nx-1,1,1:Nz-1,n-2] +coeff1*(Ez[1:Nx-1,1,1:Nz-1,n]+Ez[1:Nx-1,0,1:Nz-1,n-2]) +coeff2*(Ez[1:Nx-1,1,1:Nz-1,n-1]+Ez[1:Nx-1,0,1:Nz-1,n-1]) +coeff3*(Ez[0:Nx-2,0,1:Nz-1,n-1]-4*Ez[1:Nx-1,0,1:Nz-1,n-1] +Ez[2:Nx,0,1:Nz-1,n-1] +Ez[0:Nx-2,1,1:Nz-1,n-1]-4*Ez[1:Nx-1,1,1:Nz-1,n-1] +Ez[2:Nx,1,1:Nz-1,n-1] +Ez[1:Nx-1,0,2:Nz,n-1]+Ez[1:Nx-1,0,0:Nz-2,n-1] +Ez[1:Nx-1,1,2:Nz,n-1]+Ez[1:Nx-1,1,0:Nz-2,n-1])) #mur abc for y=Ny-1 Ez[1:Nx-1,Ny-1,1:Nz-1,n]=( -Ez[1:Nx-1,Ny-2,1:Nz-1,n-2] +coeff1*(Ez[1:Nx-1,Ny-2,1:Nz-1,n]+Ez[1:Nx-1,Ny-1,1:Nz-1,n-2]) +coeff2*(Ez[1:Nx-1,Ny-2,1:Nz-1,n-1]+Ez[1:Nx-1,Ny-1,1:Nz-1,n-1]) +coeff3*(Ez[0:Nx-2,Ny-1,1:Nz-1,n-1]-4*Ez[1:Nx-1,Ny-1,1:Nz-1,n-1] +Ez[2:Nx,Ny-1,1:Nz-1,n-1] +Ez[0:Nx-2,Ny-2,1:Nz-1,n-1]-4*Ez[1:Ny-1,Ny-2,1:Nz-1,n-1] +Ez[2:Nx,Ny-2,1:Nz-1,n-1] +Ez[1:Nx-1,Ny-1,2:Nz,n-1]+Ez[1:Nx-1,Ny-1,0:Nz-2,n-1] +Ez[1:Nx-1,Ny-2,2:Nz,n-1]+Ez[1:Nx-1,Ny-2,0:Nz-2,n-1])) #mur abc for z=0 Ez[1:Nx-1,1:Ny-1,0,n]=( -Ez[1:Nx-1,1:Ny-1,1,n-2] +coeff1*(Ez[1:Nx-1,1:Ny-1,1,n]+Ez[1:Nx-1,1:Ny-1,0,n-2]) +coeff2*(Ez[1:Nx-1,1:Ny-1,1,n-1]+Ez[1:Nx-1,1:Ny-1,0,n-1]) +coeff3*(Ez[0:Nx-2,1:Ny-1,0,n-1]-4*Ez[1:Nx-1,1:Ny-1,0,n-1] +Ez[2:Nx,1:Ny-1,0,n-1] +Ez[0:Nx-2,1:Ny-1,1,n-1]-4*Ez[1:Nx-1,1:Ny-1,1,n-1] +Ez[2:Nx,1:Ny-1,1,n-1] +Ez[1:Nx-1,2:Ny,0,n-1]+Ez[1:Nx-1,0:Ny-2,0,n-1] +Ez[1:Nx-1,2:Ny,1,n-1]+Ez[1:Nx-1,0:Ny-2,1,n-1])) #mur abc for z=Nz-1 Ez[1:Nx-1,1:Ny-1,Nz-1,n]=( -Ez[1:Nx-1,1:Ny-1,Nz-2,n-2] +coeff1*(Ez[1:Nx-1,1:Ny-1,Nz-2,n]+Ez[1:Nx-1,1:Ny-1,Nz-1,n-2]) +coeff2*(Ez[1:Nx-1,1:Ny-1,Nz-2,n-1]+Ez[1:Nx-1,1:Ny-1,Nz-1,n-1]) +coeff3*(Ez[0:Nx-2,1:Ny-1,Nz-1,n-1]-4*Ez[1:Nx-1,1:Ny-1,Nz-1,n-1] +Ez[2:Nx,1:Ny-1,Nz-1,n-1] +Ez[0:Nx-2,1:Ny-1,Nz-2,n-1]-4*Ez[1:Ny-1,1:Ny-1,Nz-2,n-1] +Ez[2:Nx,1:Ny-1,Nz-2,n-1] +Ez[1:Nx-1,2:Ny,Nz-1,n-1]+Ez[1:Nx-1,0:Ny-2,Nz-1,n-1] +Ez[1:Nx-1,2:Ny,Nz-2,n-1]+Ez[1:Nx-1,0:Ny-2,Nz-2,n-1])) #condition for edges 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] Ez[:,0,0,n]=Ez[:,1,1,n-2] Ez[:,0,Nz-1,n]=Ez[:,1,Nz-2,n-2] Ez[:,Ny-1,0,n]=Ez[:,Ny-2,1,n-2] Ez[:,Ny-1,Nz-1,n]=Ez[:,Ny-2,Nz-2,n-2] Ez[0,:,0,n]=Ez[1,:,1,n-2] Ez[0,:,Nz-1,n]=Ez[1,:,Nz-2,n-2] Ez[Nx-1,:,0,n]=Ez[Nx-2,:,1,n-2] Ez[Nx-1,:,Nz-1,n]=Ez[Nx-2,:,Nz-2,n-2] #source if n==3: Ez[48:52,48:52,48:52,3]=puls print("done!") ############################################################################ #Animation ############################################################################ #animation function def ani(i): time="Time: "+str(round(i*numpy.max(t_0)*1e9/len(t_0),1))+" ns" wave1=Ez_x_1[:,:,i+1] wave2=Ez_x_2[:,:,i+1] wave3=Ez_x_3[:,:,i+1] wave4=Ez_x_4[:,:,i+1] wave5=Ez_x_5[:,:,i+1] wave6=Ez_x_6[:,:,i+1] textz.set_text(time) pic1.set_data(wave1) pic2.set_data(wave2) pic3.set_data(wave3) pic4.set_data(wave4) pic5.set_data(wave5) pic6.set_data(wave6) return pic1,pic2,pic3,pic4,pic5,pic6,textz #init parameters time="Time: 0 ns" Ez_x_1=Ez[:,:,10,:] Ez_x_2=Ez[:,:,30,:] Ez_x_3=Ez[:,:,50,:] Ez_x_4=Ez[:,:,70,:] Ez_x_5=Ez[:,:,90,:] Ez_x_6=Ez[:,50,:,:] #create plot fig=plt.figure("3D Wellengleichung Animation") plt.suptitle("$u_{tt}+a u_{t}-c^2 (u_{xx}+ u_{yy}+u_{zz})=0$", fontsize=14) texti=fig.text(0.12, 0.06, "Grid size: "+str(Nx)+"*"+str(Ny)+"*"+str(Nz)+" cells", fontsize=11) textc=fig.text(0.45, 0.06, "Cell size: "+str(1000*delta_x)+"*"+str(1000*delta_y)+"*" +str(1000*delta_z)+" mm", fontsize=11) textz=fig.text(0.78, 0.06, time) ax=fig.add_subplot(231) plt.title("xy-plane, 0.1Z", fontsize=11) wave1=Ez_x_1[:,:,0] pic1=ax.imshow(wave1,interpolation='nearest', cmap=cm.jet, vmin=-1, vmax=1) sigma_eps[0,0,10]=sigma_obj #adjust color device=ax.imshow(-sigma_eps[:,:,10],interpolation='nearest', cmap=cm.gray,alpha=0.2) ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) ax=fig.add_subplot(232) plt.title("xy-plane, 0.3Z", fontsize=11) wave2=Ez_x_2[:,:,0] pic2=ax.imshow(wave2,interpolation='nearest', cmap=cm.jet, vmin=-1, vmax=1) sigma_eps[0,0,30]=sigma_obj #adjust color device=ax.imshow(-sigma_eps[:,:,30],interpolation='nearest', cmap=cm.gray,alpha=0.2) ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) ax=fig.add_subplot(233) plt.title("xy-plane, 0.5Z", fontsize=11) wave3=Ez_x_3[:,:,0] pic3=ax.imshow(wave3,interpolation='nearest', cmap=cm.jet, vmin=-1, vmax=1) sigma_eps[0,0,50]=sigma_obj #adjust color device=ax.imshow(-sigma_eps[:,:,50],interpolation='nearest', cmap=cm.gray,alpha=0.2) ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) ax=fig.add_subplot(234) plt.title("xy-plane, 0.7Z", fontsize=11) wave4=Ez_x_4[:,:,0] pic4=ax.imshow(wave4,interpolation='nearest', cmap=cm.jet, vmin=-1, vmax=1) device=ax.imshow(-sigma_eps[:,:,70],interpolation='nearest', cmap=cm.gray,alpha=0.2) ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) ax=fig.add_subplot(235) plt.title("xy-plane, 0.9Z", fontsize=11) wave5=Ez_x_5[:,:,0] pic5=ax.imshow(wave3,interpolation='nearest', cmap=cm.jet, vmin=-1, vmax=1) sigma_eps[0,0,90]=sigma_obj #adjust color device=ax.imshow(-sigma_eps[:,:,90],interpolation='nearest', cmap=cm.gray,alpha=0.2) ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) ax=fig.add_subplot(236) plt.title("xz-plane, 0.5Y", fontsize=11) wave6=Ez_x_6[:,:,0] pic6=ax.imshow(wave6,interpolation='nearest', cmap=cm.jet, vmin=-1, vmax=1) device=ax.imshow(-sigma_eps[:,50,:],interpolation='nearest', cmap=cm.gray,alpha=0.2) ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) #animate wave_ani = animation.FuncAnimation(fig, ani, frames=int((len(t_0)-1)),interval=1, blit=False) plt.show()