!****************************************************************************
! 電流源
!****************************************************************************
subroutine ecur_line_source
   use fdtd
   implicit none

   integer :: i,j,k,id
   real(8) :: sin_warming_up

   i=nx/2
   j=ny/2
   do k=2,nz-1
      id=media_id(i,j,k)

      ez(i,j,k)=ez(i,j,k) &
               -(dt/eps(id))/(1.0d0+(sig(id)*dt)/(2.0d0*eps(id)))*sin_warming_up(time-dt/2.0d0)
   end do

   return
end subroutine

subroutine ecur_infinitesimal_dipole_source
   use fdtd
   implicit none

   integer :: i,j,k,id
   real(8) :: sin_warming_up

   i=nx/2
   j=ny/2
   k=nz/2

   id=media_id(i,j,k)
   ez(i,j,k)=ez(i,j,k) &
            -(dt/eps(id))/(1.0d0+(sig(id)*dt)/(2.0d0*eps(id)))*sin_warming_up(time-dt/2.0d0)

   return
end subroutine

!----------------------------------------------------------------------------
! 平面波
!----------------------------------------------------------------------------
subroutine plane_wave_source
   use consts
   use fdtd
   implicit none

   integer :: i,j,k
   real(8) :: sin_warming_up

   j=3
   do k=1,nz-1
      do i=1,nx-1
         ez(i,j,k)=sin_warming_up(time-dt/2.0d0)
         hx(i,j,k)=ez(i,j,k)/z0
      end do
   end do

   return
end subroutine

!----------------------------------------------------------------------------
! 波源の関数形
! だんだん振幅が大きくなる正弦波(p.48)
!----------------------------------------------------------------------------
real(8) function sin_warming_up(t)
   use consts
   use fdtd
   implicit none

   real(8) :: t
   real(8) :: tp,t1

   tp=1.0d0/freq
   t1=4.0d0*tp
   if(t < t1) then
      ! ウォーミングアップ
      sin_warming_up=0.5d0*(1.0d0-dcos(pi*t/t1))*dsin(2.0d0*pi*freq*t)
   else
      ! 正弦波
      sin_warming_up=dsin(2.0d0*pi*freq*t)
   end if

end function

!
! End of file
!