maxwell_function.F90

Go to the documentation of this file.

!! Copyright (C) 2019 R. Jestaedt, H. Appel, F. Bonafe, M. Oliveira, N. Tancogne-Dejean
!!
!! 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 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module maxwell_function_oct_m
  use iso_c_binding
  use debug_oct_m
  use fft_oct_m
  use global_oct_m
  use io_oct_m
  use loct_math_oct_m
  use math_oct_m
  use mpi_oct_m
  use parser_oct_m
  use profiling_oct_m
  use splines_oct_m
  use string_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use namespace_oct_m
 
  implicit none
 
  public ::                      &
    mxf_t,                       &
    mxf_init,                    &
    mxf_init_const_wave,         &
    mxf_init_const_phase,        &
    mxf_init_gaussian_wave,      &
    mxf_init_cosinoidal_wave,    &
    mxf_init_fromexpr,           &
    mxf,                         &
    mxf_read,                    &
    mxf_is_empty
 
  integer, public, parameter ::     &
    MXF_EMPTY            =  10001,  &
    mxf_const_wave       =  10002,  &
    mxf_const_phase      =  10004,  &
    mxf_gaussian_wave    =  10005,  &
    mxf_cosinoidal_wave  =  10006,  &
    mxf_logistic_wave    =  10007,  &
    mxf_trapezoidal_wave =  10008,  &
    mxf_from_file        =  10009,  &
    mxf_numerical        =  10010,  &
    mxf_from_expr        =  10011,  &
    mxf_fourier_series   =  10012,  &
    mxf_zero_fourier     =  10013
 
  type mxf_t
    integer :: mode              = mxf_empty
    real(real64)   :: k_vector(3) = m_zero
    real(real64)   :: r0(3)       = m_zero  
    real(real64)   :: width             = m_zero  
    real(real64)   :: a0                = m_zero
    real(real64)   :: dx                = m_zero 
    real(real64)   :: init_x            = m_zero
    real(real64)   :: final_x           = m_zero
    real(real64)   :: growth            = m_zero
    integer :: niter             = 0
    integer :: nfreqs            = 0
 
    type(spline_t)         :: amplitude
    character(len=200)     :: expression
    type(fft_t) :: fft_handler
  end type mxf_t
 
  interface mxf
    module procedure mxf_eval
  end interface mxf
 
