test.f90

Sfdtd Simple finite-difference time-domain

! test.f90
! 
!    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: 22-05-2007 12:29:38 PM CEST

MODULE test

USE plot
USE pml
USE fdtd_gitter

IMPLICIT NONE

CHARACTER(10), PARAMETER                                :: espace = '          '

DOUBLE PRECISION, PARAMETER                             :: PREC   =  1.0d-17

  
CONTAINS

SUBROUTINE info_limits(device_id)

    INTEGER, INTENT(IN)          :: device_id
    DOUBLE PRECISION             :: pi_approx

    pi_approx = 3.14159265358979323846

    WRITE(device_id,*)  "------------------------------------------"
    WRITE(device_id,*)  "numerical limits: \n"
    WRITE(device_id,*)  "CHARACTER is KIND(' ')  = ", KIND(' '), " Byte"
    WRITE(device_id,*)  "BOOLEAN is KIND(.TRUE.) = ", KIND(.TRUE.), " Byte"

    WRITE(device_id,*)  "MIN INTEGER             = ", -HUGE(0)
    WRITE(device_id,*)  "MAX INTEGER             = ", HUGE(0)
    WRITE(device_id,*)  "INTEGER is KIND(0)      = ", KIND(0), " Byte"
    WRITE(device_id,*)  "MAX REAL                = ", HUGE(0.0)
    WRITE(device_id,*)  "MIN REAL                = ", TINY(0.0)
    WRITE(device_id,*)  "REAL  is KIND(0.0)      = ", KIND(0.0), " Byte"
    WRITE(device_id,*)  "MAX DOUBLE              = ", HUGE(0.0d0)
    WRITE(device_id,*)  "MIN DOUBLE              = ", TINY(0.0d0)
    WRITE(device_id,*)  "DOUBLE is KIND(0.0d0)   = ", KIND(0.0d0), " Byte"
    WRITE(device_id,*)  "epsilon(1.0)              = ", EPSILON(1.0)
    WRITE(device_id,*)  "-------------------------------------------"
    WRITE(device_id,*)  "\n cos() and sin(): \n"
    WRITE(device_id,*)  "pi approx     = ", pi_approx
    WRITE(device_id,*)  "cos(2*pi)     = ", DCOS(2.0d0 * pi_approx)
    WRITE(device_id,*)  "cos(2.0E6*pi) = ", DCOS(2.0E6 * pi_approx)
    WRITE(device_id,*)  "cos(2.0E9*pi) = ", DCOS(2.0E9 * pi_approx)

    WRITE(device_id,*)  "sin(2*pi)     = ", DSIN(2.0d0 * pi_approx)
    WRITE(device_id,*)  "sin(2.0E6*pi) = ", DSIN(2.0E6 * pi_approx)
    WRITE(device_id,*)  "sin(2.0E9*pi) = ", DSIN(2.0E9 * pi_approx)
    WRITE(device_id,*)  "------------------------------------------"

END SUBROUTINE info_limits

! debug
SUBROUTINE clip(what, wert, lim, wo, i1, i2, i3, re_val)

  CHARACTER(len=6)                               :: what
  DOUBLE PRECISION, INTENT(IN)                   :: wert, lim
  INTEGER, INTENT(IN)                            :: wo, i1, i2, i3
  INTEGER, INTENT(INOUT)                          :: re_val

  IF (wert > ABS(lim)) THEN
    WRITE(*,*) '* ', what, ':', wert, ' lim:', lim
    WRITE(*,*) '  - Where: ', i1, i2, i3
    WRITE(*,*) '  - Flag: ', wo
    re_val = 1
  ELSE
    re_val = 0
  ENDIF

END SUBROUTINE clip

SUBROUTINE print_cell(g, x_pos, y_pos, z_pos, timestep)

  TYPE(gitter), INTENT(IN)       :: g
  INTEGER, INTENT(IN)            :: x_pos
  INTEGER, INTENT(IN)            :: y_pos
  INTEGER, INTENT(IN)            :: z_pos
  INTEGER, INTENT(IN)            :: timestep

  WRITE(*,*)
  WRITE(*,*) 'cell:(', x_pos, y_pos, z_pos, ')', 't=', timestep 
  WRITE(*,*) '    E(', g%E(x_pos, y_pos, z_pos, :), ')'
  WRITE(*,*) '    H(', g%H(x_pos, y_pos, z_pos, :), ')'
  WRITE(*,*)

END SUBROUTINE print_cell

