poynting.f90

Sfdtd Simple finite-difference time-domain

! poynting.f90
!
! Save poynting vector of a cube faces with q^3 cells
!
!    Copyright (C) 2007  Paul Panserrieu, < [e-mail] >
!
!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License.
! 
! last modified: 17-09-2007 02:35:13 PM CEST

MODULE poynting

USE fdtd_gitter
USE plot, ONLY: toch

IMPLICIT NONE

TYPE surface

    DOUBLE PRECISION, POINTER, DIMENSION(:,:)           :: bas_x_ey, top_x_ey
    DOUBLE PRECISION, POINTER, DIMENSION(:,:)           :: bas_x_ez, top_x_ez
    DOUBLE PRECISION, POINTER, DIMENSION(:,:)           :: bas_y_ez, top_y_ez
    DOUBLE PRECISION, POINTER, DIMENSION(:,:)           :: bas_y_ex, top_y_ex
    DOUBLE PRECISION, POINTER, DIMENSION(:,:)           :: bas_z_ex, top_z_ex
    DOUBLE PRECISION, POINTER, DIMENSION(:,:)           :: bas_z_ey, top_z_ey

    DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)         :: bas_x_hz, top_x_hz
    DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)         :: bas_x_hy, top_x_hy
    DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)         :: bas_y_hx, top_y_hx
    DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)         :: bas_y_hz, top_y_hz
    DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)         :: bas_z_hy, top_z_hy
    DOUBLE PRECISION, POINTER, DIMENSION(:,:,:)         :: bas_z_hx, top_z_hx

    INTEGER                                             :: hb

END TYPE surface


CONTAINS 

