messung.f90

Sfdtd Simple finite-difference time-domain

! messung.f90
!
! Punktmessung 
!
!    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: 12-09-2007 03:34:37 PM CEST

MODULE messung

USE analytic
USE plot

IMPLICIT NONE

! (for debugging propose)
! save analytic and numeric curves separately 
INTEGER, PARAMETER                                           :: separate_plot = 0

TYPE charac
  DOUBLE PRECISION, POINTER, DIMENSION(:,:)                  :: E, H
END TYPE charac

TYPE error
  DOUBLE PRECISION, POINTER, DIMENSION(:,:)                  :: e
END TYPE error

CHARACTER(10), PARAMETER        :: empty        = "           "
CHARACTER(40), PARAMETER        :: verlauf_ex_f = "E_x fdtd   "//empty//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_ey_f = "E_y fdtd   "//empty//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_ez_f = "E_z fdtd   "//empty//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_hx_f = "H_x fdtd   "//empty//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_hy_f = "H_y fdtd   "//empty//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_hz_f = "H_z fdtd   "//empty//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_ex_a = "E_x analytic         "//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_ey_a = "E_y analytic         "//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_ez_a = "E_z analytic         "//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_hx_a = "H_x analytic         "//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_hy_a = "H_y analytic         "//empty//empty
CHARACTER(40), PARAMETER        :: verlauf_hz_a = "H_z analytic         "//empty//empty

CHARACTER(40), PARAMETER        :: error_ex     = "error: (E_x fdtd - E_x analytic) [V/m]   "
CHARACTER(40), PARAMETER        :: error_ey     = "error: (E_y fdtd - E_y analytic) [V/m]   "
CHARACTER(40), PARAMETER        :: error_ez     = "error: (E_z fdtd - E_z analytic) [V/m]   "
CHARACTER(40), PARAMETER        :: error_hx     = "error: (H_x fdtd - H_x analytic) [V/m]   "
CHARACTER(40), PARAMETER        :: error_hy     = "error: (H_y fdtd - H_y analytic) [V/m]   "

CHARACTER(40), PARAMETER        :: ver_ex       = "E_x fdtd und E_x analytic [V/m]          "
CHARACTER(40), PARAMETER        :: ver_ey       = "E_y fdtd und E_y analytic [V/m]          "
CHARACTER(40), PARAMETER        :: ver_ez       = "E_z fdtd und E_z analytic [V/m]          "
CHARACTER(40), PARAMETER        :: ver_hx       = "H_x fdtd und H_x analytic [V/m]          "
CHARACTER(40), PARAMETER        :: ver_hy       = "H_y fdtd und H_y analytic [V/m]          "
CHARACTER(40), PARAMETER        :: ver_hz       = "H_z fdtd und H_z analytic [V/m]          "

CONTAINS

SUBROUTINE messung_check(messzelle, n)
  
  INTEGER, INTENT(IN), DIMENSION(1:3)                        :: messzelle
  INTEGER, INTENT(IN)                                        :: n
  INTEGER                                                    :: i

  DO i = 1, 3, 1
    IF (messzelle(i) < (-n/2) .OR. messzelle(i) > (n/2)) THEN
      WRITE(*,*) "Messpunkt ausserhalb der Simulationsbox"
    ENDIF
  ENDDO

END SUBROUTINE messung_check

SUBROUTINE init_messung(messzelle, ch, r_err, nt, do_messung)

  INTEGER, INTENT(IN), DIMENSION(1:3)                        :: messzelle
  TYPE(charac), INTENT(INOUT)                                :: ch
  TYPE(error), INTENT(INOUT)                                 :: r_err
  INTEGER, INTENT(IN)                                        :: nt
  INTEGER, INTENT(IN)                                        :: do_messung

  IF (do_messung .EQ. 1) THEN
    ALLOCATE(ch%E(1:6, 0:nt), ch%H(1:6, 0:nt))
    ch%E = 0.0d0; ch%H = 0.0d0;
    ALLOCATE(r_err%e(0:nt, 1:5))
    r_err%e = 0.0d0
  ENDIF
END SUBROUTINE init_messung

