analytic.f90

Sfdtd Simple finite-difference time-domain

! analytic.f90
!
! Analytische Dipol Feldberechnung
!
!    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: 25-10-2007 04:47:54 PM CEST

MODULE analytic

USE fdtd_gitter

IMPLICIT NONE

CONTAINS

! return 1 if dipol is in this direction else 0
INTEGER FUNCTION direction(d, komponent)

  TYPE(dipol), INTENT(IN)                               :: d
  INTEGER, INTENT(IN)                                   :: komponent

  IF (d%E(komponent) .NE. 0.0d0) THEN
    direction = 1
  ELSE
    direction = 0
  ENDIF

END FUNCTION direction


! distanz dipol <-> messzelle, |r|
DOUBLE PRECISION FUNCTION r(g, d, messzelle) 
 
  TYPE(gitter), INTENT(IN)                              :: g
  TYPE(dipol), INTENT(IN)                               :: d
  INTEGER, DIMENSION(1:3), INTENT(IN)                   :: messzelle

  r = SQRT(( g%dx * (messzelle(1) - d%px - 0.5d0 * direction(d, 1))) ** 2.0 &
           +(g%dy * (messzelle(2) - d%py - 0.5d0 * direction(d, 2))) ** 2.0 &
           +(g%dz * (messzelle(3) - d%pz - 0.5d0 * direction(d, 3))) ** 2.0)

END FUNCTION r

! components of vector r
DOUBLE PRECISION FUNCTION rk(g, d, messzelle, k) 
 
  TYPE(gitter), INTENT(IN)                              :: g
  TYPE(dipol), INTENT(IN)                               :: d
  INTEGER, DIMENSION(1:3), INTENT(IN)                   :: messzelle
  INTEGER, INTENT(IN)                                   :: k

  IF (k .EQ. 1) THEN
    rk = g%dx * (messzelle(1) - d%px - 0.5d0 * direction(d, 1))
  ELSEIF (k .EQ. 2) THEN 
    rk = g%dy * (messzelle(2) - d%py - 0.5d0 * direction(d, 2))
  ELSE
    rk = g%dz * (messzelle(3) - d%pz - 0.5d0 * direction(d, 3))
  ENDIF 

END FUNCTION rk

SUBROUTINE stimulus(derivs, zelle, g, d, dipol_type, zt_E, zt_H)

  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3)   :: derivs
  INTEGER, INTENT(IN), DIMENSION(1:3)                   :: zelle
  TYPE(gitter), INTENT(IN)                              :: g
  TYPE(dipol), INTENT(IN)                               :: d
  INTEGER, INTENT(IN)                                   :: dipol_type
  DOUBLE PRECISION, INTENT(IN)                          :: zt_E, zt_H

  SELECT CASE (dipol_type)
    CASE(1)
      CALL cosinusoidal(derivs, zelle, zt_E, zt_H, d, g) 
    CASE(2)
      CALL sinusoidal(derivs, zelle, zt_E, zt_H, d, g)
    CASE(3)
      CALL mix(derivs, zelle, zt_E, zt_H, d, g)
    CASE DEFAULT
      WRITE(*,*) "Falscher 'dipol_type' in stimulus()"
  END SELECT

END SUBROUTINE stimulus

DOUBLE PRECISION FUNCTION cos_deriv_0 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  cos_deriv_0 = cos(d%omega * (zeit - (distanz/C)) + d%phi)
  
END FUNCTION cos_deriv_0

DOUBLE PRECISION FUNCTION cos_deriv_1 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  cos_deriv_1 =  - d%omega *  sin(d%omega * (zeit - (distanz/C)) + d%phi)
  
END FUNCTION cos_deriv_1

DOUBLE PRECISION FUNCTION cos_deriv_2 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  cos_deriv_2 = - (d%omega ** 2.0) * cos(d%omega * (zeit - (distanz/C)) + d%phi)
  
END FUNCTION cos_deriv_2

SUBROUTINE cosinusoidal(derivs, zelle, zt_E, zt_H, d, g)

  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3)   :: derivs
  INTEGER, INTENT(IN), DIMENSION(1:3)                   :: zelle
  DOUBLE PRECISION, INTENT(IN)                          :: zt_E, zt_H
  TYPE(dipol), INTENT(IN)                               :: d
  TYPE(gitter), INTENT(IN)                              :: g 
  DOUBLE PRECISION                                      :: ampl, fac

  IF(zt_E >= r(g, d, zelle)/C) THEN
    CALL find_ampl(d%E, ampl)
    fac = EPS * 3.0d0 *(g%dx * g%dx) * ampl

    derivs(1,1) = - fac * cos_deriv_0(d, zt_E, r(g, d, zelle))   
    derivs(1,2) = - fac * cos_deriv_1(d, zt_E, r(g, d, zelle))
    derivs(1,3) = - fac * cos_deriv_2(d, zt_E, r(g, d, zelle))

    derivs(2,1) = - fac * cos_deriv_0(d, zt_H, r(g, d, zelle)) 
    derivs(2,2) = - fac * cos_deriv_1(d, zt_H, r(g, d, zelle))
    derivs(2,3) = - fac * cos_deriv_2(d, zt_H, r(g, d, zelle))
  ELSE
    derivs = 0.0d0
  ENDIF

