| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module coordinate_scaling | ||
| 2 | !> Bijective coordinate scaling transforms. | ||
| 3 | !> Single responsibility: transform coordinates before/after operations | ||
| 4 | !> like splining, with explicit forward and inverse transforms. | ||
| 5 | |||
| 6 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 7 | |||
| 8 | implicit none | ||
| 9 | |||
| 10 | type, abstract :: coordinate_scaling_t | ||
| 11 | contains | ||
| 12 | procedure(transform_iface), deferred :: transform | ||
| 13 | procedure(transform_iface), deferred :: inverse | ||
| 14 | end type coordinate_scaling_t | ||
| 15 | |||
| 16 | abstract interface | ||
| 17 | pure subroutine transform_iface(self, x_in, x_out, jacobian) | ||
| 18 | !> Transform coordinates with optional Jacobian diagonal. | ||
| 19 | !> For diagonal scaling, jacobian(i) = dx_out(i)/dx_in(i). | ||
| 20 | import :: coordinate_scaling_t, dp | ||
| 21 | class(coordinate_scaling_t), intent(in) :: self | ||
| 22 | real(dp), intent(in) :: x_in(3) | ||
| 23 | real(dp), intent(out) :: x_out(3) | ||
| 24 | real(dp), intent(out), optional :: jacobian(3) | ||
| 25 | end subroutine transform_iface | ||
| 26 | end interface | ||
| 27 | |||
| 28 | type, extends(coordinate_scaling_t) :: sqrt_s_scaling_t | ||
| 29 | !> Default scaling: r = sqrt(s), angles unchanged. | ||
| 30 | !> Improves spline accuracy near magnetic axis. | ||
| 31 | contains | ||
| 32 | procedure :: transform => sqrt_s_transform | ||
| 33 | procedure :: inverse => sqrt_s_inverse | ||
| 34 | end type sqrt_s_scaling_t | ||
| 35 | |||
| 36 | type, extends(coordinate_scaling_t) :: identity_scaling_t | ||
| 37 | !> Identity scaling: no transformation. | ||
| 38 | !> Use when coordinate system already uses rho = sqrt(s) (RHO_TOR). | ||
| 39 | contains | ||
| 40 | procedure :: transform => identity_transform | ||
| 41 | procedure :: inverse => identity_inverse | ||
| 42 | end type identity_scaling_t | ||
| 43 | |||
| 44 | contains | ||
| 45 | |||
| 46 | 65 | subroutine coordinate_scaling_clone(source, dest) | |
| 47 | class(coordinate_scaling_t), intent(in) :: source | ||
| 48 | class(coordinate_scaling_t), allocatable, intent(out) :: dest | ||
| 49 | |||
| 50 | select type (source) | ||
| 51 | type is (sqrt_s_scaling_t) | ||
| 52 | ✗ | allocate (sqrt_s_scaling_t :: dest) | |
| 53 | ✗ | dest = source | |
| 54 | type is (identity_scaling_t) | ||
| 55 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 65 times.
|
65 | allocate (identity_scaling_t :: dest) |
| 56 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 65 times.
|
130 | dest = source |
| 57 | class default | ||
| 58 | ✗ | error stop 'coordinate_scaling_clone: Unsupported scaling type' | |
| 59 | end select | ||
| 60 | 65 | end subroutine coordinate_scaling_clone | |
| 61 | |||
| 62 | 22 | pure subroutine sqrt_s_transform(self, x_in, x_out, jacobian) | |
| 63 | !> Transform (s, theta, phi) -> (r=sqrt(s), theta, phi). | ||
| 64 | !> Jacobian: dr/ds = 1/(2*sqrt(s)) = 1/(2r). | ||
| 65 | class(sqrt_s_scaling_t), intent(in) :: self | ||
| 66 | real(dp), intent(in) :: x_in(3) | ||
| 67 | real(dp), intent(out) :: x_out(3) | ||
| 68 | real(dp), intent(out), optional :: jacobian(3) | ||
| 69 | |||
| 70 | real(dp) :: r | ||
| 71 | |||
| 72 | 22 | r = sqrt(x_in(1)) | |
| 73 | 22 | x_out(1) = r | |
| 74 | 22 | x_out(2) = x_in(2) | |
| 75 | 22 | x_out(3) = x_in(3) | |
| 76 | |||
| 77 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 20 times.
|
22 | if (present(jacobian)) then |
| 78 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (r > 0d0) then |
| 79 | 2 | jacobian(1) = 0.5d0 / r | |
| 80 | else | ||
| 81 | ✗ | jacobian(1) = 0d0 | |
| 82 | end if | ||
| 83 | 2 | jacobian(2) = 1d0 | |
| 84 | 2 | jacobian(3) = 1d0 | |
| 85 | end if | ||
| 86 | 22 | end subroutine sqrt_s_transform | |
| 87 | |||
| 88 | 23090562 | pure subroutine sqrt_s_inverse(self, x_in, x_out, jacobian) | |
| 89 | !> Transform (r, theta, phi) -> (s=r^2, theta, phi). | ||
| 90 | !> Jacobian: ds/dr = 2r. | ||
| 91 | class(sqrt_s_scaling_t), intent(in) :: self | ||
| 92 | real(dp), intent(in) :: x_in(3) | ||
| 93 | real(dp), intent(out) :: x_out(3) | ||
| 94 | real(dp), intent(out), optional :: jacobian(3) | ||
| 95 | |||
| 96 | 23090562 | x_out(1) = x_in(1)**2 | |
| 97 | 23090562 | x_out(2) = x_in(2) | |
| 98 | 23090562 | x_out(3) = x_in(3) | |
| 99 | |||
| 100 |
2/2✓ Branch 0 taken 23090547 times.
✓ Branch 1 taken 15 times.
|
23090562 | if (present(jacobian)) then |
| 101 | 23090547 | jacobian(1) = 2d0 * x_in(1) | |
| 102 | 23090547 | jacobian(2) = 1d0 | |
| 103 | 23090547 | jacobian(3) = 1d0 | |
| 104 | end if | ||
| 105 | 23090562 | end subroutine sqrt_s_inverse | |
| 106 | |||
| 107 | 3 | pure subroutine identity_transform(self, x_in, x_out, jacobian) | |
| 108 | !> Identity transform: x_out = x_in. | ||
| 109 | class(identity_scaling_t), intent(in) :: self | ||
| 110 | real(dp), intent(in) :: x_in(3) | ||
| 111 | real(dp), intent(out) :: x_out(3) | ||
| 112 | real(dp), intent(out), optional :: jacobian(3) | ||
| 113 | |||
| 114 | 3 | x_out = x_in | |
| 115 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
3 | if (present(jacobian)) jacobian = 1d0 |
| 116 | 3 | end subroutine identity_transform | |
| 117 | |||
| 118 | 3 | pure subroutine identity_inverse(self, x_in, x_out, jacobian) | |
| 119 | !> Identity inverse: x_out = x_in. | ||
| 120 | class(identity_scaling_t), intent(in) :: self | ||
| 121 | real(dp), intent(in) :: x_in(3) | ||
| 122 | real(dp), intent(out) :: x_out(3) | ||
| 123 | real(dp), intent(out), optional :: jacobian(3) | ||
| 124 | |||
| 125 | 3 | x_out = x_in | |
| 126 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
3 | if (present(jacobian)) jacobian = 1d0 |
| 127 | 3 | end subroutine identity_inverse | |
| 128 | |||
| 129 | 65 | end module coordinate_scaling | |
| 130 |