SUBROUTINE init_poynting(hb, f)
 
  INTEGER, INTENT(IN)                :: hb ! halb breite
  TYPE(surface), INTENT(OUT)         :: f  
  
  f%hb = hb
  
  ! x plane E_y, E_z (y,z)
  ALLOCATE(  f%bas_x_ey(-hb:hb-1,-hb:hb-1)     &
            ,f%top_x_ey(-hb:hb-1,-hb:hb-1))
  f%bas_x_ey(-hb:hb-1,-hb:hb-1) = 0.0d0
  f%top_x_ey(-hb:hb-1,-hb:hb-1) = 0.0d0

  ALLOCATE(  f%bas_x_ez(-hb:hb-1,-hb:hb-1)     &
            ,f%top_x_ez(-hb:hb-1,-hb:hb-1))
  f%bas_x_ez(-hb:hb-1,-hb:hb-1) = 0.0d0
  f%top_x_ez(-hb:hb-1,-hb:hb-1) = 0.0d0

  ! y plane E_z, E_x (z,x)
  ALLOCATE(  f%bas_y_ez(-hb:hb-1,-hb:hb-1)     &
            ,f%top_y_ez(-hb:hb-1,-hb:hb-1))
  f%bas_y_ez(-hb:hb-1,-hb:hb-1) = 0.0d0
  f%top_y_ez(-hb:hb-1,-hb:hb-1) = 0.0d0
  ALLOCATE(  f%bas_y_ex(-hb:hb-1,-hb:hb-1)     &
            ,f%top_y_ex(-hb:hb-1,-hb:hb-1))
  f%bas_y_ex(-hb:hb-1,-hb:hb-1) = 0.0d0
  f%top_y_ex(-hb:hb-1,-hb:hb-1) = 0.0d0

  ! z plane E_x, E_y (x,y)
  ALLOCATE(  f%bas_z_ex(-hb:hb-1,-hb:hb-1)     &
            ,f%top_z_ex(-hb:hb-1,-hb:hb-1))
  f%bas_z_ex(-hb:hb-1,-hb:hb-1) = 0.0d0
  f%top_z_ex(-hb:hb-1,-hb:hb-1) = 0.0d0
  ALLOCATE(  f%bas_z_ey(-hb:hb-1,-hb:hb-1)     &
            ,f%top_z_ey(-hb:hb-1,-hb:hb-1))
  f%bas_z_ey(-hb:hb-1,-hb:hb-1) = 0.0d0
  f%top_z_ey(-hb:hb-1,-hb:hb-1) = 0.0d0

  ! x plane H_z, H_y (y,z)
  ALLOCATE(  f%bas_x_hz(-hb:hb-1,-hb:hb-1, 1:2)     &
            ,f%top_x_hz(-hb:hb-1,-hb:hb-1, 1:2))
  f%bas_x_hz(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  f%top_x_hz(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  ALLOCATE(  f%bas_x_hy(-hb:hb-1,-hb:hb-1, 1:2)     &
            ,f%top_x_hy(-hb:hb-1,-hb:hb-1, 1:2))
  f%bas_x_hy(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  f%top_x_hy(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0

  ! y plane H_x, H_z (x,z)
  ALLOCATE(  f%bas_y_hx(-hb:hb-1,-hb:hb-1, 1:2)     &
            ,f%top_y_hx(-hb:hb-1,-hb:hb-1, 1:2))
  f%bas_y_hx(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  f%top_y_hx(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  ALLOCATE(  f%bas_y_hz(-hb:hb-1,-hb:hb-1, 1:2)     &
            ,f%top_y_hz(-hb:hb-1,-hb:hb-1, 1:2))
  f%bas_y_hz(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  f%top_y_hz(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0

  ! z plane H_y, H_x (x,y)
  ALLOCATE(  f%bas_z_hy(-hb:hb-1,-hb:hb-1, 1:2)     &
            ,f%top_z_hy(-hb:hb-1,-hb:hb-1, 1:2))
  f%bas_z_hy(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  f%top_z_hy(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  ALLOCATE(  f%bas_z_hx(-hb:hb-1,-hb:hb-1, 1:2)     &
            ,f%top_z_hx(-hb:hb-1,-hb:hb-1, 1:2))
  f%bas_z_hx(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  f%top_z_hx(-hb:hb-1,-hb:hb-1, 1:2) = 0.0d0
  
END SUBROUTINE init_poynting

SUBROUTINE store_e_poynting(g,f)

  TYPE(gitter), INTENT(IN)                    :: g
  TYPE(surface), INTENT(INOUT)                :: f
  INTEGER                                     :: i,j

  DO i = -f%hb, f%hb-1, 1
    DO j = -f%hb, f%hb-1, 1
      f%bas_x_ey(i,j) = g%E(-f%hb, i, j, 2)
      f%top_x_ey(i,j) = g%E(f%hb,  i, j, 2)
      f%bas_x_ez(i,j) = g%E(-f%hb, i, j, 3)
      f%top_x_ez(i,j) = g%E(f%hb,  i, j, 3)

      f%bas_y_ez(i,j) = g%E(i, -f%hb, j, 3)
      f%top_y_ez(i,j) = g%E(i,  f%hb, j, 3)
      f%bas_y_ex(i,j) = g%E(i, -f%hb, j, 1)
      f%top_y_ex(i,j) = g%E(i,  f%hb, j, 1)

      f%bas_z_ex(i,j) = g%E(i, j, -f%hb, 1)
      f%bas_z_ex(i,j) = g%E(i, j,  f%hb, 1)
      f%bas_z_ey(i,j) = g%E(i, j, -f%hb, 2)
      f%bas_z_ey(i,j) = g%E(i, j,  f%hb, 2)
    ENDDO
  ENDDO

END SUBROUTINE store_e_poynting 

SUBROUTINE store_h_poynting(g, f)

  TYPE(gitter), INTENT(IN)                    :: g
  TYPE(surface), INTENT(INOUT)                :: f
  INTEGER                                     :: i,j

  DO i = -f%hb, f%hb-1, 1
    DO j = -f%hb, f%hb-1, 1
      f%bas_x_hy(i,j,2) = f%bas_x_hy(i,j,1)
      f%top_x_hy(i,j,2) = f%top_x_hy(i,j,1)
      f%bas_x_hz(i,j,2) = f%bas_x_hz(i,j,1)
      f%top_x_hz(i,j,2) = f%top_x_hz(i,j,1)

      f%bas_y_hz(i,j,2) = f%bas_y_hz(i,j,1)
      f%top_y_hz(i,j,2) = f%top_y_hz(i,j,1)
      f%bas_y_hx(i,j,2) = f%bas_y_hx(i,j,1)
      f%top_y_hx(i,j,2) = f%top_y_hx(i,j,1)

      f%bas_z_hx(i,j,2) = f%bas_z_hx(i,j,1)
      f%bas_z_hx(i,j,2) = f%bas_z_hx(i,j,1)
      f%bas_z_hy(i,j,2) = f%bas_z_hy(i,j,1)
      f%bas_z_hy(i,j,2) = f%bas_z_hy(i,j,1)      
    ENDDO
  ENDDO

  DO i = -f%hb, f%hb-1, 1
    DO j = -f%hb, f%hb-1, 1
      f%bas_x_hy(i,j,1) = g%H(-f%hb, i, j, 2) + g%H(-f%hb-1, i, j, 2)
      f%top_x_hy(i,j,1) = g%H( f%hb, i, j, 2) + g%H( f%hb-1, i, j, 2)
      f%bas_x_hz(i,j,1) = g%H(-f%hb, i, j, 3) + g%H(-f%hb-1, i, j, 3)
      f%top_x_hz(i,j,1) = g%H( f%hb, i, j, 3) + g%H( f%hb-1, i, j, 3)

      f%bas_y_hz(i,j,1) = g%H(i, -f%hb, j, 3) + g%H(i, -f%hb-1, j, 3)
      f%top_y_hz(i,j,1) = g%H(i,  f%hb, j, 3) + g%H(i,  f%hb-1, j, 3)
      f%bas_y_hx(i,j,1) = g%H(i, -f%hb, j, 1) + g%H(i, -f%hb-1, j, 1)
      f%top_y_hx(i,j,1) = g%H(i,  f%hb, j, 1) + g%H(i,  f%hb-1, j, 1)

      f%bas_z_hx(i,j,1) = g%H(i, j, -f%hb, 1) + g%H(i, j, -f%hb-1, 1)
      f%bas_z_hx(i,j,1) = g%H(i, j,  f%hb, 1) + g%H(i, j,  f%hb-1, 1)
      f%bas_z_hy(i,j,1) = g%H(i, j, -f%hb, 2) + g%H(i, j, -f%hb-1, 2)
      f%bas_z_hy(i,j,1) = g%H(i, j,  f%hb, 2) + g%H(i, j,  f%hb-1, 2)     
    ENDDO
  ENDDO

END SUBROUTINE store_h_poynting

SUBROUTINE output_poynting(do_poynting, g, f, time_step)

  INTEGER, INTENT(IN)                         :: do_poynting
  TYPE(gitter), INTENT(IN)                    :: g
  TYPE(surface), INTENT(IN)                   :: f
  INTEGER, INTENT(IN)                         :: time_step
  CHARACTER(5)                                :: zahl
  CHARACTER, DIMENSION(1:5)                   :: x
  INTEGER                                     :: i, j 
  DOUBLE PRECISION                            :: value

  IF (do_poynting .EQ. 1) THEN
    CALL toch(x, time_step)
    zahl = x(1)//x(2)//x(3)//x(4)//x(5)
    OPEN(10,FILE='poynting_'//g%id//'_'//zahl//'.data', ACTION='WRITE')
    value = 0.0d0
    WRITE(10,*)

    DO i = -f%hb, f%hb-1, 1
      DO j = -f%hb, f%hb-1, 1
        ! x faces
        value = value + (  f%bas_x_ey(i,j) * 0.25d0 * (f%bas_x_hz(i,j,1) + f%bas_x_hz(i,j,2)) &
                         - f%bas_x_ez(i,j) * 0.25d0 * (f%bas_x_hy(i,j,1) + f%bas_x_hy(i,j,2)) &
                        )
        value = value + (  f%top_x_ey(i,j) * 0.25d0 * (f%top_x_hz(i,j,1) + f%top_x_hz(i,j,2)) &
                         - f%top_x_ez(i,j) * 0.25d0 * (f%top_x_hy(i,j,1) + f%top_x_hy(i,j,2)) &
                        )
        value = value + (  f%bas_y_ez(i,j) * 0.25d0 * (f%bas_y_hx(i,j,1) + f%bas_y_hx(i,j,2)) &
                         - f%bas_y_ex(i,j) * 0.25d0 * (f%bas_y_hz(i,j,1) + f%bas_y_hz(i,j,2)) &
                        )
        value = value + (  f%top_y_ez(i,j) * 0.25d0 * (f%top_y_hx(i,j,1) + f%top_y_hx(i,j,2)) &
                         - f%top_y_ex(i,j) * 0.25d0 * (f%top_y_hz(i,j,1) + f%top_y_hz(i,j,2)) &
                        )
        value = value + (  f%bas_x_ez(i,j) * 0.25d0 * (f%bas_z_hy(i,j,1) + f%bas_z_hy(i,j,2)) &
                         - f%bas_y_ex(i,j) * 0.25d0 * (f%bas_z_hx(i,j,1) + f%bas_z_hx(i,j,2)) &
                        )
        value = value + (  f%top_x_ez(i,j) * 0.25d0 * (f%top_z_hy(i,j,1) + f%top_z_hy(i,j,2)) &
                         - f%top_y_ex(i,j) * 0.25d0 * (f%top_z_hx(i,j,1) + f%top_z_hx(i,j,2)) &
                        )
      ENDDO
    ENDDO
    WRITE(10,*) time_step, value * g%dx * g%dx
    CLOSE(10)
  ENDIF

END SUBROUTINE output_poynting

END MODULE poynting