END SUBROUTINE cosinusoidal

DOUBLE PRECISION FUNCTION sin_deriv_0 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  sin_deriv_0 = sin(d%omega * (zeit - (distanz/C)) + d%phi)
  
END FUNCTION sin_deriv_0

DOUBLE PRECISION FUNCTION sin_deriv_1 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  sin_deriv_1 = d%omega * cos(d%omega * (zeit - (distanz/C)) + d%phi) 
  
END FUNCTION sin_deriv_1

DOUBLE PRECISION FUNCTION sin_deriv_2 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  sin_deriv_2 = - d%omega * d%omega * sin(d%omega * (zeit - (distanz/C)) + d%phi) 
  
END FUNCTION sin_deriv_2

SUBROUTINE sinusoidal(derivs, zelle, zt_E, zt_H, d, g)

  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3)   :: derivs
  INTEGER, INTENT(IN), DIMENSION(1:3)                   :: zelle
  DOUBLE PRECISION, INTENT(IN)                          :: zt_E, zt_H
  TYPE(dipol), INTENT(IN)                               :: d
  TYPE(gitter), INTENT(IN)                              :: g
  DOUBLE PRECISION                                      :: ampl, fac

  IF(zt_E >= r(g, d, zelle)/C) THEN
    CALL find_ampl(d%E, ampl)
    fac = EPS * 3.0d0 * (g%dx * g%dx) * ampl 
    derivs(1,1) = - fac * sin_deriv_0(d, zt_E, r(g, d, zelle))
    derivs(1,2) = - fac * sin_deriv_1(d, zt_E, r(g, d, zelle))
    derivs(1,3) = - fac * sin_deriv_2(d, zt_E, r(g, d, zelle))

    derivs(2,1) = - fac * sin_deriv_0(d, zt_H, r(g, d, zelle))   
    derivs(2,2) = - fac * sin_deriv_1(d, zt_H, r(g, d, zelle))
    derivs(2,3) = - fac * sin_deriv_2(d, zt_H, r(g, d, zelle))
  ELSE
    derivs = 0.0d0
  ENDIF

END SUBROUTINE sinusoidal

DOUBLE PRECISION FUNCTION mix_deriv_0 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  mix_deriv_0 =            SIN(d%omega * (zeit - (distanz/C)) + d%phi)       &
                * (1.0d0 - COS(d%omega * (zeit - (distanz/C)) + d%phi))
  
END FUNCTION mix_deriv_0

DOUBLE PRECISION FUNCTION mix_deriv_1 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  mix_deriv_1 = d%omega * (COS(d%omega * (zeit - (distanz/C)) + d%phi)                            &
                         + SIN(d%omega * (zeit - (distanz/C)) + d%phi) ** 2.0                     &
                         - COS(d%omega * (zeit - (distanz/C)) + d%phi) ** 2.0)
  
END FUNCTION mix_deriv_1

DOUBLE PRECISION FUNCTION mix_deriv_2 (d, zeit, distanz)

  TYPE(dipol), INTENT(IN)                               :: d
  DOUBLE PRECISION, INTENT(IN)                          :: zeit
  DOUBLE PRECISION, INTENT(IN)                          :: distanz

  mix_deriv_2 = (d%omega ** 2.0) * SIN(d%omega * (zeit - (distanz/C)) + d%phi)          &
                        * (4.0d0 * COS(d%omega * (zeit - (distanz/C)) + d%phi) - 1.0d0)
  
END FUNCTION mix_deriv_2

SUBROUTINE mix(derivs, zelle, zt_E, zt_H, d, g)

  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:2,1:3)   :: derivs
  INTEGER, INTENT(IN), DIMENSION(1:3)                   :: zelle
  DOUBLE PRECISION, INTENT(IN)                          :: zt_E, zt_H
  TYPE(dipol), INTENT(IN)                               :: d   
  TYPE(gitter), INTENT(IN)                              :: g  
  DOUBLE PRECISION                                      :: ampl, fac

  IF(zt_E >= r(g, d, zelle)/C) THEN
    CALL find_ampl(d%E, ampl)
    fac = ampl * EPS * 3.0d0 *(g%dx * g%dx)
    derivs(1,1) = - fac * mix_deriv_0(d, zt_E, r(g, d, zelle))
    derivs(1,2) = - fac * mix_deriv_1(d, zt_E, r(g, d, zelle))
    derivs(1,3) = - fac * mix_deriv_2(d, zt_E, r(g, d, zelle))

    derivs(2,1) = - fac * mix_deriv_0(d, zt_H, r(g, d, zelle)) 
    derivs(2,2) = - fac * mix_deriv_1(d, zt_H, r(g, d, zelle))
    derivs(2,3) = - fac * mix_deriv_2(d, zt_H, r(g, d, zelle))
  ELSE
    derivs = 0.0d0
  ENDIF

