sll_m_common_array_initializers.F90

Go to the documentation of this file.

module sll_m_common_array_initializers
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_errors.h"
#include "sll_working_precision.h"
 
   use sll_m_constants, only: &
      sll_p_pi
 
   use sll_m_gaussian, only: &
      sll_f_gaussian_deviate
 
   implicit none
 
   public :: &
      sll_f_hmf_initializer_2d, &
      sll_f_beam_initializer_2d, &
      sll_f_bump_on_tail_initializer_2d, &
      sll_f_constant_time_initializer_1d, &
      sll_f_cos_bell0_initializer_2d, &
      sll_f_cos_bell_initializer_2d, &
      sll_f_cos_sin_initializer_2d, &
      sll_f_diocotron_initializer_2d, &
      sll_f_diocotron_initializer_2d2, &
      sll_f_dsg_2d, &
      sll_f_galaxy_1d_xvx_initializer_v1v2x1x2, &
      sll_f_galaxy_2d_initializer_v1v2x1x2, &
      sll_f_gaussian_beam_initializer_4d, &
      sll_f_gaussian_initializer_2d, &
      sll_f_khp1_2d, &
      sll_f_landau_1d_xvx_initializer_v1v2x1x2, &
      sll_f_landau_1d_yvy_initializer_v1v2x1x2, &
      sll_f_landau_2d_initializer_v1v2x1x2, &
      sll_f_landau_initializer_2d, &
      sll_f_landau_initializer_4d, &
      sll_f_landau_mode_initializer_4d, &
      sll_f_landau_mode_initializer_cos_sum_4d, &
      sll_f_one_initializer_2d, &
      sll_f_periodic_gaussian_initializer_4d, &
      sll_f_rotation_a1_exact_charac_2d, &
      sll_f_rotation_a1_initializer_2d, &
      sll_f_rotation_a2_exact_charac_2d, &
      sll_f_rotation_a2_initializer_2d, &
      sll_f_rotation_phi_initializer_2d, &
      sll_i_scalar_initializer_1d, &
      sll_i_scalar_initializer_2d, &
      sll_i_scalar_initializer_4d, &
      sll_f_sdf_a1_exact_charac_2d, &
      sll_f_sdf_a1_initializer_2d, &
      sll_f_sdf_a2_exact_charac_2d, &
      sll_f_sdf_a2_initializer_2d, &
      sll_f_sdf_phi_initializer_2d, &
      sll_f_sdf_time_initializer_1d, &
      sll_f_test_vx_transport_initializer_v1v2x1x2, &
      sll_f_test_vy_transport_initializer_v1v2x1x2, &
      sll_f_test_x_transport_initializer_v1v2x1x2, &
      sll_f_test_xvx_transport_initializer_v1v2x1x2, &
      sll_f_test_y_transport_initializer_v1v2x1x2, &
      sll_f_test_yvy_transport_initializer_v1v2x1x2, &
      sll_f_translation_a1_exact_charac_2d, &
      sll_f_translation_a1_initializer_2d, &
      sll_f_translation_a2_exact_charac_2d, &
      sll_f_translation_a2_initializer_2d, &
      sll_f_translation_phi_initializer_2d, &
      sll_f_two_stream_instability_initializer_2d, &
      sll_f_twostream_1d_xvx_initializer_v1v2x1x2, &
      sll_f_bump_on_tail_initializer_4d, &
      sll_f_langmuir_cos_initializer_2d, &
      sll_f_langmuir_gaussxv_initializer_2d, &
      sll_f_langmuir_parameter_pert_initializer_2d
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   ! The functions specified here are meant to have the specific signature
   ! described in the sll_m_parallel_array_initializer. Else, they could
   ! not be used with this module.
 
   abstract interface
      function sll_i_scalar_initializer_1d(x1, params)
         use sll_m_working_precision
         sll_real64                                     :: sll_i_scalar_initializer_1d
         sll_real64, intent(in)                         :: x1
         sll_real64, dimension(:), intent(in) :: params
      end function sll_i_scalar_initializer_1d
   end interface
 
   abstract interface
      function sll_i_scalar_initializer_2d(x1, x2, params)
         use sll_m_working_precision
         sll_real64               :: sll_i_scalar_initializer_2d
         sll_real64, intent(in)                         :: x1
         sll_real64, intent(in)                         :: x2
         sll_real64, dimension(:), intent(in) :: params
      end function sll_i_scalar_initializer_2d
   end interface
 
   abstract interface
      function sll_i_scalar_initializer_4d(x1, x2, x3, x4, params)
         use sll_m_working_precision
         sll_real64                                  :: sll_i_scalar_initializer_4d
         sll_real64, intent(in)                         :: x1
         sll_real64, intent(in)                         :: x2
         sll_real64, intent(in)                         :: x3
         sll_real64, intent(in)                         :: x4
         sll_real64, dimension(:), intent(in) :: params
      end function sll_i_scalar_initializer_4d
   end interface
 
contains
 
!phi(i1,i2) =0.5_f64*(x1**2 + x2**2 - 0.1_f64*atan(x2/(x1+0.00001_f64))**2 )
 
   function sll_f_gaussian_initializer_2d(x_1, x_2, params)
      sll_real64 :: sll_f_gaussian_initializer_2d
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: xc_1
      sll_real64 :: xc_2
      sll_real64 :: sigma_1
      sll_real64 :: sigma_2
      !sll_real64 :: kx
      !sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_gaussian_initializer_2d, error: ', &
!!$         'the params array must be passed.', &
!!$         ' params(1) = xc_1', &
!!$         ' params(2) = xc_2', &
!!$         ' params(3) = sigma_1', &
!!$         ' params(4) = sigma_2'
!!$    end if
      sll_assert(size(params) >= 4)
      xc_1 = params(1)
      xc_2 = params(2)
      sigma_1 = params(3)
      sigma_2 = params(4)
      sll_f_gaussian_initializer_2d = &
         exp(-0.5_f64*(x_1 - xc_1)**2/sigma_1**2 &
             - 0.5_f64*(x_2 - xc_2)**2/sigma_2**2)
   end function sll_f_gaussian_initializer_2d
 
   function sll_f_cos_bell_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: xc_1
      sll_real64 :: xc_2
      sll_real64 :: r
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_cos_bell_initializer_2d, error: ', &
!!$         'the params array must be passed.', &
!!$         ' params(1) = xc_1', &
!!$         ' params(2) = xc_2'
!!$    end if
      sll_assert(size(params) >= 2)
      xc_1 = params(1)
      xc_2 = params(2)
 
      r = sqrt((x_1 - xc_1)**2 + (x_2 - xc_2)**2)
 
      if (r < 0.5_f64*sll_p_pi) then
         res = cos(r)**6
      else
         res = 0._f64
      end if
 
   end function sll_f_cos_bell_initializer_2d
 
   ! cossin flow
   function sll_f_cos_sin_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: a1
      sll_real64 :: a2
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_cos_sin_initializer_2d, error: ', &
!!$         'the params array must be passed.', &
!!$         ' params(1) = A1', &
!!$         ' params(2) = A2'
!!$    end if
      sll_assert(size(params) >= 2)
      a1 = params(1)
      a2 = params(2)
      res = (sin(x_1) - a1)*(cos(x_2) - a2)
 
   end function sll_f_cos_sin_initializer_2d
 
   function sll_f_cos_bell0_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: xc_1
      sll_real64 :: xc_2
      sll_real64 :: cutx_value
      sll_real64 :: cutf_value
      sll_real64 :: r
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_cos_bell_initializer_2d, error: ', &
!!$         'the params array must be passed.', &
!!$         ' params(1) = xc_1', &
!!$         ' params(2) = xc_2'
!!$    end if
      sll_assert(size(params) >= 2)
      xc_1 = params(1)
      xc_2 = params(2)
      if (size(params) >= 3) then
         cutx_value = params(3)
      else
         cutx_value = 0.25_f64
      end if
      if (size(params) >= 4) then
         cutf_value = params(4)
      else
         cutf_value = 1._f64
      end if
      r = sqrt((x_1 - xc_1)**2 + (x_2 - xc_2)**2)
 
      if (r < cutx_value*sll_p_pi) then
         res = cutf_value
      else
         res = 0._f64
      end if
 
   end function sll_f_cos_bell0_initializer_2d
 
   function sll_f_one_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x_1 + x_2
