1. ############################################################################
    #
    #   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()