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