plot.f90

Sfdtd Simple finite-difference time-domain

! plot.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: 26-09-2007 11:06:03 AM CEST

MODULE plot

USE fdtd_gitter
USE analytic, ONLY:find_ampl

IMPLICIT NONE

TYPE planeplot

  INTEGER                 :: sorte
  CHARACTER(len=1)        :: ebene1, ebene2
  INTEGER                 :: nr1, nr2

END TYPE planeplot

CONTAINS

SUBROUTINE set_scale_factor(factor, g, d, model, cf)
  
  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2)          :: factor
  TYPE(gitter), INTENT(IN)                                 :: g
  TYPE(dipol), INTENT(IN)                                  :: d
  INTEGER, INTENT(IN)                                      :: model, cf
  DOUBLE PRECISION                                         :: ampl
  
    CALL find_ampl(d%E, ampl)
  
    SELECT CASE (model)
      CASE(1)
        IF (cf .EQ. 0 .OR. cf .EQ. 2) THEN
          factor(1) = g%dx * 10 * (-g%nxl + g%nxyh) / ampl
        ENDIF
        IF (cf .EQ. 1 .OR. cf .EQ. 2) THEN
          factor(2) = g%dx * 1000 * (-g%nxl + g%nxyh) / ampl
        ENDIF
      CASE(2)
          factor(1) = 20000.0 / ampl
          factor(2) = 6000000.0 / ampl
      CASE(3)
          factor(1) = 20000.0 / ampl
          factor(2) = 6000000.0 / ampl
      CASE DEFAULT
        factor = 1.0 / ampl
    END SELECT

END SUBROUTINE set_scale_factor


SUBROUTINE run_plot(g, p1, p2, it, n, factor, do_vect_plot, do_density_plot, component)

  TYPE(gitter), INTENT(IN)                                 :: g
  TYPE(planeplot), INTENT(IN)                              :: p1, p2
  INTEGER, INTENT(IN)                                      :: it
  INTEGER, INTENT(IN)                                      :: n
  DOUBLE PRECISION, INTENT(IN), DIMENSION(1:2)             :: factor
  INTEGER, INTENT(IN)                                      :: do_vect_plot
  INTEGER, INTENT(IN)                                      :: do_density_plot
  INTEGER, INTENT(IN)                                      :: component  ! 1=x, 2=y , 3=z  oder 4 fuer alle
  INTEGER                                                  :: i


  IF (do_vect_plot .EQ. 1) THEN
    IF (p1%sorte .EQ. 0) THEN
      CALL ausgabe_feldschnitt(g, it, p1%nr1, p1%ebene1, 'e', n, factor(1))
    ELSEIF (p1%sorte .EQ. 1) THEN
      CALL ausgabe_feldschnitt(g, it, p1%nr1, p1%ebene1, 'h', n, factor(2))
    ELSE
      CALL ausgabe_feldschnitt(g, it, p1%nr1, p1%ebene1, 'e', n, factor(1))
      CALL ausgabe_feldschnitt(g, it, p1%nr2, p1%ebene2, 'h', n, factor(2))
    ENDIF
  ENDIF
  IF (do_density_plot .EQ. 1) THEN
    IF (p2%sorte .EQ. 0) THEN
      IF (component .NE. 4) THEN
        CALL density_plot(g, it, p2%nr1, p2%ebene1, 'e', n, component)
      ELSE
        DO i = 1, 3, 1
          CALL density_plot(g, it, p2%nr1, p2%ebene1, 'e', n, i)
        ENDDO
      ENDIF
    ELSEIF (p2%sorte .EQ. 1) THEN
      IF (component .NE. 4) THEN
        CALL density_plot(g, it, p2%nr1, p2%ebene1, 'h', n, component)
      ELSE
        DO i = 1, 3, 1
          CALL density_plot(g, it, p2%nr1, p2%ebene1, 'h', n, i)
        ENDDO
      ENDIF
    ELSE
      IF (component .NE. 4) THEN
        CALL density_plot(g, it, p2%nr1, p2%ebene1, 'e', n, component)
      ELSE
        DO i = 1, 3, 1
          CALL density_plot(g, it, p2%nr1, p2%ebene1, 'e', n, i)
        ENDDO
      ENDIF
      IF (component .NE. 4) THEN
        CALL density_plot(g, it, p2%nr2, p2%ebene2, 'h', n, component)
      ELSE
        DO i = 1, 3, 1
          CALL density_plot(g, it, p2%nr2, p2%ebene2, 'h', n, i)
        ENDDO
      ENDIF
    ENDIF
  ENDIF

