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 |