fdtd09.f90

FDTD algorithm.

program test

  !-- version 0.9, May 1999
  !-- by J. Bergervoet
  
    implicit none
    integer, parameter  :: wp = selected_real_kind(5)
    
    real(wp), parameter  :: c = 299792458, nu=2
    real(wp)             :: Pi, mu0, Z0, eps0, smax
    real(wp), allocatable, dimension(:)        :: sxe,sye,sze, sxh,syh,szh
    real(wp), allocatable, dimension(:,:,:,:)  :: exs,eys,ezs, hxs,hys,hzs
    character(len=20) :: s

    integer nx,ny,nz, nb, nt, n
    include 'metric.inc'
    real(wp) :: fun, g, t, tpulse

    g(t) = sin(2*t) * sin(t)**5
    ! fun(t) = max(0.0_wp,2-t)
    ! fun(t) = sin(min(t/tpulse,1.0_wp)*Pi)**2
    fun(t) = g(min(t/tpulse,1.0_wp)*Pi)
    
    dx = 1
    dy = 1
    dz = 1

    read *, nb, nx,ny,nz, nt, smax, tpulse
    nsizex = nx + 1 + 2*nb
    nsizey = ny + 1 + 2*nb
    nsizez = nz + 1 + 2*nb

    allocate( exs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &
              eys(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &
              ezs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &
              hxs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &
              hys(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &
              hzs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2) )         ! assume 0 ???
    allocate( sxe(-nb:nx+nb),  &
              sye(-nb:ny+nb),  &
              sze(-nb:nz+nb),  &
              sxh(-nb:nx+nb),  &
              syh(-nb:ny+nb),  &
              szh(-nb:nz+nb) )         ! assume 0 ???

    call setup()
    
    n = 0
    do
      n=n+1
      ! ezs(nx/2,ny/2,nz/2-2:nz/2+2,:) = fun(n+0.0_wp)
      if(n<=tpulse) then
        t = fun(n+0.0_wp)
        ezs(nx/2,ny/2,nz/2,:) = ezs(nx/2,ny/2,nz/2,:) + t
      end if

      ! if(n<=tpulse) ezs(nx/2,ny/2,nz/2,:) = fun(n+0.0_wp)
!     write(s,*) n
!     s = trim(adjustl(s)) // ".gif"
!     call plot(s, ezs(:,:,nz/2,1)+ezs(:,:,nz/2,2))

      if(n>nt) exit

      call Berenger(exs,eys,ezs, hxs,hys,hzs, nx+2*nb,ny+2*nb,nz+2*nb)

      ! combine internal
      ! call combine( Exs(0:nx-1,0:ny,0:nz,:) )
      ! call combine( Eys(0:nx,0:ny-1,0:nz,:) )
      ! call combine( Ezs(0:nx,0:ny,0:nz-1,:) )
      ! call combine( Hxs(0:nx,0:ny-1,0:nz-1,:) )
      ! call combine( Hys(0:nx-1,0:ny,0:nz-1,:) )
      ! call combine( Hzs(0:nx-1,0:ny-1,0:nz,:) )
    end do


    call plot("ez_int.gif", ezs(0:nx,0:ny,nz/2,1)+ezs(0:nx,0:ny,nz/2,2))
    call plot("ez.gif", ezs(:,:,nz/2,1)+ezs(:,:,nz/2,2))
    call plot("hx.gif", hxs(:,:,nz/2,1)+hxs(:,:,nz/2,2))
!    call plot("ez1.gif", ezs(:,:,nz/2,1))
!    call plot("hx1.gif", hxs(:,:,nz/2,1))
!    call plot("hx2.gif", hxs(:,:,nz/2,1))
    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)

    dt   =  1 / ( c*sqrt(1/dx**2+1/dy**2+1/dz**2) )     ! * 0.9 ???

    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)

    exs = 0
    eys = 0
    ezs = 0
    hxs = 0
    hys = 0
    hzs = 0

  end subroutine
        

  subroutine plot(Name, v)
    use gif_util
    character(len=*) :: Name
    real(wp)         :: range, 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( maxval(abs(v)), 1e-30_wp )
    
    Pixel = nint( (v+range)/range/2 * 255 )

    if(size(v)<130)write(*,*) v
    if(size(v)<130)write(*,*) Pixel

    write(*,*) "Max(|",trim(Name),"|) =", maxval(abs(v)) 
    call writegif (Name, Pixel, ColMap)
  end subroutine
  

  subroutine Berenger(Ex,Ey,Ez,Hx,Hy,Hz, nx,ny,nz)
  !--
  !-- Updates only the fields not on the edge of the domain.
  !-- Uses only internal E and H fields and E-fields on border.
    integer  :: nx,ny,nz
    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,1) )
    call add_diff_z(szh,  mu0,  Hx(1:,0:,0:,2:), Ey(1,0,0,1) )

    call add_diff_z(szh, -mu0,  Hy(0:,1:,0:,1:), Ex(0,1,0,1) )
    call add_diff_x(sxh,  mu0,  Hy(0:,1:,0:,2:), Ez(0,1,0,1) )
 
    call add_diff_x(sxh, -mu0,  Hz(0:,0:,1:,1:), Ey(0,0,1,1) )
    call add_diff_y(syh,  mu0,  Hz(0:,0:,1:,2:), Ex(0,0,1,1) )

    call add_diff_y(sye,  eps0, Ex(0:,1:,1:,1:), Hz(0,0,1,1) )
    call add_diff_z(sze, -eps0, Ex(0:,1:,1:,2:), Hy(0,1,0,1) )
 
    call add_diff_z(sze,  eps0, Ey(1:,0:,1:,1:), Hx(1,0,0,1) )
    call add_diff_x(sxe, -eps0, Ey(1:,0:,1:,2:), Hz(0,0,1,1) )
 
    call add_diff_x(sxe,  eps0, Ez(1:,1:,0:,1:), Hy(0,1,0,1) )
    call add_diff_y(sye, -eps0, Ez(1:,1:,0:,2:), Hx(1,0,0,1) )
 
  end subroutine


  subroutine Combine(v)
    real(wp), dimension(:,:,:,:) :: v

    v(:,:,:,1) = (v(:,:,:,1)+v(:,:,:,2)) / 2
    v(:,:,:,2) = v(:,:,:,1)
  end subroutine


      subroutine add_diff_x(s, pfac, v, w)
        implicit none
        integer i,j,k
        real v(:,:,:,:)
        real w(size(v,1),size(v,2),size(v,3),2)
        real pfac, alpha(size(v,1)), beta(size(v,1)), s(*)

     ! (x,y,z  e,h, ...  instead of  _x(s,pfac
     !   if h,e; if x,y,z
     !     s=
     !     pfac=
     !     d=
     !     size=

        call coeff(pfac, dx, s, alpha, beta, size(v,1)) 
        do k=1,size(v,3)-1
          do j=1,size(v,2)-1
            do i=1,size(v,1)-1
          v(i,j,k,1) = alpha(i) * v(i,j,k,1) + beta(i) *  &
              ( w(i+1,j,k,1) + w(i+1,j,k,2)  &
              - w(i,  j,k,1) - w(i,  j,k,2) )
            end do
          end do
        end do

      end subroutine

      subroutine add_diff_y(s, pfac, v, w)
        implicit none
        integer i,j,k
        real v(:,:,:,:)
        real w(size(v,1),size(v,2),size(v,3),2)
        real pfac, alpha(size(v,2)), beta(size(v,2)), s(*)

        call coeff(pfac, dy, s, alpha, beta, size(v,2)) 
        do k=1,size(v,3)-1
          do j=1,size(v,2)-1
            do i=1,size(v,1)-1
          v(i,j,k,1) = alpha(j) * v(i,j,k,1) + beta(j) *  &
              ( w(i,j+1,k,1) + w(i,j+1,k,2)  &
              - w(i,j,  k,1) - w(i,j,  k,2) )
            end do
          end do
        end do

      end subroutine

      subroutine add_diff_z(s, pfac, v, w)
        implicit none
        integer i,j,k
        real v(:,:,:,:)
        real w(size(v,1),size(v,2),size(v,3),2)
        real pfac, alpha(size(v,3)), beta(size(v,3)), s(*)

        call coeff(pfac, dz, s, alpha, beta, size(v,3)) 
        do k=1,size(v,3)-1
          do j=1,size(v,2)-1
            do i=1,size(v,1)-1
          v(i,j,k,1) = alpha(k) * v(i,j,k,1) + beta(k) *  &
              ( w(i,j,k+1,1) + w(i,j,k+1,2)  &
              - w(i,j,k  ,1) - w(i,j,k  ,2) )
            end do
          end do
        end do

      end subroutine

      subroutine coeff(pfac, d, s, alpha, beta, n)
        implicit none
        integer n, i
        real    :: pfac, d, s(n), alpha(n), beta(n)
        real*8  :: t

        do i=1,n
          t = exp(-s(i)*dt)
          alpha(i) = t
          if(t .gt. 1d0 - 1d-7) then
            beta(i) = dt/(pfac*d)
          else
            beta(i) = (1-t)/(pfac*s(i)*d)
          end if
        end do

      end subroutine

end program test