dip_fdtd.f90.rar_3d fdtd matlab pml_FDTD 偶极子_dip_天线 FDTD_天线 f

program test_fdtd !-- Fortran code for FDTD with Berenger PMLs, version 1.0, May 1999 !-- by Jos Bergervoet. !-- Plot field and/or Poynting vector S around radiating linear dipole implicit none integer, parameter :: wp = selected_real_kind(5) ! Choose working precision ! 5 or 15 (single/double) real(wp),parameter :: c = 299792458, & nu = 2 ! nu=order of Berenger increase real(wp), allocatable, dimension(:) :: sxe,sye,sze, sxh,syh,szh real(wp), allocatable, dimension(:,:,:,:) :: exs,eys,ezs, hxs,hys,hzs character(len=20) :: name integer :: nx,ny,nz, nb, nt, ns, n real(wp) :: Pi, mu0, Z0, eps0, & fun, g, t, period, smax, F_in, dx, dy, dz, dt ! g(t) = sin(t)**6 ! smooth, 5 times differentiable ! g(t) = sin(2*t) * sin(t)**5 ! `double-delta' ! fun(t) = g(min(t/period,1.0_wp)*Pi) ! scaled g(t) fun(t) = sin(2*Pi*t/period) ! CW print "(1x,a)", & "Give: nx,ny,nz, nt, nb, smax, period, n_s", & "For instance: 60 60 30 310 10 4e8 110 30", & "", & "nb and smax are Berenger parameters, n_s = nr. of segm. in antenna", & "" read *, nx,ny,nz, nt, nb,smax, period, ns call setup() n = 0 Fdtd_loop: do n=n+1 !-- divide source signal over both E-field species F_in = fun(n+0.0_wp) exs(nx/2-1:nx/2,ny/2,nz/2,:) = - F_in/2 if (mod(n, max(nint(period/30),1))==0) then write(name,*) n name = trim(adjustl(name)) // ".gif" !-- Plot S and some field component call plot_log("S"//name, compute_S_z(k_z=nz/2+ns/10)) call plot("E"//name, exs(:,:,nz/2+ns/10,1)+exs(:,:,nz/2+ns/10,2) ) end if if (n>nt) exit Fdtd_loop !-- update the fields including the Berenger PMLs call Berenger(exs,eys,ezs, hxs,hys,hzs) !-- apply PEC condition exs((nx-ns)/2:(nx+ns)/2-1, ny/2, nz/2, :) = 0 end do Fdtd_loop stop contains subroutine setup() real(wp) :: rho, x integer :: nmax, i rho(x,nmax) = max(-x,x-nmax,0.0_wp)/max(nb,1) Pi = acos(-1d0) mu0 = 4d-7*Pi Z0 = mu0*c eps0= 1/(Z0*c) dx = 1; dy = 1; dz = 1 dt = 1 / ( c*sqrt(1/dx**2+1/dy**2+1/dz**2) ) ! acc. to Courant allocate( exs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2), & ! All arrays include eys(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2), & ! the Berenger PML ezs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2), & ! region (nb cells) hxs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2), & ! All fields have two hys(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2), & ! species (last index) hzs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2) ) allocate( sxe(-nb:nx+nb), & sye(-nb:ny+nb), & ! these are the sigma sze(-nb:nz+nb), & ! arrays (scaled by sxh(-nb:nx+nb), & ! eps0 and mu0) syh(-nb:ny+nb), & szh(-nb:nz+nb) ) exs = 0; eys = 0; ezs = 0; hxs = 0; hys = 0; hzs = 0 sxe = (/( smax * rho(i+1.0_wp,nx)**nu, i=-nb,nx+nb )/) ! 1.0 is clumsy ? sxh = (/( smax * rho(i+0.5_wp,nx)**nu, i=-nb,nx+nb )/) sye = (/( smax * rho(i+1.0_wp,ny)**nu, i=-nb,ny+nb )/) syh = (/( smax * rho(i+0.5_wp,ny)**nu, i=-nb,ny+nb )/) sze = (/( smax * rho(i+1.0_wp,nz)**nu, i=-nb,nz+nb )/) szh = (/( smax * rho(i+0.5_wp,nz)**nu, i=-nb,nz+nb )/) write(*,*) "R_normal = ", exp(-2.0_wp/(nu+1) * smax*nb/c) end subroutine function compute_S(k_z) result(S) integer, intent(in) :: k_z real(wp), dimension(-nb+1:nx+nb-1, -nb+1:ny+nb-1) :: S real(wp), dimension(-nb:nx+nb, -nb:ny+nb, 3) :: E, H real(wp), dimension(-nb:nx+nb, -nb:ny+nb) :: temp integer :: i do i=1,3 E(:,:,1) = exs(:,:,k_z,1) + exs(:,:,k_z,2) ! problem: E and H are not E(:,:,2) = eys(:,:,k_z,1) + eys(:,:,k_z,2) ! available for exactly the E(:,:,3) = ezs(:,:,k_z,1) + ezs(:,:,k_z,2) ! same positions H(:,:,1) = hxs(:,:,k_z,1) + hxs(:,:,k_z,2) H(:,:,2) = hys(:,:,k_z,1) + hys(:,:,k_z,2) H(:,:,3) = hzs(:,:,k_z,1) + hzs(:,:,k_z,2) end do !-- Use: S^2 = (E x H)^2 = E^2 H^2 - (E dot H)^2 temp = sqrt( sum(E*E, dim=3)*sum(H*H, dim=3) - sum(E*H, dim=3)**2 ) temp( (nx-ns)/2:(nx+ns)/2-1, ny/2 ) = 0 ! Just to show the dipole-wire S = temp(-nb+1:nx+nb-1, -nb+1:ny+nb-1) end function compute_S function compute_S_z(k_z) result(S_z) integer, intent(in) :: k_z real(wp), dimension(-nb+1:nx+nb-1, -nb+1:ny+nb-1) :: S_z real(wp), dimension(-nb:nx+nb, -nb:ny+nb, 3) :: E, H real(wp), dimension(-nb:nx+nb, -nb:ny+nb) :: temp integer :: i do i=1,3 E(:,:,1) = exs(:,:,k_z,1) + exs(:,:,k_z,2) ! problem: E and H are not E(:,:,2) = eys(:,:,k_z,1) + eys(:,:,k_z,2) ! available for exactly the E(:,:,3) = ezs(:,:,k_z,1) + ezs(:,:,k_z,2) ! same positions H(:,:,1) = hxs(:,:,k_z,1) + hxs(:,:,k_z,2) H(:,:,2) = hys(:,:,k_z,1) + hys(:,:,k_z,2) H(:,:,3) = hzs(:,:,k_z,1) + hzs(:,:,k_z,2) end do !-- Use: S_z = E_x H_y - E_y H_x temp = E(:,:,1)*H(:,:,2) - E(:,:,2)*H(:,:,1) temp( (nx-ns)/2:(nx+ns)/2-1, ny/2 ) = 0 ! Just to show the dipole-wire S_z = temp(-nb+1:nx+nb-1, -nb+1:ny+nb-1) end function compute_S_z subroutine plot_log(Name, v) !-- logscale picture, peak-hold range between calls use gif_util character(len=*) :: Name real(wp) :: range=0, ndecades=3, v(:,:) integer :: i, Pixel(size(v,1), size(v,2)) type(color) :: ColMap(256) ColMap = (/ ( color(2*i-2,0,0), i=1,128 ), & ! black-to-red ( color(255,5*i/2,0), i=1,101 ), & ! red-to-yellow ( color(255,255,9*i), i=1,27 ) /) ! yellow-to-white range = max(range, maxval(abs(v)), 1e-30_wp ) Pixel = nint( max(0.0,1+log10(abs(v)/range)/ndecades) * 255 ) write(*,*) "Max(|",trim(Name),"|) =", maxval(abs(v)) call writegif (Name, Pixel, ColMap) end subroutine plot_log subroutine plot(Name, v) !-- linear scale picture, peak-hold range between calls use gif_util character(len=*) :: Name real(wp) :: range=0, v(:,:) integer :: i, Pixel(size(v,1), size(v,2)) type(color) :: ColMap(256) ColMap = (/( color(i-1,i-1,256-i), i=1,256 )/) range = max(range, maxval(abs(v)), 1e-30_wp ) Pixel = nint( (v+range)/range/2 * 255 ) write(*,*) "Max(|",trim(Name),"|) =", maxval(abs(v)) call writegif (Name, Pixel, ColMap) end subroutine plot subroutine Berenger(Ex,Ey,Ez,Hx,Hy,Hz) !-- Updates all fields, except the outer PEC faces of the domain. !-- Uses internal fields and E-fields on faces (not H on faces) real(wp), dimension(0:,0:,0:,:) :: Ex,Ey,Ez,Hx,Hy,Hz call add_diff_y(syh, -mu0, Hx(1:,0:,0:,1), Ez(1:,0:,0:,:) ) call add_diff_z(szh, mu0, Hx(1:,0:,0:,2), Ey(1:,0:,0:,:) ) call add_diff_z(szh, -mu0, Hy(0:,1:,0:,1), Ex(0:,1:,0:,:) ) call add_diff_x(sxh, mu0, Hy(0:,1:,0:,2), Ez(0:,1:,0:,:) ) call add_diff_x(sxh, -mu0, Hz(0:,0:,1:,1), Ey(0:,0:,1:,:) ) call add_diff_y(syh, mu0, Hz(0:,0:,1:,2), Ex(0:,0:,1:,:) ) call add_diff_y(sye, eps0, Ex(0:,1:,1:,1), Hz(0:,0:,1:,:) ) call add_diff_z(sze, -eps0, Ex(0:,1:,1:,2), Hy(0:,1:,0:,:) ) call add_diff_z(sze, eps0, Ey(1:,0:,1:,1), Hx(1:,0:,0:,:) ) call add_diff_x(sxe, -eps0, Ey(1:,0:,1:,2), Hz(0:,0:,1:,:) ) call add_diff_x(sxe, eps0, Ez(1:,1:,0:,1), Hy(0:,1:,0:,:) ) call add_diff_y(sye, -eps0, Ez(1:,1:,0:,2), Hx(1:,0:,0:,:) ) end subroutine su