#endif
 
      if (size(params) >= 1) then
         res = params(1)
      else
         res = 0._f64
      end if
 
   end function sll_f_one_initializer_2d
 
   !rotation flow
   function sll_f_rotation_phi_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
 
      res = -0.5_f64*(x_1**2 + x_2**2)
   end function sll_f_rotation_phi_initializer_2d
 
   function sll_f_rotation_a1_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x_1 + x_2
#endif
 
      res = -x_2
   end function sll_f_rotation_a1_initializer_2d
 
   function sll_f_rotation_a2_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x_1 + x_2
#endif
 
      res = x_1
   end function sll_f_rotation_a2_initializer_2d
 
   function sll_f_rotation_a1_exact_charac_2d(t, t0, x0_1, x0_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, intent(in) :: t0
      sll_real64, intent(in) :: x0_1
      sll_real64, intent(in) :: x0_2
      sll_real64, dimension(:), intent(in) :: params
 
      res = x0_1*cos(t - t0) - x0_2*sin(t - t0)
   end function sll_f_rotation_a1_exact_charac_2d
 
   function sll_f_rotation_a2_exact_charac_2d(t, t0, x0_1, x0_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, intent(in) :: t0
      sll_real64, intent(in) :: x0_1
      sll_real64, intent(in) :: x0_2
      sll_real64, dimension(:), intent(in) :: params
 
      res = x0_1*sin(t - t0) + x0_2*cos(t - t0)
   end function sll_f_rotation_a2_exact_charac_2d
 
   function sll_f_translation_phi_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: a1
      sll_real64 :: a2
 
      if (size(params) < 2) then
         sll_error("sll_f_translation_phi_initializer_2d", "params should be of size >=2")
         print *, 'size(params)', size(params)
      end if
      a1 = params(1)
      a2 = params(2)
 
      res = -a2*x_1 + a1*x_2
   end function sll_f_translation_phi_initializer_2d
 
   function sll_f_translation_a1_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: a1
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x_1 + x_2
#endif
 
      if (size(params) < 2) then
         sll_error("sll_f_translation_a1_initializer_2d", "params should be of size >=2")
         print *, 'size(params)', size(params)
      end if
      a1 = params(1)
 
      res = a1
   end function sll_f_translation_a1_initializer_2d
 
   function sll_f_translation_a2_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: a2
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x_1 + x_2
#endif
 
      if (size(params) < 2) then
         sll_error("sll_f_translation_a2_initializer_2d", "params should be of size >=2")
         print *, 'size(params)', size(params)
      end if
      a2 = params(2)
 
      res = a2
   end function sll_f_translation_a2_initializer_2d
 
   function sll_f_translation_a1_exact_charac_2d(t, t0, x0_1, x0_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, intent(in) :: t0
      sll_real64, intent(in) :: x0_1
      sll_real64, intent(in) :: x0_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: a1
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x0_1 + x0_2
#endif
 
      if (size(params) < 2) then
         sll_error("sll_f_translation_a1_exact_charac_2d", "params should be of size >=2")
         print *, 'size(params)', size(params)
      end if
      a1 = params(1)
 
      res = x0_1 + (t - t0)*a1
   end function sll_f_translation_a1_exact_charac_2d
 
   function sll_f_translation_a2_exact_charac_2d(t, t0, x0_1, x0_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, intent(in) :: t0
      sll_real64, intent(in) :: x0_1
      sll_real64, intent(in) :: x0_2
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: a2
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x0_1 + x0_2
#endif
 
      if (size(params) < 2) then
         sll_error("sll_f_translation_a2_exact_charac_2d", "params should be of size >=2")
         print *, 'size(params)', size(params)
      end if
      a2 = params(2)
 
      res = x0_2 + (t - t0)*a2
   end function sll_f_translation_a2_exact_charac_2d
 
   function sll_f_constant_time_initializer_1d(t, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, dimension(:), intent(in) :: params
 
      res = params(1)
 
   end function sll_f_constant_time_initializer_1d
 
   !swirling deformation flow
   function sll_f_sdf_phi_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
 
      res = 2._f64*cos(0.5_f64*x_1)**2*cos(0.5_f64*x_2)**2
   end function sll_f_sdf_phi_initializer_2d
 
   function sll_f_sdf_a1_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
 
      res = -cos(0.5_f64*x_1)**2*sin(x_2)
   end function sll_f_sdf_a1_initializer_2d
 
   function sll_f_sdf_a2_initializer_2d(x_1, x_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x_1
      sll_real64, intent(in) :: x_2
      sll_real64, dimension(:), intent(in) :: params
 
      res = cos(0.5_f64*x_2)**2*sin(x_1)
   end function sll_f_sdf_a2_initializer_2d
 
   function sll_f_sdf_a1_exact_charac_2d(t, t0, x0_1, x0_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, intent(in) :: t0
      sll_real64, intent(in) :: x0_1
      sll_real64, intent(in) :: x0_2
      sll_real64, dimension(:), intent(in) :: params
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x0_1 + x0_2 + t + t0
#endif
 
      res = x0_1 !for the moment we take the same point
   end function sll_f_sdf_a1_exact_charac_2d
 
   function sll_f_sdf_a2_exact_charac_2d(t, t0, x0_1, x0_2, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, intent(in) :: t0
      sll_real64, intent(in) :: x0_1
      sll_real64, intent(in) :: x0_2
      sll_real64, dimension(:), intent(in) :: params
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x0_1 + x0_2 + t + t0
#endif
 
      res = x0_2 !for the moment we take the same point
   end function sll_f_sdf_a2_exact_charac_2d
 
   function sll_f_sdf_time_initializer_1d(t, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: t
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: time_period
 
!!$    if( .not. present(params) ) then
!!$      print *, 'sll_f_sdf_time_initializer_1d, error: ', &
!!$        'the params array must be passed.', &
!!$        ' params(1) = time_period'
!!$      stop
!!$    end if
      sll_assert(size(params) >= 1)
      time_period = params(1)
      res = sll_p_pi*cos(sll_p_pi*t/time_period)
   end function sll_f_sdf_time_initializer_1d
 
   ! -------------------------------------------------------------------------
   !
   !             Landau damping 2d initialization function
   !
   ! -------------------------------------------------------------------------
   !
   ! The params array is declared optional to conform with the expected
   ! function signature of the initializer subroutines, but in the particular
   ! case of the landau initializer, the params array must be passed.
 
   function sll_f_landau_initializer_2d(x, vx, params)
      sll_real64 :: sll_f_landau_initializer_2d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
      sll_real64 :: v0
      sll_real64 :: sigma
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_landau_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx.'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 2)
      kx = params(1)
      eps = params(2)
      v0 = 0._f64
      if (size(params) >= 3) then
         v0 = params(3)
      end if
      sigma = 1._f64
      if (size(params) >= 4) then
         sigma = params(4)
      end if
 
      !factor1 = 1.0_f64/sqrt(2.0_f64*sll_p_pi)
      factor1 = 1.0_f64/(sigma*sqrt(2.0_f64*sll_p_pi))
      sll_f_landau_initializer_2d = factor1* &
                                    (1.0_f64 + eps*cos(kx*x))*exp(-0.5_f64*(vx - v0)**2/sigma**2)
   end function sll_f_landau_initializer_2d
 
   function sll_f_bump_on_tail_initializer_2d(x, vx, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_bump_on_tail_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx.'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 2)
      kx = params(1)
      eps = params(2)
      factor1 = 1.0_f64/sqrt(2.0_f64*sll_p_pi)
      res = factor1*(1._f64 + eps*cos(kx*x)) &
            *(0.9_f64*exp(-0.5_f64*vx**2) + 0.2_f64*exp(-0.5_f64*(vx - 4.5_f64)**2/0.5**2))
      !nbox = 3
   end function sll_f_bump_on_tail_initializer_2d
 
   function sll_f_two_stream_instability_initializer_2d(x, vx, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
      sll_real64 :: sigma
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_two_stream_instability_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx.'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 2)
      kx = params(1)
      eps = params(2)
      factor1 = 1.0_f64/sqrt(2.0_f64*sll_p_pi)
      sigma = 1._f64
      if (size(params) >= 3) then
         sigma = params(3)
      end if
      if (size(params) >= 4) then
         factor1 = params(4)
      end if
 
      res = factor1*(1._f64 + eps*cos(kx*x))*vx**2*exp(-0.5_f64*vx**2/sigma**2)
   end function sll_f_two_stream_instability_initializer_2d
 
   function sll_f_diocotron_initializer_2d(r, theta, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: r
      sll_real64, intent(in) :: theta
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: k_mode
      sll_real64 :: r_minus
      sll_real64 :: r_plus
 
!!$    if( .not. present(params) ) then
!!$       print *, '#sll_f_diocotron_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = r_minus'
!!$       print *,'#params(2)= r_plus params(3)=epsilon param(4)=k_mode'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 4)
      r_minus = params(1)
      r_plus = params(2)
      eps = params(3)
      k_mode = params(4)
 
      if ((r >= r_minus) .and. (r <= r_plus)) then
         res = (1.0_f64 + eps*cos(k_mode*theta))
      else
         res = 0._f64
      end if
 
   end function sll_f_diocotron_initializer_2d
 
   function sll_f_diocotron_initializer_2d2(x, y, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: k_mode
      sll_real64 :: r_minus
      sll_real64 :: r_plus
      sll_real64 :: r
      sll_real64 :: theta
 
!!$    if( .not. present(params) ) then
!!$       print *, '#sll_f_diocotron_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = r_minus'
!!$       print *,'#params(2)= r_plus params(3)=epsilon param(4)=k_mode'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 4)
      r_minus = params(1)
      r_plus = params(2)
      eps = params(3)
      k_mode = params(4)
 
      r = sqrt(x**2 + y**2)
 
      if (y >= 0) then
         theta = acos(x/r)
      else
         theta = 2._f64*sll_p_pi - acos(x/r)
      end if
      if ((r >= r_minus) .and. (r <= r_plus)) then
         res = (1.0_f64 + eps*cos(k_mode*theta))
      else
         res = 0._f64
      end if
 
   end function sll_f_diocotron_initializer_2d2
 
   function sll_f_beam_initializer_2d(x, vx, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: alpha
      sll_real64 :: x_plus
      sll_real64 :: x_minus
 
!!$    if( .not. present(params) ) then
!!$       print *, '#sll_f_beam_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = alpha'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 1)
      alpha = params(1)
 
      x_plus = (x + 1.2_f64)/0.3_f64
      x_minus = (x - 1.2_f64)/0.3_f64
      res = 4._f64/sqrt(2._f64*sll_p_pi*alpha)
      res = res*(0.5_f64*erf(x_plus) - 0.5_f64*erf(x_minus))
      res = res*exp(-vx**2/(2*alpha))
 
   end function sll_f_beam_initializer_2d
 
   function sll_f_hmf_initializer_2d(x, vx, params) result(res)
 
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: alpha, beta
      sll_real64 :: m, eps
 
!!$    if( .not. present(params) ) then
!!$       print *, '#sll_f_hmf_initializer_2d, error: the params array must ', &
!!$            'be passed. params = [alpha, beta, m, epsilon]'
!!$       stop
!!$    end if
      sll_assert(size(params) == 4)
      alpha = params(1)
      beta = params(2)
      m = params(3)
      eps = params(4)
 
      res = alpha
      res = res*exp(-beta*((vx*vx)*0.5_f64 - m*cos(x)))
      res = res*(1.0_f64 + eps*cos(x))
 
   end function sll_f_hmf_initializer_2d
 
   function sll_f_khp1_2d(x, y, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: k_mode_x
      sll_real64 :: k_mode_y
 
!!$   if( .not. present(params) ) then
!!$      print *, '#sll_f_khp1_2d, error: the params array must ', &
!!$           'be passed.'
!!$      print *,'#params(1)= eps  param(2)=k_mode'
!!$      stop
!!$   end if
      sll_assert(size(params) >= 3)
      eps = params(1)
      k_mode_x = params(2)
      k_mode_y = params(3)
 
      res = sin(k_mode_y*y) + eps*cos(k_mode_x*x)
 
   end function sll_f_khp1_2d
 
   function sll_f_dsg_2d(eta1, eta2, params) result(res)
      sll_real64  :: res
      sll_real64, intent(in)   :: eta1
      sll_real64, intent(in)   :: eta2
      sll_real64, dimension(:), intent(in) :: params
 
      sll_real64  :: eta1n
      sll_real64  :: eta2n
      sll_real64  :: eta1_min
      sll_real64  :: eta2_min
      sll_real64  :: eta1_max
      sll_real64  :: eta2_max
!!$    if( .not. present(params) ) then
!!$       print *, '#sll_D_sharped_Geo_2d, error: the params array must ', &
!!$            'be passed. params(1) = eta1_min, params(2) = eta2_min', &
!!$            'be passed. params(3) = eta1_max, params(4) = eta2_max'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 4)
      eta1_min = params(1)
      eta2_min = params(2)
      eta1_max = params(3)
      eta2_max = params(4)
      eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
      eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
      res = 4._f64*eta1n*(1._f64 - eta1n)*(1._f64 + 0.1_f64*sin(8.*sll_p_pi*eta2n))
   end function sll_f_dsg_2d
 
   ! This is a simplistic initializer aimed at a 4d cartesian distribution
   ! function, periodic in x and y, and compact-ish in vx and vy.
   !
   ! Basically:
   !
   ! f(x,y,vx,vy) = alpha*exp(-0.5*((x -xc )^2+(y - yc)^2)) +
   !                beta* exp(-0.5*((vx-vxc)^2+(vy-vyc)^2))
   !
   !
   ! It is meant to be used in the intervals:
   ! x:  [ 0,2*pi/kx]
   ! y:  [ 0,2*pi/ky]
   ! vx: [-6,6]
   ! vy: [-6,6]
 
   ! convention for the params array:
   ! params(1) = eta1_min
   ! params(2) = eta1_max
   ! params(3) = eta2_min
   ! params(4) = eta2_max
   ! params(5) = epsilon
 
   function sll_gaussian_initializer_4d(x, y, vx, vy, params)
      sll_real64 :: sll_gaussian_initializer_4d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: xc
      sll_real64 :: yc
      sll_real64 :: vxc
      sll_real64 :: vyc
      sll_real64 :: alpha
      sll_real64 :: beta
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_gaussian_initializer_4d, error: the params array must ', &
!!$            'be passed: ', &
!!$            'params(1) = xc, params(2) = yc, params(3) = vxc, params(4) = vyc',&
!!$            'params(5) = alpha, params(6) = beta'
!!$       stop
!!$    end if
 
      xc = params(1)
      yc = params(2)
      vxc = params(3)
      vyc = params(4)
      alpha = params(5)
      beta = params(6)
 
      sll_gaussian_initializer_4d = alpha*exp(-0.5_f64*((x - xc)**2 + (y - yc)**2)) + &
                                    beta*exp(-0.5_f64*((vx - vxc)**2 + (vy - vyc)**2))
 
   end function sll_gaussian_initializer_4d
 
   ! -------------------------------------------------------------------------
   !
   !             Landau damping 4d initialization function
   !
   ! -------------------------------------------------------------------------
   !
   ! The params array is declared optional to conform with the expected
   ! function signature of the initializer subroutines, but in the particular
   ! case of the landau initializer, the params array must be passed.
 
   function sll_f_landau_initializer_4d(x, y, vx, vy, params)
      sll_real64 :: sll_f_landau_initializer_4d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_landau_initializer_4d, error: the params array must ', &
!!$            'be passed. params(1) = eta1_min, params(2) = eta1_max, ', &
!!$            'params(3) = eta2_min, params(4) = eta2_max, params(5) = epsilon.'
!!$       print *,'does not depend on y',y
!!$       stop
!!$    end if
 
      sll_assert(size(params) >= 5)
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
      eps = params(5)
      kx = 2.*sll_p_pi/(eta1_max - eta1_min)
 
      !Normalization
      !sagemath command
      !sage : var('u v epsilon a b c d x y')
      !sage : f(a,b,c,d,epsilon) =integral(integral(integral(integral((1+epsilon*cos(2*pi/(b-a)*x))*exp(-(u*u+v*v)/2),u,-oo,oo),v,-oo,oo),x,a,b),y,c,d)
 
!!$    factor1 =  1./( (eta2_min - eta2_max) &
!!$               *(((eta1_min - eta1_max)* &
!!$               sin(2*sll_p_pi*eta1_min/(eta1_min - eta1_max)) &
!!$                - (eta1_min - eta1_max)* &
!!$               sin(2*sll_p_pi*eta1_max/(eta1_min - eta1_max)))*eps  &
!!$               + 2*sll_p_pi*eta1_min - 2*sll_p_pi*eta1_max))
      factor1 = 1.0_f64/(2.0_f64*sll_p_pi)
!!$    sll_f_landau_initializer_4d = factor1 * &
!!$         (1.0_f64/((eta2_max-eta2_min)*(eta1_max-eta1_min))+eps*cos(kx*x))*exp(-0.5_f64*(vx**2+vy**2))
      sll_f_landau_initializer_4d = factor1* &
                                    (1.0_f64 + eps*cos(kx*x))*exp(-0.5_f64*(vx**2 + vy**2))
   end function sll_f_landau_initializer_4d
 
   function sll_f_landau_mode_initializer_cos_sum_4d(x, y, vx, vy, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: ky
      sll_real64 :: ellx
      sll_real64 :: elly
      sll_real64 :: eps_ell
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_landau_initializer_4d, error: the params array must ', &
!!$            'be passed. params(1) = kx, params(2) = ky, ', &
!!$            'params(3) = eps.'
!!$       stop
!!$    end if
 
      sll_assert(size(params) >= 3)
 
      kx = params(1)
      ky = params(2)
      eps = params(3)
 
      !factor1 = eps*cos(kx*x)*cos(ky*y)
      factor1 = eps*cos(kx*x + ky*y)
      if (size(params) >= 6) then
         ellx = params(4)
         elly = params(5)
         eps_ell = params(6)
         !factor1 = factor1+eps_ell*cos(ellx*x)*cos(elly*y)
         factor1 = factor1 + eps_ell*cos(ellx*x + elly*y)
      end if
      factor1 = 1._f64 + factor1
      res = (1.0_f64/(2.0*sll_p_pi))*factor1*exp(-0.5_f64*(vx**2 + vy**2))
 
      !print *,'k=',kx,ky,eps
      !print *,'ell=',ellx,elly,eps_ell
      !print *,'size(params)=',size(params)
      !stop
 
   end function sll_f_landau_mode_initializer_cos_sum_4d
 
   function sll_f_landau_mode_initializer_4d(x, y, vx, vy, params) result(res)
      sll_real64 :: res
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: ky
      sll_real64 :: ellx
      sll_real64 :: elly
      sll_real64 :: eps_ell
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_landau_initializer_4d, error: the params array must ', &
!!$            'be passed. params(1) = kx, params(2) = ky, ', &
!!$            'params(3) = eps.'
!!$       stop
!!$    end if
 
      sll_assert(size(params) >= 3)
 
      kx = params(1)
      ky = params(2)
      eps = params(3)
 
      factor1 = eps*cos(kx*x)*cos(ky*y)
      !factor1 = eps*cos(kx*x+ky*y)
      if (size(params) >= 6) then
         ellx = params(4)
         elly = params(5)
         eps_ell = params(6)
         factor1 = factor1 + eps_ell*cos(ellx*x)*cos(elly*y)
         !factor1 = factor1+eps_ell*cos(ellx*x+elly*y)
      end if
      factor1 = 1.0_f64 + factor1
      res = (1.0_f64/(2.0_f64*sll_p_pi))*factor1*exp(-0.5_f64*(vx**2 + vy**2))
 
      !print *,'k=',kx,ky,eps
      !print *,'ell=',ellx,elly,eps_ell
      !print *,'size(params)=',size(params)
      !stop
 
   end function sll_f_landau_mode_initializer_4d
 
   ! -------------------------------------------------------------------------
   !
   !             Bump-on-tail 4d initialization function
   !
   ! -------------------------------------------------------------------------
   !
   !
 
   function sll_f_bump_on_tail_initializer_4d(x, y, vx, vy, params)
      sll_real64 :: sll_f_bump_on_tail_initializer_4d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64, dimension(:), intent(in) :: params
 
      sll_real64 :: kx
      sll_real64 :: ky
      sll_real64 :: eps
      sll_real64 :: sigmax1
      sll_real64 :: sigmax2
      sll_real64 :: sigmay
      sll_real64 :: delta
      sll_real64 :: vx2
      sll_real64 :: factor
 
      sll_assert(size(params) >= 8)
 
      kx = params(1)
      ky = params(2)
      eps = params(3)
      sigmax1 = params(4)
      sigmax2 = params(5)
      sigmay = params(6)
      delta = params(7)
      vx2 = params(8)
 
      factor = 1.0_f64/(2.0_f64*sll_p_pi*sigmay)
 
      sll_f_bump_on_tail_initializer_4d = factor* &
                                          (1.0_f64 + eps*(cos(kx*x) + cos(ky*y)))*exp(-0.5_f64*(vy/sigmay)**2)* &
                                          ((1 - delta)/sigmax1*exp(-0.5_f64*(vx/sigmax1)**2) + &
                                           delta/sigmax2*exp(-0.5_f64*((vx - vx2)/sigmax2)**2))
 
   end function sll_f_bump_on_tail_initializer_4d
 
   function sll_f_test_x_transport_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_test_x_transport_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: kx
 
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_x_transport_initializer, error: the params array', &
!!$            'must be passed. params(1) = epsilon, params(2)=kx, params(3) = ky.'
!!$       print *,'does not depend on y,vy',y,vy
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      t = params(11)
      kx = 2.0_f64*sll_p_pi/(eta1_max - eta1_min)
 
      !sll_f_test_x_transport_initializer_v1v2x1x2 =sin(kx*(x-t))
      sll_f_test_x_transport_initializer_v1v2x1x2 = sin(kx*(x - vx*t))
      !sll_f_test_x_transport_initializer_v1v2x1x2 = exp(-4*x**2)
 
   end function sll_f_test_x_transport_initializer_v1v2x1x2
 
   function sll_f_test_y_transport_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_test_y_transport_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
      sll_real64 :: kx
 
      eta2_min = params(3)
      eta2_max = params(4)
      t = params(11)
      kx = 2.0_f64*sll_p_pi/(eta2_max - eta2_min)
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_y_transport_initializer, error: the params array', &
!!$            'must be passed. params(1) = epsilon, params(2)=kx, params(3) = ky.'
!!$       print *,'does not depend on x vx',x,vx
!!$       stop
!!$    end if
 
      sll_f_test_y_transport_initializer_v1v2x1x2 = sin(kx*(y - vy*t))
      !sll_f_test_y_transport_initializer_v1v2x1x2 =exp(-4*y**2)
      !sll_f_test_x_transport_initializer_v1v2x1x2 = 2_f64
 
   end function sll_f_test_y_transport_initializer_v1v2x1x2
 
   function sll_f_test_vx_transport_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_test_vx_transport_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: kx, v, hh
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_vx_transport_initializer, error: the params array', &
!!$            ' mustbe passed. params(1)=epsilon,params(2)=kx, params(3) = ky.'
!!$       print *,'does not depend on x,y,vy',x,y,vy
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      t = params(11)
      kx = 2.*sll_p_pi/(eta1_max - eta1_min)
 
      !sll_f_test_vx_transport_initializer_v1v2x1x2 =exp(-4*(vx-t)**2)
      v = vx - t
      hh = 1.0_f64
      if ((v .gt. 1) .or. (v .lt. -1)) hh = 0.0_f64
      sll_f_test_vx_transport_initializer_v1v2x1x2 = (1 + (-3 + (3 - v**2)*v**2)*v**2)*hh
      !write(*,*) 'tini=',t
 
!!$    sll_landau_initializer_v1v2x1x2 = 1
!!$    if (x < 0) sll_landau_initializer_v1v2x1x2 = 0
 
   end function sll_f_test_vx_transport_initializer_v1v2x1x2
   function sll_f_test_vy_transport_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_test_vy_transport_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: t
 
      sll_real64, dimension(:), intent(in) :: params
      !sll_real64 :: eta1_min
      !sll_real64 :: eta1_max
      !sll_real64 :: eta2_min
      !sll_real64 :: eta2_max
 
      !sll_real64 :: eps
      !sll_real64 :: kx
      !sll_real64 :: factor1
      sll_real64 :: v, hh
      t = params(11)
 
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_vy_transport_initializer, error: the params array', &
!!$            'must be passed. params(1) = epsilon, params(2)=kx, params(3) = ky.'
!!$       print *, 'does not depend on x,y,vx',x,y,vx
!!$       stop
!!$    end if
      v = vy - t
      hh = 1.0_f64
      if ((v .gt. 1) .or. (v .lt. -1)) hh = 0.0_f64
      sll_f_test_vy_transport_initializer_v1v2x1x2 = (1 + (-3 + (3 - v**2)*v**2)*v**2)*hh
      !sll_f_test_vy_transport_initializer_v1v2x1x2 =exp(-4*vy**2)
      !sll_f_test_vy_transport_initializer_v1v2x1x2 = 2_f64
 
   end function sll_f_test_vy_transport_initializer_v1v2x1x2
 
   function sll_f_test_xvx_transport_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_test_xvx_transport_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: kx, v, hh, xx
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_xvx_transport_initializer, error: the params array', &
!!$            ' mustbe passed. params(1)=epsilon,params(2)=kx, params(3) = ky.'
!!$       print *, 'does not depend on y and vy',y,vy
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      t = params(11)
      kx = 2.*sll_p_pi/(eta1_max - eta1_min)
 
      !sll_f_test_vx_transport_initializer_v1v2x1x2 =exp(-4*(vx-t)**2)
      v = vx - t
      xx = x - t*vx + t*t/2
      hh = 1.0_f64
      if ((v .gt. 1) .or. (v .lt. -1)) hh = 0.0_f64
      sll_f_test_xvx_transport_initializer_v1v2x1x2 = (1 + (-3 + (3 - v**2)*v**2)*v**2)*hh*sin(kx*xx)
 
   end function sll_f_test_xvx_transport_initializer_v1v2x1x2
 
   function sll_f_test_yvy_transport_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_test_yvy_transport_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
      sll_real64 :: ky, v, hh, yy
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_yvy_transport_initializer, error: the params array', &
!!$            ' mustbe passed. params(1)=epsilon,params(2)=kx, params(3) = ky.'
!!$       print *,'does not depend on x and vx',x,vx
!!$       stop
!!$    end if
 
      eta2_min = params(3)
      eta2_max = params(4)
      t = params(11)
      ky = 2.*sll_p_pi/(eta2_max - eta2_min)
 
      !sll_f_test_vx_transport_initializer_v1v2x1x2 =exp(-4*(vx-t)**2)
      v = vy - t
      yy = y - t*vy + t*t/2
      hh = 1.0_f64
      if ((v .gt. 1) .or. (v .lt. -1)) hh = 0.0_f64
      sll_f_test_yvy_transport_initializer_v1v2x1x2 = (1 + (-3 + (3 - v**2)*v**2)*v**2)*hh*sin(ky*yy)
 
   end function sll_f_test_yvy_transport_initializer_v1v2x1x2
 
!!$function sll_test_transport_2d_initializer_v1v2x1x2( vx, vy, x, y, params )
!!$    sll_real64 :: sll_test_transport_2d_initializer_v1v2x1x2
!!$    sll_real64, intent(in) :: x
!!$    sll_real64, intent(in) :: y
!!$    sll_real64, intent(in) :: vx
!!$    sll_real64, intent(in) :: vy
!!$    sll_real64 :: t
!!$
!!$    sll_real64, dimension(:), intent(in) :: params
!!$    sll_real64 :: eta1_min
!!$    sll_real64 :: eta1_max
!!$    sll_real64 :: kx,v,hh,xx
!!$    if( .not. present(params) ) then
!!$       print *, ' sll_test_xvx_transport_initializer, error: the params array', &
!!$            ' mustbe passed. params(1)=epsilon,params(2)=kx, params(3) = ky.'
!!$       stop
!!$    end if
!!$
!!$    eta1_min = params(1)
!!$    eta1_max = params(2)
!!$    eta2_min = params(3)
!!$    eta2_max = params(4)
!!$    t=params(11)
!!$    kx  =  2. * sll_p_pi / (eta1_max - eta1_min)
!!$    ky  =  2. * sll_p_pi / (eta2_max - eta2_min)
!!$
!!$    !sll_f_test_vx_transport_initializer_v1v2x1x2 =exp(-4*(vx-t)**2)
!!$    vvx=vx-t
!!$    xx=x-t*vx+t*t/2
!!$    vvy=vx-t
!!$    yy=x-t*vy+t*t/2
!!$    hh=1
!!$    if ((v.gt.1).or.(v.lt.-1)) hh=0
!!$    sll_test_transport_2d_initializer_v1v2x1x2=(1+(-3+(3-v**2)*v**2)*v**2)*hh*sin(kx*(xx-v*t))
!!$
!!$  end function sll_test_transport_2d_initializer_v1v2x1x2
 
   function sll_f_landau_1d_xvx_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_landau_1d_xvx_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      !sll_real64  :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_landau_1d_initializer_v1v2x1x2, error: the params array', &
!!$            'must be passed params(1)= epsilon, params(2) = kx, params(3) = ky.'
!!$       print *,'does not depend on y and vy',y,vy
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
 
      eps = params(5)
      !kx  =  2. * sll_p_pi / (eta1_max - eta1_min)
      kx = 0.2_f64
      factor1 = 1.0_f64/sqrt((2.0_f64*sll_p_pi))
 
      sll_f_landau_1d_xvx_initializer_v1v2x1x2 = factor1* &
                                                 (1.0_f64 + eps*cos(kx*x))*exp(-0.5_f64*(vx**2))
   end function sll_f_landau_1d_xvx_initializer_v1v2x1x2
 
   function sll_f_twostream_1d_xvx_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_twostream_1d_xvx_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64  :: v0
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
#ifdef DEBUG
      sll_real64 :: dummy
      dummy = x + y + vx + vy
#endif
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_twostream_1d_initializer_v1v2x1x2, error: the params array', &
!!$            'must be passed params(1)= epsilon, params(2) = kx, params(3) = ky.'
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
 
      eps = params(5)
      v0 = 3.0_f64
      !kx  =  2. * sll_p_pi / (eta1_max - eta1_min)
      kx = 0.2_f64
      factor1 = 1.0_f64/sqrt((2.0*sll_p_pi))
 
      sll_f_twostream_1d_xvx_initializer_v1v2x1x2 = factor1/2* &
                                                    (1.0_f64 + eps*cos(kx*x))*(exp(-0.5_f64*((vx - v0)**2)) + &
                                                                               exp(-0.5_f64*((vx + v0)**2)))
   end function sll_f_twostream_1d_xvx_initializer_v1v2x1x2
 
   function sll_f_galaxy_1d_xvx_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_galaxy_1d_xvx_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      !sll_real64  :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_landau_1d_initializer_v1v2x1x2, error: the params array', &
!!$            'must be passed params(1)= epsilon, params(2) = kx, params(3) = ky.'
!!$       print *,'does not depend on y and vy',y,vy
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
 
      eps = params(5)
      !kx  =  2. * sll_p_pi / (eta1_max - eta1_min)
      kx = 0.2_f64
      factor1 = 1.0_f64/sqrt((2.0*sll_p_pi))
      sll_f_galaxy_1d_xvx_initializer_v1v2x1x2 = 0.0_f64
      if ((x .lt. 1.0_f64) .and. (x .gt. -1.0_f64)) then
         sll_f_galaxy_1d_xvx_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*(vx**2))
      end if
   end function sll_f_galaxy_1d_xvx_initializer_v1v2x1x2
 
   function sll_f_galaxy_2d_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_galaxy_2d_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      sll_real64 :: v0x
      sll_real64 :: v0y
      !sll_real64  :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
      sll_real64 :: a1_xmin, a1_xmax, a1_ymin, a1_ymax, x1mil
      sll_real64 :: a2_xmin, a2_xmax, a2_ymin, a2_ymax, x2mil
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_landau_1d_initializer_v1v2x1x2, error: the params array', &
!!$            'must be passed params(1)= epsilon, params(2) = kx, params(3) = ky.'
!!$       print *,'does not depend on y and vy',y,vy
!!$       stop
!!$    end if
      !lmin=0.0_f64
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
      !amas 1
      a1_xmin = -2.0_f64
      a1_xmax = 2.0_f64
      x1mil = 0.0_f64
      a1_ymin = -a1_xmax
      a1_ymax = -a1_xmin
      !amas 2
      a2_xmin = -6.0_f64
      a2_xmax = -2.0_f64
      x2mil = -4.0_f64
      a2_ymin = -a2_xmax
      a2_ymax = -a2_xmin
 
      eps = params(5)
      !kx  =  2. * sll_p_pi / (eta1_max - eta1_min)
      kx = 0.2_f64
      factor1 = 1.0_f64/sqrt((2.0*sll_p_pi))**2
      sll_f_galaxy_2d_initializer_v1v2x1x2 = 0.0_f64
      !un bloc de particles
!!$    if ((x.le.1.0_f64).and.(x.ge.-1.0_f64).and.(y.le.1.0_f64).and.(y.ge.-1.0_f64)) then
!!$       sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*(vx**2+vy**2))
!!$    end if
      !2 blocs de particles
      !amas 1
      if ((x .le. a1_xmax) .and. (x .ge. a1_xmin) .and. (y .le. a1_ymax) .and. (y .ge. a1_ymin)) then
         v0x = 0.0_f64
         v0y = 0.0_f64
 sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))*exp(-0.5_f64*((x-x1mil)**2+(y+x1mil)**2)*16)
      end if
      !amas 2
      if ((x .le. a2_xmax) .and. (x .ge. a2_xmin) .and. (y .le. a2_ymax) .and. (y .ge. a2_ymin)) then
         v0x = 0.0_f64 ! 1.0_f64 !0.3109793990   !0.2433894805
         v0y = -0.0_f64 !.5
  sll_f_galaxy_2d_initializer_v1v2x1x2 =factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))*exp(-0.5_f64*((x-x2mil)**2+(y+x2mil)**2)*16)
      end if
!!$ if ((x.le.lmax).and.(x.ge.lmin).and.(y.le.-lmin).and.(y.ge.-lmax)) then
!!$       v0x=-1.0
!!$       v0y=0
!!$       sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))*exp(-((x-lmil)**2+(y+lmil)**2))
!!$    end if
!!$    if ((y.le.lmax).and.(y.ge.lmin).and.(x.le.-lmin).and.(x.ge.-lmax)) then
!!$       v0x=1.0
!!$       v0y=0
!!$        sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))*exp(-((x+lmil)**2+(y-lmil)**2))
!!$     end if
      !  end if
!!$    if ((x.le.lmax).and.(x.ge.lmin).and.(y.le.-lmin).and.(y.ge.-lmax)) then
!!$       v0x=-1.0
!!$       v0y=0
!!$       sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))*exp(-((x-lmil)**2+(y+lmil)**2))
!!$    end if
!!$ if ((y.le.lmax).and.(y.ge.lmin).and.(x.le.-lmin).and.(x.ge.-lmax)) then
!!$       v0x=1.0
!!$       v0y=0
!!$        sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))*exp((x+lmil)**2+(y-lmil)**2)
!!$    end if
!!$
!!$    if ((x.le.3.0_f64).and.(x.ge.1.5_f64).and.(y.le.-1.5_f64).and.(y.ge.-3.0_f64)) then
!!$       v0x=-1.0_f64
!!$       v0y=0.0
!!$       sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))
!!$    end if
!!$ if ((y.le.3.0_f64).and.(y.ge.1.5_f64).and.(x.le.-1.5_f64).and.(x.ge.-3.0_f64)) then
!!$       v0x=1.0_f64
!!$       v0y=-0.0
!!$        sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*((vx-v0x)**2+(vy-v0y)**2))
!!$    end if
 
