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. In the following, we will use</br>
    !% <math>
    !% \mathbf{r} = \begin{pmatrix} {\rm x} \\\\ {\rm y} \\\\ {\rm z} \end{pmatrix},
    !% \mathbf{r}_0 = \begin{pmatrix} {\rm x0} \\\\ {\rm y0} \\\\ {\rm z0} \end{pmatrix},
    !% \mathbf{k} = \begin{pmatrix} {\rm kx} \\\\ {\rm ky} \\\\ {\rm kz} \end{pmatrix},
    !% w = {\rm width},
    !% g = {\rm growth}
    !% </math>
    !%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(\mathbf{r}) = \cos( \mathbf{k} (\mathbf{r} - \mathbf{r}_0) ) </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(\mathbf{r}) = (\mathbf{k} \mathbf{r}_0) </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(\mathbf{r}) = \exp( -( \mathbf{k} (\mathbf{r}-\mathbf{r}_0) )^2 / (2 {w}^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(\mathbf{r}) =
    !%   \begin{cases}
    !%   \cos( \frac{\pi}{2} \frac{ \mathbf{k}(\mathbf{r}-\mathbf{r}_0) - 2 w}{w} + \pi )
    !%   & \text{if } | \mathbf{k}(\mathbf{r}-\mathbf{r}_0) | \le w \\\\
    !%   0 & \text{otherwise}
    !%   \end{cases}
    !% </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, <br>
    !% <math>
    !%   f(\mathbf{r}) =
    !%   \left( \frac{1}{1 + e^{g \left( \frac{\mathbf{k} \cdot (\mathbf{r-r_0})}{|\mathbf{k}|} + \frac{w}{2} \right)}} \right)
    !%   \cdot
    !%   \left( \frac{1}{1 + e^{-g \left( \frac{\mathbf{k} \cdot (\mathbf{r-r_0})}{|\mathbf{k}|} - \frac{w}{2} \right)}} \right)
    !%   \cdot
    !%   {\left(1 + e^{\frac{g \cdot w}{2}}\right)}^2
    !% </math></br>
    !%
    !%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(\mathbf{r}) =
    !% \begin{cases}
    !% 0 & \text{if  }
    !% \mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_0) / |\mathbf{k}| \leq -\frac{w}{2} \\\\
    !% 1 & \text{if  }
    !% -\frac{w}{2} + \frac{1}{g} < \mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_0) / |\mathbf{k}| < \frac{w}{2} - \frac{1}{g} \\\\
    !% 1 - g\left[\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_0) / |\mathbf{k}| - \left(\frac{w}{2} - \frac{1}{g}\right)\right] &
    !% \text{if  } \frac{w}{2} - \frac{1}{g} \leq \mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_0) / |\mathbf{k}| < \frac{w}{2} \\\\
    !% 0 &
    !% \text{if  } \mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_0) / |\mathbf{k}| \geq \frac{w}{2}
    !% \end{cases}
    !% </math></br>
    !%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 = cmplx(r_re, r_im, real64)
 
    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: