| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_line_spectre | ||
| 2 | !> Symplectic field-line tracer for a single SPECTRE/SPEC volume. | ||
| 3 | !> | ||
| 4 | !> The magnetic field line is the flow of a one-degree-of-freedom canonical | ||
| 5 | !> system with the toroidal angle zeta as time, the poloidal angle theta as | ||
| 6 | !> coordinate, and the covariant vector-potential component psi = A_theta as | ||
| 7 | !> the conjugate momentum. The radial coordinate s is recovered implicitly | ||
| 8 | !> from psi at fixed (theta, zeta). The derived flow is | ||
| 9 | !> dtheta/dzeta = -dA_zeta/ds / dA_theta/ds | ||
| 10 | !> dpsi/dzeta = dA_zeta/dtheta - dA_theta/dtheta * dA_zeta/ds/dA_theta/ds | ||
| 11 | !> (the explicit d/dzeta terms cancel exactly). One step advances the pair | ||
| 12 | !> (theta, psi) by an implicit midpoint map. Both update signs | ||
| 13 | !> set the field-line orientation, so a sign flip reverses iota; keep them | ||
| 14 | !> as written. | ||
| 15 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 16 | use spectre_reader, only: spectre_data_t | ||
| 17 | use spectre_basis, only: spectre_vecpot_t, eval_spectre_vector_potential | ||
| 18 | use interpolate, only: BatchSplineData3D, construct_batch_splines_3d, & | ||
| 19 | destroy_batch_splines_3d, & | ||
| 20 | evaluate_batch_splines_3d_der | ||
| 21 | |||
| 22 | implicit none | ||
| 23 | private | ||
| 24 | |||
| 25 | public :: field_line_psi, field_line_step, spectre_rz | ||
| 26 | public :: prepare_field_line_spline, cleanup_field_line_spline | ||
| 27 | |||
| 28 | integer, parameter :: max_newton = 20 | ||
| 29 | real(dp), parameter :: newton_tol = 1.0e-13_dp | ||
| 30 | real(dp), parameter :: fd_scale = epsilon(1.0_dp)**(1.0_dp/3.0_dp) | ||
| 31 | real(dp), parameter :: twopi = 6.283185307179586476925286766559_dp | ||
| 32 | type(BatchSplineData3D) :: hotloop_spline | ||
| 33 | integer :: hotloop_lvol = 0 | ||
| 34 | |||
| 35 | contains | ||
| 36 | |||
| 37 | !> Canonical momentum psi = A_theta at (lvol, s, theta, zeta). | ||
| 38 | 27 | pure function field_line_psi(data, lvol, s, theta, zeta) result(psi) | |
| 39 | type(spectre_data_t), intent(in) :: data | ||
| 40 | integer, intent(in) :: lvol | ||
| 41 | real(dp), intent(in) :: s, theta, zeta | ||
| 42 | real(dp) :: psi | ||
| 43 | |||
| 44 | type(spectre_vecpot_t) :: av | ||
| 45 | |||
| 46 | call eval_spectre_vector_potential(data, lvol, s, theta, zeta, av, & | ||
| 47 | 27 | second_derivatives=.false.) | |
| 48 | 27 | psi = av%Ath | |
| 49 |
8/8✓ Branch 0 taken 27 times.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 27 times.
✓ Branch 3 taken 81 times.
✓ Branch 4 taken 27 times.
✓ Branch 5 taken 162 times.
✓ Branch 6 taken 162 times.
✓ Branch 7 taken 27 times.
|
594 | end function field_line_psi |
| 50 | |||
| 51 | 7 | subroutine prepare_field_line_spline(data, lvol, n_s, n_theta, n_zeta, order) | |
| 52 | type(spectre_data_t), intent(in) :: data | ||
| 53 | integer, intent(in) :: lvol, n_s, n_theta, n_zeta, order | ||
| 54 | |||
| 55 | type(spectre_vecpot_t) :: av | ||
| 56 | real(dp), allocatable :: values(:, :, :, :) | ||
| 57 | real(dp) :: xmin(3), xmax(3), x(3) | ||
| 58 | integer :: i, j, k | ||
| 59 | |||
| 60 | 7 | call cleanup_field_line_spline | |
| 61 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
|
7 | if (min(n_s, n_theta, n_zeta) <= order) then |
| 62 | ✗ | error stop 'field-line spline grid must exceed its order' | |
| 63 | end if | ||
| 64 | 7 | xmin = [-1.0_dp, 0.0_dp, 0.0_dp] | |
| 65 |
2/2✓ Branch 0 taken 21 times.
✓ Branch 1 taken 7 times.
|
28 | xmax = [1.0_dp, twopi, twopi/real(data%Nfp, dp)] |
| 66 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 7 times.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 7 times.
✓ Branch 10 taken 7 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 7 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 7 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 7 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 7 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 7 times.
|
42 | allocate (values(n_s, n_theta, n_zeta, 2)) |
| 67 |
2/2✓ Branch 0 taken 224 times.
✓ Branch 1 taken 7 times.
|
231 | do k = 1, n_zeta |
| 68 | 224 | x(3) = grid_coordinate(k, n_zeta, xmin(3), xmax(3)) | |
| 69 |
2/2✓ Branch 0 taken 7168 times.
✓ Branch 1 taken 224 times.
|
7399 | do j = 1, n_theta |
| 70 | 7168 | x(2) = grid_coordinate(j, n_theta, xmin(2), xmax(2)) | |
| 71 |
2/2✓ Branch 0 taken 229376 times.
✓ Branch 1 taken 7168 times.
|
236768 | do i = 1, n_s |
| 72 | 229376 | x(1) = grid_coordinate(i, n_s, xmin(1), xmax(1)) | |
| 73 | call eval_spectre_vector_potential(data, lvol, x(1), x(2), & | ||
| 74 | x(3), av, & | ||
| 75 | 229376 | second_derivatives=.false.) | |
| 76 |
2/2✓ Branch 0 taken 458752 times.
✓ Branch 1 taken 229376 times.
|
695296 | values(i, j, k, :) = [av%Ath, av%Azt] |
| 77 | end do | ||
| 78 | end do | ||
| 79 | end do | ||
| 80 | call construct_batch_splines_3d(xmin, xmax, values, [order, order, order], & | ||
| 81 |
2/2✓ Branch 0 taken 21 times.
✓ Branch 1 taken 7 times.
|
28 | [.false., .true., .true.], hotloop_spline) |
| 82 | 7 | hotloop_lvol = lvol | |
| 83 |
9/10✓ Branch 0 taken 7 times.
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 21 times.
✓ Branch 4 taken 7 times.
✓ Branch 5 taken 42 times.
✓ Branch 6 taken 42 times.
✓ Branch 7 taken 7 times.
✓ Branch 9 taken 7 times.
✗ Branch 10 not taken.
|
161 | end subroutine prepare_field_line_spline |
| 84 | |||
| 85 | 14 | subroutine cleanup_field_line_spline | |
| 86 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 7 times.
|
14 | if (hotloop_lvol /= 0) call destroy_batch_splines_3d(hotloop_spline) |
| 87 | 14 | hotloop_lvol = 0 | |
| 88 | 14 | end subroutine cleanup_field_line_spline | |
| 89 | |||
| 90 | 236768 | pure function grid_coordinate(i, n, xmin, xmax) result(x) | |
| 91 | integer, intent(in) :: i, n | ||
| 92 | real(dp), intent(in) :: xmin, xmax | ||
| 93 | real(dp) :: x | ||
| 94 | |||
| 95 | 236768 | x = xmin + (xmax - xmin)*real(i - 1, dp)/real(n - 1, dp) | |
| 96 | end function grid_coordinate | ||
| 97 | |||
| 98 | !> Advance the canonical pair (theta, psi) by implicit midpoint. The two | ||
| 99 | !> midpoint residuals are solved in (s_mid, theta_mid); the endpoint radius | ||
| 100 | !> is then recovered from A_theta(s1,theta1,zeta1) = psi1. Numerical | ||
| 101 | !> Jacobians use only values and first derivatives of the vector potential. | ||
| 102 | 1029762 | subroutine field_line_step(data, lvol, h, s, theta, zeta, psi) | |
| 103 | type(spectre_data_t), intent(in) :: data | ||
| 104 | integer, intent(in) :: lvol | ||
| 105 | real(dp), intent(in) :: h | ||
| 106 | real(dp), intent(inout) :: s, theta, zeta, psi | ||
| 107 | |||
| 108 | type(spectre_vecpot_t) :: av | ||
| 109 | real(dp) :: s0, theta0, zeta0, psi0, smid, theta_delta | ||
| 110 | logical :: midpoint_ok, radius_ok | ||
| 111 | |||
| 112 | 1029762 | s0 = s | |
| 113 | 1029762 | theta0 = theta | |
| 114 | 1029762 | zeta0 = zeta | |
| 115 | 1029762 | psi0 = psi | |
| 116 | 1029762 | smid = s0 | |
| 117 | 1029762 | theta_delta = 0.0_dp | |
| 118 | call solve_midpoint(data, lvol, h, s0, theta0, zeta0, psi0, smid, & | ||
| 119 | 1029762 | theta_delta, av, midpoint_ok) | |
| 120 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1029762 times.
|
1029762 | if (.not. midpoint_ok) then |
| 121 | write (*, '(A,I0,A,2(ES13.5,1X))') & | ||
| 122 | ✗ | 'field_line_step: midpoint Newton failed in volume ', lvol, & | |
| 123 | ✗ | ', (s, theta) = ', smid, theta0 + theta_delta | |
| 124 | ✗ | error stop 'field_line_step: midpoint Newton did not converge' | |
| 125 | end if | ||
| 126 | |||
| 127 | 1029762 | theta = theta0 + 2.0_dp*theta_delta | |
| 128 | 1029762 | psi = 2.0_dp*av%Ath - psi0 | |
| 129 | 1029762 | zeta = zeta0 + h | |
| 130 | 1029762 | s = 2.0_dp*smid - s0 | |
| 131 | 1029762 | call recover_radius(data, lvol, theta, zeta, psi, s, radius_ok) | |
| 132 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1029762 times.
|
1029762 | if (.not. radius_ok) error stop 'field_line_step: endpoint radius failed' |
| 133 |
8/8✓ Branch 0 taken 1029762 times.
✓ Branch 1 taken 3089286 times.
✓ Branch 2 taken 1029762 times.
✓ Branch 3 taken 3089286 times.
✓ Branch 4 taken 1029762 times.
✓ Branch 5 taken 6178572 times.
✓ Branch 6 taken 6178572 times.
✓ Branch 7 taken 1029762 times.
|
22654764 | end subroutine field_line_step |
| 134 | |||
| 135 |
8/8✓ Branch 0 taken 3089286 times.
✓ Branch 1 taken 1029762 times.
✓ Branch 2 taken 3089286 times.
✓ Branch 3 taken 1029762 times.
✓ Branch 4 taken 6178572 times.
✓ Branch 5 taken 1029762 times.
✓ Branch 6 taken 6178572 times.
✓ Branch 7 taken 1029762 times.
|
19565478 | subroutine solve_midpoint(data, lvol, h, s0, theta0, zeta0, psi0, & |
| 136 | s, theta_delta, av, converged) | ||
| 137 | type(spectre_data_t), intent(in) :: data | ||
| 138 | integer, intent(in) :: lvol | ||
| 139 | real(dp), intent(in) :: h, s0, theta0, zeta0, psi0 | ||
| 140 | real(dp), intent(inout) :: s, theta_delta | ||
| 141 | type(spectre_vecpot_t), intent(out) :: av | ||
| 142 | logical, intent(out) :: converged | ||
| 143 | |||
| 144 | real(dp) :: residual(2), jac(2, 2), det, ds, dtheta | ||
| 145 | integer :: iter | ||
| 146 | |||
| 147 | 1029762 | converged = .false. | |
| 148 |
1/2✓ Branch 0 taken 3089273 times.
✗ Branch 1 not taken.
|
3089273 | do iter = 1, max_newton |
| 149 | call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & | ||
| 150 | 3089273 | s, theta_delta, residual, av) | |
| 151 |
8/10✗ Branch 0 not taken.
✓ Branch 1 taken 3089273 times.
✓ Branch 2 taken 3089273 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6178546 times.
✓ Branch 5 taken 3089273 times.
✓ Branch 6 taken 4656957 times.
✓ Branch 7 taken 1521589 times.
✓ Branch 8 taken 1029762 times.
✓ Branch 9 taken 2059511 times.
|
12357092 | if (maxval(abs(residual)) < newton_tol) then |
| 152 | 1029762 | converged = .true. | |
| 153 | 1029762 | return | |
| 154 | end if | ||
| 155 | call midpoint_jacobian(data, lvol, h, s0, theta0, zeta0, psi0, & | ||
| 156 | 2059511 | s, theta_delta, jac) | |
| 157 | 2059511 | det = jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1) | |
| 158 |
1/2✓ Branch 0 taken 2059511 times.
✗ Branch 1 not taken.
|
2059511 | if (abs(det) <= tiny(det)) return |
| 159 | 2059511 | ds = (-residual(1)*jac(2, 2) + residual(2)*jac(1, 2))/det | |
| 160 | 2059511 | dtheta = (residual(1)*jac(2, 1) - residual(2)*jac(1, 1))/det | |
| 161 | 2059511 | s = s + ds | |
| 162 | 2059511 | theta_delta = theta_delta + dtheta | |
| 163 | end do | ||
| 164 | end subroutine solve_midpoint | ||
| 165 | |||
| 166 |
8/8✓ Branch 0 taken 11327317 times.
✓ Branch 1 taken 33981951 times.
✓ Branch 2 taken 11327317 times.
✓ Branch 3 taken 33981951 times.
✓ Branch 4 taken 11327317 times.
✓ Branch 5 taken 67963902 times.
✓ Branch 6 taken 67963902 times.
✓ Branch 7 taken 11327317 times.
|
249200974 | subroutine midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & |
| 167 | s, theta_delta, residual, av) | ||
| 168 | type(spectre_data_t), intent(in) :: data | ||
| 169 | integer, intent(in) :: lvol | ||
| 170 | real(dp), intent(in) :: h, s0, theta0, zeta0, psi0, s, theta_delta | ||
| 171 | real(dp), intent(out) :: residual(2) | ||
| 172 | type(spectre_vecpot_t), intent(out) :: av | ||
| 173 | |||
| 174 | real(dp) :: ratio, flow | ||
| 175 | |||
| 176 | call evaluate_hotloop_av(data, lvol, s, theta0 + theta_delta, & | ||
| 177 | 11327317 | zeta0 + 0.5_dp*h, av) | |
| 178 | 11327317 | ratio = av%dAzt(1)/av%dAth(1) | |
| 179 | 11327317 | flow = av%dAzt(2) - av%dAth(2)*ratio | |
| 180 | 11327317 | residual(1) = theta_delta + 0.5_dp*h*ratio | |
| 181 | 11327317 | residual(2) = av%Ath - psi0 - 0.5_dp*h*flow | |
| 182 | 11327317 | end subroutine midpoint_residual | |
| 183 | |||
| 184 | 2059511 | subroutine midpoint_jacobian(data, lvol, h, s0, theta0, zeta0, psi0, & | |
| 185 | s, theta_delta, jac) | ||
| 186 | type(spectre_data_t), intent(in) :: data | ||
| 187 | integer, intent(in) :: lvol | ||
| 188 | real(dp), intent(in) :: h, s0, theta0, zeta0, psi0, s, theta_delta | ||
| 189 | real(dp), intent(out) :: jac(2, 2) | ||
| 190 | |||
| 191 | type(spectre_vecpot_t) :: av | ||
| 192 | real(dp) :: plus(2), minus(2), delta | ||
| 193 | |||
| 194 | 2059511 | delta = fd_scale*max(1.0_dp, abs(s)) | |
| 195 | call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & | ||
| 196 | 2059511 | s + delta, theta_delta, plus, av) | |
| 197 | call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & | ||
| 198 | 2059511 | s - delta, theta_delta, minus, av) | |
| 199 |
2/2✓ Branch 0 taken 4119022 times.
✓ Branch 1 taken 2059511 times.
|
6178533 | jac(:, 1) = (plus - minus)/(2.0_dp*delta) |
| 200 | |||
| 201 | 2059511 | delta = fd_scale | |
| 202 | call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & | ||
| 203 | 2059511 | s, theta_delta + delta, plus, av) | |
| 204 | call midpoint_residual(data, lvol, h, s0, theta0, zeta0, psi0, & | ||
| 205 | 2059511 | s, theta_delta - delta, minus, av) | |
| 206 |
2/2✓ Branch 0 taken 4119022 times.
✓ Branch 1 taken 2059511 times.
|
6178533 | jac(:, 2) = (plus - minus)/(2.0_dp*delta) |
| 207 |
8/8✓ Branch 0 taken 2059511 times.
✓ Branch 1 taken 6178533 times.
✓ Branch 2 taken 2059511 times.
✓ Branch 3 taken 6178533 times.
✓ Branch 4 taken 2059511 times.
✓ Branch 5 taken 12357066 times.
✓ Branch 6 taken 12357066 times.
✓ Branch 7 taken 2059511 times.
|
45309242 | end subroutine midpoint_jacobian |
| 208 | |||
| 209 | 1029762 | subroutine recover_radius(data, lvol, theta, zeta, psi, s, converged) | |
| 210 | type(spectre_data_t), intent(in) :: data | ||
| 211 | integer, intent(in) :: lvol | ||
| 212 | real(dp), intent(in) :: theta, zeta, psi | ||
| 213 | real(dp), intent(inout) :: s | ||
| 214 | logical, intent(out) :: converged | ||
| 215 | |||
| 216 | type(spectre_vecpot_t) :: av | ||
| 217 | real(dp) :: residual | ||
| 218 | integer :: iter | ||
| 219 | |||
| 220 | 1029762 | converged = .false. | |
| 221 |
1/2✓ Branch 0 taken 2070305 times.
✗ Branch 1 not taken.
|
2070305 | do iter = 1, max_newton |
| 222 | 2070305 | call evaluate_hotloop_av(data, lvol, s, theta, zeta, av) | |
| 223 | 2070305 | residual = av%Ath - psi | |
| 224 |
2/2✓ Branch 0 taken 1029762 times.
✓ Branch 1 taken 1040543 times.
|
2070305 | if (abs(residual) < newton_tol) then |
| 225 | 1029762 | converged = .true. | |
| 226 | 1029762 | return | |
| 227 | end if | ||
| 228 | 1040543 | s = s - residual/av%dAth(1) | |
| 229 | end do | ||
| 230 |
8/8✓ Branch 0 taken 3089286 times.
✓ Branch 1 taken 1029762 times.
✓ Branch 2 taken 3089286 times.
✓ Branch 3 taken 1029762 times.
✓ Branch 4 taken 6178572 times.
✓ Branch 5 taken 1029762 times.
✓ Branch 6 taken 6178572 times.
✓ Branch 7 taken 1029762 times.
|
19565478 | end subroutine recover_radius |
| 231 | |||
| 232 |
8/8✓ Branch 0 taken 40192866 times.
✓ Branch 1 taken 13397622 times.
✓ Branch 2 taken 40192866 times.
✓ Branch 3 taken 13397622 times.
✓ Branch 4 taken 80385732 times.
✓ Branch 5 taken 13397622 times.
✓ Branch 6 taken 80385732 times.
✓ Branch 7 taken 13397622 times.
|
254554818 | subroutine evaluate_hotloop_av(data, lvol, s, theta, zeta, av) |
| 233 | type(spectre_data_t), intent(in) :: data | ||
| 234 | integer, intent(in) :: lvol | ||
| 235 | real(dp), intent(in) :: s, theta, zeta | ||
| 236 | type(spectre_vecpot_t), intent(out) :: av | ||
| 237 | |||
| 238 | real(dp) :: x(3), values(2), derivatives(3, 2) | ||
| 239 | |||
| 240 |
2/2✓ Branch 0 taken 13397594 times.
✓ Branch 1 taken 28 times.
|
13397622 | if (hotloop_lvol == lvol) then |
| 241 | x = [s, modulo(theta, twopi), & | ||
| 242 |
6/10✓ Branch 0 taken 13397594 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 13397594 times.
✓ Branch 4 taken 13397594 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 13397594 times.
✓ Branch 8 taken 40192782 times.
✓ Branch 9 taken 13397594 times.
|
53590376 | modulo(zeta, twopi/real(data%Nfp, dp))] |
| 243 | 13397594 | call evaluate_batch_splines_3d_der(hotloop_spline, x, values, derivatives) | |
| 244 | 13397594 | av%Ath = values(1) | |
| 245 | 13397594 | av%Azt = values(2) | |
| 246 |
2/2✓ Branch 0 taken 40192782 times.
✓ Branch 1 taken 13397594 times.
|
53590376 | av%dAth = derivatives(:, 1) |
| 247 |
2/2✓ Branch 0 taken 40192782 times.
✓ Branch 1 taken 13397594 times.
|
53590376 | av%dAzt = derivatives(:, 2) |
| 248 | else | ||
| 249 | call eval_spectre_vector_potential(data, lvol, s, theta, zeta, av, & | ||
| 250 | 28 | second_derivatives=.false.) | |
| 251 | end if | ||
| 252 | 13397622 | end subroutine evaluate_hotloop_av | |
| 253 | |||
| 254 | !> Cylindrical R, Z of a point in volume lvol from the SPEC interface Fourier | ||
| 255 | !> blend, mirroring libneo_coordinates_spectre for output only: interior | ||
| 256 | !> volumes blend interfaces lvol-1 and lvol linearly in s; the axis volume | ||
| 257 | !> (lvol == 1) uses the sbar**m power law (sbar**2 for m = 0) with | ||
| 258 | !> sbar = (s + 1)/2. | ||
| 259 | 3816 | pure subroutine spectre_rz(data, lvol, s, theta, zeta, R, Z) | |
| 260 | type(spectre_data_t), intent(in) :: data | ||
| 261 | integer, intent(in) :: lvol | ||
| 262 | real(dp), intent(in) :: s, theta, zeta | ||
| 263 | real(dp), intent(out) :: R, Z | ||
| 264 | |||
| 265 | integer :: ii, m, n | ||
| 266 | real(dp) :: arg, carg, sarg, sbar, fj, alss, blss | ||
| 267 | real(dp) :: rc, rs, zc, zs | ||
| 268 | |||
| 269 | 3816 | R = 0.0_dp | |
| 270 | 3816 | Z = 0.0_dp | |
| 271 |
2/2✓ Branch 0 taken 87768 times.
✓ Branch 1 taken 3816 times.
|
91584 | do ii = 1, data%mn |
| 272 | 87768 | m = data%im(ii) | |
| 273 | 87768 | n = data%in(ii) | |
| 274 | 87768 | arg = real(m, dp)*theta - real(n, dp)*zeta | |
| 275 | 87768 | carg = cos(arg) | |
| 276 | 87768 | sarg = sin(arg) | |
| 277 |
2/2✓ Branch 0 taken 27692 times.
✓ Branch 1 taken 60076 times.
|
87768 | if (lvol == 1) then |
| 278 | 27692 | sbar = 0.5_dp*(s + 1.0_dp) | |
| 279 |
2/2✓ Branch 0 taken 3612 times.
✓ Branch 1 taken 24080 times.
|
27692 | if (m == 0) then |
| 280 | 3612 | fj = sbar*sbar | |
| 281 | else | ||
| 282 | 24080 | fj = sbar**m | |
| 283 | end if | ||
| 284 | 27692 | rc = data%Rbc(ii, 0) + (data%Rbc(ii, 1) - data%Rbc(ii, 0))*fj | |
| 285 | 27692 | rs = data%Rbs(ii, 0) + (data%Rbs(ii, 1) - data%Rbs(ii, 0))*fj | |
| 286 | 27692 | zc = data%Zbc(ii, 0) + (data%Zbc(ii, 1) - data%Zbc(ii, 0))*fj | |
| 287 | 27692 | zs = data%Zbs(ii, 0) + (data%Zbs(ii, 1) - data%Zbs(ii, 0))*fj | |
| 288 | else | ||
| 289 | 60076 | alss = 0.5_dp*(1.0_dp - s) | |
| 290 | 60076 | blss = 0.5_dp*(1.0_dp + s) | |
| 291 | 60076 | rc = alss*data%Rbc(ii, lvol - 1) + blss*data%Rbc(ii, lvol) | |
| 292 | 60076 | rs = alss*data%Rbs(ii, lvol - 1) + blss*data%Rbs(ii, lvol) | |
| 293 | 60076 | zc = alss*data%Zbc(ii, lvol - 1) + blss*data%Zbc(ii, lvol) | |
| 294 | 60076 | zs = alss*data%Zbs(ii, lvol - 1) + blss*data%Zbs(ii, lvol) | |
| 295 | end if | ||
| 296 | 87768 | R = R + rc*carg + rs*sarg | |
| 297 | 91584 | Z = Z + zc*carg + zs*sarg | |
| 298 | end do | ||
| 299 | 3816 | end subroutine spectre_rz | |
| 300 | |||
| 301 | end module field_line_spectre | ||
| 302 |