END SUBROUTINE run_plot

SUBROUTINE verlauf(array, length, absis, ord, nbr, filename, title)

  DOUBLE PRECISION, INTENT(IN), DIMENSION(1)         :: array
  INTEGER, INTENT(IN)                                :: length, nbr
  CHARACTER(len=6), INTENT(IN)                       :: filename
  CHARACTER(len=40), INTENT(IN)                      :: title
  CHARACTER(2), INTENT(IN)                           :: absis, ord
  CHARACTER, DIMENSION(1:5)                          :: x
  INTEGER                                            :: i
  
  CALL toch(x,nbr)
  OPEN(20,FILE='verlauf2d_'//filename//x(1)//x(2)//x(3)//x(4)//x(5)//'.plt',ACTION='WRITE')

  WRITE(20,*) '$ DATA = CURVE2D'
  WRITE(20,*)
  WRITE(20,*)
  WRITE(20,*) '% toplabel= "'//title//' "'
  WRITE(20,*) '% xlabel= "', absis, '"' 
  WRITE(20,*) '% ylabel= "', ord, '"'
  WRITE(20,*) '% lc = 0 linetype = 1'
    
  DO i=1, length, 1
    WRITE(20,*) i-1, array(i)
  ENDDO

  WRITE(20,*) '$ END'
  CLOSE(20)

END SUBROUTINE verlauf

SUBROUTINE vergleich(array1, array2, length, absis, ord, nbr, filename, title)

  DOUBLE PRECISION, INTENT(IN), DIMENSION(1)         :: array1, array2
  INTEGER, INTENT(IN)                                :: length, nbr
  CHARACTER(2), INTENT(IN)                           :: absis, ord
  CHARACTER(len=6), INTENT(IN)                       :: filename
  CHARACTER(len=40), INTENT(IN)                      :: title
  CHARACTER, DIMENSION(1:5)                          :: x
  INTEGER                                            :: i
  
  CALL toch(x,nbr)
  OPEN(20,FILE = 'vergleich_'//filename//x(1)//x(2)//x(3)//x(4)//x(5)//'.plt',ACTION='WRITE')

  WRITE(20,*) '$ DATA = CURVE2D'
  WRITE(20,*)
  WRITE(20,*) '% toplabel= "'//title//' "'
  WRITE(20,*) '% xlabel= "', absis , '"'
  WRITE(20,*) '% ylabel= "', ord, '"'
  WRITE(20,*) '% linelabel  = analytisch'
  WRITE(20,*) '% linecolor = 4 linetype = 2'
    
  DO i=1, length, 1
    WRITE(20,*) i-1, array1(i)
  ENDDO

  WRITE(20,*)
  WRITE(20,*) '% linelabel  = numerisch'
  WRITE(20,*) '% lc = 0 linetype = 1'

  DO i=1, length, 1
    WRITE(20,*) i-1, array2(i)
  ENDDO

  WRITE(20,*) '$ END'
  CLOSE(20)

END SUBROUTINE vergleich

SUBROUTINE ausgabe_feldschnitt(g, t, ne, co, cf, n, factor)

  TYPE(gitter), INTENT(IN)                         :: g
  INTEGER, INTENT(IN)                              :: t    ! Zeitschritt
  INTEGER, INTENT(IN)                              :: ne   ! Gitterebene
  INTEGER, INTENT(IN)                              :: n  
  CHARACTER, INTENT(IN)                            :: co   ! Orientierung der Schnittflaeche
  CHARACTER, INTENT(IN)                            :: cf   ! E- oder H-Feld
  DOUBLE PRECISION, INTENT(IN)                     :: factor
  CHARACTER, DIMENSION(1:5)                        :: x
  INTEGER                                          :: i1, i2

  CALL toch(x,t)
  OPEN(10,FILE='feldbild_'//g%id//'_'//co//'_'//cf//x(1)//x(2)//x(3)//x(4)//x(5)//'.plt',ACTION='WRITE')
  WRITE(10,*) '$ DATA=VECTOR'
  WRITE(10,*) '% toplabel = " Feldschnitt von '//cf//co//'"'   
  WRITE(10,*) '% subtitle   = "Zeitschritte: ', t, '"'
  WRITE(10,*) '% vscale=', factor

  DO i1=g%nyl, g%nyyh, 1
    DO i2=g%nzl, g%nzyh, 1
      
      IF( co .EQ. 'x' ) THEN
        IF(i1 .EQ. g%nxl .AND.  i2 .EQ. g%nxl) THEN
          WRITE(10,*) '% xlabel="y"'
          WRITE(10,*) '% ylabel="z"'
        ENDIF
        IF( cf .EQ. 'E'  .OR. cf .EQ. 'e') THEN
          WRITE(10,*) i1*g%dy, i2*g%dz, 0.0, g%E(ne,i1,i2,2), g%E(ne,i1,i2,3), 0.0
        ELSE
          WRITE(10,*) i1*g%dy, i2*g%dz, 0.0, g%H(ne,i1,i2,2), g%H(ne,i1,i2,3), 0.0
        ENDIF
      ENDIF      
      IF( co .EQ. 'y' ) THEN
        IF(i1 .EQ. g%nxl .AND.  i2 .EQ. g%nxl) THEN
          WRITE(10,*) '% xlabel="x"'
          WRITE(10,*) '% ylabel="z"'
        ENDIF
        IF( cf .EQ. 'E'  .OR. cf .EQ. 'e') THEN
          WRITE(10,*) i1*g%dx, i2*g%dz, 0.0,g%E(i1,ne,i2,1), g%E(i1,ne,i2,3), 0.0
        ELSE
          WRITE(10,*) i1*g%dx, i2*g%dz, 0.0, g%H(i1,ne,i2,1), g%H(i1,ne,i2,3), 0.0
        ENDIF
      ENDIF
      IF( co .EQ. 'z' ) THEN
        IF(i1 .EQ. g%nxl .AND.  i2 .EQ. g%nxl) THEN
          WRITE(10,*) '% xlabel="x"'
          WRITE(10,*) '% ylabel="y"'
        ENDIF
        IF( cf .EQ. 'E'  .OR. cf .EQ. 'e') THEN
          WRITE(10,*) i1*g%dx, i2*g%dy, 0.0, g%E(i1,i2,ne,1), g%E(i1,i2,ne,2), 0.0
        ELSE
          WRITE(10,*) i1*g%dx, i2*g%dy, 0.0, g%H(i1,i2,ne,1), g%H(i1,i2,ne,2), 0.0
        ENDIF  
      ENDIF             
    ENDDO
  ENDDO

  WRITE(10,*) '$ END'
  CLOSE(10)

END SUBROUTINE ausgabe_feldschnitt

SUBROUTINE grid_plot(n, fd)

  INTEGER,        INTENT(IN)   :: n    ! Zellenanzahl pro Kante
  INTEGER,        INTENT(IN)   :: fd   
  INTEGER                      :: i
  
  WRITE(fd,*)  '$ DATA = CURVE2D'

  DO i = -n/2, n/2, 1
    WRITE(fd,*) '% lc = 0 linetype = 1'
    WRITE(fd,*)    i, -n/2 
    WRITE(fd,*)    i, n/2
    WRITE(fd,*)
  ENDDO

  DO i = -n/2, n/2, 1
    WRITE(fd,*) '% lc = 0 linetype = 1'
    WRITE(fd,*)    -n/2, i 
    WRITE(fd,*)     n/2, i
    WRITE(fd,*)
  ENDDO

END SUBROUTINE grid_plot 

SUBROUTINE density_plot(g, t, ne, co, cf, n, component)
    
  TYPE(gitter),   INTENT(IN)   :: g
  INTEGER,        INTENT(IN)   :: t    ! Zeitschritt
  INTEGER,        INTENT(IN)   :: ne   ! Gitterebene
  INTEGER,        INTENT(IN)   :: n    ! Zellenanzahl pro Kante
  CHARACTER,      INTENT(IN)   :: co   ! Orientierung der Schnittflaeche
  CHARACTER,      INTENT(IN)   :: cf   ! E- oder H-Feld
  INTEGER,        INTENT(IN)   :: component ! (1:3) for x, y or z
  INTEGER                      :: i1, i2
  CHARACTER, DIMENSION(1:5)    :: x
  CHARACTER                    :: komp

  CALL toch(x, t)
  IF (component .EQ. 1) THEN
    komp = 'x'
  ELSEIF (component .EQ. 2) THEN
    komp = 'y'
  ELSE
    komp = 'z'
  ENDIF
  OPEN(10,FILE='density_'//g%id//'_'//co//'_'//cf//komp//'_'//x(1)//x(2)//x(3)//x(4)//x(5)//'.plt',ACTION='WRITE')

  WRITE(10,*) '$ DATA=CONTOUR'
  WRITE(10,*) '% toplabel   = " "'
  WRITE(10,*) '% subtitle   = "Zeitschritte: ', t, '"'
  WRITE(10,*) '% interp     = 2'
  WRITE(10,*) '% contfill   = on'
  WRITE(10,*) '% meshplot   = on'
  WRITE(10,*) '% nsteps     = 50'
  WRITE(10,*) '% nx   =',  n + 1
  WRITE(10,*) '% ny   =',  n + 1
  WRITE(10,*) '% xmin = ', g%nxl, 'xmax =', g%nxgh
  WRITE(10,*) '% ymin = ', g%nxl, 'ymax =', g%nxgh

  DO i2=g%nxl, g%nxgh, 1
    DO i1=g%nxl, g%nxgh, 1     
      IF( co .EQ. 'x' ) THEN
        IF(i1 .EQ. g%nxl .AND.  i2 .EQ. g%nxl) THEN
          WRITE(10,*) '% xlabel="y"'
          WRITE(10,*) '% ylabel="z"'
        ENDIF
        IF( cf .EQ. 'E'  .OR. cf .EQ. 'e') THEN
          WRITE(10,*) g%E(ne, i1, i2, component)
        ELSE
          WRITE(10,*) g%H(ne, i1, i2, component)
        ENDIF
      ENDIF
      IF( co .EQ. 'y' ) THEN
        IF(i1 .EQ. g%nxl .AND.  i2 .EQ. g%nxl) THEN
          WRITE(10,*) '% xlabel="x"'
          WRITE(10,*) '% ylabel="z"'
        ENDIF
        IF( cf .EQ. 'E'  .OR. cf .EQ. 'e') THEN
          WRITE(10,*) g%E(i1, ne, i2, component)
        ELSE
          WRITE(10,*) g%H(i1, ne, i2, component)          
        ENDIF
      ENDIF
      IF( co .EQ. 'z' ) THEN
        IF(i1 .EQ. g%nxl .AND.  i2 .EQ. g%nxl) THEN
          WRITE(10,*) '% xlabel="x"'
          WRITE(10,*) '% ylabel="y"'
        ENDIF
        IF( cf .EQ. 'E'  .OR. cf .EQ. 'e') THEN
          WRITE(10,*) g%E(i1, i2, ne, component)
        ELSE
          WRITE(10,*) g%H(i1, i2, ne, component)
        ENDIF  
      ENDIF     
    ENDDO
  ENDDO

  !CALL grid_plot(n, 10)
  WRITE(10,*)
  WRITE(10,*) '$ END'
  CLOSE(10)

END SUBROUTINE density_plot

SUBROUTINE toch(chaine, i)

  CHARACTER, DIMENSION(1:5), INTENT(INOUT)   :: chaine
  INTEGER, INTENT(IN)                        :: i
  INTEGER                                    :: ii, m
  CHARACTER, DIMENSION(0:9)                  :: cListe

  cListe=(/ '0','1','2','3','4','5','6','7','8','9' /)

  m=i
  DO ii=5, 1, -1
    chaine(ii)= cListe(MOD(m,10))
    m=m/10
  ENDDO

END SUBROUTINE toch

END MODULE plot