| 63 |
|
logical, save :: useSpline = .true. |
| 64 |
|
integer, save :: functionType = CUBIC |
| 65 |
|
|
| 66 |
+ |
! spline structure |
| 67 |
+ |
type(cubicSpline), save :: splineSqrt |
| 68 |
|
|
| 67 |
– |
! spline variables |
| 68 |
– |
type(cubicSpline), save :: scoef |
| 69 |
– |
real ( kind = dp ), dimension(SPLINE_SEGMENTS) :: xValues |
| 70 |
– |
real ( kind = dp ), dimension(SPLINE_SEGMENTS) :: yValues |
| 71 |
– |
real ( kind = dp ), save :: dSqrt1, dSqrtN, range, dX |
| 72 |
– |
real ( kind = dp ), save :: lowerBound |
| 73 |
– |
logical, save :: uniformSpline = .true. |
| 74 |
– |
|
| 69 |
|
public::set_switch |
| 70 |
|
public::set_function_type |
| 71 |
|
public::get_switch |
| 72 |
|
|
| 73 |
|
contains |
| 74 |
|
|
| 75 |
+ |
subroutine setupSqrtSpline(rmin, rmax, np) |
| 76 |
+ |
real( kind = dp ), intent(in) :: rmin, rmax |
| 77 |
+ |
integer, intent(in) :: np |
| 78 |
+ |
real( kind = dp ) :: rvals(np), yvals(np) |
| 79 |
+ |
real( kind = dp ) :: dr, r |
| 80 |
+ |
real( kind = dp ) :: yprime1, yprimen |
| 81 |
+ |
integer :: i |
| 82 |
+ |
|
| 83 |
+ |
dr = (rmax-rmin) / float(np-1) |
| 84 |
+ |
|
| 85 |
+ |
do i = 1, np |
| 86 |
+ |
r = rmin + float(i-1)*dr |
| 87 |
+ |
rvals(i) = r |
| 88 |
+ |
yvals(i) = dsqrt(r) |
| 89 |
+ |
enddo |
| 90 |
+ |
|
| 91 |
+ |
yprime1 = 1.0d0 / ( 2.0d0 * dsqrt( rmin ) ) |
| 92 |
+ |
yprimeN = 1.0d0 / ( 2.0d0 * dsqrt( rmax ) ) |
| 93 |
+ |
|
| 94 |
+ |
call newSpline(splineSqrt, rvals, yvals, yprime1, yprimen, .true.) |
| 95 |
+ |
|
| 96 |
+ |
return |
| 97 |
+ |
end subroutine setupSqrtSpline |
| 98 |
+ |
|
| 99 |
|
subroutine set_switch(SwitchType, rinner, router) |
| 100 |
|
|
| 101 |
|
real ( kind = dp ), intent(in):: rinner, router |
| 130 |
|
|
| 131 |
|
if (.not.haveSqrtSpline) then |
| 132 |
|
! fill arrays for building the spline |
| 133 |
< |
lowerBound = 1.0d0 ! the smallest value expected for r2 |
| 134 |
< |
range = rout2(SwitchType) - lowerBound |
| 117 |
< |
dX = range / (SPLINE_SEGMENTS - 1) |
| 118 |
< |
|
| 119 |
< |
! the spline is bracketed by lowerBound and rout2 endpoints |
| 120 |
< |
xValues(1) = lowerBound |
| 121 |
< |
yValues(1) = dsqrt(lowerBound) |
| 122 |
< |
do i = 1, SPLINE_SEGMENTS-1 |
| 123 |
< |
xValues(i+1) = i * dX |
| 124 |
< |
yValues(i+1) = dsqrt( i * dX ) |
| 125 |
< |
enddo |
| 126 |
< |
|
| 127 |
< |
! set the endpoint derivatives |
| 128 |
< |
dSqrt1 = 1 / ( 2.0d0 * dsqrt( xValues(1) ) ) |
| 129 |
< |
dSqrtN = 1 / ( 2.0d0 * dsqrt( xValues(SPLINE_SEGMENTS) ) ) |
| 130 |
< |
|
| 131 |
< |
! call newSpline to fill the coefficient array |
| 132 |
< |
call newSpline(scoef, xValues, yValues, dSqrt1, dSqrtN, uniformSpline) |
| 133 |
< |
|
| 133 |
> |
call setupSqrtSpline(1.0d0, rout2(switchType), SPLINE_SEGMENTS) |
| 134 |
> |
haveSqrtSpline = .true. |
| 135 |
|
endif |
| 136 |
|
|
| 137 |
|
end subroutine set_switch |
| 154 |
|
|
| 155 |
|
real( kind = dp ), intent(in) :: r2 |
| 156 |
|
real( kind = dp ), intent(inout) :: sw, dswdr, r |
| 157 |
< |
real( kind = dp ) :: ron, roff |
| 157 |
> |
real( kind = dp ) :: ron, roff, a, b, c, d, dx |
| 158 |
|
real( kind = dp ) :: rval, rval2, rval3, rval4, rval5 |
| 159 |
|
real( kind = dp ) :: rvaldi, rvaldi2, rvaldi3, rvaldi4, rvaldi5 |
| 160 |
|
integer, intent(in) :: SwitchType |
| 161 |
|
logical, intent(inout) :: in_switching_region |
| 162 |
+ |
integer :: j |
| 163 |
|
|
| 164 |
|
sw = 0.0d0 |
| 165 |
|
dswdr = 0.0d0 |
| 180 |
|
|
| 181 |
|
else |
| 182 |
|
if (useSpline) then |
| 183 |
< |
call lookup_uniform_spline(scoef, r2, r) |
| 183 |
> |
j = MAX(1, MIN(splineSqrt%np, idint((r2-splineSqrt%x(1)) * splineSqrt%dx_i) + 1)) |
| 184 |
> |
|
| 185 |
> |
dx = r2 - splineSqrt%x(j) |
| 186 |
> |
|
| 187 |
> |
a = splineSqrt%c(1,j) |
| 188 |
> |
b = splineSqrt%c(2,j) |
| 189 |
> |
c = splineSqrt%c(3,j) |
| 190 |
> |
d = splineSqrt%c(4,j) |
| 191 |
> |
|
| 192 |
> |
r = c + dx * d |
| 193 |
> |
r = b + dx * r |
| 194 |
> |
r = a + dx * r |
| 195 |
|
else |
| 196 |
|
r = dsqrt(r2) |
| 197 |
|
endif |
