module deviation use constants, only: dp, pi use field_base, only: field_t use fieldline_mod, only: flock_of_fieldlines_t implicit none private public :: calc_deviation contains subroutine calc_deviation(flock, deviation_A, deviation_B) use fieldline_labels, only: fieldline_modes_t use fieldline_labels, only: allocate_modes use fieldline_labels, only: fourier_transform_over_label use surface_average_mod, only: surface_average_t use surface_average_mod, only: calc_surface_averages use fit_functions, only: S_A, S_B use error_handling, only: failed_sanity_check type(flock_of_fieldlines_t), intent(in) :: flock real(dp), intent(out) :: deviation_A, deviation_B real(dp), parameter :: tol = 1e-12 type(fieldline_modes_t) :: modes real(dp) :: iota_p, eta_b type(surface_average_t) :: average real(dp) :: symmetric_remainder real(dp) :: B_squared_sqrtg logical :: any_has_sin_part call fourier_transform_over_label(flock, modes) any_has_sin_part = .false. if (has_sin_modes(modes%delta_aspect_ratio)) then print *, "error: non-vanishing sin part of delta aspect ratio: " print *, "sin part: ", sum(abs(modes%delta_aspect_ratio%sin_coeffs)) print *, "cos part: ", sum(abs(modes%delta_aspect_ratio%cos_coeffs)) any_has_sin_part = .true. end if if (has_sin_modes(modes%delta_eta)) then print *, "error: non-vanishing sin part of delta eta: " print *, "sin part: ", sum(abs(modes%delta_eta%sin_coeffs)) print *, "cos part: ", sum(abs(modes%delta_eta%cos_coeffs)) any_has_sin_part = .true. end if if (any_has_sin_part) call failed_sanity_check() call calc_surface_averages(flock, average) iota_p = flock%iota_p eta_b = flock%eta_b deviation_A = pi*sum(modes%radial_drift%sin_coeffs* & modes%delta_aspect_ratio%cos_coeffs* & S_A(iota_p*modes%delta_aspect_ratio%mode_numbers)) symmetric_remainder = pi*sum(modes%radial_drift%cos_coeffs* & modes%delta_aspect_ratio%sin_coeffs* & S_A(iota_p*modes%delta_aspect_ratio%mode_numbers)) if (abs(symmetric_remainder/deviation_A) > tol) then print *, "warning: non-vanishing symmetric part of deviation A: " print *, "symmetric: ", symmetric_remainder print *, "antisymmetric: ", deviation_A print *, "ratio: ", symmetric_remainder/deviation_A end if deviation_A = deviation_A*average%B_squared/average%lambda_b* & sqrt(eta_b)*sqrt(flock%I_ref)/average%normalization deviation_B = pi*sum(modes%radial_drift%sin_coeffs* & modes%delta_eta%cos_coeffs* & S_B(iota_p*modes%delta_eta%mode_numbers)) symmetric_remainder = pi*sum(modes%radial_drift%cos_coeffs* & modes%delta_eta%sin_coeffs* & S_B(iota_p*modes%delta_eta%mode_numbers)) if (abs(symmetric_remainder/deviation_B) > tol) then print *, "warning: non-vanishing symmetric part of deviation B: " print *, "symmetric: ", symmetric_remainder print *, "antisymmetric: ", deviation_B print *, "ratio: ", symmetric_remainder/deviation_B end if deviation_B = deviation_B*average%B_squared/average%lambda_b*0.5_dp/ & average%normalization end subroutine calc_deviation function has_sin_modes(modes) use fieldline_labels, only: modes_t type(modes_t), intent(in) :: modes logical :: has_sin_modes real(dp), parameter :: tol = 1e-3, numerical_zero = 1e-8 real(dp) :: sum_sin, sum_cos sum_sin = sum(abs(modes%sin_coeffs)) sum_cos = sum(abs(modes%cos_coeffs)) if (sum_cos > numerical_zero) then has_sin_modes = sum_sin/sum_cos > tol else has_sin_modes = sum_sin > numerical_zero end if end function has_sin_modes end module deviation