diff --git a/src/IO.f90 b/src/IO.f90
index 3d80feefa..d067a84c0 100644
--- a/src/IO.f90
+++ b/src/IO.f90
@@ -1557,13 +1557,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
  case (400_pInt)
    msg = 'matrix inversion error'
  case (401_pInt)
-   msg = 'math_check: quat -> axisAngle -> quat failed'
- case (402_pInt)
-   msg = 'math_check: quat -> R -> quat failed'
- case (403_pInt)
-   msg = 'math_check: quat -> euler -> quat failed'
- case (404_pInt)
-   msg = 'math_check: R -> euler -> R failed'
+   msg = 'math_check failed'
  case (405_pInt)
    msg = 'I_TO_HALTON-error: an input base BASE is <= 1'
  case (406_pInt)
diff --git a/src/math.f90 b/src/math.f90
index b5507d868..1bf903ced 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -180,15 +180,11 @@ subroutine math_init
 
  implicit none
  integer(pInt) :: i
- real(pReal), dimension(3,3) :: R,R2
- real(pReal), dimension(3) ::   Eulers,v
- real(pReal), dimension(4) ::   q,q2,axisangle,randTest
+ real(pReal), dimension(4) ::   randTest
 ! the following variables are system dependend and shound NOT be pInt
  integer :: randSize                                                                                ! gfortran requires a variable length to compile
  integer, dimension(:), allocatable :: randInit                                                     ! if recalculations of former randomness (with given seed) is necessary
                                                                                                     ! comment the first random_seed call out, set randSize to 1, and use ifort
- character(len=64) :: error_msg
-
  write(6,'(/,a)')   ' <<<+-  math init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
@@ -211,7 +207,7 @@ subroutine math_init
  enddo
 
  write(6,'(a,I2)') ' size  of random seed:    ', randSize
- do i =1, randSize
+ do i = 1_pInt,randSize
    write(6,'(a,I2,I14)') ' value of random seed:    ', i, randInit(i)
  enddo
  write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
@@ -221,6 +217,27 @@ subroutine math_init
  call halton_seed_set(int(randInit(1), pInt))
  call halton_ndim_set(3_pInt)
 
+ call math_check()
+
+end subroutine math_init
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief check correctness of (some) math functions
+!--------------------------------------------------------------------------------------------------
+subroutine math_check
+
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use prec,     only: tol_math_check
+ use numerics, only: &
+   fixedSeed
+ use IO,       only: IO_error, IO_timeStamp
+
+ character(len=64) :: error_msg
+
+ real(pReal), dimension(3,3) :: R,R2
+ real(pReal), dimension(3) ::   Eulers,v
+ real(pReal), dimension(4) ::   q,q2,axisangle
+
  ! --- check rotation dictionary ---
 
  q = math_qRand()          ! random quaternion
@@ -230,7 +247,8 @@ subroutine math_init
  q2 = math_axisAngleToQ(axisangle(1:3),axisangle(4))
  if ( any(abs( q-q2) > tol_math_check) .and. &
       any(abs(-q-q2) > tol_math_check) ) then
-   write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
+   write (error_msg, '(a,e14.6)' ) &
+          'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
    call IO_error(401_pInt,ext_msg=error_msg)
  endif
 
@@ -239,8 +257,9 @@ subroutine math_init
  q2 = math_RtoQ(R)
  if ( any(abs( q-q2) > tol_math_check) .and. &
       any(abs(-q-q2) > tol_math_check) ) then
-   write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
-   call IO_error(402_pInt,ext_msg=error_msg)
+   write (error_msg, '(a,e14.6)' ) &
+          'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
+   call IO_error(401_pInt,ext_msg=error_msg)
  endif
 
  ! +++ q -> euler -> q  +++
@@ -248,28 +267,46 @@ subroutine math_init
  q2 = math_EulerToQ(Eulers)
  if ( any(abs( q-q2) > tol_math_check) .and. &
       any(abs(-q-q2) > tol_math_check) ) then
-   write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
-   call IO_error(403_pInt,ext_msg=error_msg)
+   write (error_msg, '(a,e14.6)' ) &
+          'quat -> euler -> quatmaximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
+   call IO_error(401_pInt,ext_msg=error_msg)
  endif
 
  ! +++ R -> euler -> R  +++
  Eulers = math_RtoEuler(R)
  R2 = math_EulerToR(Eulers)
  if ( any(abs( R-R2) > tol_math_check) ) then
-   write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2))
-   call IO_error(404_pInt,ext_msg=error_msg)
+   write (error_msg, '(a,e14.6)' ) &
+          'R -> euler -> R maximum deviation ',maxval(abs( R-R2))
+   call IO_error(401_pInt,ext_msg=error_msg)
  endif
 
  ! +++ check rotation sense of q and R +++
- q = math_qRand()          ! random quaternion
  call halton(3_pInt,v)     ! random vector
  R = math_qToR(q)
  if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
-   write(6,'(a,4(f8.3,1x))') 'q',q
-   call IO_error(409_pInt)
+   write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
+   call IO_error(401_pInt,ext_msg=error_msg)
  endif
 
-end subroutine math_init
+ ! +++ check vector expansion +++
+ if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
+             math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt,3_pInt,0_pInt])) > tol_math_check)) then
+   write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'
+   call IO_error(401_pInt,ext_msg=error_msg)
+ endif
+ if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - &
+             math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt])) > tol_math_check)) then
+   write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]'
+   call IO_error(401_pInt,ext_msg=error_msg)
+ endif
+ if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - &
+             math_expand([1.0_pReal,2.0_pReal],[1_pInt,2_pInt,3_pInt])) > tol_math_check)) then
+   write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
+   call IO_error(401_pInt,ext_msg=error_msg)
+ endif
+ 
+end subroutine math_check
  
 
 !--------------------------------------------------------------------------------------------------
@@ -334,6 +371,30 @@ recursive subroutine math_qsort(a, istart, iend)
 end subroutine math_qsort
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief vector expansion
+!> @details takes a set of numbers (a,b,c,...) and corresponding multipliers (x,y,z,...)
+!> to return a vector of x times a, y times b, z times c, ...
+!--------------------------------------------------------------------------------------------------
+pure function math_expand(what,how)
+
+ implicit none
+ real(pReal), dimension(:), intent(in) ::  what
+ integer(pInt), dimension(:), intent(in) ::  how
+ real(pReal), dimension(sum(how)) ::  math_expand
+ integer(pInt) :: i,j,o
+
+ o = 0_pInt
+ do i = 1, size(how)
+   do j = 1, how(i)
+     o = o + 1_pInt
+     math_expand(o) = what(1+mod(i-1,size(what)))
+   enddo 
+ enddo
+
+end function math_expand
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief range of integers starting at one
 !--------------------------------------------------------------------------------------------------