source: trunk/LMDZ.GENERIC/utilities/aerosols_utilities/math.f90 @ 3436

module math
    ! """
    ! Contains useful math functions.
    ! """
    implicit none

    integer, parameter :: &
        dp = selected_real_kind(15,300)

    real(kind=dp), parameter :: &
        pi = 3.141592653589793238460_dp

    doubleprecision, parameter :: &
        sqrtpi = 0.5641895835477563d0  ! = 1 / sqrt(pi)

    doubleprecision, parameter :: &
        prec_high = 10d0 ** (-precision(0d0)), prec_low = 10d0 ** (-precision(0.))

    save

    contains
        doubleprecision function chi2(observed_data, calculated_data)
            ! """
            ! Evaluate the goodeness of fit of calculated data in respect to observations.
            !
            ! input:
            !   observed_data: array of observed data
            !   calculated_data: array of calculated data
            ! """
            implicit None

            doubleprecision, dimension(:), intent(in) :: observed_data, calculated_data

            chi2 = sum((observed_data(:) - calculated_data(:))**2 / calculated_data(:))
        end function chi2


        doubleprecision function chi2_reduced(observed_data, calculated_data, input_deviation)
            ! """
            ! Evaluate the goodeness of fit of calculated data in respect to observations and observational errors.
            !
            ! input:
            !   observed_data: array of observed data
            !   calculated_data: array of calculated data
            !   input_deviation: array of observational error
            ! """
            implicit None

            doubleprecision, dimension(:), intent(in) :: observed_data, calculated_data, input_deviation

            chi2_reduced = sum((observed_data(:) - calculated_data(:))**2 / input_deviation(:)**2)
        end function chi2_reduced


        doubleprecision function deg2rad(angle)
            ! """
            ! Return an angle in radians from an angle in degrees.
            !
            ! input:
            !   angle: (degree) an angle
            ! """
            implicit none

            doubleprecision, intent(in) :: angle

            deg2rad = angle * (pi / 180.0_dp)
        end function deg2rad


        doubleprecision function ellipse_polar_form(semi_major_axis, semi_minor_axis, angle)
            ! """
            ! Return the ellipse polar form of an ellipse relative to its center.
            !
            ! inputs:
            !   angle: (deg) angle between the semi-major axis, the center, and the point of the ellipse
            !   semi_major_axis: semi-major axis of the ellipse
            !   semi_minor_axis: semi-minor axis of the ellipse
            ! """
            implicit none

            doubleprecision, intent(in) :: angle, semi_major_axis, semi_minor_axis

            doubleprecision :: &
                theta

            theta = deg2rad(angle)

            ellipse_polar_form = semi_major_axis * semi_minor_axis / &
                sqrt((semi_minor_axis * cos(theta))**2 + (semi_major_axis * sin(theta))**2)

            return

        end function ellipse_polar_form


        doubleprecision function gaussian(x, fwhm)
            ! """
            ! Return the gaussian of a value, for a given Full Width Half Maximum.
            !
            ! inputs:
            !   x: value for which to calculate the gaussian
            !   fwhm: Full Width Half Maximum of the gaussian function
            ! """
            implicit none

            doubleprecision, intent(in) :: x, fwhm

            doubleprecision :: &
                sigma

            sigma = fwhm / (2D0 * sqrt(2D0 * log(2D0)))
            gaussian = 1D0 / (sigma * sqrt(2D0 * pi)) * exp(-1D0/2D0 * (x / sigma)**2D0)
        end function gaussian


        doubleprecision function gaussian_noise() result(n)
            ! """
            ! Return a random value following a standard normal distriubtion PDF.
            !
            ! output:
            !   n: random gaussian noise
            ! """
            implicit none

            doubleprecision ::&
                r  ! random number between 0 and 1

            call init_random_seed
            call random_number(r)

            n = sqrt(2D0) * erfinv(2 * r - 1)

            return

        end function gaussian_noise


        doubleprecision function sec(angle)
            ! """
            ! Return the secant of an angle in degrees.
            !
            ! input:
            !   angle: (degree) an angle
            ! """
            implicit none

            doubleprecision, intent(in) :: angle

            sec = 1D0 / cos(deg2rad(angle))
        end function sec


        doubleprecision function sgn(value)
            ! """
            ! Return the sign of a value.
            ! Not to be confused with Fortran built-in function SIGN(A, B).
            ! """
            implicit none

            doubleprecision, intent(in) :: value

            if (value >= 0D0) then
                sgn = 1D0
            else
                sgn = -1D0
            end if

            return
        end function sgn


        doubleprecision function sinc_fwhm(x, fwhm)
            ! """
            ! Return the sine cardinal of a value, for a given Full Width Half Maximum.
            !
            ! inputs:
            !   x: value for which to calculate sinc
            !   fwhm: Full Width Half Maximum of the sinc function
            ! """
            implicit none

            doubleprecision,intent(in) :: x, fwhm

            doubleprecision, parameter :: &
                k = 1D0 / (2D0 * 1.89549D0)  ! sinc fwhm constant

            if (x > 0D0 - tiny(0.) .and. x < 0D0 - tiny(0.)) then
                sinc_fwhm = 1D0  ! avoid NaN at x = 0
            else
                sinc_fwhm = sin(x / (k * fwhm)) / (x / (k * fwhm))
            end if
        end function sinc_fwhm


        function arange(start, stop, step) result(array)
            ! """
            ! Return evenly spaced values within a given interval.
            !
            ! inputs:
            !   start: start of interval; the array starts with this value
            !   stop: end of interval; the array does not inlude this value
            !   step: spacing between values
            !
            ! outputs:
            !   array: array of evenly spaced values
            ! """
            implicit none

            doubleprecision, intent(in) :: start, stop, step
            doubleprecision, dimension(:), allocatable :: array

            integer :: &
                i, &  ! index
                n     ! number of elements in array

            n = ceiling((stop - start) / step)

            allocate(array(n))

            do i = 1, n
                array(i) = start + (i - 1) * step
            end do

            return
        end function arange


        function arange_include(start, stop, step) result(array)
            ! """
            ! Return evenly spaced values including a given interval.
            !
            ! inputs:
            !   start: start of interval; the array starts with this value
            !   stop: end of interval; the last element of array is always greater than this value by a max of 1 step
            !   step: spacing between values
            !
            ! outputs:
            !   array: array of evenly spaced values
            ! """
            implicit none

            doubleprecision, intent(in) :: start, stop, step
            doubleprecision, dimension(:), allocatable :: array

            integer :: &
                i, &  ! index
                n     ! number of elements in array

            n = ceiling((stop - start) / step) + 1  ! operation +1 ensures that stop is included in array

            allocate(array(n))

            do i = 1, n
                array(i) = start + (i - 1) * step
            end do

            return
        end function arange_include


        function convolve(signal, filter) result(convolved_signal)
            ! """
            ! Convolve the signal by the filter using classical convolution.
            !
            ! inputs:
            !   signal: the signal array
            !   filter: the filter array
            !
            ! output:
            !   convolved_signal: signal convolved by filter
            !
            ! source: https://fortrandev.wordpress.com/2013/04/01/fortran-convolution-algorithm/
            ! """
            implicit none

            doubleprecision, dimension(:), intent(in) :: signal, filter
            doubleprecision, dimension(:), allocatable :: convolved_signal, convolution

            integer :: &
                i, &
                j, &
                k, &
                size_filter, &
                size_signal, &
                start_index

            size_signal = size(signal)
            size_filter = size(filter)

            start_index = floor(size_filter / 2d0)

            allocate(convolved_signal(size_signal), convolution(size_signal + size_filter))

            ! Last part
            do i = size_signal, size_signal + size_filter
                convolution(i) = 0D0
                j = size_signal

                do k = 1, size_filter
                    convolution(i) = convolution(i) + signal(j) * filter(k)
                    j = j - 1
                end do
            end do

            ! Middle part
            do i = size_filter, size_signal
                convolution(i) = 0D0
                j = i

                do k = 1, size_filter
                    convolution(i) = convolution(i) + signal(j) * filter(k)
                    j = j - 1
                end do
            end do

            ! First part
            do i = 1, size_filter
                convolution(i) = 0D0
                j = i
                k = 1

                do while (j > 0)
                    convolution(i) = convolution(i) + signal(j) * filter(k)
                    j = j - 1
                    k = k + 1
                end do
            end do

            convolved_signal(:) = convolution(start_index:start_index + size_signal - 1)
        end function convolve


        function slide_convolve(signal, filter) result(convolved_signal)
            ! """
            ! Slide convolve the signal by a sliding filter using classical convolution.
            !
            ! inputs:
            !   signal: the signal array
            !   filter: the filter 2D-array
            !
            ! output:
            !   convolved_signal: signal convolved by filter
            ! """
            implicit none

            doubleprecision, dimension(:), intent(in) :: signal
            doubleprecision, dimension(:, :), intent(in) :: filter
            doubleprecision, dimension(:), allocatable :: convolved_signal, convolution

            integer :: &
                i, &  ! index
                j, &  ! index
                k, &  ! index
                size_filter, &  ! size of the filter
                size_signal, &  ! size of the signal
                start_index     ! start index of the result

            size_signal = size(signal)
            size_filter = size(filter, 1)

            ! Check sizes
            ! TODO [low] remove the size constraint
            if (size_signal < size_filter) then
                print '("ERROR: slide_convolve: filter size must be lower than signal size, &
                        &but sizes are ", I10, "and ", I10)', size_filter, size_signal
                print '("Catched error: ", ES15.8)', signal(size_filter)
                stop  ! be sure the program stops here
            end if

            start_index = floor(size_filter / 2D0) + 1

            allocate(convolved_signal(size_signal), convolution(size_signal + size_filter))

            convolution(:) = 0D0

            ! Last part
            do i = size_signal + 1, start_index + size_signal - 1
                j = size_signal

                do k = 1, size_filter
                    convolution(i) = convolution(i) + signal(j) * filter(k, j)
                    j = j - 1
                end do
            end do

            ! Middle part
            do i = size_filter, size_signal
                j = i

                do k = 1, size_filter
                    convolution(i) = convolution(i) + signal(j) * filter(k, j)
                    j = j - 1
                end do
            end do

            ! First part
            do i = start_index + 1, size_filter - 1
                j = i

                do k = 1, i
                    convolution(i) = convolution(i) + signal(j) * filter(k, j)
                    j = j - 1
                end do
            end do

            ! First index
            i = start_index
            convolution(i) = 0D0
            j = i

            do k = 1, i
                convolution(i) = convolution(i) + signal(j) * filter(k, j)
                j = j - 1
            end do

            k = i + 1
            convolution(i) = convolution(i) + signal(1) * filter(k, 1)

            convolved_signal(:) = convolution(start_index:start_index + size_signal - 1)
        end function slide_convolve


        function search_sorted(array, value) result(index)
            ! """
            ! Find the index into a sorted array such that the corresponding value is the closest to 'value'.
            ! """
            implicit none

            doubleprecision, intent(in) :: value
            doubleprecision, dimension(:), intent(in) :: array

            integer :: index

            integer :: &
                i_low, &
                i_high, &
                i_mid

            i_low = 1
            i_mid = 0
            i_high = size(array)

            if(value < array(1)) then
                index = 1
                return
            elseif(value > array(i_high)) then
                index = i_high
                return
            end if

            do while(i_low <= i_high)
                i_mid = i_low + (i_high - i_low) / 2

                if(value < array(i_mid)) then
                    i_high = i_mid - 1
                elseif(value > array(i_mid)) then
                    i_low = i_mid + 1
                else
                    index = i_mid
                    return
                end if
            end do

            if(array(i_low) - value < value - array(i_high)) then
                index = i_low
            else
                index = i_high
            end if

            return
        end function search_sorted


        function erfinv(x) result(r)
            ! """
            ! Calculate the inverse of the erf function.
            !
            ! input:
            !   x: a number between -1 and 1
            !
            ! output:
            !   r: x = erf(r)
            ! """
            implicit none

            doubleprecision, parameter :: &  ! erfinv-approximation factors
                a0 = 0.886226899, &
                a1 = -1.645349621, &
                a2 = 0.914624893, &
                a3 = -0.140543331, &
                b0 = 1, &
                b1 = -2.118377725, &
                b2 = 1.442710462, &
                b3 = -0.329097515, &
                b4 = 0.012229801, &
                c0 = -1.970840454, &
                c1 = -1.62490649, &
                c2 = 3.429567803, &
                c3 = 1.641345311, &
                d0 = 1, &
                d1 = 3.543889200, &
                d2 = 1.637067800

            doubleprecision :: x, r

            integer ::&
                sign_x  ! sign of elements of array x

            doubleprecision :: &
                x2, &  ! square of x
                y  ! intermediate value

            if(x < -1 .or. x > 1) then
                print '("ERROR: erfinv(x): x must be in [-1;1]")'
                stop
            end if

            if (x > -tiny(0.) .and. x < tiny(0.)) then
                r = 0
                return
            end if

            if (x > 0) then
                sign_x = 1
            else
                sign_x = -1
                x = -x
            end if

            if (x <= 0.7) then
                x2 = x * x
                r = x * (((a3 * x2 + a2) * x2 + a1) * x2 + a0)
                r = r / (((b4 * x2 + b3) * x2 + b2) * x2 + b1) * x2 + b0
            else
                y = sqrt(-log((1 - x) / 2))
                r = (((c3 * y + c2) * y + c1) * y + c0)
                r = r / ((d2 * y + d1) * y + d0)
            end if

            r = r * sign_x
            x = x * sign_x

            r = r - (erf(r) - x) / (2 / sqrt(PI) * exp (-r * r))
            r = r - (erf(r) - x) / (2 / sqrt(PI) * exp (-r * r))

            return

        end function erfinv


        function interp(x_new, x, y) result(y_new)
            ! """
            ! Interpolate array y of abscisse x to array y_new of abscisse x_new.
            !
            ! inputs:
            !   x_new: abscisses on which y will be interpolated
            !   x: abscisses of y
            !   y: array to interpole
            !
            ! output:
            !   y_new: interpolation of y(x) on abscisses x_new
            ! """

            implicit none

            doubleprecision, dimension(:), intent(in) :: x, y, x_new
            doubleprecision, dimension(size(x_new)) :: y_new

            integer ::&
                i, &  ! index
                j     ! index

            i = 1
            j = 1
            y_new(:) = 0D0

            if(size(x) /= size(y)) then
                print '("ERROR: interp: x and y must have the same size")'
                stop
            end if

            if(x(1) < x(size(x))) then  ! ascending numerical order x array
                if(x_new(1) < x(1) .or. x_new(size(x_new)) > x(size(x))) then
                    print '("ERROR: interp: x_new is outside x boundaries")'
                    print *, x(1), ' < ' , x_new(1),  '--', x_new(size(x_new)), ' < ', x(size(x))
                    stop
                end if

                do while (i <= size(x_new))
                    if(x_new(i) <= x(j + 1)) then
                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)

                        i = i + 1
                    else
                        j = j + 1
                    end if
                end do
            else  ! descending numerical order x array
                if(x_new(1) > x(1) .or. x_new(size(x_new)) < x(size(x))) then
                    print '("ERROR: interp: x_new is outside x boundaries")'
                    print *,x(1), ' > ', x_new(1), '--', x_new(size(x_new)), ' > ', x(size(x))
                    stop
                end if

                do while (i <= size(x_new))
                    if(x_new(i) >= x(j + 1)) then
                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)

                        i = i + 1
                    else
                        j = j + 1
                    end if
                end do
            end if

            return
        end function interp


        function interp_fast(x_new, x, y) result(y_new)
            ! """
            ! Interpolate array y of abscisse x to array y_new of abscisse x_new.
            ! Assumptions (no check performed):
            !   - x_new is within x boundaries
            !   - x and y have the same size
            !
            ! inputs:
            !   x_new: abscisses on which y will be interpolated
            !   x: abscisses of y
            !   y: array to interpole
            !
            ! output:
            !   y_new: interpolation of y(x) on abscisses x_new
            ! """

            implicit none

            doubleprecision, dimension(:), intent(in) :: x, y, x_new
            doubleprecision, dimension(size(x_new)) :: y_new

            integer ::&
                i, &  ! index
                j     ! index

            i = 1
            j = 1
            y_new(:) = 0D0

            if(x(1) < x(size(x))) then  ! ascending numerical order x array
                do while(i <= size(x_new))
                    if(x_new(i) <= x(j + 1)) then
                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)

                        i = i + 1
                    else
                        j = j + 1
                    end if
                end do
            else  ! descending numerical order x array
                do while(i <= size(x_new))
                    if(x_new(i) >= x(j + 1)) then
                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)

                        i = i + 1
                    else
                        j = j + 1
                    end if
                end do
            end if

            return
        end function interp_fast


        function interp_ex(x_new, x, y) result(y_new)
            ! """
            ! Interpolate array y of abscisse x to array y_new of abscisse x_new.
            ! The points outside x range are linearly extrapolated.
            ! Data must be in strict increasing order.
            !
            ! inputs:
            !   x_new: abscisses on which y will be interpolated
            !   x: abscisses of y
            !   y: array to interpole
            !
            ! output:
            !   y_new: interpolation of y(x) on abscisses x_new
            ! """

            implicit none

            doubleprecision, dimension(:), intent(in) :: x, y, x_new
            doubleprecision, dimension(size(x_new)) :: y_new

            integer ::&
                i, &  ! index
                i_max, &
                i_min, &
                n     ! index

            n = size(x)
            i_min = 0
            i_max = 0

            ! Check order
            if(x(1) > x(n)) then
                if(x_new(1) < x_new(size(x_new))) then
                    write(*, '("Error: interp_ex: x_new numerical order must be the same than x (descending)")')

                    stop
                end if

                ! Find where x_new is in data range
                do i = 1, size(x_new)
                    if(x_new(i) <= x(1) .and. i_min == 0) then
                        i_min = i
                    end if

                    if(x_new(i) >= x(n)) then
                        i_max = i
                    end if

                    if(x_new(i) < x(n)) then
                        exit
                    end if
                end do
            else
                if(x_new(1) > x_new(size(x_new))) then
                    write(*, '("Error: interp_ex: x_new numerical order must be the same than x (ascending)")')

                    stop
                end if

                ! Find where x_new is in data range
                do i = 1, size(x_new)
                    if(x_new(i) >= x(1) .and. i_min == 0) then
                        i_min = i
                    end if

                    if(x_new(i) <= x(n)) then
                        i_max = i
                    end if

                    if(x_new(i) > x(n)) then
                        exit
                    end if
                end do
            end if

            ! Interpolation within data range
            if(i_min > 0 .and. i_max > 0) then
                y_new(i_min:i_max) = interp_fast(x_new(i_min:i_max), x, y)
            end if

            if(i_min == 0 .or. i_max == 0) then
                if(i_min > 0) then
                    ! Higher bound extrapolation
                    do i = 1, size(x_new)
                        y_new(i) = (x_new(i) - x(n - 1)) * (y(n) - y(n - 1)) / (x(n) - x(n - 1)) + y(n - 1)
                    end do
                else if (i_max > 0) then
                    ! Lower bound extrapolation
                    do i = 1, size(x_new)
                        y_new(i) = (x_new(i) - x(1)) * (y(2) - y(1)) / (x(2) - x(1)) + y(1)
                    end do
                else
                    write(*, '("Error: interp_ex: unable to determine boundary")')

                    stop
                end if
            else
                ! Lower bound extrapolation
                if(i_min > 1) then
                    do i = 1, i_min - 1
                        y_new(i) = (x_new(i) - x(1)) * (y(2) - y(1)) / (x(2) - x(1)) + y(1)
                    end do
                end if

                ! Higher bound extrapolation
                if(i_max < size(x_new)) then
                    do i = i_max + 1, size(x_new)
                        y_new(i) = (x_new(i) - x(n - 1)) * (y(n) - y(n - 1)) / (x(n) - x(n - 1)) + y(n - 1)
                    end do
                end if
            end if
        end function interp_ex


        function interp_ex_0d(x_new, x, y) result(y_new)
            ! """
            ! Interpolate array y of abscisse x to value y_new at x_new.
            ! This function only takes a double as x_new, not an array. Return a double, not an array.
            ! The points outside x range are linearly extrapolated.
            ! Data must be in strict increasing order.
            !
            ! inputs:
            !   x_new: value at which y will be interpolated
            !   x: abscisses of y
            !   y: array to interpole
            !
            ! output:
            !   y_new: interpolation of y(x) at x_new
            ! """
            implicit none

            doubleprecision, dimension(:), intent(in) :: x, y
            doubleprecision :: x_new, y_new

            doubleprecision, dimension(1) :: arr_tmp

            arr_tmp = interp_ex([x_new], x, y)
            y_new = arr_tmp(1)
        end function interp_ex_0d


        function mean_restep(x_new, x, y) result(y_new)
            ! """
            ! Change the values of y of abscisse x so that the new values y_new of abcisse x_new are the mean of y
            ! between a step of x_new centered on x_new.
            ! Both x and x_new must be regularly spaced.
            ! Example:
            ! >>> x_new = [0, 3, 6]
            ! >>> x = [0, 1, 2, 3, 4, 5, 6]
            ! >>> y = [1, 2, 4, 3, 6, 2, 0]
            ! >>> mean_restep(x_new, x, y)
            ! >>> [1.50, 4.33, 1.00]
            !
            ! inputs:
            !   x_new: abscisses on which y will be interpolated
            !   x: abscisses of y
            !   y: array to interpole
            !
            ! output:
            !   y_new: interpolation of y(x) on abscisses x_new
            ! """
            implicit none

            doubleprecision, dimension(:), intent(in) :: x, y, x_new
            doubleprecision, dimension(size(x_new)) :: y_new

            integer :: &
                i, &  ! index
                j, &  ! index
                n     ! number of elements

            doubleprecision :: &
                sum_x, &
                min_x, &
                max_x

            call check_inputs()

            j = 1

            ! First value
            i = 1
            n = 0
            sum_x = 0D0
            max_x = x_new(i) + (x_new(i + 1) - x_new(i)) / 2D0
            min_x = x_new(i)

            do while (x(j) < max_x .and. j < size(x))
                if (x(j) >= min_x) then
                    sum_x = sum_x + y(j)
                    n = n + 1
                end if

                j = j + 1
            end do

            y_new(i) = sum_x / n

            ! Intermediate values
            do i = 2, size(x_new) - 1
                n = 0
                sum_x = 0D0
                max_x = x_new(i) + (x_new(i + 1) - x_new(i)) / 2D0
                min_x = x_new(i) - (x_new(i) - x_new(i - 1)) / 2D0

                do while (x(j) < max_x .and. j < size(x))
                    if (x(j) >= min_x) then
                        sum_x = sum_x + y(j)
                        n = n + 1
                    end if

                    j = j + 1
                end do

                y_new(i) = sum_x / n
            end do

            ! Last value
            i = size(x_new)
            n = 0
            sum_x = 0D0
            max_x = x_new(i)
            min_x = x_new(i) - (x_new(i) - x_new(i - 1)) / 2D0

            do while (x(j) < max_x .and. j < size(x))
                if (x(j) >= min_x) then
                    sum_x = sum_x + y(j)
                    n = n + 1
                end if

                j = j + 1
            end do

            if(n == 0) then
                y_new(i) = y(j)
            else
                y_new(i) = sum_x / n
            end if

            return

            contains
                subroutine check_inputs()
                    implicit none

                    if (size(x) /= size(y)) then
                        print '("ERROR: mean_restep: x and y must have the same size")'
                        stop
                    end if

                    if (x(1) > x(size(x))) then
                        print '("ERROR: mean_restep: x must be in increasing order")'
                        stop
                    end if

                    if(x_new(1) < x(1) .or. x_new(size(x_new)) > x(size(x))) then
                        print '("ERROR: interp: x_new is outside x boundaries")'
                        print *,x_new(1),x(1),x_new(size(x_new)),x(size(x))
                        stop
                    end if
                end subroutine check_inputs
        end function mean_restep


        function reverse_array(array) result(reversed_array)
            ! """
            ! Reverse a 1-D double precision array.
            ! """
            implicit none

            doubleprecision, dimension(:), intent(in) :: array
            integer :: head, tail

            doubleprecision, dimension(size(array)) :: reversed_array

            head = 1
            tail = size(array)
            reversed_array(:) = 0d0

            do while(head < tail)
                reversed_array(tail) = array(head)
                reversed_array(head) = array(tail)
                head = head + 1
                tail = tail - 1
            end do

            return
        end function reverse_array


        function voigt(x, y)
            ! """
            ! Calculate the voigt function using an algorithm written by Humlicek JQSRT, 27, 437 (1982).
            ! Calculate the complex probability function W(z) = exp(-z^2) * erfc(-z^2) in the superior complex plan
            ! (i.e. for y >= 0). The real part of this function is the Voigt function.
            !
            ! The article shows that the rational function W can be written as (Eq. 11):
            !   W_n(z) = (-iz) * sum_{k=1}^{n/2}(c_k((-iz)^2)^{k-1}) / ((-iz)^n + sum_{k=1}^{n/2}(d_k((-iz)^2)^{k-1}))
            ! With z = x + iy, (-iz = t). The coefficent c and d are given in Table 3 of the article.
            ! The maximal relative error on the imaginary and real parts is < 1e-4.
            !
            ! inputs:
            !   x: real part coordinate
            !   y: imaginary part coordinate
            !
            ! output:
            !   voigt: value of the voigt function at z = x + iy
            ! """
            implicit none

            doubleprecision, dimension(:), intent(in) :: x
            doubleprecision, intent(in) :: y
            doubleprecision, dimension(size(x)) :: voigt

            integer :: &
                i  ! index

            double precision, dimension(size(x)) :: &
                s  ! intermediate value

            complex (kind=8), dimension(size(x)) :: &
                t, &  ! intermediate value
                u, &  ! intermediate value
                w     ! intermediate value

            do i=1,size(x)
                t(i) = cmplx(y, -x(i), kind=8)
            end do

            s(:) = abs(x(:)) + y

            where (s >= 15D0)  ! n = 2 (region i)
                w = t * sqrtpi / (0.5D0 + t**2)
                voigt = dble(w * sqrtpi)
            elsewhere (s >= 5.5D0)  ! n = 4 (region II)
                u = t**2
                w = t * (1.410474D0 + u * sqrtpi) / (0.75D0 + u * (3D0 + u))
                voigt = dble(w * sqrtpi)
            elsewhere (0.195D0 * abs(x) - 0.176D0 <= y)  ! n = 6 (region III)
                w = (16.4955D0 + t * (20.20933D0 + t * (11.96482D0 + t * (3.778987D0 + t * 0.5642236D0)))) / &
                         (16.4955D0 + t * (38.82363D0 + t * (39.27121D0 + t * (21.69274D0 + t * (6.699398D0 + t)))))
                voigt = dble(w * sqrtpi)
            elsewhere  ! n = 8 (region IV)
                u = t**2
                w = exp(u) - &
                         t * (36183.31D0 - &
                              u * (3321.9905D0 - &
                                   u * (1540.787D0 - &
                                        u * (219.0313D0 - u *(35.76683D0 - u * (1.320522D0 - u * sqrtpi)))))) / &
                         (32066.6D0 - &
                          u * (24322.84D0 - &
                               u * (9022.228D0 - &
                                    u * (2186.181D0 - u * (364.2191D0 - u * (61.57037D0 - u * (1.841439D0 - u)))))))
                voigt = dble(w * sqrtpi)
            end where

            return
        end function voigt


        doubleprecision function voigt_from_data(x, y, v) result(voigt)
            ! """
            ! """
            implicit none

            doubleprecision, intent(in) :: &
                    x, &
                    y, &
                    v(400, 100)

            doubleprecision, parameter :: &
                    dp = 0.0025d0, &
                    ds = 0.01d0, &
                    yl = 1d0 / 99d0, &
                    xl = 399d0

            integer :: &
                    ip, &
                    is

            doubleprecision :: &
                    a, &
                    b, &
                    pressure_space, &
                    x2, &
                    x2y2, &
                    s, &
                    y2, &
                    f

            if(y > 5.4d0) then
                x2 = x * x
                y2 = y * y
                x2y2 = x2 + y2
                voigt = sqrtpi * y * (1d0 + (3d0 * x2 - y2) / (2d0 * x2y2 * x2y2) + &
                        0.75d0 * (5d0 * x2 * x2 + y2 * y2 - 10d0 * x2 * y2) / (x2y2 * x2y2 * x2y2 * x2y2)) / x2y2
            else
                if(y <= yl) then
                    b = exp(-x * x)
                    pressure_space = 1d0 / (1d0 + x)
                    ip = int(pressure_space / dp)
                    f = pressure_space / dp - ip
                    if(ip >= 1) then
                        if(ip < 400) then
                            a = (1d0 - f) * v(ip, 1) + f * v(ip + 1, 1)
                        else
                            a = v(ip, 1)
                        end if
                    else
                        a = f * v(ip + 1, 1)
                    end if
                    a = a * a
                    f = y / yl
                    voigt = (1d0 - f) * b + f * a
                else
                    if(x > xl) then
                        s = y / (1d0 + y)
                        is = int(s / ds)
                        b = v(1, is) * xl / x
                        b = b * b
                        a = v(1, is + 1) * xl / x
                        a = a * a
                        f = s / ds - is
                        voigt = (1d0 - f) * b + f * a
                    else
                        s = y / (1d0 + y)
                        pressure_space = 1d0 / (1d0 + x)
                        is = int(s / ds)
                        ip = int(pressure_space / dp)
                        f = pressure_space / dp - ip
                        if(ip < 400) then
                            b = (1d0 - f) * v(ip, is) + f * v(ip + 1, is)
                            a = (1d0 - f) * v(ip, is + 1) + f * v(ip + 1, is + 1)
                        else
                            b = v(ip, is)
                            a = v(ip, is + 1)
                        end if
                        b = b * b
                        a = a * a
                        f = s / ds - is
                        voigt = (1d0 - f) * b + f * a
                    end if
                end if
            end if

            return
        end function voigt_from_data


        recursive subroutine fft(x)
            ! """
            ! Calculate the Cooley-Tukey FFT functions.
            !
            ! input:
            !   x: double complex vector of size 2^n
            !
            ! notes:
            !   Source: https://rosettacode.org/wiki/Fast_Fourier_transform#Fortran
            ! """
            implicit none

            complex(kind=dp), dimension(:), intent(inout) :: x
            complex(kind=dp) ::&
                t

            integer ::&
                i, &  ! index
                n     ! size of array x

            complex(kind=dp), dimension(:), allocatable ::&
                even, &  ! even-number indexed values of x
                odd      ! odd-number indexed values of x

            n = size(x)

            if(n <= 1) return

            allocate(odd((n+1)/2))
            allocate(even(n/2))

            ! divide
            odd(:) = x(1:n:2)
            even(:) = x(2:n:2)

            ! conquer
            call fft(odd)
            call fft(even)

            ! combine
            do i = 1, n/2
               t = exp(cmplx(0.0_dp, -2.0_dp * pi * real(i-1, dp) / real(n, dp), kind=dp)) * even(i)
               x(i) = odd(i) + t
               x(i+n/2) = odd(i) - t
            end do

            deallocate(odd)
            deallocate(even)

        end subroutine fft


        recursive subroutine quicksort(array)
            ! """
            ! Sort double precision numbers into ascending numerical order.
            !
            ! input:
            !   array: double precision array to sort
            !
            ! output:
            !   array: sorted double precision array
            !
            ! notes:
            !   Author: t-nissie, some tweaks by 1AdAstra1
            !   Source: https://gist.github.com/t-nissie/479f0f16966925fa29ea
            ! """
            implicit none

            double precision, intent(inout), dimension(:) :: array

            double precision :: &
                x, &  ! pivot point
                tmp   ! temporary value

            integer :: &
                first = 1, & ! index of the beginning of the array
                last, &      ! index of the end of the array
                insertion_size_threshold = 32

            integer :: &
                i, &
                j

            last = size(array, 1)

            if(last < insertion_size_threshold) then  ! use insertion sort on small arrays
                do i = 2, last
                    tmp = array(i)

                    do j = i - 1, 1, -1
                        if(array(j) < tmp) then
                            exit
                        else
                            array(j + 1) = array(j)
                        end if
                    end do

                    array(j + 1) = tmp
                end do
            else
                x = array((first+last)/2)
                i = first
                j = last

                do
                    do while (array(i) < x)
                        i = i + 1
                    end do

                    do while (x < array(j))
                        j = j - 1
                    end do

                    if (i >= j) then
                        exit
                    end if

                    tmp = array(i)
                    array(i) = array(j)
                    array(j) = tmp

                    i = i + 1
                    j = j - 1
                end do

                if (first < i - 1) call quicksort(array(first:i-1))
                if (j + 1 < last)  call quicksort(array(j+1:last))
            end if
        end subroutine quicksort


        recursive subroutine quicksort_index(array, sorted_index)
            ! """
            ! Sort double precision numbers into ascending numerical order.
            !
            ! input:
            !   array: double precision array to sort
            !
            ! output:
            !   sorted_index: sorted index of the array
            ! """
            implicit none

            double precision, intent(inout), dimension(:) :: array
            integer, intent(out), dimension(:) :: sorted_index

            integer :: &
                i

            do i = 1, size(array)
                sorted_index(i) = i
            end do

            call partition(array, sorted_index)

            contains
                recursive subroutine partition(array, sorted_index)
                    implicit none

                    integer, intent(inout), dimension(:) :: sorted_index
                    double precision, intent(inout), dimension(:) :: array

                    integer :: &
                        tmp_i

                    double precision :: &
                        x, &  ! pivot point
                        tmp   ! temporary value

                    integer :: &
                        first = 1, & ! index of the beginning of the array
                        last, &      ! index of the end of the array
                        insertion_size_threshold = 128

                    integer :: &
                        i, &
                        j

                    last = size(array, 1)

                    if(last < insertion_size_threshold) then  ! use insertion sort on small arrays
                        do i = 2, last
                            tmp = array(i)
                            tmp_i = sorted_index(i)

                            do j = i - 1, 1, -1
                                if(array(j) < tmp) then
                                    exit
                                else
                                    array(j + 1) = array(j)
                                    sorted_index(j + 1) = sorted_index(j)
                                end if
                            end do

                                array(j + 1) = tmp
                                sorted_index(j + 1) = tmp_i
                        end do
                    else
                        x = array((first + last) / 2)
                        i = first
                        j = last

                        do
                            do while (array(i) < x)
                                i = i + 1
                            end do

                            do while (x < array(j))
                                j = j - 1
                            end do

                            if (i >= j) then
                                exit
                            end if

                            tmp = array(i)
                            array(i) = array(j)
                            array(j) = tmp

                            tmp_i = sorted_index(i)
                            sorted_index(i) = sorted_index(j)
                            sorted_index(j) = tmp_i

                            i = i + 1
                            j = j - 1
                        end do

                        if (first < i - 1) call partition(array(first:i-1), sorted_index(first:i-1))
                        if (j + 1 < last)  call partition(array(j+1:last), sorted_index(j+1:last))
                    end if
                end subroutine partition
        end subroutine quicksort_index


        subroutine reallocate_1Ddouble(double_array, new_size)
            ! """
            ! Change the size of a 1D double precision array.
            ! :param double_array: array to change the size
            ! :param new_size: new size of the array
            ! :return double_array: resized array
            ! """
            implicit none

            integer, intent(in) :: new_size
            doubleprecision, intent(inout), dimension(:), allocatable :: double_array

            integer :: &
                size_array  ! size of the array
            doubleprecision, dimension(:), allocatable :: &
                tmp         ! temporary array

            size_array = size(double_array)

            call move_alloc(double_array, tmp)
            allocate(double_array(new_size))

            double_array(:) = 0d0

            size_array = min(new_size, size_array)

            double_array(:size_array) = tmp(:size_array)
        end subroutine reallocate_1Ddouble


        subroutine reallocate_1Dinteger(integer_array, new_size)
            ! """
            ! Change the size of a 1D integer array.
            ! :param integer_array: array to change the size
            ! :param new_size: new size of the array
            ! :return integer_array: resized array
            ! """
            implicit none

            integer, intent(in) :: new_size
            integer, intent(inout), dimension(:), allocatable :: integer_array

            integer :: &
                size_array  ! size of the array
            integer, dimension(:), allocatable :: &
                tmp         ! temporary array

            size_array = size(integer_array)

            call move_alloc(integer_array, tmp)
            allocate(integer_array(new_size))

            integer_array(:) = 0

            size_array = min(new_size, size_array)

            integer_array(:size_array) = tmp(:size_array)
        end subroutine reallocate_1Dinteger


        subroutine reallocate_2Ddouble(double_array, new_size_1, new_size_2)
            ! """
            ! Change the size of a 2D double precision array.
            ! :param double_array: array to change the size
            ! :param new_size_1: new size of dimension 1 of the array
            ! :param new_size_2: new size of dimension 2 the array
            ! :return double_array: resized array
            ! """
            implicit none

            integer, intent(in) :: new_size_1, new_size_2
            doubleprecision, intent(inout), dimension(:, :), allocatable :: double_array

            integer :: &
                shape_array(2)  ! shape of the array

            doubleprecision, dimension(:, :), allocatable :: &
                tmp             ! temporary array

            shape_array = shape(double_array)

            call move_alloc(double_array, tmp)
            allocate(double_array(new_size_1, new_size_2))

            double_array(:, :) = 0d0

            shape_array(1) = min(new_size_1, shape_array(1))
            shape_array(2) = min(new_size_2, shape_array(2))

            double_array(:shape_array(1), :shape_array(2)) = tmp(:shape_array(1), :shape_array(2))
        end subroutine reallocate_2Ddouble


        subroutine reallocate_2Dinteger(integer_array, new_size_1, new_size_2)
            ! """
            ! Change the size of a 2D integer array.
            ! :param integer_array: array to change the size
            ! :param new_size_1: new size of dimension 1 of the array
            ! :param new_size_2: new size of dimension 2 the array
            ! :return integer_array: resized array
            ! """
            implicit none

            integer, intent(in) :: new_size_1, new_size_2
            integer, intent(inout), dimension(:, :), allocatable :: integer_array

            integer :: &
                shape_array(2)  ! shape of the array

            integer, dimension(:, :), allocatable :: &
                tmp             ! temporary array

            shape_array = shape(integer_array)

            call move_alloc(integer_array, tmp)
            allocate(integer_array(new_size_1, new_size_2))

            integer_array(:, :) = 0

            shape_array(1) = min(new_size_1, shape_array(1))
            shape_array(2) = min(new_size_2, shape_array(2))

            integer_array(:shape_array(1), :shape_array(2)) = tmp(:shape_array(1), :shape_array(2))
        end subroutine reallocate_2Dinteger


        subroutine reallocate_3Ddouble(double_array, new_size_1, new_size_2, new_size_3)
            ! """
            ! Change the size of a 2D double precision array.
            ! :param double_array: array to change the size
            ! :param new_size_1: new size of dimension 1 of the array
            ! :param new_size_2: new size of dimension 2 of the array
            ! :param new_size_3: new size of dimension 3 of the array
            ! :return double_array: resized array
            ! """
            implicit none

            integer, intent(in) :: new_size_1, new_size_2, new_size_3
            doubleprecision, intent(inout), dimension(:, :, :), allocatable :: double_array

            integer :: &
                shape_array(3)  ! shape of the array

            doubleprecision, dimension(:, :, :), allocatable :: &
                tmp             ! temporary array

            shape_array = shape(double_array)

            call move_alloc(double_array, tmp)
            allocate(double_array(new_size_1, new_size_2, new_size_3))

            double_array(:, :, :) = 0d0

            shape_array(1) = min(new_size_1, shape_array(1))
            shape_array(2) = min(new_size_2, shape_array(2))
            shape_array(3) = min(new_size_3, shape_array(3))

            double_array(:shape_array(1), :shape_array(2), :shape_array(3)) = &
                tmp(:shape_array(1), :shape_array(2), :shape_array(3))
        end subroutine reallocate_3Ddouble


        subroutine ifft(x)
            ! """
            ! Calculate the Cooley-Tukey IFFT functions.
            !
            ! input:
            !   x: double complex vector of size 2^n
            ! """
            implicit none

            complex(kind=dp), dimension(:), intent(inout)  :: x

            call ifft_r(x)

            x = x / real(size(x), dp)

            contains
                recursive subroutine ifft_r(x)
                    implicit none

                    complex(kind=dp), dimension(:), intent(inout) :: x

                    integer ::&
                        i, &  ! index
                        n     ! size of array x

                    complex(kind=dp) ::&
                        t  ! intermediate value

                    complex(kind=dp), dimension(:), allocatable ::&
                        even, &  ! even-number indexed values of x
                        odd      ! odd-number indexed values of x

                    n=size(x)

                    if(n <= 1) return

                    allocate(odd((n + 1) / 2))
                    allocate(even(n / 2))

                    ! divide
                    odd(:) = x(1:n:2)
                    even(:) = x(2:n:2)

                    ! conquer
                    call ifft_r(odd)
                    call ifft_r(even)

                    ! combine
                    do i = 1, n/2
                       t = exp(cmplx(0.0_dp, 2.0_dp * pi * real(i - 1, dp) / real(n, dp), kind=dp)) * even(i)
                       x(i) = odd(i) + t
                       x(i+n/2) = odd(i) - t
                    end do

                    deallocate(odd)
                    deallocate(even)

                end subroutine ifft_r

        end subroutine ifft


        subroutine init_random_seed()
            ! """
            ! Initialize the random seed with a varying seed in order to ensure a different random number sequence for
            ! each invocation of the program.
            !
            ! notes:
            !   Using the random number generator file takes far too long to be useful.
            !   Source: https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
            ! """
            use iso_fortran_env, only: int64

            implicit none

            integer :: &
                i, &      ! index
                n, &      ! size of random seed
                dt(8), &  ! date and time
                pid       ! pid

            integer, allocatable :: &
                seed(:)  ! random generato seed

            integer(int64) :: &
                t  ! (s) time

            call random_seed(size = n)

            allocate(seed(n))

            ! XOR:ing the current time and pid. The PID is useful in case one launches multiple instances of the same
            ! program in parallel.
            call system_clock(t)

            if (t == 0) then
                call date_and_time(values=dt)

                t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
                   + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
                   + dt(3) * 24_int64 * 60 * 60 * 1000 &
                   + dt(5) * 60 * 60 * 1000 &
                   + dt(6) * 60 * 1000 + dt(7) * 1000 &
                   + dt(8)
            end if

            pid = getpid()
            t = ieor(t, int(pid, kind(t)))

            do i = 1, n
                seed(i) = lcg(t)
            end do

            call random_seed(put=seed)

            contains
                function lcg(s)
                    integer :: lcg
                    integer(int64) :: s

                    if (s == 0) then
                        s = 104729
                    else
                        s = mod(s, 4294967296_int64)
                    end if
                        s = mod(s * 279470273_int64, 4294967291_int64)
                        lcg = int(mod(s, int(huge(0), int64)), kind(0))
                end function lcg
        end subroutine init_random_seed


        subroutine matinv(a, n)
            ! """
            ! Inverse a 2D matrix of dimension (n,n)
            ! The original matrix a(n,n) will be destroyed during the calculation
            ! Method: Based on Doolittle LU factorization for Ax=b
            ! Source: http://ww2.odu.edu/~agodunov/computing/programs/book2/Ch06/Inverse.f90
            ! by Alex G. December 2009
            !
            ! inputs:
            !   a(n,n): array of coefficients for matrix A
            !   n: dimension
            !
            ! output:
            !   c(n,n): inverse matrix of A
            ! """
            implicit none

            integer, intent(in) :: n
            double precision, intent(inout) :: a(n,n)

            double precision :: &
                l(n,n), &  ! lower triangular matrix
                U(n,n), &  ! upper triangular matrix
                b(n), &    ! intermediate value
                c(n,n), &  ! temporary matrix
                d(n), &    ! intermediate value
                x(n)       ! intermediate value

            double precision :: &
                coeff  ! intermediate value

            integer :: &
                i, &  ! index
                j, &  ! index
                k     ! index

            ! Step 0: initialization for matrices l and U and b
            l(:, :) = 0d0
            U(:, :) = 0d0
            b = 0d0
            c(:, :) = 0d0

            ! Step 1: forward elimination
            do k = 1, n - 1
                do i = k + 1, n
                    coeff = a(i, k) / a(k, k)

                    l(i, k) = coeff

                    do j = k + 1, n
                        a(i, j) = a(i, j) - coeff * a(k, j)
                    end do
                end do
            end do

            ! Step 2: prepare l and U matrices
            ! l matrix is a matrix of the elimination coefficient + the diagonal elements are 1.0
            do i = 1, n
                l(i, i) = 1.0
            end do

            ! U matrix is the upper triangular part of A
            do j = 1, n
                do i = 1, j
                    U(i, j) = a(i, j)
                end do
            end do

            ! Step 3: compute columns of the inverse matrix C
            do k = 1, n
                b(k) = 1.0
                d(1) = b(1)

                ! Step 3a: Solve Ld=b using the forward substitution
                do i = 2, n
                    d(i) = b(i)

                    do j=1, i - 1
                        d(i) = d(i) - l(i, j) * d(j)
                    end do
                end do

                ! Step 3b: Solve Ux=d using the back substitution
                x(n) = d(n) / U(n,n)

                do i = n - 1, 1, -1
                    x(i) = d(i)

                    do j = n, i + 1, -1
                        x(i) = x(i) - U(i, j) * x(j)
                    end do

                    x(i) = x(i) / u(i, i)
                end do

                ! Step 3c: fill the solutions x(n) into column k of C
                do i = 1, n
                    c(i, k) = x(i)
                end do

                b(k) = 0d0
            end do

            a(:, :) = c(:, :)
        end subroutine matinv
end module math