32 #include "../base/base_uses.f90"
38 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'splines'
56 REAL(
dp),
INTENT(in) :: x(:), y(:), xnew(:)
57 REAL(
dp) :: ynew(size(xnew))
60 REAL(
dp) :: c(0:4, size(x) - 1)
63 CALL spline3pars(x, y, [2, 2], [0._dp, 0._dp], c)
67 ip = iixmin(xnew(i), x, ip)
68 ynew(i) = poly3(xnew(i), c(:, ip))
83 REAL(
dp),
INTENT(in) :: x(:), y(:), xnew(:)
84 REAL(
dp),
INTENT(out),
OPTIONAL :: ynew(:), dynew(:), d2ynew(:)
87 REAL(
dp) :: c(0:4, size(x) - 1)
89 CALL spline3pars(x, y, [2, 2], [0._dp, 0._dp], c)
93 ip = iixmin(xnew(i), x, ip)
94 IF (
PRESENT(ynew)) ynew(i) = poly3(xnew(i), c(:, ip))
95 IF (
PRESENT(dynew)) dynew(i) = dpoly3(xnew(i), c(:, ip))
96 IF (
PRESENT(d2ynew)) d2ynew(i) = d2poly3(xnew(i), c(:, ip))
108 SUBROUTINE spline3pars(xi, yi, bctype, bcval, c)
111 REAL(
dp),
INTENT(in) :: xi(:), yi(:)
112 INTEGER,
INTENT(in) :: bctype(2)
113 REAL(
dp),
INTENT(in) :: bcval(2)
114 REAL(
dp),
INTENT(out) :: c(0:, :)
116 INTEGER :: i, i2, info, ipiv(4), j, n
117 REAL(
dp) :: ae(4, 4), be(4), bemat(4, 1), c1, c2, c3, c4, ce(4), d2p1, d2pn, x0, xe(4), &
118 ye(4), hi(size(c, 2)), cs(2*size(c, 2)), bs(2*size(c, 2)), bmat(2*size(c, 2), 1), &
120 INTEGER :: ipiv2(2*size(c, 2))
147 IF (bctype(1) < 1 .OR. bctype(1) > 2)
CALL stop_error(
"spline3pars error: bctype /= 1 or 2.")
148 IF (bctype(2) < 1 .OR. bctype(2) > 2)
CALL stop_error(
"spline3pars error: bctype /= 1 or 2.")
149 IF (
SIZE(c, 1) /= 5)
CALL stop_error(
"spline3pars error: size(c,1) /= 5.")
150 IF (
SIZE(c, 2) /=
SIZE(xi) - 1)
CALL stop_error(
"spline3pars error: size(c,2) /= size(xi)-1.")
151 IF (
SIZE(xi) /=
SIZE(yi))
CALL stop_error(
"spline3pars error: size(xi) /= size(yi)")
160 hi(i) = xi(i + 1) - xi(i)
165 IF (bctype(1) == 2)
THEN
166 IF (n < 4)
CALL stop_error(
"spline3pars error: n < 4")
172 ae(i, j) = (xe(i) - x0)**(j - 1)
176 be = ye; bemat(:, 1) = be
177 CALL dgesv(4, 1, ae, 4, ipiv, bemat, 4, info)
178 IF (info /= 0)
CALL stop_error(
"spline3pars error: dgesv error.")
183 IF (bctype(2) == 2)
THEN
184 IF (n < 4)
CALL stop_error(
"spline3pars error: n < 4")
190 ae(i, j) = (xe(i) - x0)**(j - 1)
194 be = ye; bemat(:, 1) = be
195 CALL dgesv(4, 1, ae, 4, ipiv, bemat, 4, info)
196 IF (info /= 0)
CALL stop_error(
"spline3pars error: dgesv error.")
202 IF (bctype(1) == 1) d2p1 = bcval(1)
203 IF (bctype(2) == 1) d2pn = bcval(2)
211 as(4, 1) = 6/hi(1)**2
216 as(5, i2 - 1) = 1/hi(i - 1)
217 as(4, i2) = 2/hi(i - 1)
218 as(3, i2 + 1) = 2/hi(i)
219 as(2, i2 + 2) = 1/hi(i)
220 as(5, i2) = 1/hi(i - 1)**2
221 as(4, i2 + 1) = -1/hi(i)**2
222 bs(i2) = (yi(i + 1) - yi(i))/hi(i) - (yi(i) - yi(i - 1))/hi(i - 1)
226 as(4, 2*(n - 1)) = 6/hi(n - 1)**2
231 CALL lapack_sgbsv(2*(n - 1), 1, 2, 1, as, 5, ipiv2, bmat, 2*(n - 1), info)
232 IF (info /= 0)
CALL stop_error(
"spline3pars error: dgbsv error.")
245 c(2, i) = -(c1 - c2 + 2*c3 + c4)/hi(i)
246 c(3, i) = 3*c3/hi(i)**2
247 c(4, i) = (-c3 + c4)/hi(i)**3
261 SUBROUTINE spline3valder(x, xi, c, val, der)
264 REAL(
dp),
INTENT(in) :: x, xi(:), c(0:, :)
265 REAL(
dp),
INTENT(out) :: val, der
280 IF (
SIZE(c, 1) /= 5)
CALL stop_error(
"spline3 error: size(c,1) /= 5.")
281 IF (
SIZE(c, 2) /=
SIZE(xi) - 1)
CALL stop_error(
"spline3 error: size(c,2) /= size(xi)-1.")
285 val = poly3(x, c(:, i1))
286 der = dpoly3(x, c(:, i1))
297 INTEGER FUNCTION iix(x, xi)
RESULT(i1)
303 REAL(
dp),
INTENT(in) :: x, xi(:)
314 CALL stop_error(
"error in iix: n < 2")
323 ELSEIF (x <= xi(1))
THEN
325 ELSEIF (x <= xi(2))
THEN
327 ELSEIF (x <= xi(3))
THEN
329 ELSEIF (x >= xi(n))
THEN
335 IF (i2 - i1 == 1)
EXIT
336 ic = i1 + (i2 - i1)/2
337 IF (x >= xi(ic))
THEN
353 INTEGER FUNCTION iixmin(x, xi, i_min)
RESULT(ip)
355 REAL(
dp),
INTENT(in) :: x, xi(:)
356 INTEGER,
INTENT(in) :: i_min
358 IF (i_min >= 1 .AND. i_min <=
SIZE(xi) - 1)
THEN
359 ip = iix(x, xi(i_min:)) + i_min - 1
375 FUNCTION iixun(x, n, x1, xn)
380 REAL(
dp),
INTENT(in) :: x
381 INTEGER,
INTENT(in) :: n
382 REAL(
dp),
INTENT(in) :: x1, xn
394 i = int((x - x1)/(xn - x1)*(n - 1)) + 1
397 IF (i > n - 1) i = n - 1
412 FUNCTION iixexp(x, n, x1, alpha, beta)
419 REAL(
dp),
INTENT(in) :: x
420 INTEGER,
INTENT(in) :: n
421 REAL(
dp),
INTENT(in) :: x1, alpha, beta
437 i = int(log((x - x1)/alpha + 1)/beta) + 1
440 IF (i > n - 1) i = n - 1
454 REAL(
dp),
INTENT(in) :: x, c(0:)
463 poly3 = c(1) + c(2)*dx + c(3)*dx**2 + c(4)*dx**3
474 FUNCTION dpoly3(x, c)
476 REAL(
dp),
INTENT(in) :: x, c(0:)
485 dpoly3 = c(2) + 2*c(3)*dx + 3*c(4)*dx**2
496 FUNCTION d2poly3(x, c)
498 REAL(
dp),
INTENT(in) :: x, c(0:)
507 d2poly3 = 2*c(3) + 6*c(4)*dx
516 SUBROUTINE stop_error(msg)
518 CHARACTER(LEN=*) :: msg
Defines the basic variable types.
integer, parameter, public dp
Interface to the LAPACK F77 library.
Simple splines Splines are fully specified by the interpolation points, except that at the ends,...
subroutine, public spline3ders(x, y, xnew, ynew, dynew, d2ynew)
...
real(dp) function, dimension(size(xnew)), public spline3(x, y, xnew)
...