END SUBROUTINE mix

SUBROUTINE feldberechnung(f, s, zelle, zt_E, d, g)
  
  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(1:6)    :: f               ! Feld array
  DOUBLE PRECISION, INTENT(IN), DIMENSION(1:2,1:3)   :: s               ! Stimulus
  INTEGER, INTENT(IN), DIMENSION(1:3)                :: zelle           ! zelle
  DOUBLE PRECISION, INTENT(IN)                       :: zt_E            ! Zeit_E
  TYPE(dipol), INTENT(IN)                            :: d               ! Dipolstruktur
  TYPE(gitter), INTENT(IN)                           :: g
  INTEGER                                            :: pos1, pos2, pos3, pos, i

  pos1 = 1; pos2 = 2; pos3 = 3

  IF (zt_E >= r(g, d, zelle)/C) THEN

    pos = direction(d, 1) +  2 * direction(d, 2) + 3 * direction(d, 3)

    IF (pos .EQ. pos1) THEN
      CALL permut(pos1, 1)
      CALL permut(pos2, 1)
      CALL permut(pos3, 1)
    ELSEIF(pos .EQ. pos2) THEN
      CALL permut(pos1, -1)
      CALL permut(pos2, -1)
      CALL permut(pos3, -1)
    ENDIF

    ! Ex, Ey, Ez
    f(pos1) = g%dx * rk(g, d, zelle, 1) * rk(g, d, zelle, 3)                 &
                   / (4.0d0 * PI * EPS * r(g, d, zelle) ** 5.0)              &
                  * ( 3.0d0 * (s(1,1) + (r(g, d, zelle)/C) * s(1,2))         &
                  + (r(g, d, zelle)/C)**2.0 * s(1,3) )

    f(pos2) = g%dx * rk(g, d, zelle, 2) * rk(g, d, zelle, 3)                 &
                    / (4.0d0 * PI * EPS * r(g, d, zelle) ** 5.0)             &
                  * ( 3.0d0 * (s(1,1) + (r(g, d, zelle)/C) * s(1,2))         &
                  + (r(g, d, zelle)/C)**2.0 * s(1,3) )   
 
    f(pos3) = g%dx / (4.0d0 * PI * EPS * r(g, d, zelle) ** 5)                &
                * ( ( 2.0d0 * rk(g, d, zelle, 3) ** 2.0                      &
                    - (rk(g, d, zelle, 1)**2.0 + rk(g, d, zelle, 2) ** 2.0)  &
                    ) * (s(1,1) + (r(g, d, zelle)/C) * s(1,2))               &
                    - (rk(g, d, zelle, 1)**2.0 + rk(g, d, zelle, 2) ** 2.0)  &
                      * (r(g, d, zelle)/C) ** 2.0 * s(1,3) )                                                              
    ! Hx, Hy, Hz
    f(pos1+3) = - g%dx * rk(g, d, zelle, 2)                                  &
                / (4.0d0 * PI * r(g, d, zelle) ** 3.0)                       &
                * (s(2,2) + (r(g, d, zelle)/C) * s(2,3))

    f(pos2+3) = g%dx * rk(g, d, zelle, 1)                                    &
                / (4.0d0 * PI * r(g, d, zelle) ** 3.0)                       &
                * (s(2,2) + (r(g, d, zelle)/C) * s(2,3))

    f(pos3+3) = 0.0d0
  ELSE
    f = 0.0d0
  ENDIF

END SUBROUTINE feldberechnung

SUBROUTINE permut(i, sinn)

  INTEGER, INTENT(INOUT)     :: i
  INTEGER, INTENT(IN)        :: sinn
  INTEGER, DIMENSION(1:3)    :: triplet
  
  triplet = (/1, 2, 3/); triplet = CSHIFT(triplet, SIGN(1, sinn))
  i = triplet(i)

END SUBROUTINE permut

SUBROUTINE find_ampl(array, ampl)

  DOUBLE PRECISION, INTENT(IN), DIMENSION(1:3)     :: array 
  DOUBLE PRECISION, INTENT(INOUT)                  :: ampl
  INTEGER                                          :: i

  DO i = 1, 3, 1
    IF (array(i) .NE. 0.0d0) THEN
      ampl  = array(i)
    ENDIF
  ENDDO

END SUBROUTINE find_ampl

END MODULE analytic