| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module collis_alp | ||
| 2 | use, intrinsic :: iso_fortran_env, only: wp => real64 | ||
| 3 | implicit none | ||
| 4 | |||
| 5 | integer, parameter :: nsorts = 3 | ||
| 6 | integer, parameter :: N_S_GRID = 101 | ||
| 7 | |||
| 8 | real(wp), dimension(nsorts) :: efcolf, velrat, enrat | ||
| 9 | real(wp), dimension(nsorts, N_S_GRID) :: efcolf_grid, velrat_grid, enrat_grid | ||
| 10 | real(wp) :: lorentz_nu = 0.0_wp | ||
| 11 | real(wp) :: lorentz_dhh = 0.0_wp | ||
| 12 | |||
| 13 | contains | ||
| 14 | |||
| 15 | 1 | subroutine configure_lorentz_nustar(nu_star, major_radius_cm, iota, vref) | |
| 16 | real(wp), intent(in) :: nu_star, major_radius_cm, iota, vref | ||
| 17 | |||
| 18 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (nu_star < 0.0_wp) error stop 'nu_star_standard must be nonnegative' |
| 19 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (major_radius_cm <= 0.0_wp) & |
| 20 | ✗ | error stop 'lorentz_major_radius_cm must be positive' | |
| 21 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (abs(iota) <= tiny(1.0_wp)) error stop 'lorentz_iota must be nonzero' |
| 22 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (vref <= 0.0_wp) error stop 'Lorentz reference velocity must be positive' |
| 23 | |||
| 24 | 1 | lorentz_nu = nu_star*abs(iota)*vref/major_radius_cm | |
| 25 | ! stost advances in path length d tau = vref dt. This coefficient gives | ||
| 26 | ! the Lorentz generator (nu/2) d/dlambda[(1-lambda^2)d/dlambda]. | ||
| 27 | 1 | lorentz_dhh = lorentz_nu/(2.0_wp*vref) | |
| 28 | 1 | end subroutine configure_lorentz_nustar | |
| 29 | |||
| 30 | 8000000 | subroutine stost_lorentz(z, dtauc, ierr) | |
| 31 | real(wp), intent(inout) :: z(5) | ||
| 32 | real(wp), intent(in) :: dtauc | ||
| 33 | integer, intent(out) :: ierr | ||
| 34 | |||
| 35 | real(wp) :: alam, coala, dalam | ||
| 36 | real :: ur | ||
| 37 | |||
| 38 | 8000000 | ierr = 0 | |
| 39 | 8000000 | alam = z(5) | |
| 40 | 8000000 | coala = 1.0_wp - alam**2 | |
| 41 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8000000 times.
|
8000000 | if (coala < 0.0_wp) then |
| 42 | ✗ | ierr = 1 | |
| 43 | ✗ | return | |
| 44 | end if | ||
| 45 | |||
| 46 | 8000000 | call getran(1, ur) | |
| 47 | dalam = sqrt(2.0_wp*lorentz_dhh*coala*dtauc)*real(ur, wp) & | ||
| 48 | 8000000 | - 2.0_wp*alam*lorentz_dhh*dtauc | |
| 49 | 8000000 | alam = alam + dalam | |
| 50 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8000000 times.
|
8000000 | if (alam > 1.0_wp) then |
| 51 | ✗ | alam = 2.0_wp - alam | |
| 52 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8000000 times.
|
8000000 | else if (alam < -1.0_wp) then |
| 53 | ✗ | alam = -2.0_wp - alam | |
| 54 | end if | ||
| 55 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8000000 times.
|
8000000 | if (abs(alam) > 1.0_wp) then |
| 56 | ✗ | ierr = 2 | |
| 57 | ✗ | return | |
| 58 | end if | ||
| 59 | 8000000 | z(5) = alam | |
| 60 | end subroutine stost_lorentz | ||
| 61 | |||
| 62 | 229 | subroutine coleff(p, dpp, dhh, fpeff) | |
| 63 | real(wp), intent(in) :: p | ||
| 64 | real(wp), intent(out) :: dpp, dhh, fpeff | ||
| 65 | |||
| 66 | 229 | call coleff_local(p, efcolf, velrat, enrat, dpp, dhh, fpeff) | |
| 67 | 229 | end subroutine coleff | |
| 68 | |||
| 69 | 142483808 | subroutine coleff_local(p, efcolf_loc, velrat_loc, enrat_loc, dpp, dhh, fpeff) | |
| 70 | real(wp), intent(in) :: p | ||
| 71 | real(wp), intent(in) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) | ||
| 72 | real(wp), intent(out) :: dpp, dhh, fpeff | ||
| 73 | |||
| 74 | integer :: i | ||
| 75 | real(wp) :: plim, xbeta, dp, dh, dpd | ||
| 76 | |||
| 77 | 142483808 | plim = max(p, 1.d-8) | |
| 78 | |||
| 79 | 142483808 | dpp = 0.0d0 | |
| 80 | 142483808 | dhh = 0.0d0 | |
| 81 | 142483808 | fpeff = 0.0d0 | |
| 82 | |||
| 83 |
2/2✓ Branch 0 taken 427451424 times.
✓ Branch 1 taken 142483808 times.
|
569935232 | do i = 1, nsorts |
| 84 | 427451424 | xbeta = p*velrat_loc(i) | |
| 85 | 427451424 | call onseff(xbeta, dp, dh, dpd) | |
| 86 | 427451424 | dpp = dpp + dp*efcolf_loc(i) | |
| 87 | 427451424 | dhh = dhh + dh*efcolf_loc(i) | |
| 88 | 569935232 | fpeff = fpeff + (dpd/plim - 2.0d0*dp*p*enrat_loc(i))*efcolf_loc(i) | |
| 89 | end do | ||
| 90 | |||
| 91 | 142483808 | dhh = dhh/plim**2 | |
| 92 | 142483808 | end subroutine coleff_local | |
| 93 | |||
| 94 | 427451424 | subroutine onseff(v, dp, dh, dpd) | |
| 95 | real(wp), intent(in) :: v | ||
| 96 | real(wp), intent(out) :: dp, dh, dpd | ||
| 97 | |||
| 98 | real(wp), parameter :: sqp = 1.7724538d0 | ||
| 99 | real(wp), parameter :: cons = 0.75225278d0 | ||
| 100 | real(wp) :: v2, v3, ex, er | ||
| 101 | |||
| 102 | 427451424 | v2 = v**2 | |
| 103 | 427451424 | v3 = v2*v | |
| 104 |
2/2✓ Branch 0 taken 291348640 times.
✓ Branch 1 taken 136102784 times.
|
427451424 | if (v < 0.01d0) then |
| 105 | 291348640 | dp = cons*(1.d0 - 0.6d0*v2) | |
| 106 | 291348640 | dh = cons*(1.d0 - 0.2d0*v2) | |
| 107 | 291348640 | dpd = 2.d0*cons*(1.d0 - 1.2d0*v2) | |
| 108 |
2/2✓ Branch 0 taken 51109 times.
✓ Branch 1 taken 136051675 times.
|
136102784 | else if (v > 6.d0) then |
| 109 | 51109 | dp = 1.d0/v3 | |
| 110 | 51109 | dh = (1.d0 - 0.5d0/v2)/v | |
| 111 | 51109 | dpd = -1.d0/v3 | |
| 112 | else | ||
| 113 | 136051675 | ex = exp(-v2)/sqp | |
| 114 | 136051675 | er = erf(v) | |
| 115 | 136051675 | dp = er/v3 - 2.d0*ex/v2 | |
| 116 | 136051675 | dh = er*(1.d0 - 0.5d0/v2)/v + ex/v2 | |
| 117 | 136051675 | dpd = 4.d0*ex - dp | |
| 118 | end if | ||
| 119 | 427451424 | end subroutine onseff | |
| 120 | |||
| 121 | 8 | subroutine loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, & | |
| 122 | tempe, ealpha, v0, dchichi, slowrate, & | ||
| 123 | dchichi_norm, slowrate_norm) | ||
| 124 | real(wp), intent(in) :: am1, am2, Z1, Z2, densi1, densi2 | ||
| 125 | real(wp), intent(in) :: tempi1, tempi2, tempe, ealpha | ||
| 126 | real(wp), intent(out) :: v0, dchichi, slowrate, dchichi_norm, slowrate_norm | ||
| 127 | |||
| 128 | call compute_collision_coeffs(am1, am2, Z1, Z2, densi1, densi2, & | ||
| 129 | tempi1, tempi2, tempe, ealpha, v0, & | ||
| 130 | 8 | efcolf, velrat, enrat) | |
| 131 | |||
| 132 | 8 | call compute_rates(v0, dchichi, slowrate, dchichi_norm, slowrate_norm) | |
| 133 | 8 | end subroutine loacol_alpha | |
| 134 | |||
| 135 | 1725 | subroutine compute_collision_coeffs(am1, am2, Z1, Z2, densi1, densi2, & | |
| 136 | tempi1, tempi2, tempe, ealpha, v0, & | ||
| 137 | efcolf_out, velrat_out, enrat_out) | ||
| 138 | real(wp), intent(in) :: am1, am2, Z1, Z2, densi1, densi2 | ||
| 139 | real(wp), intent(in) :: tempi1, tempi2, tempe, ealpha | ||
| 140 | real(wp), intent(out) :: v0 | ||
| 141 | real(wp), intent(out) :: efcolf_out(nsorts), velrat_out(nsorts), enrat_out(nsorts) | ||
| 142 | |||
| 143 | real(wp), parameter :: pi = 3.14159265358979d0 | ||
| 144 | real(wp), parameter :: pmass = 1.6726d-24 | ||
| 145 | real(wp), parameter :: emass = 9.1094d-28 | ||
| 146 | real(wp), parameter :: e_cgs = 4.8032d-10 | ||
| 147 | real(wp), parameter :: ev_cgs = 1.6022d-12 | ||
| 148 | |||
| 149 | real(wp) :: dense, vti1, vti2, vte | ||
| 150 | real(wp) :: alami1, alami2, alame, frecol_base | ||
| 151 | |||
| 152 | 1725 | enrat_out(1) = ealpha/tempi1 | |
| 153 | 1725 | enrat_out(2) = ealpha/tempi2 | |
| 154 | 1725 | enrat_out(3) = ealpha/tempe | |
| 155 | |||
| 156 | 1725 | v0 = sqrt(2.d0*ealpha*ev_cgs/(4.d0*pmass)) | |
| 157 | 1725 | vti1 = sqrt(2.d0*tempi1*ev_cgs/(pmass*am1)) | |
| 158 | 1725 | vti2 = sqrt(2.d0*tempi2*ev_cgs/(pmass*am2)) | |
| 159 | 1725 | vte = sqrt(2.d0*tempe*ev_cgs/emass) | |
| 160 | |||
| 161 | 1725 | velrat_out(1) = v0/vti1 | |
| 162 | 1725 | velrat_out(2) = v0/vti2 | |
| 163 | 1725 | velrat_out(3) = v0/vte | |
| 164 | |||
| 165 | 1725 | dense = densi1*Z1 + densi2*Z2 | |
| 166 |
2/2✓ Branch 0 taken 113 times.
✓ Branch 1 taken 1612 times.
|
1725 | if (dense == 0.0_wp) then |
| 167 | 113 | efcolf_out = 0.0_wp | |
| 168 | 113 | return | |
| 169 | end if | ||
| 170 | alami1 = 23.d0 - log(max(epsilon(1.d0), & | ||
| 171 | 1612 | sqrt(densi1*Z1**2/tempi1)*2.d0*Z1*(4.d0 + am1)/(4.d0*tempi1 + am1*ealpha))) | |
| 172 | alami2 = 23.d0 - log(max(epsilon(1.d0), & | ||
| 173 | 1612 | sqrt(densi2*Z2**2/tempi2)*2.d0*Z2*(4.d0 + am2)/(4.d0*tempi2 + am2*ealpha))) | |
| 174 | 1612 | alame = 24.d0 - log(sqrt(dense)/tempe) | |
| 175 | |||
| 176 | 1612 | frecol_base = 2.d0*pi*dense*e_cgs**4*2.d0**2/((4.d0*pmass)**2*v0**3) | |
| 177 | 1612 | frecol_base = frecol_base/v0 | |
| 178 | |||
| 179 | 1612 | efcolf_out(1) = frecol_base*Z1**2*alami1*densi1/dense | |
| 180 | 1612 | efcolf_out(2) = frecol_base*Z2**2*alami2*densi2/dense | |
| 181 | 1612 | efcolf_out(3) = frecol_base*alame | |
| 182 | |||
| 183 |
2/2✓ Branch 0 taken 4836 times.
✓ Branch 1 taken 1612 times.
|
6448 | efcolf_out = efcolf_out*velrat_out |
| 184 | end subroutine compute_collision_coeffs | ||
| 185 | |||
| 186 | 8 | subroutine compute_rates(v0, dchichi, slowrate, dchichi_norm, slowrate_norm) | |
| 187 | real(wp), intent(in) :: v0 | ||
| 188 | real(wp), intent(out) :: dchichi, slowrate, dchichi_norm, slowrate_norm | ||
| 189 | |||
| 190 | real(wp) :: p, dpp, dhh, fpeff | ||
| 191 | |||
| 192 | 8 | p = 1.d0 | |
| 193 | 8 | call coleff(p, dpp, dhh, fpeff) | |
| 194 | |||
| 195 | 8 | dchichi = dhh*v0 | |
| 196 | 8 | slowrate = abs(fpeff)*v0 | |
| 197 | 8 | dchichi_norm = dhh | |
| 198 | 8 | slowrate_norm = abs(fpeff) | |
| 199 | 8 | end subroutine compute_rates | |
| 200 | |||
| 201 | 17 | subroutine init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) | |
| 202 | use simple_profiles, only: get_plasma_params | ||
| 203 | |||
| 204 | real(wp), intent(in) :: am1, am2, Z1, Z2, ealpha, v0 | ||
| 205 | real(wp) :: s, Te, Ti1, Ti2, ni1, ni2 | ||
| 206 | real(wp) :: densi1_loc, densi2_loc, tempi1_loc, tempi2_loc, tempe_loc | ||
| 207 | real(wp) :: v0_loc | ||
| 208 | integer :: i | ||
| 209 | |||
| 210 |
2/2✓ Branch 0 taken 1717 times.
✓ Branch 1 taken 17 times.
|
1734 | do i = 1, N_S_GRID |
| 211 | 1717 | s = dble(i - 1)/dble(N_S_GRID - 1) | |
| 212 | |||
| 213 | 1717 | call get_plasma_params(s, Te, Ti1, Ti2, ni1, ni2) | |
| 214 | |||
| 215 | 1717 | tempe_loc = max(Te, 1.0d0) | |
| 216 | 1717 | tempi1_loc = max(Ti1, 1.0d0) | |
| 217 | 1717 | tempi2_loc = max(Ti2, 1.0d0) | |
| 218 | ! Preserve exact vacuum; retain the floor for positive densities. | ||
| 219 |
3/4✓ Branch 0 taken 112 times.
✓ Branch 1 taken 1605 times.
✓ Branch 2 taken 112 times.
✗ Branch 3 not taken.
|
1717 | if (ni1 == 0.0_wp .and. ni2 == 0.0_wp) then |
| 220 | 112 | densi1_loc = 0.0_wp | |
| 221 | 112 | densi2_loc = 0.0_wp | |
| 222 | else | ||
| 223 | 1605 | densi1_loc = max(ni1*1.0d-6, 1.0d0) | |
| 224 | 1605 | densi2_loc = max(ni2*1.0d-6, 1.0d0) | |
| 225 | end if | ||
| 226 | |||
| 227 | call compute_collision_coeffs(am1, am2, Z1, Z2, & | ||
| 228 | densi1_loc, densi2_loc, & | ||
| 229 | tempi1_loc, tempi2_loc, tempe_loc, & | ||
| 230 | ealpha, v0_loc, & | ||
| 231 | 3451 | efcolf_grid(:, i), velrat_grid(:, i), enrat_grid(:, i)) | |
| 232 | end do | ||
| 233 | 17 | end subroutine init_collision_profiles | |
| 234 | |||
| 235 | 32853756 | subroutine get_local_coeffs(s, efcolf_loc, velrat_loc, enrat_loc) | |
| 236 | real(wp), intent(in) :: s | ||
| 237 | real(wp), intent(out) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) | ||
| 238 | |||
| 239 | real(wp) :: s_clamped, s_norm, weight | ||
| 240 | integer :: idx_lo, idx_hi | ||
| 241 | |||
| 242 | 32853756 | s_clamped = max(0.0d0, min(1.0d0, s)) | |
| 243 | 32853756 | s_norm = s_clamped*dble(N_S_GRID - 1) | |
| 244 | 32853756 | idx_lo = max(1, min(N_S_GRID - 1, int(s_norm) + 1)) | |
| 245 | 32853756 | idx_hi = idx_lo + 1 | |
| 246 | 32853756 | weight = s_norm - dble(idx_lo - 1) | |
| 247 | |||
| 248 | ! Stable interpolation form: reproduces a node value exactly when both | ||
| 249 | ! bracketing nodes are equal (flat profile -> scalar coefficients bitwise). | ||
| 250 | efcolf_loc = efcolf_grid(:, idx_lo) & | ||
| 251 |
2/2✓ Branch 0 taken 32853756 times.
✓ Branch 1 taken 98561268 times.
|
131415024 | + weight*(efcolf_grid(:, idx_hi) - efcolf_grid(:, idx_lo)) |
| 252 | velrat_loc = velrat_grid(:, idx_lo) & | ||
| 253 |
2/2✓ Branch 0 taken 32853756 times.
✓ Branch 1 taken 98561268 times.
|
131415024 | + weight*(velrat_grid(:, idx_hi) - velrat_grid(:, idx_lo)) |
| 254 | enrat_loc = enrat_grid(:, idx_lo) & | ||
| 255 |
2/2✓ Branch 0 taken 98561268 times.
✓ Branch 1 taken 32853756 times.
|
131415024 | + weight*(enrat_grid(:, idx_hi) - enrat_grid(:, idx_lo)) |
| 256 | 32853756 | end subroutine get_local_coeffs | |
| 257 | |||
| 258 | 32853391 | subroutine stost(z, dtauc, iswmode, ierr, max_substep_fraction, s_profile) | |
| 259 | use finite_check, only: is_finite | ||
| 260 | |||
| 261 | integer, intent(in) :: iswmode | ||
| 262 | real(wp), intent(in) :: dtauc | ||
| 263 | real(wp), intent(inout) :: z(5) | ||
| 264 | integer, intent(out) :: ierr | ||
| 265 | real(wp), intent(in), optional :: max_substep_fraction | ||
| 266 | !> Normalised toroidal flux of the profiles at z. Defaults to z(1), | ||
| 267 | !> which is s only for integ_coords TEST, CANFLUX, VMEC and BOOZER; | ||
| 268 | !> tracers pass profile_radius%profile_s(z(1:3)). | ||
| 269 | real(wp), intent(in), optional :: s_profile | ||
| 270 | |||
| 271 | real(wp), parameter :: default_fraction = 0.1_wp | ||
| 272 | real(wp), parameter :: twopi = 8.0_wp*atan(1.0_wp) | ||
| 273 | integer, parameter :: max_substeps = 100000 | ||
| 274 | real(wp) :: initial(5), fraction, remaining, step, p, dpp, dhh, fpeff | ||
| 275 | real(wp) :: alam, angle, azimuth | ||
| 276 | real(wp) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) | ||
| 277 | real :: ur | ||
| 278 | integer :: substep | ||
| 279 | logical :: pitch, energy, diffusion | ||
| 280 | |||
| 281 | 32853391 | initial = z | |
| 282 | 32853391 | ierr = 0 | |
| 283 | 32853391 | fraction = default_fraction | |
| 284 |
2/2✓ Branch 0 taken 811009 times.
✓ Branch 1 taken 32042382 times.
|
32853391 | if (present(max_substep_fraction)) fraction = max_substep_fraction |
| 285 |
6/10✓ Branch 0 taken 164266955 times.
✓ Branch 1 taken 32853391 times.
✓ Branch 2 taken 164266955 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 32853391 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 32853391 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 32853391 times.
|
197120346 | if (.not. all(is_finite(z)) .or. .not. is_finite(dtauc) .or. & |
| 286 | .not. is_finite(fraction)) then | ||
| 287 | ✗ | ierr = 1 | |
| 288 | 32853384 | return | |
| 289 | end if | ||
| 290 | if (z(4) <= 0.0_wp .or. abs(z(5)) > 1.0_wp .or. dtauc < 0.0_wp .or. & | ||
| 291 | fraction <= 0.0_wp .or. fraction > 0.25_wp .or. & | ||
| 292 |
7/14✓ Branch 0 taken 32853391 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 32853391 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 32853391 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 32853391 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 32853391 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 32853391 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 32853391 times.
|
32853391 | iswmode < 1 .or. iswmode > 4) then |
| 293 | ✗ | ierr = 1 | |
| 294 | ✗ | return | |
| 295 | end if | ||
| 296 |
1/2✓ Branch 0 taken 32853391 times.
✗ Branch 1 not taken.
|
32853391 | if (dtauc == 0.0_wp) return |
| 297 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 32853387 times.
|
32853391 | if (present(s_profile)) then |
| 298 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (.not. is_finite(s_profile)) then |
| 299 | ✗ | ierr = 1 | |
| 300 | ✗ | return | |
| 301 | end if | ||
| 302 | 4 | call get_local_coeffs(s_profile, efcolf_loc, velrat_loc, enrat_loc) | |
| 303 | else | ||
| 304 | 32853387 | call get_local_coeffs(z(1), efcolf_loc, velrat_loc, enrat_loc) | |
| 305 | end if | ||
| 306 | ! A vacuum neither advances the state nor consumes random numbers. | ||
| 307 |
6/6✓ Branch 0 taken 32853407 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 32853383 times.
✓ Branch 4 taken 32853383 times.
✓ Branch 5 taken 8 times.
|
32853415 | if (all(efcolf_loc == 0.0_wp)) return |
| 308 | 32853383 | pitch = iswmode == 1 .or. iswmode == 4 | |
| 309 | 32853383 | energy = iswmode /= 4 | |
| 310 | 32853383 | diffusion = iswmode < 3 | |
| 311 | 32853383 | remaining = dtauc | |
| 312 | |||
| 313 |
2/2✓ Branch 0 taken 142483379 times.
✓ Branch 1 taken 1 time.
|
142483380 | do substep = 1, max_substeps |
| 314 | 142483379 | p = z(4) | |
| 315 | 142483379 | call coleff_local(p, efcolf_loc, velrat_loc, enrat_loc, dpp, dhh, fpeff) | |
| 316 | if (.not. all(is_finite([dpp, dhh, fpeff])) .or. & | ||
| 317 |
6/8✓ Branch 0 taken 427450137 times.
✓ Branch 1 taken 142483379 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 427450137 times.
✓ Branch 4 taken 142483374 times.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 142483374 times.
|
569933516 | dpp < 0.0_wp .or. dhh < 0.0_wp) then |
| 318 | 5 | ierr = 1 | |
| 319 | 5 | exit | |
| 320 | end if | ||
| 321 | ! Choose from the current state before drawing increments. Bound | ||
| 322 | ! angular diffusion and both the relative drag and largest possible | ||
| 323 | ! radial noise (getran(0) has support [-sqrt(3),sqrt(3)]). | ||
| 324 | 142483374 | step = remaining | |
| 325 |
3/4✓ Branch 0 taken 84885643 times.
✓ Branch 1 taken 57597731 times.
✓ Branch 2 taken 84885643 times.
✗ Branch 3 not taken.
|
142483374 | if (pitch .and. dhh > 0.0_wp) & |
| 326 | 84885643 | step = min(step, fraction**2/dhh) | |
| 327 |
2/2✓ Branch 0 taken 141089038 times.
✓ Branch 1 taken 1394336 times.
|
142483374 | if (energy) then |
| 328 |
1/2✓ Branch 0 taken 141089038 times.
✗ Branch 1 not taken.
|
141089038 | if (abs(fpeff) > 0.0_wp) & |
| 329 | 141089038 | step = min(step, fraction*p/abs(fpeff)) | |
| 330 |
3/4✓ Branch 0 taken 141081738 times.
✓ Branch 1 taken 7300 times.
✓ Branch 2 taken 141081738 times.
✗ Branch 3 not taken.
|
141089038 | if (diffusion .and. dpp > 0.0_wp) & |
| 331 | 141081738 | step = min(step, fraction**2*p**2/(6.0_wp*dpp)) | |
| 332 | end if | ||
| 333 |
3/4✓ Branch 0 taken 142483373 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 142483373 times.
|
142483374 | if (step <= 0.0_wp .or. remaining - step == remaining) then |
| 334 | 1 | ierr = 4 ! Cannot represent further time progress. | |
| 335 | 1 | exit | |
| 336 | end if | ||
| 337 |
2/2✓ Branch 0 taken 84885642 times.
✓ Branch 1 taken 57597731 times.
|
142483373 | if (pitch) then |
| 338 | ! Isotropic random rotation of fixed small angle. Its weak | ||
| 339 | ! generator is dhh*d/dlambda[(1-lambda**2)*d/dlambda]. | ||
| 340 | ! Unlike Euler pitch kicks, it stays on the velocity sphere. | ||
| 341 | 84885642 | angle = sqrt(4.0_wp*dhh*step) | |
| 342 | 84885642 | call random_number(ur) | |
| 343 | 84885642 | azimuth = twopi*real(ur, wp) | |
| 344 | 84885642 | alam = z(5) | |
| 345 | z(5) = alam*cos(angle) + sqrt(1.0_wp - alam**2)* & | ||
| 346 | 84885642 | sin(angle)*cos(azimuth) | |
| 347 | end if | ||
| 348 |
2/2✓ Branch 0 taken 141089037 times.
✓ Branch 1 taken 1394336 times.
|
142483373 | if (energy) then |
| 349 | 141089037 | z(4) = p + fpeff*step | |
| 350 |
2/2✓ Branch 0 taken 141081737 times.
✓ Branch 1 taken 7300 times.
|
141089037 | if (diffusion) then |
| 351 | 141081737 | call getran(0, ur) | |
| 352 | 141081737 | z(4) = z(4) + sqrt(2.0_wp*dpp*step)*real(ur, wp) | |
| 353 | end if | ||
| 354 | end if | ||
| 355 |
6/10✓ Branch 0 taken 712416865 times.
✓ Branch 1 taken 142483373 times.
✓ Branch 2 taken 712416865 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 142483373 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 142483373 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 142483373 times.
|
854900238 | if (.not. all(is_finite(z)) .or. z(4) <= 0.0_wp .or. & |
| 356 | abs(z(5)) > 1.0_wp) then | ||
| 357 | ✗ | ierr = 1 | |
| 358 | ✗ | exit | |
| 359 | end if | ||
| 360 | 142483373 | remaining = remaining - step | |
| 361 |
2/2✓ Branch 1 taken 109629997 times.
✓ Branch 2 taken 32853376 times.
|
284966753 | if (remaining == 0.0_wp) return |
| 362 | end do | ||
| 363 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 6 times.
|
7 | if (ierr == 0) ierr = 5 ! Explicit computational budget exhaustion. |
| 364 | 7 | z = initial | |
| 365 | end subroutine stost | ||
| 366 | |||
| 367 | 149081737 | subroutine getran(irand, ur) | |
| 368 | integer, intent(in) :: irand | ||
| 369 | real, intent(out) :: ur | ||
| 370 | |||
| 371 | 149081737 | call random_number(ur) | |
| 372 | |||
| 373 |
2/2✓ Branch 0 taken 141081737 times.
✓ Branch 1 taken 8000000 times.
|
149081737 | if (irand == 0) then |
| 374 | 141081737 | ur = 3.464102*(ur - 0.5) | |
| 375 | else | ||
| 376 |
2/2✓ Branch 0 taken 4001417 times.
✓ Branch 1 taken 3998583 times.
|
8000000 | if (ur > 0.5) then |
| 377 | 4001417 | ur = 1.0 | |
| 378 | else | ||
| 379 | 3998583 | ur = -1.0 | |
| 380 | end if | ||
| 381 | end if | ||
| 382 | 149081737 | end subroutine getran | |
| 383 | |||
| 384 | end module collis_alp | ||
| 385 |