!!$    if ((x.le.2.5_f64).and.(x.ge.0.5_f64).and.(y.le.-0.5_f64).and.(y.ge.-2.5_f64)) then
!!$       sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*(vx**2+vy**2))
!!$    end if
!!$ if ((y.le.2.5_f64).and.(y.ge.0.5_f64).and.(x.le.-0.5_f64).and.(x.ge.-2.5_f64)) then
!!$       sll_f_galaxy_2d_initializer_v1v2x1x2 = factor1*exp(-0.5_f64*(vx**2+vy**2))
!!$    end if
   end function sll_f_galaxy_2d_initializer_v1v2x1x2
 
   function sll_f_landau_1d_yvy_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_landau_1d_yvy_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
      !sll_real64  :: t
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      sll_real64 :: eps
      sll_real64 :: ky
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_landau_1d_initializer_v1v2x1x2, error: the params array', &
!!$            'must be passed params(1)= epsilon, params(2) = kx, params(3) = ky.'
!!$       print *,'does not depend on vx and x',vx,x
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
 
      eps = params(5)
      ky = 2.*sll_p_pi/(eta2_max - eta2_min)
      factor1 = 1.0_f64/sqrt((2.0*sll_p_pi))
 
      sll_f_landau_1d_yvy_initializer_v1v2x1x2 = factor1* &
                                                 (1.0_f64 + eps*cos(ky*y))*exp(-0.5_f64*(vy**2))
   end function sll_f_landau_1d_yvy_initializer_v1v2x1x2
 
   function sll_f_landau_2d_initializer_v1v2x1x2(vx, vy, x, y, params)
      sll_real64 :: sll_f_landau_2d_initializer_v1v2x1x2
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: ky
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_f_landau_initializer_4d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx, params(3) = ky.'
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
 
      eps = params(5)
      kx = 2.*sll_p_pi/(eta1_max - eta1_min)
      ky = kx
      factor1 = 1.0_f64/(2.0*sll_p_pi)
      sll_f_landau_2d_initializer_v1v2x1x2 = factor1* &
                                             (1.0_f64 + eps*cos(kx*x)*cos(ky*y))*exp(-0.5_f64*(vx**2 + vy**2))
 
   end function sll_f_landau_2d_initializer_v1v2x1x2
 
   ! this function is a 1D landau initializer used for debugging
   ! 4D drift kinetic simulations in variables x1,x2,x3 ,v1
   ! the function is constant with respect to x2 and x3
 
   function sll_landau_initializer_dk_test_4d(v1, x1, x2, x3, params)
      sll_real64 :: sll_landau_initializer_dk_test_4d
      sll_real64, intent(in) :: x1
      sll_real64, intent(in) :: x2
      sll_real64, intent(in) :: x3
      sll_real64, intent(in) :: v1
      sll_real64, dimension(:), intent(in) :: params
 
      sll_real64 :: epsilon
      sll_real64 :: kx
      sll_real64 :: factor1
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_landau_initializer_dk_test_4d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx'
!!$       print *,'does not depend on x2,x3',x2,x3
!!$       stop
!!$    end if
 
      epsilon = params(1)
      kx = params(2)
      factor1 = 0.5_f64/sll_p_pi
 
      !write(*,*) 'x1 v1',x1,v1
 
      sll_landau_initializer_dk_test_4d = factor1* &
                                          (1.0_f64 + epsilon*cos(kx*x1))*exp(-0.5_f64*(v1**2))
   end function sll_landau_initializer_dk_test_4d
 
   !---------------------------------------------------------------------------
   !
   !                         Periodic Maxwellian
   !
   ! 4D distribution in [0,1]X[0,1][-6,6]X[-6,6]  with the property of being
   ! periodic in the spatial directions (x1,x2) and sll_m_gaussian in velocity space.
   !
   ! f(x,y,vx,vy) = sin(kx*x)*sin(ky*y) +
   !                beta* exp(-0.5*((vx-vxc)^2+(vy-vyc)^2))
   !
   !
   ! It is meant to be used in the intervals:
   ! x:  [ 0,2*pi/kx]
   ! y:  [ 0,2*pi/ky]
   ! vx: [-6,6]
   ! vy: [-6,6]
 
   ! convention for the params array:
   ! params(1) = eta1_min
   ! params(2) = eta1_max
   ! params(3) = eta2_min
   ! params(4) = eta2_max
   ! params(5) = epsilon
   !
   !---------------------------------------------------------------------------
   function sll_f_periodic_gaussian_initializer_4d(x, y, vx, vy, params) &
      result(val)
 
      sll_real64 :: val !sll_gaussian_initializer_4d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: xc
      sll_real64 :: yc
      sll_real64 :: vxc
      sll_real64 :: vyc
      sll_real64 :: alpha
      sll_real64 :: beta
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_gaussian_initializer_4d, error: the params array must ', &
!!$            'be passed: ', &
!!$            'params(1) = xc, params(2) = yc, params(3) = vxc, params(4) = vyc',&
!!$            'params(5) = alpha, params(6) = beta'
!!$       stop
!!$    end if
 
      xc = params(1)
      yc = params(2)
      vxc = params(3)
      vyc = params(4)
      alpha = params(5)
      beta = params(6)
 
      val = alpha*exp(-0.5_f64*((x - xc)**2 + (y - yc)**2)) + &
            beta*exp(-0.5_f64*((vx - vxc)**2 + (vy - vyc)**2))
   end function sll_f_periodic_gaussian_initializer_4d
 
   !---------------------------------------------------------------------------
   !
   !                         Periodic-Periodic initializer
   !
   ! 4D distribution in [0,2*pi/kx]X[0,2*pi/ky][-6,6]X[-6,6]  with the property of being
   ! periodic in the spatial directions (x1,x2) and sll_m_gaussian in velocity space.
   !
   ! f(x,y,vx,vy) = (1 + alpha * cos(kx*x)*cos(ky*y) )
   !                * exp(-0.5*((vx-vxc)^2+(vy-vyc)^2))/ 2*pi
   !
   !  This function is described in the article of Crouseilles and al. 2009
   ! ' A parallel Vlasov solver based on local cubic spline interpolation on patches'
   !
   ! It is meant to be used in the intervals:
   ! x:  [ 0,2*pi/kx]
   ! y:  [ 0,2*pi/ky]
   ! vx: [-6,6]
   ! vy: [-6,6]
 
   ! convention for the params array:
   ! params(1) = eta1_min
   ! params(2) = eta1_max
   ! params(3) = eta2_min
   ! params(4) = eta2_max
   ! params(5) = vxc
   ! params(6) = vyc
   ! params(7) = alpha
   ! params(8) = beta
   !---------------------------------------------------------------------------
 
   function sll_periodic_periodic_gaussian2009_initializer_4d(x, y, vx, vy, params) &
      result(val)
 
      sll_real64 :: val !sll_gaussian_initializer_4d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min, eta1_max
      sll_real64 :: eta2_min, eta2_max
      sll_real64 :: vxc
      sll_real64 :: vyc
      sll_real64 :: alpha
      sll_real64 :: beta
      sll_real64 :: kx, ky
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_periodic_periodic_gaussian2009_initializer_4d, error: the params array must ', &
!!$            'be passed: ', &
!!$            'params(1) = eta1_min, params(2) = eta1_max, params(3) = eta2_min, params(4) = eta2_max',&
!!$            'params(5) = vxc, params(6) = vyx, params(7) = alpha, params(8) = beta'
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
      vxc = params(5)
      vyc = params(6)
      alpha = params(7)
      beta = params(8)
      kx = 2.*sll_p_pi/(eta1_max - eta1_min)
      ky = 2.*sll_p_pi/(eta2_max - eta2_min)
 
      val = beta*(1.0_f64 + alpha*cos(kx*x)*cos(ky*y)) &
            *exp(-0.5*((vx - vxc)**2 + (vy - vyc)**2)) &
            /(2*sll_p_pi)
 
   end function sll_periodic_periodic_gaussian2009_initializer_4d
 
   !---------------------------------------------------------------------------
   !
   !                         Periodic-Periodic initializer another
   !
   ! 4D distribution in [0,2*pi/kx]X[0,2*pi/ky][-6,6]X[-6,6]  with the property of being
   ! periodic in the spatial directions (x1,x2) and sll_m_gaussian in velocity space.
   !
   ! f(x,y,vx,vy) = (1 + alpha * (cos(kx*x) + cos(ky*y)) )
   !                * exp(-0.5*((vx-vxc)^2+(vy-vyc)^2))/ (2*sll_p_pi)
   !
   !  This function is described in the article of Filbet and Sonnendrucker 2002
   ! 'Comparison of Eulerian Vlasov solvers'
   !
   ! It is meant to be used in the intervals:
   ! x:  [ 0,2*pi/kx]
   ! y:  [ 0,2*pi/ky]
   ! vx: [-6,6]
   ! vy: [-6,6]
 
   ! convention for the params array:
   ! params(1) = eta1_min
   ! params(2) = eta1_max
   ! params(3) = eta2_min
   ! params(4) = eta2_max
   ! params(5) = vxc
   ! params(6) = vyc
   ! params(7) = alpha
   ! params(8) = beta
   !---------------------------------------------------------------------------
 
   function sll_periodic_periodic_gaussian2002_initializer_4d(x, y, vx, vy, params) &
      result(val)
 
      sll_real64 :: val !sll_gaussian_initializer_4d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eta1_min, eta1_max
      sll_real64 :: eta2_min, eta2_max
      sll_real64 :: vxc
      sll_real64 :: vyc
      sll_real64 :: alpha
      sll_real64 :: beta
      sll_real64 :: kx, ky
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_periodic_periodic_gaussian2002_initializer_4d, error: the params array must ', &
!!$            'be passed: ', &
!!$            'params(1) = eta1_min, params(2) = eta1_max, params(3) = eta2_min, params(4) = eta2_max',&
!!$            'params(5) = vxc, params(6) = vyx, params(7) = alpha, params(8) = beta'
!!$       stop
!!$    end if
 
      eta1_min = params(1)
      eta1_max = params(2)
      eta2_min = params(3)
      eta2_max = params(4)
      vxc = params(5)
      vyc = params(6)
      alpha = params(7)
      beta = params(8)
      kx = 2.*sll_p_pi/(eta1_max - eta1_min)
      ky = 2.*sll_p_pi/(eta2_max - eta2_min)
 
      val = beta*(1.0_f64 + alpha*(cos(kx*x) + cos(ky*y))) &
            *exp(-0.5*((vx - vxc)**2 + (vy - vyc)**2)) &
            /(2*sll_p_pi)
 
   end function sll_periodic_periodic_gaussian2002_initializer_4d
 
   !---------------------------------------------------------------------------
   !
   !                         Gaussian beam 4d initializer
   !
   ! convention for the params array:
   ! params(1) = vth
   ! params(2) = xth
   ! params(3) = vxc
   ! params(4) = vyc
   ! params(5) = xc
   ! params(6) = yc
   ! params(7) = n0
   ! params(8) = radius
   !---------------------------------------------------------------------------
 
   function sll_f_gaussian_beam_initializer_4d(x, y, vx, vy, params) &
      result(val)
 
      sll_real64 :: val
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, intent(in) :: vx
      sll_real64, intent(in) :: vy
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: vxc
      sll_real64 :: vyc
      sll_real64 :: xc, yc, radius
      sll_real64 :: vt, xt, n0
 