contains
 
  !------------------------------------------------------------
  subroutine mxf_read(f, namespace, function_name, ierr)
    type(mxf_t),       intent(inout) :: f
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: function_name
    integer,           intent(out)   :: ierr
 
    type(block_t) :: blk
    integer :: nrows, ncols, i, function_type, idim
    character(len=100) :: row_name, function_expression
    real(real64) :: a0, r0(3), growth, width, k_vector(3)
 
    push_sub(mxf_read)
 
    !%Variable MaxwellFunctions
    !%Type block
    !%Section Maxwell
    !%Description
    !% This block specifies the shape of a "spatial-dependent function", such as the
    !% envelope needed when using the <tt>MaxwellFunctions</tt> block. Each line in the block
    !% specifies one function. The first element of each line will be a string
    !% that defines the name of the function. The second element specifies which type
    !% of function we are using; in the following we provide an example for each of the
    !% possible types:
    !%
    !%Option mxf_const_wave 10002
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_const_wave | kx | ky | kz | x0 | y0 | z0
    !% <br>%</tt>
    !%
    !% The function is constant plane wave <math> f(x,y,z) = a0 * \cos( kx*(x-x0) + ky*(y-y0) + kz*(z-z0) ) </math>
    !%
    !%Option mxf_const_phase 10004
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_const_phase | kx | ky | kz | x0 | y0 | z0
    !% <br>%</tt>
    !%
    !% The function is a constant phase of <math> f(x,y,z) = a0 * (kx * x0 + ky * y0 + kz * z0) </math>
    !%
    !%Option mxf_gaussian_wave 10005
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_gaussian_wave | kx | ky | kz | x0 | y0 | z0 | width
    !% <br>%</tt>
    !%
    !% The function is a Gaussian, <math> f(x,y,z) = a0 * \exp( -( kx*(x-x0) + ky*(y-y0) + kz*(z-z0) )^2 / (2 width^2) ) </math>
    !%
    !%Option mxf_cosinoidal_wave 10006
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_cosinoidal_wave | kx | ky | kz | x0 | y0 | z0 | width
    !% <br>%</tt>
    !%
    !% <math> f(x,y,z) =  \cos( \frac{\pi}{2} \frac{kx*(x-x0)+ky*(y-y0)+kz*(z-z0)-2 width}{width} + \pi )  </math>
    !%
    !% If <math> | kx*x + ky*y + kz*z - x0 | > \xi\_0 </math>, then <math> f(x,y,z) = 0 </math>.
    !%
    !%Option mxf_logistic_wave 10007
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_logistic_wave | kx | ky | kz | x0 | y0 | z0 | growth | width
    !% <br>%</tt>
    !%
    !% The function is a logistic function, <math> f(x,y,z) = a0 * 1/(1+\exp(growth*(kx*(x-x0)+ky*(y-y0)+kz*(kz*(z-z0))+width/2))) * 1/(1+\exp(-growth*(kx*(x-x0)+ky*(y-y0)+kz*(kz*(z-z0))-width/2)))  </math>
    !%
    !%Option mxf_trapezoidal_wave 10008
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_trapezoidal_wave | kx | ky | kz | x0 | y0 | z0 | growth | width
    !% <br>%</tt>
    !%
    !% The function is a logistic function,
    !%   <br> <math> f(x,y,z)      = a0 * ( ( 1-growth*(k*(r-r0)-width/2)*\Theta(k*(r-r0)-width/2))*\Theta(-(k*(r-r0)+width/2+1/growth)) </math>
    !%   <br> <math> \qquad \qquad \qquad + (-1+growth*(k*(r-r0)+width/2)*\Theta(k*(r-r0)+width/2))*\Theta(-(k*(r-r0)-width/2+1/growth)) ) </math>
    !%
    !%Option mxf_from_expr 10011
    !%
    !% <tt>%MaxwellFunctions
    !% <br>&nbsp;&nbsp; "function-name" | mxf_from_expr | kx | ky | kz | "expression"
    !% <br>%</tt>
    !%
    !% The temporal shape of the field is given as an expression (e.g., <tt>cos(2.0*x-3*y+4*z)</tt>. The
    !% letter <i>x</i>, <i>y</i>, <i>z</i> means spatial coordinates, obviously.
    !% The expression is used to construct the function <i>f</i>
    !% that defines the field.
    !%End
    ierr = -3
    if (parse_block(namespace, 'MaxwellFunctions', blk) /= 0) then
      ierr = -1
      pop_sub(mxf_read)
      return
    end if
 
    nrows = parse_block_n(blk)
    row_loop: do i = 1, nrows
      call parse_block_string(blk, i-1, 0, row_name)
      if (trim(row_name) == trim(function_name)) then
 
        ncols = parse_block_cols(blk, i-1)
        call parse_block_integer(blk, i-1, 1, function_type)
 
        a0 = m_one
        r0 = m_zero
        width = m_zero
        k_vector = m_zero
        select case (function_type)
        case (mxf_const_wave)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          do idim = 1, 3
            call parse_block_float(blk, i-1, 4+idim, r0(idim), units_inp%length)
          end do
          call mxf_init_const_wave(f, a0, k_vector, r0)
        case (mxf_const_phase)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          do idim = 1, 3
            call parse_block_float(blk, i-1, 4+idim, r0(idim), units_inp%length)
          end do
          call mxf_init_const_phase(f, a0, k_vector, r0)
        case (mxf_gaussian_wave)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          do idim = 1, 3
            call parse_block_float(blk, i-1, 4+idim, r0(idim), units_inp%length)
          end do
          call parse_block_float(blk, i-1, 8, width, units_inp%length)
          call mxf_init_gaussian_wave(f, a0, k_vector, r0, width)
        case (mxf_cosinoidal_wave)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          do idim = 1, 3
            call parse_block_float(blk, i-1, 4+idim, r0(idim), units_inp%length)
          end do
          call parse_block_float(blk, i-1, 8, width, units_inp%length)
          call mxf_init_cosinoidal_wave(f, a0, k_vector, r0, width)
        case (mxf_logistic_wave)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          do idim = 1, 3
            call parse_block_float(blk, i-1, 4+idim, r0(idim), units_inp%length)
          end do
          call parse_block_float(blk, i-1, 8, growth, units_inp%length)
          call parse_block_float(blk, i-1, 9, width, units_inp%length)
          call mxf_init_logistic_wave(f, a0, k_vector, r0, growth, width)
        case (mxf_trapezoidal_wave)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          do idim = 1, 3
            call parse_block_float(blk, i-1, 4+idim, r0(idim), units_inp%length)
          end do
          call parse_block_float(blk, i-1, 8, growth, units_inp%length)
          call parse_block_float(blk, i-1, 9, width, units_inp%length)
          call mxf_init_trapezoidal_wave(f, a0, k_vector, r0, growth, width)
        case (mxf_from_expr)
          do idim = 1, 3
            call parse_block_float(blk, i-1, 1+idim, k_vector(idim), unit_one/units_inp%length)
          end do
          call parse_block_string(blk, i-1, 5, function_expression)
          call conv_to_c_string(function_expression)
          call mxf_init_fromexpr(f, k_vector, trim(function_expression))
        case default
          ierr = -2
          call parse_block_end(blk)
          pop_sub(mxf_read)
          return
        end select
 
        ierr = 0
        exit row_loop
      end if
    end do row_loop
 
    call parse_block_end(blk)
 
    pop_sub(mxf_read)
  end subroutine mxf_read
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init(f)
    type(mxf_t), intent(inout) :: f
 
    push_sub(mxf_init)
 
    f%mode = mxf_empty
    f%niter = 0
    f%dx = m_zero
 
    pop_sub(mxf_init)
  end subroutine mxf_init
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  logical function mxf_is_empty(f)
    type(mxf_t), intent(in) :: f
 
    push_sub(mxf_is_empty)
    mxf_is_empty = (f%mode == mxf_empty)
 
    pop_sub(mxf_is_empty)
  end function mxf_is_empty
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_const_wave(f, a0, k_vector, r0)
    type(mxf_t), intent(inout) :: f
    real(real64),       intent(in)    :: a0, k_vector(3), r0(3)
 
    push_sub(mxf_init_const_wave)
 
    f%mode = mxf_const_wave
    f%a0 = a0
    f%k_vector = k_vector
    f%r0 = r0
 
    pop_sub(mxf_init_const_wave)
  end subroutine mxf_init_const_wave
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_const_phase(f, a0, k_vector, r0)
    type(mxf_t), intent(inout) :: f
    real(real64),       intent(in)    :: a0, k_vector(3), r0(3)
 
    push_sub(mxf_init_const_phase)
 
    f%mode = mxf_const_phase
    f%a0 = a0
    f%k_vector = k_vector
    f%r0 = r0
 
    pop_sub(mxf_init_const_phase)
  end subroutine mxf_init_const_phase
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_gaussian_wave(f, a0, k_vector, r0, width)
    type(mxf_t), intent(inout) :: f
    real(real64),       intent(in)    :: a0, k_vector(3), r0(3), width
 
    push_sub(mxf_init_gaussian_wave)
 
    f%mode = mxf_gaussian_wave
    f%a0 = a0
    f%k_vector = k_vector
    f%r0 = r0
    f%width = width
 
    pop_sub(mxf_init_gaussian_wave)
  end subroutine mxf_init_gaussian_wave
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_cosinoidal_wave(f, a0, k_vector, r0, width)
    type(mxf_t), intent(inout) :: f
    real(real64),       intent(in)    :: a0, k_vector(3), r0(3), width
 
    push_sub(mxf_init_cosinoidal_wave)
 
    f%mode = mxf_cosinoidal_wave
    f%a0 = a0
    f%k_vector = k_vector
    f%r0 = r0
    f%width = width
 
    pop_sub(mxf_init_cosinoidal_wave)
  end subroutine mxf_init_cosinoidal_wave
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_logistic_wave(f, a0, k_vector, r0, growth, width)
    type(mxf_t), intent(inout) :: f
    real(real64),       intent(in)    :: a0, k_vector(3), r0(3), growth, width
 
    push_sub(mxf_init_logistic_wave)
 
    f%mode = mxf_logistic_wave
    f%a0 = a0
    f%k_vector = k_vector
    f%r0 = r0
    f%growth = growth
    f%width = width
 
    pop_sub(mxf_init_logistic_wave)
  end subroutine
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_trapezoidal_wave(f, a0, k_vector, r0, growth, width)
    type(mxf_t), intent(inout) :: f
    real(real64),       intent(in)    :: a0, k_vector(3), r0(3), growth, width
 
    push_sub(mxf_init_trapezoidal_wave)
 
    f%mode = mxf_trapezoidal_wave
    f%a0 = a0
    f%k_vector = k_vector
    f%r0 = r0
    f%growth = growth
    f%width = width
 
    pop_sub(mxf_init_trapezoidal_wave)
  end subroutine
  !------------------------------------------------------------
 
 
  !------------------------------------------------------------
  subroutine mxf_init_fromexpr(f, k_vector, expression)
    type(mxf_t),      intent(inout) :: f
    real(real64),      intent(in)   :: k_vector(3)
    character(len=*), intent(in)    :: expression
 
    push_sub(mxf_init_fromexpr)
 
    f%mode       = mxf_from_expr
    f%k_vector   = k_vector
    f%expression = trim(expression)
 
    pop_sub(mxf_init_fromexpr)
  end subroutine mxf_init_fromexpr
  !------------------------------------------------------------
 
  !------------------------------------------------------------
  complex(real64) function mxf_envelope_eval(f, x) result(env)
    type(mxf_t), intent(in) :: f
    real(real64),       intent(in) :: x(:)
    real(real64) :: xx, limit_1, limit_2, limit_3, limit_4, r_re, r_im, rr
    complex(real64) :: r
    integer :: xdim
 
    xdim = size(x)
 
    select case (f%mode)
    case (mxf_const_wave)
 
      env = f%a0
 
    case (mxf_const_phase)
 
      env = f%a0 * sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim)))
 
    case (mxf_gaussian_wave)
 
      r = exp(-((sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))) &
        / norm2(f%k_vector(1:xdim)) )**2 / (m_two*f%width**2)))
      env = f%a0 * r
 
    case (mxf_cosinoidal_wave)
 
      r = m_zero
      if (abs( sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))/norm2(f%k_vector(1:xdim)))) <= f%width) then
        r = - cos((m_pi/m_two) * ((sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))) &
          / norm2(f%k_vector(1:xdim)) - m_two*f%width) / f%width))
      end if
      env = f%a0 * r
 
    case (mxf_logistic_wave)
 
      r = m_one/(m_one + exp(f%growth*(sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))) &
        / norm2(f%k_vector(1:xdim)) - f%width/m_two))) &
        + m_one/(m_one + exp(-f%growth*(sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))) &
        / norm2(f%k_vector(1:xdim)) + f%width/m_two))) - m_one
      env = f%a0 * r
 
    case (mxf_trapezoidal_wave)
 
      xx = sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))) / norm2(f%k_vector(1:xdim))
      limit_1 = - f%width/m_two - m_one/f%growth
      limit_2 = - f%width/m_two
      limit_3 =   f%width/m_two
      limit_4 =   f%width/m_two + m_one/f%growth
      if ((xx > limit_1) .and. (xx <= limit_2)) then
        r = m_one + f%growth * (sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim)))/norm2(f%k_vector(1:xdim)) + f%width/m_two)
      else if ((xx > limit_2) .and. (xx <= limit_3)) then
        r = m_one
      else if ((xx > limit_3) .and. (xx <= limit_4)) then
        r = m_one - f%growth * (sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim)))/norm2(f%k_vector(1:xdim)) - f%width/m_two)
      else
        r = m_zero
      end if
      env = f%a0 * r
 
    case (mxf_from_expr)
      ! here x(:) is actually x(:) - v(:)*t
      rr = norm2(x)
      call parse_expression(r_re, r_im, 3, x(:), rr, m_zero, trim(f%expression))
      env = r_re + m_zi*r_im
 
    case default
 
      env = m_zero
 
    end select
 
  end function mxf_envelope_eval
  !------------------------------------------------------------
 
  !------------------------------------------------------------
  complex(real64) function mxf_eval(f, x, phi) result(y)
    type(mxf_t), intent(in) :: f
    real(real64),       intent(in) :: x(:)
    real(real64), optional, intent(in) :: phi
 
    real(real64) :: phi_
    integer :: xdim
 
    ! no push_sub because it is called too frequently
 
    phi_ = optional_default(phi, m_zero)
 
    xdim = size(x)
 
    y = mxf_envelope_eval(f, x)
    if (f%mode /= mxf_const_phase) then
      y = y * exp(m_zi * (sum(f%k_vector(1:xdim)*(x(:) - f%r0(1:xdim))) + phi_))
    endif
 
  end function mxf_eval
  !------------------------------------------------------------
 
end module maxwell_function_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: