############################################################################ # # Name: Parabolic.py # Python: 3.3.2 on win32 # Autor: Adrian Haas # ############################################################################ #Parabolic antenna plane wave reflection ############################################################################ # # Parabolic reflection of a gaussian pulse. Change parameter ef for # different kind of parabolics. 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=1e-3 delta_y=delta_x S=1/2**0.5 delta_t=delta_x/c*S Nx=400 Ny=400 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 sigma_eps=numpy.zeros(Nx*Ny).reshape(Nx,Ny) #add parabolic sigma_obj=50 s=0 ef=1.2 #5,0.5 while s<Ny-2: t=(s-(Ny-2)/2)**2/(Ny-2)**2 sigma_eps[55+int(220*t*ef):60+int(220*t*ef),s]=sigma_obj/eps0 sigma_eps[:,:50]=0 sigma_eps[:,len(y_0)-2-50:]=0 s+=1 #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]=((1/(1/delta_t**2 +sigma_eps[1:Nx-1,1:Ny-1]/delta_t)) *(cz[1:Nx-1,1:Ny-1]**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_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/delta_t**2 +sigma_eps[1:Nx-1,1:Ny-1]/delta_t) *Ez[1:Nx-1,1:Ny-1,n-1]-1/delta_t**2 *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[190:210,190:210,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("Parabol") plt.title("$u_{tt}+au_{t}-c^2 (u_{xx}+ u_{yy})=0$") ax=fig.add_subplot(111) pic=ax.imshow(wave,interpolation='nearest', cmap=cm.Blues, vmin=-0.5, vmax=0.5) device=ax.imshow(-sigma_eps,interpolation='nearest', cmap=cm.gray,alpha=0.2) texti=fig.text(0.3, 0.85, "Parabol reflects plane wave", color='white') texti=fig.text(0.5, 0.2, "Grid size: "+str(Nx)+"*"+str(Ny)+" cells", color='white') textc=fig.text(0.5, 0.17, "Cell size: "+str(1000*delta_x)+"*"+str(1000*delta_y)+" mm", color='white') textz=fig.text(0.5, 0.14, zeit, color='white') 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()