!!$    if( .not. present(params) ) then
!!$       print *, 'sll_gaussian_initializer_4d, error: the params array must ', &
!!$            'be passed: ', &
!!$            'params(1) = vt, params(2) = xt, params(3) = sigma_x, params(4) = sigma_v',&
!!$            ' params(5) = vxc, params(6) = vyc, params(7) = xc, params(8) = yc, params(9) = n0'
!!$       stop
!!$    end if
 
      vt = params(1)
      xt = params(2)
      vxc = params(3)
      vyc = params(4)
      xc = params(5)
      yc = params(6)
      n0 = params(7)
      radius = params(8)
 
      val = n0*exp(-0.5*(radius**2*(x - xc)**2 + radius**2*(y - yc)**2)/(xt*xt)) &
            /(2*sll_p_pi*xt**2) &
            *exp(-0.5*((vx - vxc)**2 + (vy - vyc)**2)/(vt*vt))/(2*sll_p_pi*vt**2)
 
   end function sll_f_gaussian_beam_initializer_4d
 
   function sll_f_langmuir_cos_initializer_2d(x, vx, params)
      sll_real64 :: sll_f_langmuir_cos_initializer_2d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
      sll_real64 :: sum
      sll_int32  :: ns
      sll_real64 :: vb2
      sll_real64 :: vb1
      sll_real64 :: dvb
      sll_real64 :: vj
      sll_real64 :: vtb
      sll_real64 :: db
      sll_real64 :: ve
      sll_int32 :: j
      sll_real64 :: sigma
 
!!$    if( .not. present(params) ) then
!!$       print *, 'langmuir_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx', &
!!$            'params(3) = ns, params(4) = vb2, params(5) = vb1', &
!!$            ' params(6) = vtb, params(7) = db', &
!!$            'params(8) = ve, params(9) = sigma.'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 9)
      eps = params(1)
      kx = params(2)
      ns = int(params(3), i32)
      vb2 = params(4)
      vb1 = params(5)
      vtb = params(6)
      db = params(7)
      ve = params(8)
      sigma = params(9)
 
      sum = 0.0_f64
      do j = 0, ns - 1
         dvb = (vb2 - vb1)/(ns - 1)
         vj = vb2 - j*dvb
         sum = sum + exp(-(vx - vj)**2/(2.0_f64*vtb**2*sigma**2))
      end do
 
      factor1 = 1.0_f64/(sigma*sqrt(2.0_f64*sll_p_pi))
      !langmuir initializer with cosine function
      sll_f_langmuir_cos_initializer_2d = factor1*(1.0_f64 + eps*cos(kx*x)) &
                                          *((1.0_f64 - db)*exp(-(vx - ve)**2/(2.0_f64*sigma**2)) + db*sum/(ns*vtb))
 
   end function sll_f_langmuir_cos_initializer_2d
 
   function sll_f_langmuir_gaussxv_initializer_2d(x, vx, params)
      sll_real64 :: sll_f_langmuir_gaussxv_initializer_2d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
      sll_real64 :: sum
      sll_int32  :: ns
      sll_real64 :: vb2
      sll_real64 :: vb1
      sll_real64 :: dvb
      sll_real64 :: vj
      sll_real64 :: vtb
      sll_real64 :: db
      sll_real64 :: ve
      sll_int32 :: j
      sll_real64 :: sigma
 