SUBROUTINE print_pml_cell(b, typ, x_pos, y_pos, z_pos, timestep)

  TYPE(pml_boundary), INTENT(IN)       :: b
  INTEGER, INTENT(IN)                  :: typ ! 1:bas_x, 2:top_x, 3:bas_y, 4:top_y, 5:bas_z and 6:top_z
  INTEGER, INTENT(IN)                  :: x_pos
  INTEGER, INTENT(IN)                  :: y_pos
  INTEGER, INTENT(IN)                  :: z_pos
  INTEGER, INTENT(IN)                  :: timestep

  WRITE(*,*)
  WRITE(*,*) 'pml cell:(', x_pos, y_pos, z_pos, ')', 't=', timestep
  SELECT CASE (typ)
    CASE(1)
      WRITE(*,*) '    bas_x E(', b%bas_x%E(x_pos, y_pos, z_pos, :), ')'
      WRITE(*,*) '    bas_x H(', b%bas_x%H(x_pos, y_pos, z_pos, :), ')'
    CASE(2)
      WRITE(*,*) '    top_x E(', b%top_x%E(x_pos, y_pos, z_pos, :), ')'
      WRITE(*,*) '    top_x H(', b%top_x%H(x_pos, y_pos, z_pos, :), ')'
    CASE(3)
      WRITE(*,*) '    bas_y E(', b%bas_y%E(x_pos, y_pos, z_pos, :), ')'
      WRITE(*,*) '    bas_y H(', b%bas_y%H(x_pos, y_pos, z_pos, :), ')'
    CASE(4)
      WRITE(*,*) '    top_y E(', b%top_y%E(x_pos, y_pos, z_pos, :), ')'
      WRITE(*,*) '    top_y H(', b%top_y%H(x_pos, y_pos, z_pos, :), ')'
    CASE(5)
      WRITE(*,*) '    bas_z E(', b%bas_z%E(x_pos, y_pos, z_pos, :), ')'
      WRITE(*,*) '    bas_z H(', b%bas_z%H(x_pos, y_pos, z_pos, :), ')'
    CASE(6)
      WRITE(*,*) '    top_z E(', b%top_z%E(x_pos, y_pos, z_pos, :), ')'
      WRITE(*,*) '    top_z H(', b%top_z%H(x_pos, y_pos, z_pos, :), ')'
    CASE DEFAULT
      WRITE(*,*) "??"
  END SELECT
  WRITE(*,*)

END SUBROUTINE print_pml_cell

SUBROUTINE store_max(wert, maxi, x, y, z, ix,  iy, iz)

  DOUBLE PRECISION, INTENT(IN)                   :: wert
  DOUBLE PRECISION, INTENT(INOUT)                :: maxi
  INTEGER, INTENT(IN)                            :: ix,  iy, iz
  INTEGER, INTENT(INOUT)                         :: x, y, z

  IF (ABS(wert) > ABS(maxi)) THEN
    maxi = wert
    x = ix; y = iy; z = iz;
  ENDIF

END SUBROUTINE store_max 

SUBROUTINE print_max(maxi, comment, x, y, z)

  DOUBLE PRECISION, INTENT(IN)                :: maxi
  CHARACTER(10), INTENT(IN)                   :: comment
  INTEGER, INTENT(IN)                         :: x, y, z

  WRITE(*,*) maxi, comment, x, y, z

END SUBROUTINE print_max

SUBROUTINE plot_anregung()

  INTEGER                                            :: i
  DOUBLE PRECISION, DIMENSION(1:3, 0:100)            :: curves

  DO i = 0, 100, 1 
    curves(1,i) = COS( i * 2.0d0 * PI / 100.0d0)
    curves(2,i) = SIN( i * 2.0d0 * PI / 100.0d0)
    curves(3,i) = SIN( i * 2.0d0 * PI / 100.0d0) * (1.0d0 - COS( i * 2.0d0 * PI / 100.0d0))
  ENDDO
  
  CALL verlauf(curves(1,:), 101, '  ', '  ', 987, 'anregungdip_____', "Cosinus Anregung    "//espace//espace)
  CALL verlauf(curves(2,:), 101, '  ', '  ', 986, 'anregungdip_____', "Sinus Anregung      "//espace//espace)
  CALL verlauf(curves(3,:), 101, '  ', '  ', 985, 'anregungdip_____', "sin(x)(1-cos(x))    "//espace//espace)

END SUBROUTINE plot_anregung

SUBROUTINE init_debug_log(logname, device_id)

  INTEGER, INTENT(IN)                      :: device_id
  CHARACTER(len=11), INTENT(INOUT)         :: logname      

  logname = 'simu__sfdtd'
  OPEN(device_id, FILE= logname//'.debug.log', ACTION='WRITE')

END SUBROUTINE init_debug_log

SUBROUTINE end_debug_log(device_id)

  INTEGER, INTENT(IN)         :: device_id
  CLOSE(device_id)

END SUBROUTINE end_debug_log

END MODULE test