From a2869260c41e920d903a961ae6071d1c6103a2c3 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 13 Jul 2026 14:55:00 +0200 Subject: [PATCH] fix(KIM): verify off-diagonal FP kernel --- KIM/src/asymptotics/flr2_fourier_kernel.f90 | 318 ++++++++++++------ KIM/tests/CMakeLists.txt | 4 + KIM/tests/test_flr2_fourier_kernel.f90 | 220 +++++++++--- verification/FORMULA_INDEX.md | 18 +- .../04_charge_potential_kernels.wl | 23 +- verification/oracles/rho_phi_kernels.dat | 28 +- 6 files changed, 450 insertions(+), 161 deletions(-) diff --git a/KIM/src/asymptotics/flr2_fourier_kernel.f90 b/KIM/src/asymptotics/flr2_fourier_kernel.f90 index e5b63e8f..1a430129 100644 --- a/KIM/src/asymptotics/flr2_fourier_kernel.f90 +++ b/KIM/src/asymptotics/flr2_fourier_kernel.f90 @@ -1,7 +1,8 @@ module flr2_fourier_kernel_m - ! Fused Fokker-Planck off-diagonal Fourier-space kernel integrand - ! G(k_r, k'_r, r_g) for the forced-periodicity solver. Collapses to the - ! diagonal source-of-truth calc_hatK_Phi_in_Fourier at k_r = k'_r. + ! Oracle-verified Fokker-Planck off-diagonal Fourier-space kernel integrand + ! G(k_r, k'_r, r_g) for the forced-periodicity solver. Production uses the + ! full finite-radius k_perp model; the historical radial-only diagonal is + ! retained as an explicitly named approximation for compatibility checks. use KIM_kinds_m, only: dp use species_m, only: plasma_t implicit none @@ -9,132 +10,258 @@ module flr2_fourier_kernel_m public :: hatG_rho_phi, hatG_rho_B public :: hatG_rho_phi_diag_sp, hatG_rho_B_diag_sp public :: core_rho_phi_sp, core_rho_B_sp - public :: scaled_bessel_pair - - ! FLR-parameter convention for the two-wavenumber kernel. When .false. - ! (default) the perpendicular wavenumber k_s is dropped from b_+ / b_x, so - ! the kernel depends only on the radial wavenumbers k_r, k'_r and collapses - ! exactly to the diagonal source-of-truth calc_hatK_Phi_in_Fourier. The - ! .true. path reintroduces k_s^2 via kperp = sqrt(k_s^2 + k_r^2) (see - ! hatG_rho_phi) — the canonical form of the recovered kernel_mod.f90. - logical, save, public :: kern_include_ks2 = .false. + public :: scaled_bessel_pair, flr_arg_pair + public :: fp_rho_phi_harmonic, fp_rho_B_harmonic, fp_rho_phi_debye + + integer, parameter, public :: FP_KPERP_FULL = 1 + integer, parameter, public :: FP_KR_ONLY_APPROX = 2 + integer, save, public :: fp_ks_model = FP_KPERP_FULL contains + pure subroutine flr_arg_pair(ks, kr, krp, rho_L, model, bplus, bcross, status) + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite + real(dp), intent(in) :: ks, kr, krp, rho_L + integer, intent(in) :: model + real(dp), intent(out) :: bplus, bcross + integer, intent(out) :: status + real(dp) :: rho_L2 + + status = 0 + bplus = 0.0_dp + bcross = 0.0_dp + if (.not. ieee_is_finite(kr) .or. .not. ieee_is_finite(krp) .or. & + .not. ieee_is_finite(rho_L) .or. rho_L < 0.0_dp) then + status = 1 + return + end if + + rho_L2 = rho_L*rho_L + select case (model) + case (FP_KPERP_FULL) + if (.not. ieee_is_finite(ks)) then + status = 2 + return + end if + bplus = rho_L2*(2.0_dp*ks*ks + kr*kr + krp*krp)/2.0_dp + bcross = rho_L2*sqrt(ks*ks + kr*kr)*sqrt(ks*ks + krp*krp) + case (FP_KR_ONLY_APPROX) + bplus = rho_L2*(kr*kr + krp*krp)/2.0_dp + bcross = rho_L2*abs(kr*krp) + case default + status = 3 + end select + end subroutine flr_arg_pair + subroutine flr_arg_pair_sp(plasma_in, sp, kr, krp, j, bplus, bcross) - ! Shared FLR argument pair (b_+, b_x) for the two-wavenumber kernels. - ! Single home for the kern_include_ks2 branch used by both hatG_rho_phi - ! and hatG_rho_B. When .false. (default) the perpendicular wavenumber - ! k_s is dropped, so the pair depends only on the radial wavenumbers - ! k_r, k'_r and collapses exactly to the diagonal source-of-truth. The - ! .true. path reintroduces k_s^2 via kperp = sqrt(k_s^2 + k_r^2), giving - ! b_+ = rho_L^2/2 (kperp^2 + kperpp^2) and b_x = rho_L^2 kperp kperpp. type(plasma_t), intent(in) :: plasma_in integer, intent(in) :: sp, j real(dp), intent(in) :: kr, krp real(dp), intent(out) :: bplus, bcross - real(dp) :: rho_L2, ks + integer :: status - rho_L2 = plasma_in%spec(sp)%rho_L(j)**2.0d0 - if (kern_include_ks2) then - ks = plasma_in%ks(j) - bplus = rho_L2 * (2.0d0 * ks**2 + kr**2 + krp**2) / 2.0d0 - bcross = rho_L2 * sqrt(ks**2 + kr**2) * sqrt(ks**2 + krp**2) - else - bplus = rho_L2 * (kr**2 + krp**2) / 2.0d0 - bcross = rho_L2 * abs(kr * krp) - end if + call flr_arg_pair(plasma_in%ks(j), kr, krp, plasma_in%spec(sp)%rho_L(j), & + fp_ks_model, bplus, bcross, status) + if (status /= 0) error stop 'invalid finite-radius FLR arguments' end subroutine flr_arg_pair_sp - subroutine scaled_bessel_pair(bplus, bcross, sI0, sIm1) - ! Overflow-safe scaled modified Bessel products: - ! sI0 = exp(-bplus) * I_0(bcross) - ! sIm1 = exp(-bplus) * I_{-1}(bcross) - ! For large bcross the unscaled I_n(bcross) would overflow before the - ! exp(-bplus) damping (b_+ >= b_x by construction), so use the large- - ! argument asymptotics (lifted from the recovered kernel_mod.f90). + subroutine scaled_bessel_harmonic(bplus, bcross, order, value, status) use fortnum_special, only: bessel_in real(dp), intent(in) :: bplus, bcross - real(dp), intent(out) :: sI0, sIm1 - real(dp), parameter :: asymptotic_threshold = 500.0d0 ! safely below the ~710 I_0 overflow - real(dp), parameter :: twopi = 6.283185307179586d0 + integer, intent(in) :: order + real(dp), intent(out) :: value + integer, intent(out) :: status + real(dp), parameter :: asymptotic_threshold = 500.0_dp + real(dp), parameter :: twopi = 6.283185307179586_dp + real(dp) :: mu, term, series + integer :: k, n - if (bcross > asymptotic_threshold) then - sI0 = exp(bcross - bplus) / sqrt(twopi * bcross) - sIm1 = exp(-bplus + asinh(-1.0d0/bcross) + bcross*sqrt(1.0d0 + 1.0d0/bcross**2)) & - / sqrt(twopi * bcross * sqrt(1.0d0 + 1.0d0/bcross**2)) + status = 0 + value = 0.0_dp + if (bcross < 0.0_dp .or. bplus < 0.0_dp .or. & + bplus + 128.0_dp*epsilon(1.0_dp)*max(1.0_dp, bplus) < bcross) then + status = 1 + return + end if + n = abs(order) + if (bcross == 0.0_dp) then + if (n == 0) value = exp(-bplus) + else if (bcross > asymptotic_threshold) then + mu = 4.0_dp*real(n*n, dp) + series = 1.0_dp + term = 1.0_dp + do k = 1, 5 + term = term * (-(mu - real((2*k - 1)**2, dp))) / & + (real(k, dp)*8.0_dp*bcross) + series = series + term + end do + value = exp(bcross - bplus)*series/sqrt(twopi*bcross) else - sI0 = exp(-bplus) * bessel_in(0, bcross) - sIm1 = exp(-bplus) * bessel_in(-1, bcross) + value = exp(-bplus)*bessel_in(n, bcross) end if + end subroutine scaled_bessel_harmonic + + subroutine scaled_bessel_pair(bplus, bcross, sI0, sIm1) + real(dp), intent(in) :: bplus, bcross + real(dp), intent(out) :: sI0, sIm1 + integer :: status + + call scaled_bessel_harmonic(bplus, bcross, 0, sI0, status) + if (status /= 0) error stop 'invalid scaled Bessel arguments' + call scaled_bessel_harmonic(bplus, bcross, -1, sIm1, status) + if (status /= 0) error stop 'invalid scaled Bessel arguments' end subroutine scaled_bessel_pair + subroutine fp_rho_phi_harmonic(lambda_D, vT, omega_c, nu, ks, kr, krp, rg, & + mphi, A1, A2, I00, I20, model, value, status) + use constants_m, only: pi, com_unit + real(dp), intent(in) :: lambda_D, vT, omega_c, nu, ks, kr, krp, rg + integer, intent(in) :: mphi, model + real(dp), intent(in) :: A1, A2 + complex(dp), intent(in) :: I00, I20 + complex(dp), intent(out) :: value + integer, intent(out) :: status + real(dp) :: bplus, bcross, sIm, sImm1, energy_moment, alpha, alphap + complex(dp) :: phase + + value = (0.0_dp, 0.0_dp) + if (lambda_D <= 0.0_dp .or. omega_c == 0.0_dp .or. nu <= 0.0_dp) then + status = 4 + return + end if + call flr_arg_pair(ks, kr, krp, abs(vT/omega_c), model, bplus, bcross, status) + if (status /= 0) return + call scaled_bessel_harmonic(bplus, bcross, mphi, sIm, status) + if (status /= 0) return + call scaled_bessel_harmonic(bplus, bcross, mphi - 1, sImm1, status) + if (status /= 0) return + + energy_moment = (1.0_dp - bplus - real(mphi, dp))*sIm + bcross*sImm1 + alpha = atan2(kr, ks) + alphap = atan2(krp, ks) + phase = exp(com_unit*((kr - krp)*rg + real(mphi, dp)*(alphap - alpha))) + value = phase/(4.0_dp*pi*lambda_D**2) * com_unit*vT**2*ks/(omega_c*nu) * & + (I00*(A1*sIm + A2*energy_moment) + 0.5_dp*A2*I20*sIm) + end subroutine fp_rho_phi_harmonic + + subroutine fp_rho_B_harmonic(lambda_D, vT, omega_c, nu, ks, kr, krp, rg, & + mphi, A1, A2, I01, I21, model, value, status) + use constants_m, only: pi, com_unit, sol + real(dp), intent(in) :: lambda_D, vT, omega_c, nu, ks, kr, krp, rg + integer, intent(in) :: mphi, model + real(dp), intent(in) :: A1, A2 + complex(dp), intent(in) :: I01, I21 + complex(dp), intent(out) :: value + integer, intent(out) :: status + real(dp) :: bplus, bcross, sIm, sImm1, energy_moment, alpha, alphap + complex(dp) :: phase + + value = (0.0_dp, 0.0_dp) + if (lambda_D <= 0.0_dp .or. omega_c == 0.0_dp .or. nu <= 0.0_dp) then + status = 4 + return + end if + call flr_arg_pair(ks, kr, krp, abs(vT/omega_c), model, bplus, bcross, status) + if (status /= 0) return + call scaled_bessel_harmonic(bplus, bcross, mphi, sIm, status) + if (status /= 0) return + call scaled_bessel_harmonic(bplus, bcross, mphi - 1, sImm1, status) + if (status /= 0) return + + energy_moment = (1.0_dp - bplus - real(mphi, dp))*sIm + bcross*sImm1 + alpha = atan2(kr, ks) + alphap = atan2(krp, ks) + phase = exp(com_unit*((kr - krp)*rg + real(mphi, dp)*(alphap - alpha))) + value = -phase/(4.0_dp*pi*lambda_D**2) * vT**3/(omega_c*nu*sol) * & + (I01*(A1*sIm + A2*energy_moment) + 0.5_dp*A2*I21*sIm) + end subroutine fp_rho_B_harmonic + + pure subroutine fp_rho_phi_debye(lambda_D, kr, krp, rg, include_debye, value) + use constants_m, only: pi, com_unit + real(dp), intent(in) :: lambda_D, kr, krp, rg + logical, intent(in) :: include_debye + complex(dp), intent(out) :: value + + value = (0.0_dp, 0.0_dp) + if (include_debye) then + value = -exp(com_unit*(kr - krp)*rg)/(4.0_dp*pi*lambda_D**2) + end if + end subroutine fp_rho_phi_debye + complex(dp) function hatG_rho_phi(plasma_in, kr, krp, j) result(G) - ! Fused two-wavenumber rho-Phi kernel Ghat^{rho Phi}(k_r, k'_r, r_g): - ! sum the per-species core over the non-turned-off species, then apply - ! the Fourier phase exp(i (k_r - k'_r) r_g) and the /(4*pi) source-of- - ! truth normalization. On the diagonal k_r = k'_r the FLR pair collapses - ! to a single b and the phase is unity, so this equals the diagonal - ! integrand summed over species and divided by 4*pi (collapse by - ! construction, since the diagonal shares this same core). + ! Full two-wavenumber rho-Phi kernel Ghat^{rho Phi}(k_r,k'_r,r_g). + ! Each retained harmonic carries its gyro-angle/radial Fourier phase, + ! finite-radius moments, and 1/(4*pi) normalization. The Debye term is + ! independent of m_phi and is therefore included exactly once/species. use grid_m, only: rg_grid - use constants_m, only: pi, com_unit use config_m, only: turn_off_ions, turn_off_electrons + use config_m, only: artificial_debye_case + use setup_m, only: mphi_max type(plasma_t), intent(in) :: plasma_in real(dp), intent(in) :: kr, krp integer, intent(in) :: j - integer :: sp - real(dp) :: bplus, bcross - complex(dp) :: acc + integer :: sp, mphi, status + complex(dp) :: term - acc = (0.0_dp, 0.0_dp) + G = (0.0_dp, 0.0_dp) do sp = 0, plasma_in%n_species - 1 if (turn_off_ions .and. sp >= 1) cycle if (turn_off_electrons .and. sp == 0) cycle - call flr_arg_pair_sp(plasma_in, sp, kr, krp, j, bplus, bcross) - acc = acc + core_rho_phi_sp(plasma_in, sp, bplus, bcross, j) - end do + call fp_rho_phi_debye(plasma_in%spec(sp)%lambda_D(j), kr, krp, & + rg_grid%xb(j), artificial_debye_case <= 1, term) + G = G + term + if (artificial_debye_case == 1) cycle - G = exp(com_unit * (kr - krp) * rg_grid%xb(j)) / (4.0d0 * pi) * acc + do mphi = -mphi_max, mphi_max + call fp_rho_phi_harmonic(plasma_in%spec(sp)%lambda_D(j), & + plasma_in%spec(sp)%vT(j), plasma_in%spec(sp)%omega_c(j), & + plasma_in%spec(sp)%nu(j), plasma_in%ks(j), kr, krp, rg_grid%xb(j), & + mphi, plasma_in%spec(sp)%A1(j), plasma_in%spec(sp)%A2(j), & + plasma_in%spec(sp)%I00(j, mphi), plasma_in%spec(sp)%I20(j, mphi), & + fp_ks_model, term, status) + if (status /= 0) error stop 'invalid rho-Phi harmonic arguments' + G = G + term + end do + end do end function hatG_rho_phi complex(dp) function hatG_rho_B(plasma_in, kr, krp, j) result(G) - ! Fused two-wavenumber rho-B kernel Ghat^{rho B}(k_r, k'_r, r_g): sum the - ! signed per-species core over the non-turned-off species, then apply the - ! Fourier phase exp(i (k_r - k'_r) r_g) and the /(4*pi) source-of-truth - ! normalization. On the diagonal k_r = k'_r the FLR pair collapses to a - ! single b and the phase is unity, so this equals the diagonal integrand - ! summed over species and divided by 4*pi (collapse by construction). - ! rho-B has no Debye term (handled inside core_rho_B_sp). + ! Full two-wavenumber rho-B kernel Ghat^{rho B}(k_r,k'_r,r_g). + ! Each retained harmonic carries its gyro-angle/radial Fourier phase, + ! signed current moment, finite-radius factors, and 1/(4*pi) + ! normalization. rho-B has no Debye contribution. use grid_m, only: rg_grid - use constants_m, only: pi, com_unit use config_m, only: turn_off_ions, turn_off_electrons + use config_m, only: artificial_debye_case + use setup_m, only: mphi_max type(plasma_t), intent(in) :: plasma_in real(dp), intent(in) :: kr, krp integer, intent(in) :: j - integer :: sp - real(dp) :: bplus, bcross - complex(dp) :: acc + integer :: sp, mphi, status + complex(dp) :: term - acc = (0.0_dp, 0.0_dp) + G = (0.0_dp, 0.0_dp) + if (artificial_debye_case == 1) return do sp = 0, plasma_in%n_species - 1 if (turn_off_ions .and. sp >= 1) cycle if (turn_off_electrons .and. sp == 0) cycle - call flr_arg_pair_sp(plasma_in, sp, kr, krp, j, bplus, bcross) - acc = acc + core_rho_B_sp(plasma_in, sp, bplus, bcross, j) + do mphi = -mphi_max, mphi_max + call fp_rho_B_harmonic(plasma_in%spec(sp)%lambda_D(j), & + plasma_in%spec(sp)%vT(j), plasma_in%spec(sp)%omega_c(j), & + plasma_in%spec(sp)%nu(j), plasma_in%ks(j), kr, krp, rg_grid%xb(j), & + mphi, plasma_in%spec(sp)%A1(j), plasma_in%spec(sp)%A2(j), & + plasma_in%spec(sp)%I01(j, mphi), plasma_in%spec(sp)%I21(j, mphi), & + fp_ks_model, term, status) + if (status /= 0) error stop 'invalid rho-B harmonic arguments' + G = G + term + end do end do - - G = exp(com_unit * (kr - krp) * rg_grid%xb(j)) / (4.0d0 * pi) * acc end function hatG_rho_B complex(dp) function core_rho_phi_sp(plasma_in, sp, bplus, bcross, j) result(G) - ! Per-species rho-Phi core (Debye + FP) parameterized by the FLR - ! argument pair (b_+, b_x). This is the single home for the rho-Phi - ! integrand math: the diagonal integrand calls it with bplus=bcross=b, - ! and the fused two-wavenumber kernel calls it with the off-diagonal - ! pair. b_+ scales exp(-b) and the (1 - b) polynomial factor; b_x is the - ! argument (and the leading multiplier) of the modified Bessel terms. - ! No phase, NO /(4*pi); species selection (turn_off_*) stays in callers. + ! Historical m_phi=0 per-species rho-Phi core used by diagonal + ! compatibility tests and diagnostics. The production off-diagonal + ! path uses fp_rho_phi_harmonic instead. No phase or 1/(4*pi). use config_m, only: artificial_debye_case use constants_m, only: com_unit type(plasma_t), intent(in) :: plasma_in @@ -165,14 +292,9 @@ complex(dp) function core_rho_phi_sp(plasma_in, sp, bplus, bcross, j) result(G) end function core_rho_phi_sp complex(dp) function core_rho_B_sp(plasma_in, sp, bplus, bcross, j) result(G) - ! Per-species rho-B core (SIGNED FP term) parameterized by the FLR - ! argument pair (b_+, b_x). Single home for the rho-B integrand math: - ! the diagonal integrand calls it with bplus=bcross=b, and the fused - ! two-wavenumber kernel calls it with the off-diagonal pair. b_+ scales - ! exp(-b) and the (1 - b) polynomial factor; b_x is the argument (and the - ! leading multiplier) of the modified Bessel terms. rho-B carries the - ! leading minus, so callers accumulate with '+'. NO Debye term, no phase, - ! NO /(4*pi); species selection (turn_off_*) stays in callers. + ! Historical m_phi=0 per-species signed rho-B core used by diagonal + ! compatibility tests and diagnostics. The production off-diagonal + ! path uses fp_rho_B_harmonic instead. No phase or 1/(4*pi). use config_m, only: artificial_debye_case use constants_m, only: sol type(plasma_t), intent(in) :: plasma_in @@ -206,10 +328,10 @@ complex(dp) function hatG_rho_phi_diag_sp(plasma_in, sp, kr, j) result(G) type(plasma_t), intent(in) :: plasma_in integer, intent(in) :: sp, j real(dp), intent(in) :: kr - real(dp) :: b + real(dp) :: b, bcross - b = (kr**2.0d0) * plasma_in%spec(sp)%rho_L(j)**2.0d0 - G = core_rho_phi_sp(plasma_in, sp, b, b, j) + call flr_arg_pair_sp(plasma_in, sp, kr, kr, j, b, bcross) + G = core_rho_phi_sp(plasma_in, sp, b, bcross, j) end function hatG_rho_phi_diag_sp complex(dp) function hatG_rho_B_diag_sp(plasma_in, sp, kr, j) result(G) @@ -221,9 +343,9 @@ complex(dp) function hatG_rho_B_diag_sp(plasma_in, sp, kr, j) result(G) type(plasma_t), intent(in) :: plasma_in integer, intent(in) :: sp, j real(dp), intent(in) :: kr - real(dp) :: b + real(dp) :: b, bcross - b = (kr**2.0d0) * plasma_in%spec(sp)%rho_L(j)**2.0d0 - G = core_rho_B_sp(plasma_in, sp, b, b, j) + call flr_arg_pair_sp(plasma_in, sp, kr, kr, j, b, bcross) + G = core_rho_B_sp(plasma_in, sp, b, bcross, j) end function hatG_rho_B_diag_sp end module flr2_fourier_kernel_m diff --git a/KIM/tests/CMakeLists.txt b/KIM/tests/CMakeLists.txt index 83dac984..cbb3d7de 100644 --- a/KIM/tests/CMakeLists.txt +++ b/KIM/tests/CMakeLists.txt @@ -68,8 +68,12 @@ set_target_properties(test_flr2_fourier_kernel PROPERTIES OUTPUT_NAME test_flr2_fourier_kernel.x RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") target_link_libraries(test_flr2_fourier_kernel KIM_lib kilca_lib lapack cerf ddeabm slatec) +configure_file(${CMAKE_SOURCE_DIR}/verification/oracles/rho_phi_kernels.dat + ${CMAKE_BINARY_DIR}/tests/rho_phi_kernels.dat COPYONLY) add_test(NAME test_flr2_fourier_kernel COMMAND ${CMAKE_BINARY_DIR}/tests/test_flr2_fourier_kernel.x) +set_tests_properties(test_flr2_fourier_kernel PROPERTIES + WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/tests/) add_executable(test_ampere_matrices ${CMAKE_SOURCE_DIR}/KIM/tests/test_ampere_matrices.f90) set_target_properties(test_ampere_matrices PROPERTIES diff --git a/KIM/tests/test_flr2_fourier_kernel.f90 b/KIM/tests/test_flr2_fourier_kernel.f90 index bb9075ed..c6d43f00 100644 --- a/KIM/tests/test_flr2_fourier_kernel.f90 +++ b/KIM/tests/test_flr2_fourier_kernel.f90 @@ -20,9 +20,11 @@ program test_flr2_fourier_kernel use KIM_kinds_m, only: dp use flr2_fourier_kernel_m, only: hatG_rho_phi, hatG_rho_B use flr2_fourier_kernel_m, only: hatG_rho_phi_diag_sp, hatG_rho_B_diag_sp + use flr2_fourier_kernel_m, only: fp_ks_model, FP_KPERP_FULL, FP_KR_ONLY_APPROX implicit none + call test_independent_rho_phi_oracle() call test_populated_plasma_and_kernel_stub() call test_diagonal_matches_inline() call test_collapse_rho_phi() @@ -40,6 +42,147 @@ program test_flr2_fourier_kernel contains + subroutine test_independent_rho_phi_oracle() + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use flr2_fourier_kernel_m, only: FP_KPERP_FULL, FP_KR_ONLY_APPROX + use flr2_fourier_kernel_m, only: flr_arg_pair, fp_rho_phi_harmonic, fp_rho_B_harmonic + use flr2_fourier_kernel_m, only: fp_rho_phi_debye + + integer :: unit, ios, case_id, charge_sign, mphi, include_debye + integer :: status, nrows + real(dp) :: lambda_D, vT, omega_c, nu, ks, kr, krp, rg, A1, A2 + real(dp) :: kpar, omega_E, Vpar, omega, x1, x2, rho_L + real(dp) :: bplus_ref, bcross_ref, phase_re, phase_im + real(dp) :: I00_re, I00_im, I20_re, I20_im, I01_re, I01_im, I21_re, I21_im + real(dp) :: F0_ref, F2_ref + real(dp) :: dynamic_re, dynamic_im, debye_re, debye_im + real(dp) :: total_re, total_im, bplus, bcross, scaled_error + real(dp) :: rho_B_re, rho_B_im + real(dp) :: nan_value, max_scaled_error + complex(dp) :: I00, I20, I01, I21, actual_dynamic, actual_debye, actual_rho_B, expected + complex(dp) :: representative(-2:2), representative_total + character(len=2048) :: line + logical :: seen_charge(-1:1), seen_large, seen_zero_flr, seen_debye + + nrows = 0 + max_scaled_error = 0.0_dp + seen_charge = .false. + seen_large = .false. + seen_zero_flr = .false. + seen_debye = .false. + representative = (0.0_dp, 0.0_dp) + open(newunit=unit, file='rho_phi_kernels.dat', status='old', action='read', iostat=ios) + if (ios /= 0) error stop 'cannot open rho_phi_kernels.dat' + + do + read(unit, '(A)', iostat=ios) line + if (ios < 0) exit + if (ios > 0) error stop 'cannot read rho-Phi oracle line' + if (len_trim(line) == 0 .or. line(1:1) == '#') cycle + read(line, *, iostat=ios) case_id, charge_sign, lambda_D, vT, omega_c, nu, & + ks, kr, krp, rg, mphi, A1, A2, kpar, omega_E, Vpar, omega, & + include_debye, x1, x2, rho_L, bplus_ref, bcross_ref, phase_re, & + phase_im, I00_re, I00_im, I20_re, I20_im, I01_re, I01_im, & + I21_re, I21_im, F0_ref, F2_ref, & + dynamic_re, dynamic_im, debye_re, debye_im, total_re, total_im, & + rho_B_re, rho_B_im + if (ios /= 0) error stop 'malformed rho-Phi oracle row' + + call flr_arg_pair(ks, kr, krp, rho_L, FP_KPERP_FULL, bplus, bcross, status) + if (status /= 0) error stop 'full finite-radius FLR arguments rejected oracle row' + call assert_oracle_real('bplus', bplus, bplus_ref, case_id, max_scaled_error) + call assert_oracle_real('bcross', bcross, bcross_ref, case_id, max_scaled_error) + + I00 = cmplx(I00_re, I00_im, kind=dp) + I20 = cmplx(I20_re, I20_im, kind=dp) + call fp_rho_phi_harmonic(lambda_D, vT, omega_c, nu, ks, kr, krp, rg, & + mphi, A1, A2, I00, I20, FP_KPERP_FULL, actual_dynamic, status) + if (status /= 0) error stop 'rho-Phi harmonic rejected oracle row' + expected = cmplx(dynamic_re, dynamic_im, kind=dp) + call assert_oracle_complex('dynamic', actual_dynamic, expected, case_id, max_scaled_error) + + call fp_rho_phi_debye(lambda_D, kr, krp, rg, include_debye == 1 .and. mphi == 0, actual_debye) + expected = cmplx(debye_re, debye_im, kind=dp) + call assert_oracle_complex('Debye', actual_debye, expected, case_id, max_scaled_error) + expected = cmplx(total_re, total_im, kind=dp) + call assert_oracle_complex('total', actual_dynamic + actual_debye, expected, & + case_id, max_scaled_error) + + I01 = cmplx(I01_re, I01_im, kind=dp) + I21 = cmplx(I21_re, I21_im, kind=dp) + call fp_rho_B_harmonic(lambda_D, vT, omega_c, nu, ks, kr, krp, rg, & + mphi, A1, A2, I01, I21, FP_KPERP_FULL, actual_rho_B, status) + if (status /= 0) error stop 'rho-B harmonic rejected oracle row' + expected = cmplx(rho_B_re, rho_B_im, kind=dp) + call assert_oracle_complex('rho-B', actual_rho_B, expected, case_id, max_scaled_error) + + nrows = nrows + 1 + if (charge_sign >= -1 .and. charge_sign <= 1) seen_charge(charge_sign) = .true. + if (bcross_ref >= 500.0_dp) seen_large = .true. + if (rho_L <= 1.0e-8_dp) seen_zero_flr = .true. + if (include_debye == 1) seen_debye = .true. + if (case_id >= 9 .and. case_id <= 13) representative(mphi) = actual_dynamic + end do + close(unit) + + if (nrows /= 13) error stop 'rho-Phi oracle row coverage changed' + if (.not. seen_charge(-1) .or. .not. seen_charge(1)) error stop 'oracle lacks both charge signs' + if (.not. seen_large) error stop 'oracle lacks large-Bessel case' + if (.not. seen_zero_flr) error stop 'oracle lacks zero-FLR case' + if (.not. seen_debye) error stop 'oracle lacks Debye case' + + representative_total = sum(representative) + scaled_error = abs(representative_total - representative(0))/abs(representative_total) + if (scaled_error >= 5.0e-3_dp) error stop 'mphi=0 truncation exceeds verified bound' + scaled_error = (abs(representative(-2)) + abs(representative(2)))/abs(representative_total) + if (scaled_error >= 5.0e-5_dp) error stop '|mphi|=2 shell exceeds verified bound' + + nan_value = ieee_value(0.0_dp, ieee_quiet_nan) + call flr_arg_pair(nan_value, 1.0_dp, 2.0_dp, 0.5_dp, FP_KPERP_FULL, & + bplus, bcross, status) + if (status == 0) error stop 'full finite-radius model accepted undefined ks' + call flr_arg_pair(nan_value, 1.0_dp, 2.0_dp, 0.5_dp, FP_KR_ONLY_APPROX, & + bplus, bcross, status) + if (status /= 0) error stop 'named radial-only approximation used undefined ks' + + print *, 'PASS: ', nrows, ' independent finite-radius rho-Phi oracle rows' + print *, 'Maximum scaled kernel-oracle error: ', max_scaled_error + end subroutine test_independent_rho_phi_oracle + + subroutine assert_oracle_real(label, actual, expected, case_id, max_error) + character(len=*), intent(in) :: label + real(dp), intent(in) :: actual, expected + integer, intent(in) :: case_id + real(dp), intent(inout) :: max_error + real(dp) :: error + + error = abs(actual - expected)/max(1.0_dp, abs(expected)) + max_error = max(max_error, error) + if (error > 1.0e-12_dp) then + print *, 'FAIL: ', trim(label), ' oracle case ', case_id + print *, ' actual/expected/error = ', actual, expected, error + error stop + end if + end subroutine assert_oracle_real + + subroutine assert_oracle_complex(label, actual, expected, case_id, max_error) + character(len=*), intent(in) :: label + complex(dp), intent(in) :: actual, expected + integer, intent(in) :: case_id + real(dp), intent(inout) :: max_error + real(dp) :: error + + error = abs(actual - expected)/max(1.0_dp, abs(expected)) + max_error = max(max_error, error) + if (error > 1.0e-12_dp) then + print *, 'FAIL: ', trim(label), ' oracle case ', case_id + print *, ' actual = ', actual + print *, ' expected = ', expected + print *, ' scaled error = ', error + error stop + end if + end subroutine assert_oracle_complex + subroutine test_diagonal_matches_inline() ! Characterization test: the shared per-species diagonal integrand ! functions must reproduce, bit-for-bit within tolerance, the original @@ -80,6 +223,10 @@ subroutine test_diagonal_matches_inline() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() + ! This characterization block intentionally exercises the named + ! historical radial-only approximation copied below. + fp_ks_model = FP_KR_ONLY_APPROX + kr_arr = [0.1_dp, 1.0_dp, 5.0_dp] ! A few interior guiding-centre boundary-grid indices. jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] @@ -168,7 +315,7 @@ subroutine test_collapse_rho_phi() ! diagonal and the fused kernel share the same per-species core. use config_m, only: profiles_in_memory, nml_config_path use config_m, only: turn_off_ions, turn_off_electrons - use flr2_fourier_kernel_m, only: hatG_rho_phi, hatG_rho_phi_diag_sp, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_phi, hatG_rho_phi_diag_sp use constants_m, only: pi use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -198,7 +345,7 @@ subroutine test_collapse_rho_phi() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - kern_include_ks2 = .false. + fp_ks_model = FP_KR_ONLY_APPROX kr_arr = [0.1_dp, 1.0_dp, 5.0_dp] jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] @@ -218,7 +365,7 @@ subroutine test_collapse_rho_phi() gfused = hatG_rho_phi(plasma, kr, kr, j) - tol = 1.0e-12_dp * (1.0_dp + abs(dref)) + tol = 5.0e-12_dp * (1.0_dp + abs(dref)) if (abs(gfused - dref) >= tol) then print *, 'FAIL: rho-Phi collapse mismatch at kr=', kr, ' j=', j print *, ' diag/4pi = ', dref @@ -242,7 +389,7 @@ subroutine test_phase_rho_phi() ! test: this compares two independent complex-exp/Bessel evaluations, ! not a construction identity. use config_m, only: profiles_in_memory, nml_config_path - use flr2_fourier_kernel_m, only: hatG_rho_phi, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_phi use constants_m, only: com_unit use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -272,7 +419,7 @@ subroutine test_phase_rho_phi() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - kern_include_ks2 = .false. + fp_ks_model = FP_KR_ONLY_APPROX ci = com_unit kr_arr = [0.1_dp, 1.0_dp, 2.0_dp] @@ -307,8 +454,8 @@ subroutine test_phase_rho_phi() end subroutine subroutine test_collapse_rho_phi_ks2() - ! Structural collapse for the k_s^2-inclusive path (kern_include_ks2 = - ! .true.). This branch keeps k_s in the FLR arguments, so it does NOT + ! Structural collapse for the full finite-radius model. This model + ! keeps k_s in the FLR arguments, so it does NOT ! collapse to the k_s-dropped diagonal source-of-truth. Instead, on the ! diagonal k_r = k'_r the two FLR arguments both reduce to the full ! perpendicular argument bf = (k_s^2 + k_r^2) * rho_L^2, so the fused @@ -317,7 +464,7 @@ subroutine test_collapse_rho_phi_ks2() ! sqrt/2 algebra of the .true. branch without needing a k_s = 0 setup. use config_m, only: profiles_in_memory, nml_config_path use config_m, only: turn_off_ions, turn_off_electrons - use flr2_fourier_kernel_m, only: hatG_rho_phi, core_rho_phi_sp, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_phi, core_rho_phi_sp use constants_m, only: pi use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -347,7 +494,7 @@ subroutine test_collapse_rho_phi_ks2() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - kern_include_ks2 = .true. + fp_ks_model = FP_KPERP_FULL kr_arr = [0.1_dp, 1.0_dp, 5.0_dp] jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] @@ -368,7 +515,7 @@ subroutine test_collapse_rho_phi_ks2() gfused = hatG_rho_phi(plasma, kr, kr, j) - tol = 1.0e-12_dp * (1.0_dp + abs(ref)) + tol = 5.0e-12_dp * (1.0_dp + abs(ref)) if (abs(gfused - ref) >= tol) then print *, 'FAIL: rho-Phi ks2 collapse mismatch at kr=', kr, ' j=', j print *, ' core(bf,bf)/4pi = ', ref @@ -379,15 +526,14 @@ subroutine test_collapse_rho_phi_ks2() end do end do - ! Defensive: restore the default so any later test sees the drop-k_s^2 - ! path in the shared module state. - kern_include_ks2 = .false. + fp_ks_model = FP_KR_ONLY_APPROX - print *, 'PASS: fused rho-Phi kernel (ks2 path) collapses to core(bf,bf) at kr=krp' + print *, 'PASS: full rho-Phi model collapses to core(bf,bf) at kr=krp' end subroutine subroutine test_zero_wavenumber_limit() - ! Analytic b -> 0 anchor: at kr = krp = 0 in the default drop-k_s^2 path, + ! Analytic b -> 0 anchor in the named radial-only approximation: at + ! kr = krp = 0, ! both FLR arguments vanish (b_+ = b_x = 0). The Bessel/exp assembly then ! reduces to a closed form INDEPENDENT of the Bessel routine, because ! exp(-b_+) = 1, (1 - b_+) = 1, I_0(0) = 1, I_{-1}(0) = 0, and the @@ -398,10 +544,10 @@ subroutine test_zero_wavenumber_limit() ! I_{-1}(0)=0 here rather than calling bessel_in, so this test validates ! the kernel's argument assembly against math, not against the Bessel ! routine it shares with the collapse/phase tests. Only valid for the - ! .false. path (in the .true. path the arguments are rho_L^2*ks^2 /= 0). + ! radial-only model (the full model has arguments rho_L^2*ks^2 /= 0). use config_m, only: profiles_in_memory, nml_config_path use config_m, only: artificial_debye_case, turn_off_ions, turn_off_electrons - use flr2_fourier_kernel_m, only: hatG_rho_phi, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_phi use constants_m, only: com_unit, pi use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -431,9 +577,8 @@ subroutine test_zero_wavenumber_limit() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - ! Closed form derived above only holds for the drop-k_s^2 path; the - ! toggle is deliberately fixed here (do NOT sweep it). - kern_include_ks2 = .false. + ! This closed form only holds for the named radial-only approximation. + fp_ks_model = FP_KR_ONLY_APPROX jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] @@ -482,14 +627,14 @@ subroutine test_zero_wavenumber_limit() end subroutine subroutine test_collapse_rho_B() - ! Collapse-by-construction for the rho-B kernel (drop-k_s^2 path): the + ! Collapse-by-construction for the radial-only rho-B approximation: the ! fused two-wavenumber hatG_rho_B evaluated on the diagonal k_r = k'_r ! must equal the per-species signed diagonal integrand summed over the ! non-turned-off species and divided by 4*pi. Exact to machine precision ! because both share the same per-species core. use config_m, only: profiles_in_memory, nml_config_path use config_m, only: turn_off_ions, turn_off_electrons - use flr2_fourier_kernel_m, only: hatG_rho_B, hatG_rho_B_diag_sp, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_B, hatG_rho_B_diag_sp use constants_m, only: pi use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -519,7 +664,7 @@ subroutine test_collapse_rho_B() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - kern_include_ks2 = .false. + fp_ks_model = FP_KR_ONLY_APPROX kr_arr = [0.1_dp, 1.0_dp, 5.0_dp] jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] @@ -539,7 +684,7 @@ subroutine test_collapse_rho_B() gfused = hatG_rho_B(plasma, kr, kr, j) - tol = 1.0e-12_dp * (1.0_dp + abs(dref)) + tol = 5.0e-12_dp * (1.0_dp + abs(dref)) if (abs(gfused - dref) >= tol) then print *, 'FAIL: rho-B collapse mismatch at kr=', kr, ' j=', j print *, ' diag/4pi = ', dref @@ -554,15 +699,15 @@ subroutine test_collapse_rho_B() end subroutine subroutine test_collapse_rho_B_ks2() - ! Structural collapse for the rho-B k_s^2-inclusive path - ! (kern_include_ks2 = .true.). On the diagonal k_r = k'_r both FLR + ! Structural collapse for the full finite-radius rho-B model. On the + ! diagonal k_r = k'_r both FLR ! arguments reduce to the full perpendicular argument ! bf = (k_s^2 + k_r^2) * rho_L^2, so the fused kernel must equal the ! per-species core evaluated at (bf, bf), summed over the non-turned-off ! species and divided by 4*pi. use config_m, only: profiles_in_memory, nml_config_path use config_m, only: turn_off_ions, turn_off_electrons - use flr2_fourier_kernel_m, only: hatG_rho_B, core_rho_B_sp, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_B, core_rho_B_sp use constants_m, only: pi use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -592,7 +737,7 @@ subroutine test_collapse_rho_B_ks2() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - kern_include_ks2 = .true. + fp_ks_model = FP_KPERP_FULL kr_arr = [0.1_dp, 1.0_dp, 5.0_dp] jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] @@ -613,7 +758,7 @@ subroutine test_collapse_rho_B_ks2() gfused = hatG_rho_B(plasma, kr, kr, j) - tol = 1.0e-12_dp * (1.0_dp + abs(ref)) + tol = 5.0e-12_dp * (1.0_dp + abs(ref)) if (abs(gfused - ref) >= tol) then print *, 'FAIL: rho-B ks2 collapse mismatch at kr=', kr, ' j=', j print *, ' core(bf,bf)/4pi = ', ref @@ -624,15 +769,13 @@ subroutine test_collapse_rho_B_ks2() end do end do - ! Defensive: restore the default so any later test sees the drop-k_s^2 - ! path in the shared module state. - kern_include_ks2 = .false. + fp_ks_model = FP_KR_ONLY_APPROX - print *, 'PASS: fused rho-B kernel (ks2 path) collapses to core(bf,bf) at kr=krp' + print *, 'PASS: full rho-B model collapses to core(bf,bf) at kr=krp' end subroutine subroutine test_zero_wavenumber_limit_rho_B() - ! Analytic b -> 0 anchor for the rho-B kernel (drop-k_s^2 path): at + ! Analytic b -> 0 anchor for the radial-only rho-B approximation: at ! kr = krp = 0 both FLR arguments vanish (b_+ = b_x = 0). With ! I_0(0) = 1, I_{-1}(0) = 0, exp(-b_+) = 1, (1 - b_+) = 1 and the ! A2 * b_x * I_{-1}(b_x) term vanishing, and the Fourier phase = 1, the @@ -640,10 +783,11 @@ subroutine test_zero_wavenumber_limit_rho_B() ! -1/lambda_D^2 * vT^3/(omega_c*nu*sol) * ( I01*(A1+A2) + 0.5*I21*A2 ). ! We fold I_0(0)=1 and I_{-1}(0)=0 in analytically (no bessel_in call), ! so this validates the kernel's argument assembly against math, not the - ! Bessel routine. rho-B has NO Debye term. Only valid for the .false. path. + ! Bessel routine. rho-B has NO Debye term. Only valid for this named + ! approximation. use config_m, only: profiles_in_memory, nml_config_path use config_m, only: artificial_debye_case, turn_off_ions, turn_off_electrons - use flr2_fourier_kernel_m, only: hatG_rho_B, kern_include_ks2 + use flr2_fourier_kernel_m, only: hatG_rho_B use constants_m, only: pi, sol use species_m, only: plasma, set_profiles_from_arrays use grid_m, only: rg_grid @@ -673,8 +817,8 @@ subroutine test_zero_wavenumber_limit_rho_B() call from_kim_factory_get_kim('electrostatic', kim_instance) call kim_instance%init() - ! Closed form derived above only holds for the drop-k_s^2 path. - kern_include_ks2 = .false. + ! Closed form derived above only holds for the radial-only approximation. + fp_ks_model = FP_KR_ONLY_APPROX jj_arr = [rg_grid%npts_b / 4, rg_grid%npts_b / 2, (3 * rg_grid%npts_b) / 4] diff --git a/verification/FORMULA_INDEX.md b/verification/FORMULA_INDEX.md index 3154671f..68603a36 100644 --- a/verification/FORMULA_INDEX.md +++ b/verification/FORMULA_INDEX.md @@ -561,9 +561,9 @@ path exists and that each code family has a source-to-check chain. | KIM-01 | `KIM/src/background_equilibrium/species_mod.f90::calculate_plasma_backs` | vT, gyrofrequency, Larmor/Debye lengths, collision scale | (4.14)–(4.16), (5.8) | conventions script; #198 pending | | KIM-02 | `KIM/src/background_equilibrium/species_mod.f90::calculate_thermodynamic_forces_and_susc` | A1/A2, x1/x2, FP susceptibilities | (5.8), (5.12)–(5.20) | flow test; independent #194 oracle | | KIM-03 | `KIM/src/background_equilibrium/calculate_equil.f90::calculate_equil` | h_theta/h_z, k_s, k_parallel, omega_E | (3.5), (5.18)–(5.20) | conventions script; #198 pending | -| KIM-04 | `KIM/src/asymptotics/flr2_fourier_kernel.f90` | b_plus, b_cross, phases, charge/current kernels | Ch. 14 | independent #196 oracle; production #187 pending | -| KIM-05 | `KIM/src/kernels/FP_kernel_plasma_prefacs.f90` | finite-radius FP kernel prefactors | Ch. 14 | independent #196 derivation; production #187 pending | -| KIM-06 | `KIM/src/kernels/kernel.f90` | real-space charge/current kernels and 1/(4 pi) normalization | Ch. 13–14 | independent #196 derivation; #187/#197 pending | +| KIM-04 | `KIM/src/asymptotics/flr2_fourier_kernel.f90` | b_plus, b_cross, phases, charge/current kernels | Ch. 14 | independent #196 oracle; production #187 verified | +| KIM-05 | `KIM/src/kernels/FP_kernel_plasma_prefacs.f90` | finite-radius FP kernel prefactors | Ch. 14 | independent #196 derivation; production Fourier consumption #187 verified | +| KIM-06 | `KIM/src/kernels/kernel.f90` | real-space charge/current kernels and 1/(4 pi) normalization | Ch. 13–14 | independent #196 derivation; off-diagonal normalization #187 verified; #197 pending | | KIM-07 | `KIM/src/electrostatic_poisson/poisson.f90` | global Gaussian-CGS weak Poisson problem | (12.25), Ch. 15 | #197 pending | | KIM-08 | `KIM/src/electrostatic_poisson/periodic_solve.f90` | periodic Fourier Poisson matrix/RHS | (12.25), Ch. 15 | #197/#188 pending | | KIM-09 | `KIM/src/asymptotics/FLR2_asymptotics.f90` | aligned potential and FLR2 limits | Ch. 5 | #197 pending | @@ -600,12 +600,18 @@ upstream Weber integral proves that the energy moment is plasma-dispersion integral, homogeneous harmonic sum (13.70)–(13.73), and Horton small-FLR limit so Krook assumptions are never substituted for FP. -`verification/oracles/rho_phi_kernels.dat` records signed electron/ion, +`verification/oracles/rho_phi_kernels.dat` records rho-Phi and rho-B values for signed electron/ion, diagonal/off-diagonal, finite-`k_s`, `m_phi=-2..2`, zero-FLR, Debye-only, and `b_x=600` cases at 34 committed digits. Six negative fixtures independently mutate `4 pi`, Fourier measure, phase, Bessel order, the `b_+` sign, and hidden -removal of `k_s^2`. Production consumption and the finite-radius default are -completed by #187. +removal of `k_s^2`. The #187 production path defaults to the full +`k_perp^2=k_s^2+k_r^2` model, sums every configured harmonic, includes the +Debye source once rather than once per harmonic, rejects undefined `k_s`, and +uses a stabilized scaled-Bessel expansion for large arguments. The explicitly +named radial-only model remains available solely as a compatibility +approximation. `test_flr2_fourier_kernel` consumes all 13 rows and currently +agrees within `1.90e-16` scaled error; the independent script passes 52 checks +and rejects all six mutations. ### Fokker–Planck susceptibility oracle (#194) diff --git a/verification/mathematica/04_charge_potential_kernels.wl b/verification/mathematica/04_charge_potential_kernels.wl index 9ac0e8b1..58ded120 100644 --- a/verification/mathematica/04_charge_potential_kernels.wl +++ b/verification/mathematica/04_charge_potential_kernels.wl @@ -213,7 +213,8 @@ x1FromCase[c_] := c["kpar"] c["vT"]/c["nu"]; x2FromCase[c_] := -(c["omegaE"] + c["kpar"] c["Vpar"] + c["mphi"] c["omegaC"] - c["omega"])/c["nu"]; kernelParts[c_] := Module[ - {rho, bp, bx, ph, x1, x2, i00, i20, f0, f2, dynamic, debye, total}, + {rho, bp, bx, ph, x1, x2, i00, i20, i01, i21, f0, f2, dynamic, + rhoBDynamic, debye, total, speedOfLight = 29979245800}, rho = Abs[c["vT"]/c["omegaC"]]; bp = bPlus[c["ks"], c["kr"], c["krp"], rho]; bx = bCross[c["ks"], c["kr"], c["krp"], rho]; @@ -221,18 +222,24 @@ kernelParts[c_] := Module[ x1 = x1FromCase[c]; x2 = x2FromCase[c]; i00 = fpConserving[0, 0, x1, x2]; i20 = fpConserving[2, 0, x1, x2]; + i01 = fpConserving[0, 1, x1, x2]; + i21 = fpConserving[2, 1, x1, x2]; f0 = flrDensity[bp, bx, c["mphi"]]; f2 = flrEnergy[bp, bx, c["mphi"]]; dynamic = ph kernelNormalization I c["vT"]^2 c["ks"] / (c["lambda"]^2 c["omegaC"] c["nu"]) * (i00 (c["A1"] f0 + c["A2"] f2) + c["A2"] i20 f0/2); + rhoBDynamic = -ph kernelNormalization c["vT"]^3 / + (c["lambda"]^2 c["omegaC"] c["nu"] speedOfLight) * + (i01 (c["A1"] f0 + c["A2"] f2) + c["A2"] i21 f0/2); debye = If[c["debye"] == 1 && c["mphi"] == 0, -ph kernelNormalization/c["lambda"]^2, 0]; total = debye + dynamic; <|"rho" -> rho, "bp" -> bp, "bx" -> bx, "phase" -> ph, "x1" -> x1, "x2" -> x2, "I00" -> i00, "I20" -> i20, + "I01" -> i01, "I21" -> i21, "F0" -> f0, "F2" -> f2, "dynamic" -> dynamic, - "debyeTerm" -> debye, "total" -> total|> + "debyeTerm" -> debye, "total" -> total, "rhoBDynamic" -> rhoBDynamic|> ]; (* Geometry and limit proofs do not assume a collision model. *) @@ -283,9 +290,11 @@ oracleRows = Table[ {ToString[c["debye"]]}, formatNumber /@ {p["x1"], p["x2"], p["rho"], p["bp"], p["bx"], Re[p["phase"]], Im[p["phase"]], Re[p["I00"]], Im[p["I00"]], - Re[p["I20"]], Im[p["I20"]], p["F0"], p["F2"], + Re[p["I20"]], Im[p["I20"]], Re[p["I01"]], Im[p["I01"]], + Re[p["I21"]], Im[p["I21"]], p["F0"], p["F2"], Re[p["dynamic"]], Im[p["dynamic"]], Re[p["debyeTerm"]], - Im[p["debyeTerm"]], Re[p["total"]], Im[p["total"]]} + Im[p["debyeTerm"]], Re[p["total"]], Im[p["total"]], + Re[p["rhoBDynamic"]], Im[p["rhoBDynamic"]]} ] ], {c, cases} @@ -304,7 +313,7 @@ oracleHeader = { "# Full finite-radius b+/bx include ks^2. The upstream Gaussian moment proves 1-b+-mphi.", "# Case 6 has bx=600; case 7 is zero-FLR; case 8 is Debye-only; cases 9-13 record mphi=-2..2.", "# Values generated at 70-digit working precision and committed with 34 digits.", - "# id charge lambdaD vT omegaC nu ks kr krp rg mphi A1 A2 kpar omegaE Vpar omega debye x1 x2 rhoL bplus bcross phaseRe phaseIm I00Re I00Im I20Re I20Im F0 F2 dynamicRe dynamicIm debyeRe debyeIm totalRe totalIm" + "# id charge lambdaD vT omegaC nu ks kr krp rg mphi A1 A2 kpar omegaE Vpar omega debye x1 x2 rhoL bplus bcross phaseRe phaseIm I00Re I00Im I20Re I20Im I01Re I01Im I21Re I21Im F0 F2 rhoPhiDynamicRe rhoPhiDynamicIm debyeRe debyeIm rhoPhiTotalRe rhoPhiTotalIm rhoBDynamicRe rhoBDynamicIm" }; If[regenerate, @@ -322,6 +331,10 @@ If[regenerate, If[Length[committedRows] == Length[oracleRows], Do[ checkNumeric["kernel oracle row " <> ToString[row], + committedRows[[row, -4]] + I committedRows[[row, -3]], + parseNumber[oracleRows[[row, -4]]] + I parseNumber[oracleRows[[row, -3]]], + 10^-30]; + checkNumeric["rho-B oracle row " <> ToString[row], committedRows[[row, -2]] + I committedRows[[row, -1]], parseNumber[oracleRows[[row, -2]]] + I parseNumber[oracleRows[[row, -1]]], 10^-30], diff --git a/verification/oracles/rho_phi_kernels.dat b/verification/oracles/rho_phi_kernels.dat index f770eab6..98404c3f 100644 --- a/verification/oracles/rho_phi_kernels.dat +++ b/verification/oracles/rho_phi_kernels.dat @@ -2,17 +2,17 @@ # Full finite-radius b+/bx include ks^2. The upstream Gaussian moment proves 1-b+-mphi. # Case 6 has bx=600; case 7 is zero-FLR; case 8 is Debye-only; cases 9-13 record mphi=-2..2. # Values generated at 70-digit working precision and committed with 34 digits. -# id charge lambdaD vT omegaC nu ks kr krp rg mphi A1 A2 kpar omegaE Vpar omega debye x1 x2 rhoL bplus bcross phaseRe phaseIm I00Re I00Im I20Re I20Im F0 F2 dynamicRe dynamicIm debyeRe debyeIm totalRe totalIm -1 1 2. 1.1 3. 0.8 0.4 0.2 0.2 0.7 0 0.3 -0.2 0.5 0.1 0.25 0.7 1 0.6875 0.59375 0.3666666666666666666666666666666667 0.02688888888888888888888888888888889 0.02688888888888888888888888888888889 1. 0 0.4769623207499516668365775931719617 1.34337632825321339566065283021071 0.7839061781255760478643132032375341 0.3449420361888029943597827445343369 0.9736453643317647842445142523052613 0.9478170693156334195008871827455435 -0.0004178590961293672899336853811358227 -0.0001100163923241111610126729967799579 -0.0198943678864869169711104704215643 0 -0.02031222698261628426104415580270012 -0.0001100163923241111610126729967799579 -2 -1 1.5 0.9 -4. 1. 0.3333333333333333333333333333333333 -0.4 0.7 0.8 -1 -0.25 0.4 0.6 -0.2 0.3333333333333333333333333333333333 0.5 0 0.54 -3.5 0.225 0.022078125 0.02043694953387124259549447252480515 -0.9666062290617209856953406799442396 -0.2562662637552586252047982881698298 0.002077683956804036006494684747967483 -0.2924084728658036108361060112168703 0.006934065296765442912170263729614332 -0.3061580468275320945302813238153212 0.009995864093294536935210747909060029 0.01977208197189641732791396010577808 5.047415454043973692153457479209403e-6 1.400166683123637810553028951820593e-6 0 0 5.047415454043973692153457479209403e-6 1.400166683123637810553028951820593e-6 -3 1 2.5 1.2 2.5 0.75 0.5 -0.75 0.6666666666666666666666666666666667 0.6 1 0.2 0.3333333333333333333333333333333333 0.4 0.1666666666666666666666666666666667 -0.2 0.75 0 0.64 -2.448888888888888888888888888888889 0.48 0.1736 0.1730664612222714860697226208385838 0.4887944915408868160211973948072708 0.8723989597880581648166673990205585 0.0148238114854477958981925926578536 -0.435814834304506581458537913442345 0.05282503413832141426335510727730936 -0.4907697501249466210393516461114387 0.07301538524579436045614679393210793 0.1339013577618891883343405199639868 0.0000694820099044233620111156542140318 0.0001392248759887814818924784799184243 0 0 0.0000694820099044233620111156542140318 0.0001392248759887814818924784799184243 -4 1 1.75 1. 2. 1. 0.25 0.6 -0.8 0.9 2 -0.4 0.5 0.3333333333333333333333333333333333 0.1 0.2 0.6666666666666666666666666666666667 0 0.3333333333333333333333333333333333 -3.5 0.5 0.140625 0.1361998187407017072936212763966931 -0.8840958541814920694768016513522915 0.4673055966057949018049342890624715 0.0007212265367171929539324657251207403 -0.2881669773376818932432168415721195 0.002255439536090886009723836672798034 -0.2931278873860203081135294536674759 0.002017723452994739431459647998928217 0.00577566384872274732347176787834151 -2.150419740868720426700321048562586e-6 1.126951894879956860071538985849193e-6 0 0 -2.150419740868720426700321048562586e-6 1.126951894879956860071538985849193e-6 -5 -1 1.8 0.8 -3. 1.2 0.2857142857142857142857142857142857 -0.5 0.8333333333333333333333333333333333 0.4 -2 0.1666666666666666666666666666666667 -0.3 0.5 -0.125 0.25 0.6 0 0.3333333333333333333333333333333333 -4.5 0.2666666666666666666666666666666667 0.03938523557571176618795666414714034 0.03607601350712578529252483017641597 0.3942418351125960285331024428597845 0.9190067330803690282515360107946175 0.0002666613214760912095989950399634674 -0.2234018321923041475564299573377105 0.0008248220506996469140035738560665456 -0.2257837005672125861481190239162597 0.000156418964203493298504921476517882 0.0004631302232632412750048735340668457 1.495045516406575146033182851936704e-8 3.500709274575635256396610564120795e-8 0 0 1.495045516406575146033182851936704e-8 3.500709274575635256396610564120795e-8 -6 1 2. 1. 1. 1. 0.1 24.49469330283602469094718087269004 24.49469330283602469094718087269004 0.3333333333333333333333333333333333 0 0.25 0.2 0.5 0 0 0.3333333333333333333333333333333333 0 0.5 0.3333333333333333333333333333333333 1. 600. 600. 1. 0 0.7578298877360103424236057937223274 2.37564007535752074950740987216793 1.169293182961532486178389412097128 0.7329465504735785674569898903710722 0.01629014665630598169061102304622454 0.008141673876517436643251074954017288 -0.00002931872652008470597166766142424933 0.00001238441610829100447778937328371508 0 0 -0.00002931872652008470597166766142424933 0.00001238441610829100447778937328371508 -7 1 2. 1. 1.e8 1. 0.3333333333333333333333333333333333 0.4 -0.6 0.5 0 0.3333333333333333333333333333333333 -0.25 0.5 0 0 0.25 1 0.5 0.25 1.e-8 3.711111111111111111111111111111111e-17 3.573841865417711186815181460669049e-17 0.8775825618903727161162815826038297 0.4794255386042030002732879352155714 1.082366905762994932172263587775514 3.042045058014564051046381804906527 1.228546668426184681996684092037351 0.8428781702666359449338590390021461 0.99999999999999996288888888888889 0.999999999999999925777777777777781 -6.606711147386609432971320169072563e-12 -8.397901678580790655353527013053519e-12 -0.01745895033701274825900937419124877 -0.009537868039169149864391913156104389 -0.01745895034361945940639598362422009 -0.009537868047567051542972703811457916 -8 -1 3. 1. -5. 1. 0.4 0.3333333333333333333333333333333333 -0.25 0.75 0 0 0 0.5 0 0 0.3333333333333333333333333333333333 1 0.5 0.3333333333333333333333333333333333 0.2 0.009872222222222222222222222222222222 0.009824233077220610992864928663855168 0.9058136834259364207445166606527002 0.4236762572039380103616839880311025 0.7578298877360103424236057937223274 2.37564007535752074950740987216793 1.169293182961532486178389412097128 0.7329465504735785674569898903710722 0.9902002401987167834328944214018998 0.9804725476696356878873796996807079 0 0 -0.00800915140208416713042740975492236 -0.003746120589148886962650392576170468 -0.00800915140208416713042740975492236 -0.003746120589148886962650392576170468 -9 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 -2 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 20.9 0.2 0.0187 0.0169955876626846887032793677040316 0.9716675638238312370428044040637474 -0.2363517408708064698861197740969023 1.311762434073074816495327966437032e-6 0.04787425850143288323013032578022972 3.953075520801203478061608872986528e-6 0.04792908738307940406756942566656722 0.00003543818996847633211279510916319949 0.0001056535817962807876375553015996531 7.483725232672084528061297431773954e-10 -1.820679875767347984781833965519952e-10 0 0 7.483725232672084528061297431773954e-10 -1.820679875767347984781833965519952e-10 -10 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 -1 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 10.9 0.2 0.0187 0.0169955876626846887032793677040316 -0.9717653852472468488004535641904121 -0.2359492234258674488796405567793635 0.00001778699334156899338252558591527229 0.09193575400971078093503566930381928 0.00005421589225719964755627780385339402 0.09232324848464393967782959549298052 0.008340662807935853625112714843862475 0.01652595751422752290008315323905084 -1.902839687412108709961870079521953e-7 -4.613182649408617525486807342356357e-8 0 0 -1.902839687412108709961870079521953e-7 -4.613182649408617525486807342356357e-8 -11 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 0 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 0.9 0.2 0.0187 0.0169955876626846887032793677040316 0.7550932412115528473007442832390617 0.655617416697140275668825905436999 0.2690284498695408509328378171555757 1.126581811092725591534733630462194 0.8159576576435002274973534579201671 1.018627487470777979930753104457581 0.9815446361623944909477242149063205 0.9633315059320748815941626317668765 0.00008349620193227154873549819076117614 -0.00001780412501276963834263734267993515 0 0 0.00008349620193227154873549819076117614 -0.00001780412501276963834263734267993515 -12 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 1 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 -9.1 0.2 0.0187 0.0169955876626846887032793677040316 -0.3699837991399858731379044236481369 -0.9290382060894711065696532922191437 0.00003668408398036629216948820610267005 -0.1102208437279178642803478408665778 0.0001125242814462708135598658460174802 -0.1108875770923986972573881793473588 0.008340662807935853625112714843862475 0.01652595751422752290008315323905084 8.678039981580448341430975492881512e-8 2.182885490007452953032826565391461e-7 0 0 8.678039981580448341430975492881512e-8 2.182885490007452953032826565391461e-7 -13 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 2 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 -19.1 0.2 0.0187 0.0169955876626846887032793677040316 -0.09765726689923169068126435331717347 0.9952201054147631247835241473834049 1.88107602764291949540440505407339e-6 -0.05239187523155515129255189496374619 5.673691764095093509148800897345533e-6 -0.05246372710305089119287610344311752 0.00003543818996847633211279510916319949 0.0001056535817962807876375553015996531 8.236093430026791294796165699554261e-11 -8.38889268468837257086372476680649e-10 0 0 8.236093430026791294796165699554261e-11 -8.38889268468837257086372476680649e-10 +# id charge lambdaD vT omegaC nu ks kr krp rg mphi A1 A2 kpar omegaE Vpar omega debye x1 x2 rhoL bplus bcross phaseRe phaseIm I00Re I00Im I20Re I20Im I01Re I01Im I21Re I21Im F0 F2 rhoPhiDynamicRe rhoPhiDynamicIm debyeRe debyeIm rhoPhiTotalRe rhoPhiTotalIm rhoBDynamicRe rhoBDynamicIm +1 1 2. 1.1 3. 0.8 0.4 0.2 0.2 0.7 0 0.3 -0.2 0.5 0.1 0.25 0.7 1 0.6875 0.59375 0.3666666666666666666666666666666667 0.02688888888888888888888888888888889 0.02688888888888888888888888888888889 1. 0 0.4769623207499516668365775931719617 1.34337632825321339566065283021071 0.7839061781255760478643132032375341 0.3449420361888029943597827445343369 0.4119220042840491668134079213757851 -0.2943568074176793401112543739089323 0.6770098811084520413373614027960522 -1.156640968746033777598369447902164 0.9736453643317647842445142523052613 0.9478170693156334195008871827455435 -0.0004178590961293672899336853811358227 -0.0001100163923241111610126729967799579 -0.0198943678864869169711104704215643 0 -0.02031222698261628426104415580270012 -0.0001100163923241111610126729967799579 8.715660611107301685524669448235418e-15 -3.033819550037769938817898681825842e-14 +2 -1 1.5 0.9 -4. 1. 0.3333333333333333333333333333333333 -0.4 0.7 0.8 -1 -0.25 0.4 0.6 -0.2 0.3333333333333333333333333333333333 0.5 0 0.54 -3.5 0.225 0.022078125 0.02043694953387124259549447252480515 -0.9666062290617209856953406799442396 -0.2562662637552586252047982881698298 0.002077683956804036006494684747967483 -0.2924084728658036108361060112168703 0.006934065296765442912170263729614332 -0.3061580468275320945302813238153212 -0.0134664700903965296717248085516411 0.04338825005613451467846488751675209 -0.04494301581236861146777022787712993 0.1325058590673376497333048765807853 0.009995864093294536935210747909060029 0.01977208197189641732791396010577808 5.047415454043973692153457479209403e-6 1.400166683123637810553028951820593e-6 0 0 5.047415454043973692153457479209403e-6 1.400166683123637810553028951820593e-6 6.134282515848604783953922068925282e-17 -9.487173827847892652473852120568338e-17 +3 1 2.5 1.2 2.5 0.75 0.5 -0.75 0.6666666666666666666666666666666667 0.6 1 0.2 0.3333333333333333333333333333333333 0.4 0.1666666666666666666666666666666667 -0.2 0.75 0 0.64 -2.448888888888888888888888888888889 0.48 0.1736 0.1730664612222714860697226208385838 0.4887944915408868160211973948072708 0.8723989597880581648166673990205585 0.0148238114854477958981925926578536 -0.435814834304506581458537913442345 0.05282503413832141426335510727730936 -0.4907697501249466210393516461114387 -0.05672166755890094124933415662831483 0.1050970395957161554420443771300841 -0.2021291236820493004104768340958157 0.3153759188808721402269635903291857 0.07301538524579436045614679393210793 0.1339013577618891883343405199639868 0.0000694820099044233620111156542140318 0.0001392248759887814818924784799184243 0 0 0.0000694820099044233620111156542140318 0.0001392248759887814818924784799184243 4.549772719983424964648322420903311e-15 6.191203949457254417238813639677053e-17 +4 1 1.75 1. 2. 1. 0.25 0.6 -0.8 0.9 2 -0.4 0.5 0.3333333333333333333333333333333333 0.1 0.2 0.6666666666666666666666666666666667 0 0.3333333333333333333333333333333333 -3.5 0.5 0.140625 0.1361998187407017072936212763966931 -0.8840958541814920694768016513522915 0.4673055966057949018049342890624715 0.0007212265367171929539324657251207403 -0.2881669773376818932432168415721195 0.002255439536090886009723836672798034 -0.2931278873860203081135294536674759 -0.007572878635530526016290890113767773 0.02575326204565987905377683650725453 -0.02368211512895430310210028506437936 0.07784281755321323519205926350849654 0.002017723452994739431459647998928217 0.00577566384872274732347176787834151 -2.150419740868720426700321048562586e-6 1.126951894879956860071538985849193e-6 0 0 -2.150419740868720426700321048562586e-6 1.126951894879956860071538985849193e-6 8.189946099765136770153184382303977e-18 4.118623409312546921139243055149149e-17 +5 -1 1.8 0.8 -3. 1.2 0.2857142857142857142857142857142857 -0.5 0.8333333333333333333333333333333333 0.4 -2 0.1666666666666666666666666666666667 -0.3 0.5 -0.125 0.25 0.6 0 0.3333333333333333333333333333333333 -4.5 0.2666666666666666666666666666666667 0.03938523557571176618795666414714034 0.03607601350712578529252483017641597 0.3942418351125960285331024428597845 0.9190067330803690282515360107946175 0.0002666613214760912095989950399634674 -0.2234018321923041475564299573377105 0.0008248220506996469140035738560665456 -0.2257837005672125861481190239162597 -0.00359992783992723132958643303950681 0.01592473459610599201180442405909179 -0.01113509768444523333904824705689837 0.0480799576573699129996068228695054 0.000156418964203493298504921476517882 0.0004631302232632412750048735340668457 1.495045516406575146033182851936704e-8 3.500709274575635256396610564120795e-8 0 0 1.495045516406575146033182851936704e-8 3.500709274575635256396610564120795e-8 3.439308554134740367471013688517628e-19 -6.290138297389169401170684486857027e-20 +6 1 2. 1. 1. 1. 0.1 24.49469330283602469094718087269004 24.49469330283602469094718087269004 0.3333333333333333333333333333333333 0 0.25 0.2 0.5 0 0 0.3333333333333333333333333333333333 0 0.5 0.3333333333333333333333333333333333 1. 600. 600. 1. 0 0.7578298877360103424236057937223274 2.37564007535752074950740987216793 1.169293182961532486178389412097128 0.7329465504735785674569898903710722 0.5052199251573402282824038624815516 -0.4162399497616528336617267518880466 0.7795287886410216574522596080647518 -1.511368966350947621695340073085952 0.01629014665630598169061102304622454 0.008141673876517436643251074954017288 -0.00002931872652008470597166766142424933 0.00001238441610829100447778937328371508 0 0 -0.00002931872652008470597166766142424933 0.00001238441610829100447778937328371508 -2.753997702479672226619606573228087e-15 3.208510491537903695645312684831888e-15 +7 1 2. 1. 1.e8 1. 0.3333333333333333333333333333333333 0.4 -0.6 0.5 0 0.3333333333333333333333333333333333 -0.25 0.5 0 0 0.25 1 0.5 0.25 1.e-8 3.711111111111111111111111111111111e-17 3.573841865417711186815181460669049e-17 0.8775825618903727161162815826038297 0.4794255386042030002732879352155714 1.082366905762994932172263587775514 3.042045058014564051046381804906527 1.228546668426184681996684092037351 0.8428781702666359449338590390021461 0.5411834528814974660861317938877571 -0.4789774709927179744768090975467363 0.6142733342130923409983420460186756 -1.578560914866682027533070480498927 0.99999999999999996288888888888889 0.999999999999999925777777777777781 -6.606711147386609432971320169072563e-12 -8.397901678580790655353527013053519e-12 -0.01745895033701274825900937419124877 -0.009537868039169149864391913156104389 -0.01745895034361945940639598362422009 -0.009537868047567051542972703811457916 6.853099647972290040704898205096596e-22 -8.158708249399612065099168709414146e-22 +8 -1 3. 1. -5. 1. 0.4 0.3333333333333333333333333333333333 -0.25 0.75 0 0 0 0.5 0 0 0.3333333333333333333333333333333333 1 0.5 0.3333333333333333333333333333333333 0.2 0.009872222222222222222222222222222222 0.009824233077220610992864928663855168 0.9058136834259364207445166606527002 0.4236762572039380103616839880311025 0.7578298877360103424236057937223274 2.37564007535752074950740987216793 1.169293182961532486178389412097128 0.7329465504735785674569898903710722 0.5052199251573402282824038624815516 -0.4162399497616528336617267518880466 0.7795287886410216574522596080647518 -1.511368966350947621695340073085952 0.9902002401987167834328944214018998 0.9804725476696356878873796996807079 0 0 -0.00800915140208416713042740975492236 -0.003746120589148886962650392576170468 -0.00800915140208416713042740975492236 -0.003746120589148886962650392576170468 0 0 +9 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 -2 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 20.9 0.2 0.0187 0.0169955876626846887032793677040316 0.9716675638238312370428044040637474 -0.2363517408708064698861197740969023 1.311762434073074816495327966437032e-6 0.04787425850143288323013032578022972 3.953075520801203478061608872986528e-6 0.04792908738307940406756942566656722 0.00005483166974425452732950470899706795 0.001144005359894519019447617613602227 0.0001652385567694903053829752508908369 0.003435852612719090024401992862509665 0.00003543818996847633211279510916319949 0.0001056535817962807876375553015996531 7.483725232672084528061297431773954e-10 -1.820679875767347984781833965519952e-10 0 0 7.483725232672084528061297431773954e-10 -1.820679875767347984781833965519952e-10 1.666371547849503411079556181669591e-21 5.654688777140692179777192278395748e-21 +10 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 -1 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 10.9 0.2 0.0187 0.0169955876626846887032793677040316 -0.9717653852472468488004535641904121 -0.2359492234258674488796405567793635 0.00001778699334156899338252558591527229 0.09193575400971078093503566930381928 0.00005421589225719964755627780385339402 0.09232324848464393967782959549298052 0.0003877564548462040557390577729529359 0.004199437411695024383777590823260305 0.00118190645120695231672685612400399 0.01264681696523788497668518174697524 0.008340662807935853625112714843862475 0.01652595751422752290008315323905084 -1.902839687412108709961870079521953e-7 -4.613182649408617525486807342356357e-8 0 0 -1.902839687412108709961870079521953e-7 -4.613182649408617525486807342356357e-8 5.156808787467341959119170633504908e-19 -3.521847500248144406137878979168398e-18 +11 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 0 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 0.9 0.2 0.0187 0.0169955876626846887032793677040316 0.7550932412115528473007442832390617 0.655617416697140275668825905436999 0.2690284498695408509328378171555757 1.126581811092725591534733630462194 0.8159576576435002274973534579201671 1.018627487470777979930753104457581 0.4842512097651735316791080708800362 0.02784725996690606476252053483194908 1.468723783758300409495236224256301 -0.1664705225525996361246444119763546 0.9815446361623944909477242149063205 0.9633315059320748815941626317668765 0.00008349620193227154873549819076117614 -0.00001780412501276963834263734267993515 0 0 0.00008349620193227154873549819076117614 -0.00001780412501276963834263734267993515 2.315170268081289265638711228835109e-14 1.496501236158219097229221303127778e-14 +12 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 1 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 -9.1 0.2 0.0187 0.0169955876626846887032793677040316 -0.3699837991399858731379044236481369 -0.9290382060894711065696532922191437 0.00003668408398036629216948820610267005 -0.1102208437279178642803478408665778 0.0001125242814462708135598658460174802 -0.1108875770923986972573881793473588 -0.0006676503284426665174846853510685949 0.006019355848105129902330703771715473 -0.00204794192232212880678955839751814 0.01815390308165629008446486412193028 0.008340662807935853625112714843862475 0.01652595751422752290008315323905084 8.678039981580448341430975492881512e-8 2.182885490007452953032826565391461e-7 0 0 8.678039981580448341430975492881512e-8 2.182885490007452953032826565391461e-7 4.935030327244258009942639720131854e-18 -1.351445831832823713374145495197086e-18 +13 1 2. 1. 5. 0.5 0.15 0.8 -0.5 0.55 2 0.2 -0.1666666666666666666666666666666667 0.25 0.1 0.2 0.6 0 0.5 -19.1 0.2 0.0187 0.0169955876626846887032793677040316 -0.09765726689923169068126435331717347 0.9952201054147631247835241473834049 1.88107602764291949540440505407339e-6 -0.05239187523155515129255189496374619 5.673691764095093509148800897345533e-6 -0.05246372710305089119287610344311752 -0.00007185710425595952472444827306560348 0.001369633845406779375482387615104354 -0.0002167350253884325720494841942785993 0.004114375336544043567867151527089332 0.00003543818996847633211279510916319949 0.0001056535817962807876375553015996531 8.236093430026791294796165699554261e-11 -8.38889268468837257086372476680649e-10 0 0 8.236093430026791294796165699554261e-11 -8.38889268468837257086372476680649e-10 -6.980463824529203206251372484473233e-21 -1.057327059694032471450276019543135e-21