gouy.f90

一个关于dip 的 源代码 适合 初学FDTD 的下载学习

program test_fdtd

  !-- Test code for FDTD + PML, version 1.02, Nov 1999
  !-- by Jos Bergervoet.

    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

  !-- For the Berenger method, each of the E and H fields has two species,
  !-- denoted here by the last index. This is strictly speaking only
  !-- needed in the boundary regions, but for convenience, here we do it
  !-- on the entire grid.

    character(len=20)  :: name
    integer            :: nx,ny,nz, nb, nt, tstep_gif, n
    real(wp)           :: Pi, mu0, Z0, eps0,  &
                          fun, g, g2, t, tpulse, smax, F_in, dx, dy, dz, dt

    g(t) = sin(t)**6                      ! smooth, 5 times differentiable
    fun(t) = g(min(t/tpulse,1.0_wp)*Pi)   ! scaled g(t) to [0,tpulse]

    print "(A)",  &
     "Give:  nx,ny,nz, nt, nb, smax, tpulse, tstep_gif",  &
     "",  &
     "nx,ny,nz, nt are space/time steps for the simulation,",  &
     "nb, smax are thickness and absorption parameter of Berenger layer,",  &
     "tpulse is the duration in timesteps of the pulse,",  &
     "tstep_gif is the number of timesteps between pictures.",  &
     "",  &
     "Example input: 60 60 60  200  10 4e8  20  5"

    read *, nx,ny,nz, nt, nb,smax, tpulse, tstep_gif

    call setup()
    
    n = 0
    Fdtd_loop: do
        n=n+1

        if (n<=tpulse) then        ! divide pulse over both E-field species
          F_in = fun(n+0.0_wp)
          ezs(nx/2,:,:,:) = F_in/2
        end if

        if (mod(n,tstep_gif) == 0) then     ! put this picture in a file
          write(name,*) n
          name = trim(adjustl(name)) // ".gif"
          call plot(name, ezs(:,:,nz/2,1)+ezs(:,:,nz/2,2))
        end if

        if (n>nt) exit Fdtd_loop

        call Berenger(exs,eys,ezs, hxs,hys,hzs)

        call PEC()    ! short circuit the field on all conductors
        
    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
        

  subroutine PEC()    ! set E-field (both species) zero on parabolic mirror
    integer  :: i, j, k

    do i=1,ny
      do j=1,nz
        k = nx+1 - ((i-dble(ny+1)/2)**2 + (j-dble(nz+1)/2)**2) / (2*nx)
        exs(k-1, i-1:i, j-1:j, :) = 0
        eys(k-1:k, i-1, j-1:j, :) = 0
        ezs(k-1:k, i-1:i, j-1, :) = 0
      end do
    end do

  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 )

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

  end subroutine
  

  subroutine Berenger(Ex,Ey,Ez,Hx,Hy,Hz)

  !-- Updates the fields NOT on the PEC faces of the domain, and
  !-- 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


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

  !-- Using exponential stepping, as described by Katz, Thiele en Tavlove
  !-- in IEEE microw., Vol 4, Nr. 8, p268(3), 1994

    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) = alpha(i) * v(i,j,k) + 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(wp) :: v(:,:,:), w(:,:,:,:)
    real(wp) :: 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) = alpha(j) * v(i,j,k) + 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(wp) :: v(:,:,:), w(:,:,:,:)
    real(wp) :: 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) = alpha(k) * v(i,j,k) + 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(wp) :: pfac, d, s(n), alpha(n), beta(n)
    real*8   :: t  ! double precision here, whatever the working precision

    do i=1,n
      t = exp(dble(-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_fdtd