| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_invariants | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | use finite_check, only: is_finite | ||
| 4 | use field_can_mod, only: evaluate_field => evaluate, integ_to_ref | ||
| 5 | use parmot_mod, only: rho0 => ro0 | ||
| 6 | use simple, only: tracer_t, encode_symplectic_state | ||
| 7 | use util, only: twopi | ||
| 8 | use velo_mod, only: isw_field_type | ||
| 9 | use magfie_sub, only: TEST | ||
| 10 | |||
| 11 | implicit none | ||
| 12 | private | ||
| 13 | |||
| 14 | integer, parameter, public :: INVARIANT_SUCCESS = 0 | ||
| 15 | integer, parameter, public :: INVARIANT_INVALID_INPUT = 1 | ||
| 16 | integer, parameter, public :: INVARIANT_NOT_AXISYMMETRIC = 2 | ||
| 17 | integer, parameter, public :: INVARIANT_NO_INTERSECTION = 3 | ||
| 18 | integer, parameter, public :: INVARIANT_ROOT_FAILURE = 4 | ||
| 19 | integer, parameter, public :: INVARIANT_CAPACITY = 5 | ||
| 20 | |||
| 21 | integer, parameter, public :: SECTION_EXTREMUM_UNKNOWN = 0 | ||
| 22 | integer, parameter, public :: SECTION_B_MINIMUM = 1 | ||
| 23 | integer, parameter, public :: SECTION_B_MAXIMUM = 2 | ||
| 24 | |||
| 25 | integer, parameter :: N_THETA_SCAN = 128 | ||
| 26 | integer, parameter :: N_RADIAL_SCAN = 257 | ||
| 27 | integer, parameter :: MAX_SECTION_BRANCHES = 16 | ||
| 28 | integer, parameter :: MAX_ROOT_ITERATIONS = 80 | ||
| 29 | real(dp), parameter :: R_MIN = 1.0e-6_dp | ||
| 30 | real(dp), parameter :: UNIT_RADIAL_MAX = 1.0_dp - 1.0e-8_dp | ||
| 31 | real(dp), parameter :: TEST_RADIAL_MAX = 0.5_dp - 1.0e-8_dp | ||
| 32 | real(dp), parameter :: AXISYM_TOL = 1.0e-10_dp | ||
| 33 | real(dp), parameter :: ROOT_TOL = 1.0e-11_dp | ||
| 34 | |||
| 35 | type, public :: guiding_center_invariants_t | ||
| 36 | real(dp) :: h0 = 0.0_dp | ||
| 37 | real(dp) :: j_perp = 0.0_dp | ||
| 38 | real(dp) :: p_phi = 0.0_dp | ||
| 39 | end type guiding_center_invariants_t | ||
| 40 | |||
| 41 | type, public :: invariant_start_t | ||
| 42 | real(dp) :: state(5) = 0.0_dp | ||
| 43 | integer :: sigma = 0 | ||
| 44 | integer :: section_branch = 0 | ||
| 45 | integer :: section_kind = SECTION_EXTREMUM_UNKNOWN | ||
| 46 | real(dp) :: residual = 0.0_dp | ||
| 47 | end type invariant_start_t | ||
| 48 | |||
| 49 | type, public :: invariant_start_result_t | ||
| 50 | integer :: status = INVARIANT_INVALID_INPUT | ||
| 51 | type(invariant_start_t), allocatable :: starts(:) | ||
| 52 | end type invariant_start_result_t | ||
| 53 | |||
| 54 | public :: invariants_from_state, states_from_invariants | ||
| 55 | public :: potato_to_simple_invariants, invariant_flux_convention | ||
| 56 | |||
| 57 | contains | ||
| 58 | |||
| 59 | 3 | subroutine invariants_from_state(tracer, state, invariants, status) | |
| 60 | type(tracer_t), intent(in) :: tracer | ||
| 61 | real(dp), intent(in) :: state(5) | ||
| 62 | type(guiding_center_invariants_t), intent(out) :: invariants | ||
| 63 | integer, intent(out) :: status | ||
| 64 | |||
| 65 | type(tracer_t) :: work | ||
| 66 | real(dp) :: axis_flux, outward_sign, psi_star_native, zsympl(4) | ||
| 67 | |||
| 68 | invariants = guiding_center_invariants_t() | ||
| 69 | 3 | status = INVARIANT_INVALID_INPUT | |
| 70 |
4/6✓ Branch 0 taken 15 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
18 | if (.not. all(is_finite(state))) return |
| 71 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (rho0 <= 0.0_dp) return |
| 72 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (state(4) <= 0.0_dp) return |
| 73 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (abs(state(5)) > 1.0_dp) return |
| 74 | |||
| 75 | 3 | work = tracer | |
| 76 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
|
3 | if (.not. field_is_axisymmetric(work)) then |
| 77 | ✗ | status = INVARIANT_NOT_AXISYMMETRIC | |
| 78 | ✗ | return | |
| 79 | end if | ||
| 80 | 3 | call flux_convention(work, axis_flux, outward_sign) | |
| 81 | call encode_symplectic_state(work%f, state, zsympl, invariants%h0, & | ||
| 82 | 3 | invariants%j_perp, psi_star_native) | |
| 83 | 3 | invariants%p_phi = outward_sign*(psi_star_native - axis_flux) | |
| 84 | 3 | status = INVARIANT_SUCCESS | |
| 85 |
12/12✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 192 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 9 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 9 times.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 9 times.
✓ Branch 11 taken 3 times.
|
249 | end subroutine invariants_from_state |
| 86 | |||
| 87 | 2 | pure subroutine potato_to_simple_invariants(potato, psi_axis, psi_edge, & | |
| 88 | v0_ratio, simple_invariants) | ||
| 89 | !> Convert raw POTATO invariants to SIMPLE's portable convention. | ||
| 90 | !> | ||
| 91 | !> POTATO stores psi*=(c/q)P_phi in the equilibrium flux gauge. SIMPLE's | ||
| 92 | !> public invariant uses an axis-zero, outward-positive flux gauge. v0_ratio | ||
| 93 | !> is v0_POTATO/v0_SIMPLE; it is one when both runs use the same E_ref, | ||
| 94 | !> particle mass, and charge state. | ||
| 95 | type(guiding_center_invariants_t), intent(in) :: potato | ||
| 96 | real(dp), intent(in) :: psi_axis, psi_edge, v0_ratio | ||
| 97 | type(guiding_center_invariants_t), intent(out) :: simple_invariants | ||
| 98 | |||
| 99 | real(dp) :: direction, velocity_scale2 | ||
| 100 | |||
| 101 | 2 | direction = sign(1.0_dp, psi_edge - psi_axis) | |
| 102 | 2 | velocity_scale2 = v0_ratio**2 | |
| 103 | 2 | simple_invariants%h0 = velocity_scale2*potato%h0 | |
| 104 | 2 | simple_invariants%j_perp = velocity_scale2*potato%j_perp | |
| 105 | 2 | simple_invariants%p_phi = direction*(potato%p_phi - psi_axis) | |
| 106 | 2 | end subroutine potato_to_simple_invariants | |
| 107 | |||
| 108 | 1 | subroutine invariant_flux_convention(tracer, axis_flux, outward_sign) | |
| 109 | type(tracer_t), intent(in) :: tracer | ||
| 110 | real(dp), intent(out) :: axis_flux, outward_sign | ||
| 111 | |||
| 112 | 1 | call flux_convention(tracer, axis_flux, outward_sign) | |
| 113 | 1 | end subroutine invariant_flux_convention | |
| 114 | |||
| 115 | ✗ | subroutine states_from_invariants(tracer, invariants, result) | |
| 116 | type(tracer_t), intent(in) :: tracer | ||
| 117 | type(guiding_center_invariants_t), intent(in) :: invariants | ||
| 118 | type(invariant_start_result_t), intent(out) :: result | ||
| 119 | |||
| 120 | type(invariant_start_t) :: candidate | ||
| 121 | 1 | type(invariant_start_t), allocatable :: buffer(:) | |
| 122 | real(dp) :: residuals(N_RADIAL_SCAN, MAX_SECTION_BRANCHES) | ||
| 123 | real(dp) :: thetas(N_RADIAL_SCAN, MAX_SECTION_BRANCHES) | ||
| 124 | logical :: valid(N_RADIAL_SCAN, MAX_SECTION_BRANCHES) | ||
| 125 | integer :: kinds(N_RADIAL_SCAN, MAX_SECTION_BRANCHES) | ||
| 126 | real(dp) :: axis_flux, outward_sign, radius, radial_max, step_radius | ||
| 127 | integer :: counts(N_RADIAL_SCAN) | ||
| 128 | integer :: branch, branch_right, i, nfound, sigma, solve_status | ||
| 129 | |||
| 130 | 1 | result%status = INVARIANT_INVALID_INPUT | |
| 131 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (allocated(result%starts)) deallocate (result%starts) |
| 132 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
|
1 | allocate (result%starts(0)) |
| 133 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (.not. is_finite(invariants%h0)) return |
| 134 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (.not. is_finite(invariants%j_perp)) return |
| 135 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (.not. is_finite(invariants%p_phi)) return |
| 136 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (rho0 <= 0.0_dp) return |
| 137 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (invariants%h0 <= 0.0_dp) return |
| 138 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (invariants%j_perp < 0.0_dp) return |
| 139 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
|
1 | if (.not. field_is_axisymmetric(tracer)) then |
| 140 | ✗ | result%status = INVARIANT_NOT_AXISYMMETRIC | |
| 141 | ✗ | return | |
| 142 | end if | ||
| 143 | 1 | call flux_convention(tracer, axis_flux, outward_sign) | |
| 144 | |||
| 145 | 1 | residuals = 0.0_dp | |
| 146 | 1 | thetas = 0.0_dp | |
| 147 | 1 | valid = .false. | |
| 148 | 1 | kinds = SECTION_EXTREMUM_UNKNOWN | |
| 149 | 1 | counts = 0 | |
| 150 | 1 | radial_max = radial_upper_bound() | |
| 151 | 1 | step_radius = (radial_max - R_MIN)/real(N_RADIAL_SCAN - 1, dp) | |
| 152 | |||
| 153 |
2/2✓ Branch 0 taken 257 times.
✓ Branch 1 taken 1 time.
|
258 | do i = 1, N_RADIAL_SCAN |
| 154 | 257 | radius = R_MIN + real(i - 1, dp)*step_radius | |
| 155 | call section_geometry_at_radius(tracer, radius, thetas(i, :), & | ||
| 156 |
4/4✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 257 times.
✓ Branch 3 taken 4112 times.
✓ Branch 4 taken 257 times.
|
8482 | kinds(i, :), counts(i)) |
| 157 | end do | ||
| 158 | |||
| 159 |
5/6✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 8224 times.
✓ Branch 6 taken 1 time.
|
8230 | allocate (buffer(2*MAX_SECTION_BRANCHES*N_RADIAL_SCAN)) |
| 160 | 1 | nfound = 0 | |
| 161 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
3 | do sigma = -1, 1, 2 |
| 162 |
2/2✓ Branch 0 taken 514 times.
✓ Branch 1 taken 2 times.
|
516 | do i = 1, N_RADIAL_SCAN |
| 163 | 514 | radius = R_MIN + real(i - 1, dp)*step_radius | |
| 164 | call section_residuals_at_radius(tracer, invariants, axis_flux, & | ||
| 165 | outward_sign, radius, sigma, thetas(i, :), counts(i), & | ||
| 166 |
6/6✓ Branch 0 taken 8224 times.
✓ Branch 1 taken 514 times.
✓ Branch 3 taken 8224 times.
✓ Branch 4 taken 514 times.
✓ Branch 5 taken 8224 times.
✓ Branch 6 taken 514 times.
|
25188 | residuals(i, :), valid(i, :)) |
| 167 | end do | ||
| 168 | |||
| 169 |
2/2✓ Branch 0 taken 512 times.
✓ Branch 1 taken 2 times.
|
515 | do i = 1, N_RADIAL_SCAN - 1 |
| 170 |
2/2✓ Branch 0 taken 1024 times.
✓ Branch 1 taken 512 times.
|
1538 | do branch = 1, counts(i) |
| 171 |
2/2✓ Branch 0 taken 180 times.
✓ Branch 1 taken 844 times.
|
1024 | if (.not. valid(i, branch)) cycle |
| 172 | branch_right = matching_branch(thetas(i, branch), & | ||
| 173 | kinds(i, branch), thetas(i + 1, :), kinds(i + 1, :), & | ||
| 174 |
6/6✓ Branch 0 taken 13504 times.
✓ Branch 1 taken 844 times.
✓ Branch 2 taken 13504 times.
✓ Branch 3 taken 844 times.
✓ Branch 4 taken 13504 times.
✓ Branch 5 taken 844 times.
|
41356 | valid(i + 1, :), counts(i + 1)) |
| 175 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 842 times.
|
844 | if (branch_right == 0) cycle |
| 176 |
2/2✓ Branch 0 taken 840 times.
✓ Branch 1 taken 2 times.
|
842 | if (residuals(i, branch)*residuals(i + 1, branch_right) > & |
| 177 | 0.0_dp) cycle | ||
| 178 | |||
| 179 | call refine_invariant_root(tracer, invariants, axis_flux, & | ||
| 180 | outward_sign, sigma, branch, & | ||
| 181 | R_MIN + real(i - 1, dp)*step_radius, & | ||
| 182 | R_MIN + real(i, dp)*step_radius, thetas(i, branch), & | ||
| 183 | 2 | thetas(i + 1, branch_right), candidate, solve_status) | |
| 184 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (solve_status /= INVARIANT_SUCCESS) cycle |
| 185 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (is_duplicate(buffer, nfound, candidate)) cycle |
| 186 | 2 | nfound = nfound + 1 | |
| 187 | 1536 | buffer(nfound) = candidate | |
| 188 | end do | ||
| 189 | end do | ||
| 190 | end do | ||
| 191 | |||
| 192 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | deallocate (result%starts) |
| 193 |
10/16✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 time.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 time.
✓ Branch 14 taken 5 times.
✓ Branch 15 taken 1 time.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 1 time.
|
10 | allocate (result%starts(nfound)) |
| 194 |
4/6✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1 time.
|
4 | if (nfound > 0) result%starts = buffer(1:nfound) |
| 195 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (nfound > 0) then |
| 196 | 1 | result%status = INVARIANT_SUCCESS | |
| 197 | else | ||
| 198 | ✗ | result%status = INVARIANT_NO_INTERSECTION | |
| 199 | end if | ||
| 200 |
4/6✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 6 not taken.
|
8 | end subroutine states_from_invariants |
| 201 | |||
| 202 | 4 | logical function field_is_axisymmetric(tracer) | |
| 203 | type(tracer_t), intent(in) :: tracer | ||
| 204 | |||
| 205 | real(dp), parameter :: RADII(3) = [0.19_dp, 0.51_dp, 0.83_dp] | ||
| 206 | real(dp), parameter :: THETAS(3) = [0.17_dp, 1.31_dp, 4.73_dp] | ||
| 207 | real(dp), parameter :: PHI_TEST = 0.731_dp | ||
| 208 | type(tracer_t) :: left, right | ||
| 209 | real(dp) :: scale | ||
| 210 | integer :: i | ||
| 211 | |||
| 212 | 4 | field_is_axisymmetric = .false. | |
| 213 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
|
16 | do i = 1, size(RADII) |
| 214 | 12 | left = tracer | |
| 215 | 12 | right = tracer | |
| 216 | call evaluate_field(left%f, RADII(i)*radial_upper_bound(), & | ||
| 217 | 12 | THETAS(i), 0.0_dp, 0) | |
| 218 | call evaluate_field(right%f, RADII(i)*radial_upper_bound(), & | ||
| 219 | 12 | THETAS(i), PHI_TEST, 0) | |
| 220 | 12 | scale = max(1.0_dp, abs(left%f%Bmod), abs(right%f%Bmod)) | |
| 221 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(left%f%Bmod - right%f%Bmod) > AXISYM_TOL*scale) return |
| 222 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(left%f%dBmod(3)) > AXISYM_TOL*scale) return |
| 223 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(right%f%dBmod(3)) > AXISYM_TOL*scale) return |
| 224 | 12 | scale = max(1.0_dp, abs(left%f%Aph), abs(right%f%Aph)) | |
| 225 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(left%f%Aph - right%f%Aph) > AXISYM_TOL*scale) return |
| 226 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(left%f%dAph(3)) > AXISYM_TOL*scale) return |
| 227 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(right%f%dAph(3)) > AXISYM_TOL*scale) return |
| 228 | 12 | scale = max(1.0_dp, abs(left%f%hph), abs(right%f%hph)) | |
| 229 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(left%f%hph - right%f%hph) > AXISYM_TOL*scale) return |
| 230 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | if (abs(left%f%dhph(3)) > AXISYM_TOL*scale) return |
| 231 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
16 | if (abs(right%f%dhph(3)) > AXISYM_TOL*scale) return |
| 232 | end do | ||
| 233 | 4 | field_is_axisymmetric = .true. | |
| 234 |
24/24✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 256 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 12 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 12 times.
✓ Branch 9 taken 4 times.
✓ Branch 10 taken 12 times.
✓ Branch 11 taken 4 times.
✓ Branch 12 taken 12 times.
✓ Branch 13 taken 4 times.
✓ Branch 14 taken 24 times.
✓ Branch 15 taken 4 times.
✓ Branch 16 taken 256 times.
✓ Branch 17 taken 4 times.
✓ Branch 18 taken 12 times.
✓ Branch 19 taken 4 times.
✓ Branch 20 taken 12 times.
✓ Branch 21 taken 4 times.
✓ Branch 22 taken 12 times.
✓ Branch 23 taken 4 times.
|
668 | end function field_is_axisymmetric |
| 235 | |||
| 236 | 257 | subroutine section_geometry_at_radius(tracer, radius, thetas, kinds, count) | |
| 237 | type(tracer_t), intent(in) :: tracer | ||
| 238 | real(dp), intent(in) :: radius | ||
| 239 | real(dp), intent(out) :: thetas(MAX_SECTION_BRANCHES) | ||
| 240 | integer, intent(out) :: kinds(MAX_SECTION_BRANCHES) | ||
| 241 | integer, intent(out) :: count | ||
| 242 | |||
| 243 | integer :: branch | ||
| 244 | |||
| 245 | 257 | call find_section_extrema(tracer, radius, thetas, count) | |
| 246 | 257 | kinds = SECTION_EXTREMUM_UNKNOWN | |
| 247 |
2/2✓ Branch 0 taken 514 times.
✓ Branch 1 taken 257 times.
|
771 | do branch = 1, count |
| 248 | 771 | kinds(branch) = extremum_kind(tracer, radius, thetas(branch)) | |
| 249 | end do | ||
| 250 | 257 | end subroutine section_geometry_at_radius | |
| 251 | |||
| 252 | 514 | subroutine section_residuals_at_radius(tracer, invariants, axis_flux, & | |
| 253 | outward_sign, radius, sigma, & | ||
| 254 | thetas, count, residuals, valid) | ||
| 255 | type(tracer_t), intent(in) :: tracer | ||
| 256 | type(guiding_center_invariants_t), intent(in) :: invariants | ||
| 257 | real(dp), intent(in) :: axis_flux, outward_sign, radius | ||
| 258 | integer, intent(in) :: sigma | ||
| 259 | real(dp), intent(in) :: thetas(MAX_SECTION_BRANCHES) | ||
| 260 | integer, intent(in) :: count | ||
| 261 | real(dp), intent(out) :: residuals(MAX_SECTION_BRANCHES) | ||
| 262 | logical, intent(out) :: valid(MAX_SECTION_BRANCHES) | ||
| 263 | |||
| 264 | type(tracer_t) :: work | ||
| 265 | real(dp) :: lambda2 | ||
| 266 | integer :: branch | ||
| 267 | |||
| 268 | 514 | residuals = 0.0_dp | |
| 269 | 514 | valid = .false. | |
| 270 |
2/2✓ Branch 0 taken 1028 times.
✓ Branch 1 taken 514 times.
|
1542 | do branch = 1, count |
| 271 | 1028 | work = tracer | |
| 272 | 1028 | call evaluate_field(work%f, radius, thetas(branch), 0.0_dp, 0) | |
| 273 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1028 times.
|
1028 | if (work%f%Bmod <= 0.0_dp) cycle |
| 274 | 1028 | lambda2 = 1.0_dp - invariants%j_perp*work%f%Bmod/invariants%h0 | |
| 275 |
2/2✓ Branch 0 taken 182 times.
✓ Branch 1 taken 846 times.
|
1028 | if (lambda2 < 0.0_dp) cycle |
| 276 | residuals(branch) = outward_sign*(work%f%Aph - axis_flux + & | ||
| 277 | rho0*sqrt(invariants%h0)*real(sigma, dp)* & | ||
| 278 | 846 | sqrt(max(0.0_dp, lambda2))*work%f%hph) - invariants%p_phi | |
| 279 | 1542 | valid(branch) = .true. | |
| 280 | end do | ||
| 281 |
12/12✓ Branch 0 taken 1542 times.
✓ Branch 1 taken 514 times.
✓ Branch 2 taken 3084 times.
✓ Branch 3 taken 514 times.
✓ Branch 4 taken 32896 times.
✓ Branch 5 taken 514 times.
✓ Branch 6 taken 1542 times.
✓ Branch 7 taken 514 times.
✓ Branch 8 taken 1542 times.
✓ Branch 9 taken 514 times.
✓ Branch 10 taken 1542 times.
✓ Branch 11 taken 514 times.
|
42662 | end subroutine section_residuals_at_radius |
| 282 | |||
| 283 | 844 | integer function matching_branch(theta, kind, candidates, candidate_kinds, & | |
| 284 | candidate_valid, count) | ||
| 285 | real(dp), intent(in) :: theta | ||
| 286 | integer, intent(in) :: kind, count | ||
| 287 | real(dp), intent(in) :: candidates(MAX_SECTION_BRANCHES) | ||
| 288 | integer, intent(in) :: candidate_kinds(MAX_SECTION_BRANCHES) | ||
| 289 | logical, intent(in) :: candidate_valid(MAX_SECTION_BRANCHES) | ||
| 290 | |||
| 291 | real(dp) :: distance, best_distance | ||
| 292 | integer :: i | ||
| 293 | |||
| 294 | 844 | matching_branch = 0 | |
| 295 | 844 | best_distance = huge(1.0_dp) | |
| 296 |
2/2✓ Branch 0 taken 1688 times.
✓ Branch 1 taken 844 times.
|
2532 | do i = 1, count |
| 297 |
2/2✓ Branch 0 taken 184 times.
✓ Branch 1 taken 1504 times.
|
1688 | if (.not. candidate_valid(i)) cycle |
| 298 |
2/2✓ Branch 0 taken 662 times.
✓ Branch 1 taken 842 times.
|
1504 | if (candidate_kinds(i) /= kind) cycle |
| 299 | 842 | distance = angular_distance(candidates(i), theta) | |
| 300 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 842 times.
|
842 | if (distance >= best_distance) cycle |
| 301 | 1688 | matching_branch = i | |
| 302 | 3378 | best_distance = distance | |
| 303 | end do | ||
| 304 | 844 | end function matching_branch | |
| 305 | |||
| 306 | 303 | subroutine find_section_extrema(tracer, radius, roots, count) | |
| 307 | type(tracer_t), intent(in) :: tracer | ||
| 308 | real(dp), intent(in) :: radius | ||
| 309 | real(dp), intent(out) :: roots(MAX_SECTION_BRANCHES) | ||
| 310 | integer, intent(out) :: count | ||
| 311 | |||
| 312 | real(dp) :: dtheta, left, right, fleft, fright, root | ||
| 313 | integer :: i | ||
| 314 | |||
| 315 | 303 | roots = 0.0_dp | |
| 316 | 303 | count = 0 | |
| 317 | 303 | dtheta = twopi/real(N_THETA_SCAN, dp) | |
| 318 | 303 | left = 0.5_dp*dtheta | |
| 319 | 303 | fleft = theta_derivative(tracer, radius, left) | |
| 320 |
2/2✓ Branch 1 taken 38784 times.
✓ Branch 2 taken 303 times.
|
39087 | do i = 0, N_THETA_SCAN - 1 |
| 321 | 38784 | left = (real(i, dp) + 0.5_dp)*dtheta | |
| 322 | 38784 | right = left + dtheta | |
| 323 | 38784 | fright = theta_derivative(tracer, radius, right) | |
| 324 |
2/2✓ Branch 0 taken 606 times.
✓ Branch 1 taken 38178 times.
|
38784 | if (fleft*fright <= 0.0_dp) then |
| 325 | 606 | call bisect_theta_extremum(tracer, radius, left, right, root) | |
| 326 |
3/4✓ Branch 0 taken 303 times.
✓ Branch 1 taken 303 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 303 times.
|
606 | root = modulo(root, twopi) |
| 327 | 606 | if (.not. root_is_duplicate(roots, count, root)) then | |
| 328 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 606 times.
|
606 | if (count == MAX_SECTION_BRANCHES) exit |
| 329 | 606 | count = count + 1 | |
| 330 | 606 | roots(count) = root | |
| 331 | end if | ||
| 332 | end if | ||
| 333 | 39087 | fleft = fright | |
| 334 | end do | ||
| 335 | 303 | call sort_roots(roots, count) | |
| 336 | 303 | end subroutine find_section_extrema | |
| 337 | |||
| 338 | 606 | subroutine bisect_theta_extremum(tracer, radius, left_in, right_in, root) | |
| 339 | type(tracer_t), intent(in) :: tracer | ||
| 340 | real(dp), intent(in) :: radius, left_in, right_in | ||
| 341 | real(dp), intent(out) :: root | ||
| 342 | |||
| 343 | real(dp) :: left, right, middle, fleft, fmiddle | ||
| 344 | integer :: iteration | ||
| 345 | |||
| 346 | 606 | left = left_in | |
| 347 | 606 | right = right_in | |
| 348 | 606 | fleft = theta_derivative(tracer, radius, left) | |
| 349 |
1/2✓ Branch 1 taken 606 times.
✗ Branch 2 not taken.
|
606 | do iteration = 1, MAX_ROOT_ITERATIONS |
| 350 | 606 | middle = 0.5_dp*(left + right) | |
| 351 | 606 | fmiddle = theta_derivative(tracer, radius, middle) | |
| 352 |
1/2✓ Branch 0 taken 606 times.
✗ Branch 1 not taken.
|
606 | if (abs(fmiddle) <= ROOT_TOL) exit |
| 353 | ✗ | if (right - left <= ROOT_TOL) exit | |
| 354 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
606 | if (fleft*fmiddle <= 0.0_dp) then |
| 355 | right = middle | ||
| 356 | else | ||
| 357 | ✗ | left = middle | |
| 358 | ✗ | fleft = fmiddle | |
| 359 | end if | ||
| 360 | end do | ||
| 361 | 606 | root = 0.5_dp*(left + right) | |
| 362 | 606 | end subroutine bisect_theta_extremum | |
| 363 | |||
| 364 | 41419 | real(dp) function theta_derivative(tracer, radius, theta) | |
| 365 | type(tracer_t), intent(in) :: tracer | ||
| 366 | real(dp), intent(in) :: radius, theta | ||
| 367 | |||
| 368 | type(tracer_t) :: work | ||
| 369 | |||
| 370 | 41419 | work = tracer | |
| 371 | 41419 | call evaluate_field(work%f, radius, theta, 0.0_dp, 0) | |
| 372 | 41419 | theta_derivative = work%f%dBmod(2) | |
| 373 |
12/12✓ Branch 0 taken 41419 times.
✓ Branch 1 taken 124257 times.
✓ Branch 2 taken 41419 times.
✓ Branch 3 taken 248514 times.
✓ Branch 4 taken 2650816 times.
✓ Branch 5 taken 41419 times.
✓ Branch 6 taken 41419 times.
✓ Branch 7 taken 124257 times.
✓ Branch 8 taken 41419 times.
✓ Branch 9 taken 124257 times.
✓ Branch 10 taken 124257 times.
✓ Branch 11 taken 41419 times.
|
3644872 | end function theta_derivative |
| 374 | |||
| 375 | 2 | subroutine refine_invariant_root(tracer, invariants, axis_flux, outward_sign, & | |
| 376 | sigma, branch, & | ||
| 377 | left_in, right_in, theta_left_in, & | ||
| 378 | theta_right_in, solution, status) | ||
| 379 | type(tracer_t), intent(in) :: tracer | ||
| 380 | type(guiding_center_invariants_t), intent(in) :: invariants | ||
| 381 | real(dp), intent(in) :: axis_flux, outward_sign | ||
| 382 | integer, intent(in) :: sigma, branch | ||
| 383 | real(dp), intent(in) :: left_in, right_in, theta_left_in, theta_right_in | ||
| 384 | type(invariant_start_t), intent(out) :: solution | ||
| 385 | integer, intent(out) :: status | ||
| 386 | |||
| 387 | real(dp) :: left, right, middle, fleft, fmiddle | ||
| 388 | real(dp) :: theta_left, theta_right, theta_middle, theta_guess | ||
| 389 | real(dp) :: z_integ(5), z_ref(3) | ||
| 390 | integer :: iteration, kind | ||
| 391 | logical :: ok | ||
| 392 | |||
| 393 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 2 times.
|
12 | solution = invariant_start_t() |
| 394 | 2 | status = INVARIANT_ROOT_FAILURE | |
| 395 | 2 | left = left_in | |
| 396 | 2 | right = right_in | |
| 397 | 2 | theta_left = theta_left_in | |
| 398 | 2 | theta_right = theta_right_in | |
| 399 | 2 | theta_guess = theta_left | |
| 400 | call residual_near_theta(tracer, invariants, axis_flux, outward_sign, & | ||
| 401 | 2 | left, theta_guess, sigma, fleft, theta_left, kind, ok) | |
| 402 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (.not. ok) return |
| 403 | |||
| 404 |
1/2✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
|
44 | do iteration = 1, MAX_ROOT_ITERATIONS |
| 405 | 44 | middle = 0.5_dp*(left + right) | |
| 406 | 44 | theta_middle = angular_midpoint(theta_left, theta_right) | |
| 407 | 44 | theta_guess = theta_middle | |
| 408 | call residual_near_theta(tracer, invariants, axis_flux, outward_sign, & | ||
| 409 | 44 | middle, theta_guess, sigma, fmiddle, theta_middle, kind, ok) | |
| 410 |
1/2✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
|
44 | if (.not. ok) return |
| 411 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 42 times.
|
44 | if (abs(fmiddle) <= ROOT_TOL*max(1.0_dp, abs(invariants%p_phi))) exit |
| 412 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if (right - left <= ROOT_TOL) exit |
| 413 |
2/2✓ Branch 1 taken 20 times.
✓ Branch 2 taken 22 times.
|
88 | if (fleft*fmiddle <= 0.0_dp) then |
| 414 | 20 | right = middle | |
| 415 | 20 | theta_right = theta_middle | |
| 416 | else | ||
| 417 | 22 | left = middle | |
| 418 | 22 | theta_left = theta_middle | |
| 419 | 22 | fleft = fmiddle | |
| 420 | end if | ||
| 421 | end do | ||
| 422 | |||
| 423 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | z_integ(1:3) = [middle, theta_middle, 0.0_dp] |
| 424 | 2 | z_integ(4) = sqrt(invariants%h0) | |
| 425 | call pitch_at_position(tracer, invariants, z_integ(1:3), sigma, & | ||
| 426 | 2 | z_integ(5), ok) | |
| 427 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (.not. ok) return |
| 428 | 2 | call integ_to_ref(z_integ(1:3), z_ref) | |
| 429 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | solution%state(1:3) = z_ref |
| 430 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | solution%state(4:5) = z_integ(4:5) |
| 431 | 2 | solution%sigma = sigma | |
| 432 | 2 | solution%section_branch = branch | |
| 433 | 2 | solution%section_kind = kind | |
| 434 | 2 | solution%residual = fmiddle | |
| 435 | 2 | status = INVARIANT_SUCCESS | |
| 436 | end subroutine refine_invariant_root | ||
| 437 | |||
| 438 | 46 | subroutine residual_near_theta(tracer, invariants, axis_flux, outward_sign, & | |
| 439 | radius, theta_guess, sigma, & | ||
| 440 | residual, theta, kind, ok) | ||
| 441 | type(tracer_t), intent(in) :: tracer | ||
| 442 | type(guiding_center_invariants_t), intent(in) :: invariants | ||
| 443 | real(dp), intent(in) :: axis_flux, outward_sign, radius, theta_guess | ||
| 444 | integer, intent(in) :: sigma | ||
| 445 | real(dp), intent(out) :: residual, theta | ||
| 446 | integer, intent(out) :: kind | ||
| 447 | logical, intent(out) :: ok | ||
| 448 | |||
| 449 | type(tracer_t) :: work | ||
| 450 | real(dp) :: roots(MAX_SECTION_BRANCHES), distance, best_distance | ||
| 451 | real(dp) :: lambda2 | ||
| 452 | integer :: count, i, best | ||
| 453 | |||
| 454 | 46 | residual = 0.0_dp | |
| 455 | 46 | theta = 0.0_dp | |
| 456 | 46 | kind = SECTION_EXTREMUM_UNKNOWN | |
| 457 | 46 | ok = .false. | |
| 458 | 46 | call find_section_extrema(tracer, radius, roots, count) | |
| 459 |
1/2✓ Branch 0 taken 46 times.
✗ Branch 1 not taken.
|
46 | if (count == 0) return |
| 460 | 46 | best = 1 | |
| 461 | 46 | best_distance = angular_distance(roots(1), theta_guess) | |
| 462 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 46 times.
|
92 | do i = 2, count |
| 463 | 46 | distance = angular_distance(roots(i), theta_guess) | |
| 464 |
2/2✓ Branch 0 taken 23 times.
✓ Branch 1 taken 23 times.
|
92 | if (distance < best_distance) then |
| 465 | 23 | best = i | |
| 466 | 23 | best_distance = distance | |
| 467 | end if | ||
| 468 | end do | ||
| 469 | 46 | theta = roots(best) | |
| 470 | 46 | work = tracer | |
| 471 | 46 | call evaluate_field(work%f, radius, theta, 0.0_dp, 0) | |
| 472 | 46 | lambda2 = 1.0_dp - invariants%j_perp*work%f%Bmod/invariants%h0 | |
| 473 |
1/2✓ Branch 0 taken 46 times.
✗ Branch 1 not taken.
|
46 | if (lambda2 < 0.0_dp) return |
| 474 | residual = outward_sign*(work%f%Aph - axis_flux + & | ||
| 475 | rho0*sqrt(invariants%h0)*real(sigma, dp)*sqrt(lambda2)*work%f%hph) - & | ||
| 476 | 46 | invariants%p_phi | |
| 477 | 46 | kind = extremum_kind(tracer, radius, theta) | |
| 478 | 46 | ok = .true. | |
| 479 |
12/12✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 276 times.
✓ Branch 3 taken 46 times.
✓ Branch 4 taken 2944 times.
✓ Branch 5 taken 46 times.
✓ Branch 6 taken 138 times.
✓ Branch 7 taken 46 times.
✓ Branch 8 taken 138 times.
✓ Branch 9 taken 46 times.
✓ Branch 10 taken 138 times.
✓ Branch 11 taken 46 times.
|
3818 | end subroutine residual_near_theta |
| 480 | |||
| 481 | 560 | integer function extremum_kind(tracer, radius, theta) | |
| 482 | type(tracer_t), intent(in) :: tracer | ||
| 483 | real(dp), intent(in) :: radius, theta | ||
| 484 | |||
| 485 | real(dp), parameter :: ANGLE_STEP = 1.0e-5_dp | ||
| 486 | real(dp) :: slope | ||
| 487 | |||
| 488 | slope = theta_derivative(tracer, radius, theta + ANGLE_STEP) - & | ||
| 489 | 560 | theta_derivative(tracer, radius, theta - ANGLE_STEP) | |
| 490 | 560 | extremum_kind = SECTION_EXTREMUM_UNKNOWN | |
| 491 |
2/2✓ Branch 0 taken 280 times.
✓ Branch 1 taken 280 times.
|
560 | if (slope > 0.0_dp) extremum_kind = SECTION_B_MINIMUM |
| 492 |
2/2✓ Branch 0 taken 280 times.
✓ Branch 1 taken 280 times.
|
560 | if (slope < 0.0_dp) extremum_kind = SECTION_B_MAXIMUM |
| 493 | 560 | end function extremum_kind | |
| 494 | |||
| 495 | 2 | subroutine pitch_at_position(tracer, invariants, position, sigma, pitch, ok) | |
| 496 | type(tracer_t), intent(in) :: tracer | ||
| 497 | type(guiding_center_invariants_t), intent(in) :: invariants | ||
| 498 | real(dp), intent(in) :: position(3) | ||
| 499 | integer, intent(in) :: sigma | ||
| 500 | real(dp), intent(out) :: pitch | ||
| 501 | logical, intent(out) :: ok | ||
| 502 | |||
| 503 | type(tracer_t) :: work | ||
| 504 | real(dp) :: lambda2 | ||
| 505 | |||
| 506 | 2 | work = tracer | |
| 507 | 2 | call evaluate_field(work%f, position(1), position(2), position(3), 0) | |
| 508 | 2 | lambda2 = 1.0_dp - invariants%j_perp*work%f%Bmod/invariants%h0 | |
| 509 | 2 | ok = lambda2 >= 0.0_dp | |
| 510 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (ok) then |
| 511 | 2 | pitch = real(sigma, dp)*sqrt(max(0.0_dp, lambda2)) | |
| 512 | else | ||
| 513 | ✗ | pitch = 0.0_dp | |
| 514 | end if | ||
| 515 |
12/12✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 128 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 6 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 6 times.
✓ Branch 11 taken 2 times.
|
166 | end subroutine pitch_at_position |
| 516 | |||
| 517 | 5 | subroutine flux_convention(tracer, axis_flux, outward_sign) | |
| 518 | type(tracer_t), intent(in) :: tracer | ||
| 519 | real(dp), intent(out) :: axis_flux, outward_sign | ||
| 520 | |||
| 521 | type(tracer_t) :: axis_field, edge_field | ||
| 522 | real(dp) :: delta_flux | ||
| 523 | |||
| 524 | 5 | axis_field = tracer | |
| 525 | 5 | edge_field = tracer | |
| 526 | 5 | call evaluate_field(axis_field%f, R_MIN, 0.0_dp, 0.0_dp, 0) | |
| 527 | 5 | call evaluate_field(edge_field%f, radial_upper_bound(), 0.0_dp, 0.0_dp, 0) | |
| 528 | 5 | axis_flux = axis_field%f%Aph | |
| 529 | 5 | delta_flux = edge_field%f%Aph - axis_flux | |
| 530 | 5 | outward_sign = sign(1.0_dp, delta_flux) | |
| 531 |
24/24✓ Branch 0 taken 15 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 30 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 320 times.
✓ Branch 5 taken 5 times.
✓ Branch 6 taken 15 times.
✓ Branch 7 taken 5 times.
✓ Branch 8 taken 15 times.
✓ Branch 9 taken 5 times.
✓ Branch 10 taken 15 times.
✓ Branch 11 taken 5 times.
✓ Branch 12 taken 15 times.
✓ Branch 13 taken 5 times.
✓ Branch 14 taken 30 times.
✓ Branch 15 taken 5 times.
✓ Branch 16 taken 320 times.
✓ Branch 17 taken 5 times.
✓ Branch 18 taken 15 times.
✓ Branch 19 taken 5 times.
✓ Branch 20 taken 15 times.
✓ Branch 21 taken 5 times.
✓ Branch 22 taken 15 times.
✓ Branch 23 taken 5 times.
|
825 | end subroutine flux_convention |
| 532 | |||
| 533 | 30 | pure real(dp) function radial_upper_bound() | |
| 534 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5 times.
|
30 | if (isw_field_type == TEST) then |
| 535 | radial_upper_bound = TEST_RADIAL_MAX | ||
| 536 | else | ||
| 537 | ✗ | radial_upper_bound = UNIT_RADIAL_MAX | |
| 538 | end if | ||
| 539 | end function radial_upper_bound | ||
| 540 | |||
| 541 | 606 | logical function root_is_duplicate(roots, count, root) | |
| 542 | real(dp), intent(in) :: roots(MAX_SECTION_BRANCHES), root | ||
| 543 | integer, intent(in) :: count | ||
| 544 | |||
| 545 | integer :: i | ||
| 546 | |||
| 547 | 606 | root_is_duplicate = .false. | |
| 548 |
2/2✓ Branch 0 taken 303 times.
✓ Branch 1 taken 606 times.
|
909 | do i = 1, count |
| 549 |
1/2✓ Branch 0 taken 303 times.
✗ Branch 1 not taken.
|
909 | if (angular_distance(roots(i), root) < 1.0e-8_dp) then |
| 550 |
1/2✓ Branch 0 taken 606 times.
✗ Branch 1 not taken.
|
606 | root_is_duplicate = .true. |
| 551 | return | ||
| 552 | end if | ||
| 553 | end do | ||
| 554 | end function root_is_duplicate | ||
| 555 | |||
| 556 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | logical function is_duplicate(buffer, count, candidate) |
| 557 | type(invariant_start_t), intent(in) :: buffer(:), candidate | ||
| 558 | integer, intent(in) :: count | ||
| 559 | |||
| 560 | integer :: i | ||
| 561 | |||
| 562 | 2 | is_duplicate = .false. | |
| 563 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
3 | do i = 1, count |
| 564 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (buffer(i)%sigma /= candidate%sigma) cycle |
| 565 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (abs(buffer(i)%state(1) - candidate%state(1)) > 1.0e-7_dp) cycle |
| 566 | ✗ | if (angular_distance(buffer(i)%state(2), candidate%state(2)) > & | |
| 567 | 1.0e-7_dp) cycle | ||
| 568 | 2 | is_duplicate = .true. | |
| 569 | 3 | return | |
| 570 | end do | ||
| 571 | end function is_duplicate | ||
| 572 | |||
| 573 | 303 | subroutine sort_roots(roots, count) | |
| 574 | real(dp), intent(inout) :: roots(MAX_SECTION_BRANCHES) | ||
| 575 | integer, intent(in) :: count | ||
| 576 | |||
| 577 | real(dp) :: value | ||
| 578 | integer :: i, j | ||
| 579 | |||
| 580 |
2/2✓ Branch 0 taken 303 times.
✓ Branch 1 taken 303 times.
|
606 | do i = 2, count |
| 581 | 303 | value = roots(i) | |
| 582 | 303 | j = i - 1 | |
| 583 |
2/2✓ Branch 0 taken 303 times.
✓ Branch 1 taken 303 times.
|
606 | do while (j >= 1) |
| 584 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 303 times.
|
303 | if (roots(j) <= value) exit |
| 585 | 303 | roots(j + 1) = roots(j) | |
| 586 | 303 | j = j - 1 | |
| 587 | end do | ||
| 588 | 606 | roots(j + 1) = value | |
| 589 | end do | ||
| 590 | 303 | end subroutine sort_roots | |
| 591 | |||
| 592 | 1237 | pure real(dp) function angular_distance(left, right) | |
| 593 | real(dp), intent(in) :: left, right | ||
| 594 | |||
| 595 | angular_distance = abs(modulo(left - right + 0.5_dp*twopi, twopi) - & | ||
| 596 |
3/4✓ Branch 0 taken 888 times.
✓ Branch 1 taken 349 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 888 times.
|
1237 | 0.5_dp*twopi) |
| 597 | 1237 | end function angular_distance | |
| 598 | |||
| 599 | 44 | pure real(dp) function angular_midpoint(left, right) | |
| 600 | real(dp), intent(in) :: left, right | ||
| 601 | |||
| 602 | angular_midpoint = modulo(left + 0.5_dp* & | ||
| 603 |
5/8✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 44 times.
✓ Branch 4 taken 22 times.
✓ Branch 5 taken 22 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 22 times.
|
44 | (modulo(right - left + 0.5_dp*twopi, twopi) - 0.5_dp*twopi), twopi) |
| 604 | 44 | end function angular_midpoint | |
| 605 | |||
| 606 | ✗ | end module orbit_invariants | |
| 607 |