!!$    if( .not. present(params) ) then
!!$       print *, 'langmuir_initializer_2d, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx', &
!!$            'params(3) = ns, params(4) = vb2, params(5) = vb1', &
!!$            ' params(6) = vtb, params(7) = db', &
!!$            'params(8) = ve, params(9) = sigma.'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 9)
      eps = params(1)
      kx = params(2)
      ns = int(params(3), i32)
      vb2 = params(4)
      vb1 = params(5)
      vtb = params(6)
      db = params(7)
      ve = params(8)
      sigma = params(9)
 
      sum = 0.0_f64
      do j = 0, ns - 1
         dvb = (vb2 - vb1)/(ns - 1)
         vj = vb2 - j*dvb
         sum = sum + exp(-(vx - vj)**2/(2.0_f64*vtb**2*sigma**2))
      end do
 
      factor1 = 1.0_f64/(sigma*sqrt(2.0_f64*sll_p_pi))
 
      ! langmuir initializer with gaussian_deviate
      sll_f_langmuir_gaussxv_initializer_2d = factor1*(1.0_f64 + eps*sll_f_gaussian_deviate()) &
                                              *((1.0_f64 - db)*exp(-(vx - ve)**2/(2.0_f64*sigma**2)) + db*sum/(ns*vtb))
 
   end function sll_f_langmuir_gaussxv_initializer_2d
 
   function sll_f_langmuir_parameter_pert_initializer_2d(x, vx, params)
      sll_real64 :: sll_f_langmuir_parameter_pert_initializer_2d
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: vx
 
      sll_real64, dimension(:), intent(in) :: params
      sll_real64 :: eps
      sll_real64 :: kx
      sll_real64 :: factor1
      sll_real64 :: sum
      sll_int32  :: ns
      sll_real64 :: vb2
      sll_real64 :: vb1
      sll_real64 :: dvb
      sll_real64 :: vj
      sll_real64 :: vtb
      sll_real64 :: db
      sll_real64 :: ve
      sll_int32 :: j
      sll_real64 :: sigma
      sll_real64 :: pert
 