SUBROUTINE store_timestep(g, d, messzelle, ch, r_err, it, do_messung, dipol_type, debug)

    TYPE(charac), INTENT(INOUT)                              :: ch
    TYPE(error), INTENT(INOUT)                               :: r_err
    TYPE(gitter), INTENT(IN)                                 :: g
    TYPE(dipol), INTENT(IN)                                  :: d
    INTEGER, INTENT(IN)                                      :: it
    INTEGER, INTENT(IN)                                      :: debug
    INTEGER, INTENT(IN)                                      :: do_messung 
    INTEGER, INTENT(IN)                                      :: dipol_type
    INTEGER, INTENT(IN), DIMENSION(1:3)                      :: messzelle 
    DOUBLE PRECISION, DIMENSION(1:5,1:3)                     :: ableitungen 
    DOUBLE PRECISION, DIMENSION(1:6)                         :: feldarray
    DOUBLE PRECISION                                         :: zt_E, zt_H
    
    IF (do_messung .EQ. 1) THEN
      zt_E = g%dt * (it) 
      zt_H = g%dt * (it) - g%dt * 0.5d0

      CALL stimulus(ableitungen, messzelle, g, d, dipol_type, zt_E, zt_H)
      CALL feldberechnung(feldarray, ableitungen, messzelle, zt_E, d, g)
      ch%E(1:5:2,it) = g%E(messzelle(1), messzelle(2), messzelle(3), 1:3) 
      ch%E(2:6:2,it) = feldarray(1:3)                               
      ch%H(1:5:2,it) = g%H(messzelle(1), messzelle(2), messzelle(3), 1:3)  
      ch%H(2:6:2,it) = feldarray(4:6)
      IF (debug .NE. 0) THEN
        WRITE(*,*) '(', it, ')     E_analy = ', feldarray(1:3)
        WRITE(*,*) '(', it, ')     E_num   = ', ch%E(1:5:2,it)
        WRITE(*,*)
        WRITE(*,*) '(', it, ')     H_analy = ', feldarray(4:6)
        WRITE(*,*) '(', it, ')     H_num   = ', ch%H(1:5:2,it)
        WRITE(*,*)
        ! absolute error
        ! numeric - analytic 
        r_err%e(it,1:3) = ch%E(1:5:2,it) - feldarray(1:3)
        r_err%e(it,4) = ch%H(1,it) - feldarray(4)
        r_err%e(it,5) = ch%H(3,it) - feldarray(5)
      ENDIF
    ENDIF

END SUBROUTINE store_timestep 

SUBROUTINE save_messung(do_messung, ch, r_err, nt, logname)
  
  TYPE(charac), INTENT(IN)                                   :: ch
  INTEGER, INTENT(IN)                                        :: do_messung 
  TYPE(error), INTENT(INOUT)                                 :: r_err
  INTEGER, INTENT(IN)                                        :: nt
  CHARACTER(len=6), INTENT(IN)                               :: logname
  INTEGER                                                    :: i

  IF (do_messung .EQ. 1) THEN
    CALL vergleich(ch%E(2,:), ch%E(1,:), nt+1, 't ', 'Ex', 101, logname,  ver_ex)
    CALL vergleich(ch%E(4,:), ch%E(3,:), nt+1, 't ', 'Ey', 102 , logname, ver_ey)
    CALL vergleich(ch%E(6,:), ch%E(5,:), nt+1, 't ', 'Ez', 103 , logname, ver_ez)
    CALL vergleich(ch%H(2,:), ch%H(1,:), nt+1, 't ', 'Hx', 201 , logname, ver_hx)
    CALL vergleich(ch%H(4,:), ch%H(3,:), nt+1, 't ', 'Hy', 202 , logname, ver_hy)
    CALL vergleich(ch%H(6,:), ch%H(5,:), nt+1, 't ', 'Hz', 203 , logname, ver_hz)
    CALL verlauf(r_err%e(:,1), nt+1, 't ', 'er', 901, logname, error_ex)
    CALL verlauf(r_err%e(:,2), nt+1, 't ', 'er', 902, logname, error_ey)
    CALL verlauf(r_err%e(:,3), nt+1, 't ', 'er', 903, logname, error_ez)
    CALL verlauf(r_err%e(:,4), nt+1, 't ', 'er', 904, logname, error_hx)
    CALL verlauf(r_err%e(:,5), nt+1, 't ', 'er', 905, logname, error_hy)

    IF (separate_plot .EQ. 1) THEN
      CALL verlauf(ch%E(1,:), nt+1, 't ', 'Ex', 301, logname, verlauf_ex_f)
      CALL verlauf(ch%E(3,:), nt+1, 't ', 'Ey', 302, logname, verlauf_ey_f)
      CALL verlauf(ch%E(5,:), nt+1, 't ', 'Ez', 303, logname, verlauf_ez_f)
      CALL verlauf(ch%H(1,:), nt+1, 't ', 'Hx', 401, logname, verlauf_hx_f)
      CALL verlauf(ch%H(3,:), nt+1, 't ', 'Hy', 402, logname, verlauf_hy_f)
      CALL verlauf(ch%H(5,:), nt+1, 't ', 'Hz', 403, logname, verlauf_hz_f)
      CALL verlauf(ch%E(2,:), nt+1, 't ', 'Ex', 501, logname, verlauf_ex_a)
      CALL verlauf(ch%E(4,:), nt+1, 't ', 'Ey', 502, logname, verlauf_ey_a)
      CALL verlauf(ch%E(6,:), nt+1, 't ', 'Ez', 503, logname, verlauf_ez_a)
      CALL verlauf(ch%H(2,:), nt+1, 't ', 'Hx', 601, logname, verlauf_hx_a)
      CALL verlauf(ch%H(4,:), nt+1, 't ', 'Hy', 602, logname, verlauf_hy_a)
      CALL verlauf(ch%H(6,:), nt+1, 't ', 'Hz', 603, logname, verlauf_hz_a)
    ENDIF
  ENDIF

END SUBROUTINE save_messung

END MODULE messung