!!$    if( .not. present(params) ) then
!!$       print *, 'langmuir_initializer_2d_random, error: the params array must ', &
!!$            'be passed. params(1) = epsilon, params(2) = kx', &
!!$            'params(3) = ns, params(4) = vb2, params(5) = vb1', &
!!$            ' params(6) = vtb, params(7) = db', &
!!$            'params(8) = ve, params(9) = sigma, params(10) = perturbation.'
!!$       stop
!!$    end if
      sll_assert(size(params) >= 10)
      eps = params(1)
      kx = params(2)
      ns = int(params(3), i32)
      vb2 = params(4)
      vb1 = params(5)
      vtb = params(6)
      db = params(7)
      ve = params(8)
      sigma = params(9)
      pert = params(10)
 
      sum = 0.0_f64
      do j = 0, ns - 1
         dvb = (vb2 - vb1)/(ns - 1)
         vj = vb2 - j*dvb
         sum = sum + exp(-(vx - vj)**2/(2.0_f64*vtb**2*sigma**2))
      end do
 
      factor1 = 1.0_f64/(sigma*sqrt(2.0_f64*sll_p_pi))
 
      ! langmuir initializer with given perturbation
      sll_f_langmuir_parameter_pert_initializer_2d = factor1*(1.0_f64 + eps*pert) &
                                                     *((1.0_f64 - db)*exp(-(vx - ve)**2/(2.0_f64*sigma**2)) + db*sum/(ns*vtb))
 
   end function sll_f_langmuir_parameter_pert_initializer_2d
 
end module sll_m_common_array_initializers