12 logical:: initialized=.false.
14 procedure interpolation4
15 procedure interpolation8
16 generic:: interpolation => interpolation4, interpolation8
19 generic:: sequence => sequence4, sequence8
22 generic:: interval => interval4, interval8
25 generic:: weights => weights4, weights8
26 procedure deriv_weights4
27 procedure deriv_weights8
28 generic:: deriv_weights => deriv_weights4, deriv_weights8
29 procedure deriv2_weights4
30 procedure deriv2_weights8
31 generic:: deriv2_weights => deriv2_weights4, deriv2_weights8
32 procedure sequence_weights4
33 procedure sequence_weights8
34 generic:: sequence_weights => sequence_weights4, sequence_weights8
35 procedure deriv_sequence_weights4
36 procedure deriv_sequence_weights8
37 generic:: deriv_sequence_weights => deriv_sequence_weights4, deriv_sequence_weights8
38 procedure deriv2_sequence_weights4
39 procedure deriv2_sequence_weights8
40 generic:: deriv2_sequence_weights => deriv2_sequence_weights4, deriv2_sequence_weights8
41 procedure interpolate3d
46 procedure interpolate1d4
47 procedure interpolate1d8
48 generic:: interpolate1d => interpolate1d4, interpolate1d8
56 FUNCTION interpolation4 (self, x0, x, y)
RESULT(f)
58 real(kind=4),
dimension(:):: x, y
59 real(kind=4):: x0, f, g
67 g = merge(g,g*(x0-x(j))/(x(i)-x(j)),j==i)
71 END FUNCTION interpolation4
72 FUNCTION interpolation8 (self, x0, x, y)
RESULT(f)
74 real(kind=8),
dimension(:):: x, y
75 real(kind=8):: x0, f, g
83 g = merge(g,g*(x0-x(j))/(x(i)-x(j)),j==i)
87 END FUNCTION interpolation8
94 FUNCTION weights4 (self, x0, x)
RESULT(w)
96 real(kind=4),
dimension(:):: x
97 real(kind=4),
dimension(size(x)):: w
111 w(i) = w(i)*(x0-x(j))/(x(i)-x(j))
114 END FUNCTION weights4
115 FUNCTION weights8 (self, x0, x)
RESULT(w)
117 real(kind=8),
dimension(:):: x
118 real(kind=8),
dimension(size(x)):: w
132 w(i) = w(i)*(x0-x(j))/(x(i)-x(j))
135 END FUNCTION weights8
142 FUNCTION deriv_weights4 (self, x0, x)
RESULT(w)
144 real(kind=4),
dimension(:):: x
145 real(kind=4),
dimension(size(x)):: w
146 real(kind=4):: x0, nomin, denom, g
161 if (k==j .or. k==i) cycle
165 denom = denom*(x(i)-x(j))
169 END FUNCTION deriv_weights4
170 FUNCTION deriv_weights8 (self, x0, x)
RESULT(w)
172 real(kind=8),
dimension(:):: x
173 real(kind=8),
dimension(size(x)):: w
174 real(kind=8):: x0, nomin, denom, g
189 if (k==j .or. k==i) cycle
193 denom = denom*(x(i)-x(j))
197 END FUNCTION deriv_weights8
204 FUNCTION deriv2_weights4 (self, x0, x)
RESULT(w)
206 real(kind=4),
dimension(:):: x
207 real(kind=4),
dimension(size(x)):: w
208 real(kind=4):: x0, nomin, denom, g
209 integer:: n, i, j, k, l
222 if (k==j .or. k==i) cycle
225 if (l==k .or. l==j .or. l==i) cycle
230 denom = denom*(x(i)-x(j))
234 END FUNCTION deriv2_weights4
235 FUNCTION deriv2_weights8 (self, x0, x)
RESULT(w)
237 real(kind=8),
dimension(:):: x
238 real(kind=8),
dimension(size(x)):: w
239 real(kind=8):: x0, nomin, denom, g
240 integer:: n, i, j, k, l
253 if (k == j .or. k == i) cycle
256 if (l == k .or. l == j .or. l == i) cycle
261 denom = denom*(x(i)-x(j))
265 END FUNCTION deriv2_weights8
271 SUBROUTINE interval4 (self, x0, x, i0, i2, o, order)
273 integer,
optional:: order
274 real(kind=4),
dimension(:):: x
275 real(kind=4),
dimension(size(x)):: y
277 integer:: i0, i1, i2, o, n
282 if (
present(order)) o = order
286 do while (x(i1)==x(1) .and. i1<n)
290 if (x(i1)==x(1))
then 303 do while (x(i1)<x0 .and. i1<n)
306 if (i1 < 3 .or. i1 > n-1)
then 335 if (self%verbose>0)
then 336 print*,
'lagrange_interval: i0, i2, o =', i0, i2, o
338 END SUBROUTINE interval4
339 SUBROUTINE interval8 (self, x0, x, i0, i2, o, order)
341 integer,
optional:: order
342 real(kind=8),
dimension(:):: x
343 real(kind=8),
dimension(size(x)):: y
345 integer:: i0, i1, i2, o, n
350 if (
present(order)) o = order
354 do while (x(i1)==x(1) .and. i1<n)
358 if (x(i1)==x(1))
then 371 do while (x(i1)<x0 .and. i1<n)
374 if (i1 < 3 .or. i1 > n-1)
then 403 if (self%debug(1))
then 404 print
'("lagrange_interval: i0,i2:",2i3,3x,"o:",i2,3x,"x0:",f10.6,3x,"x:",10f10.6)', &
407 END SUBROUTINE interval8
413 FUNCTION sequence4 (self, x0, x, y, order, correct)
RESULT(f)
415 integer,
optional:: order
416 real(kind=4),
optional:: correct
417 real(kind=4),
dimension(:):: x, y
419 integer:: i0, i2, o, n
422 call self%interval (x0, x, i0, i2, o, order)
423 if (self%verbose>0)
then 424 print
'(" i0,i2,order=",3i3,4x,"x0,x,y =",1p,18g12.4)', &
425 i0, i2, o, x0, x(i0:i2), y(i0:i2)
427 f = self%interpolation (x0, x(i0:i2), y(i0:i2))
428 if (self%verbose>0)
then 429 if (
present(correct))
then 430 print
'(" i0,i2,order=",3i3,4x,"x0,x,y,f,correct=",1p,18g12.4)', &
431 i0, i2, o, x0, x(i0:i2), y(i0:i2), f, correct
433 print
'(" i0,i2,order=",3i3,4x,"x0,x,y,f=",1p,18g12.4)', &
434 i0, i2, o, x0, x(i0:i2), y(i0:i2), f
437 END FUNCTION sequence4
438 FUNCTION sequence8 (self, x0, x, y, order, correct)
RESULT(f)
440 integer,
optional:: order
441 real(kind=8),
optional:: correct
442 real(kind=8),
dimension(:):: x, y
444 integer:: i0, i2, o, n
447 call self%interval (x0, x, i0, i2, o, order)
448 if (self%verbose>0)
then 449 print
'(" i0,i2,order=",3i3,4x,"x0,x,y =",1p,18g12.4)', &
450 i0, i2, o, x0, x(i0:i2), y(i0:i2)
452 f = self%interpolation (x0, x(i0:i2), y(i0:i2))
453 if (self%verbose>0)
then 454 if (
present(correct))
then 455 print
'(" i0,i2,order=",3i3,4x,"x0,x,y,f,correct=",1p,18g12.4)', &
456 i0, i2, o, x0, x(i0:i2), y(i0:i2), f, correct
458 print
'(" i0,i2,order=",3i3,4x,"x0,x,y,f=",1p,18g12.4)', &
459 i0, i2, o, x0, x(i0:i2), y(i0:i2), f
462 END FUNCTION sequence8
468 SUBROUTINE sequence_weights4 (self, x0, x, i0, i2, w, order)
470 integer,
optional:: order
471 real(kind=4),
dimension(:):: x, w
473 integer:: i0, i2, o, n
476 call self%interval (x0, x, i0, i2, o, order)
481 w(1:i2-i0+1) = self%weights (x0, x(i0:i2))
483 END SUBROUTINE sequence_weights4
484 SUBROUTINE sequence_weights8 (self, x0, x, i0, i2, w, order)
486 integer,
optional:: order
487 real(kind=8),
dimension(:):: x, w
489 integer:: i0, i2, o, n
492 call self%interval (x0, x, i0, i2, o, order)
497 w(1:i2-i0+1) = self%weights (x0, x(i0:i2))
499 END SUBROUTINE sequence_weights8
505 SUBROUTINE deriv_sequence_weights4 (self, x0, x, i0, i2, w, order)
507 integer,
optional:: order
508 real(kind=4),
dimension(:):: x, w
510 integer:: i0, i2, o, n
513 call self%interval (x0, x, i0, i2, o, order)
515 w(1:i2-i0+1) = self%deriv_weights (x0, x(i0:i2))
516 END SUBROUTINE deriv_sequence_weights4
517 SUBROUTINE deriv_sequence_weights8 (self, x0, x, i0, i2, w, order)
519 integer,
optional:: order
520 real(kind=8),
dimension(:):: x, w
522 integer:: i0, i2, o, n
525 call self%interval (x0, x, i0, i2, o, order)
527 w(1:i2-i0+1) = self%deriv_weights (x0, x(i0:i2))
528 END SUBROUTINE deriv_sequence_weights8
534 SUBROUTINE deriv2_sequence_weights4 (self, x0, x, i0, i2, w, order)
536 integer,
optional:: order
537 real(kind=4),
dimension(:):: x, w
539 integer:: i0, i2, o, n
542 call self%interval (x0, x, i0, i2, o, order)
544 w(1:i2-i0+1) = self%deriv2_weights (x0, x(i0:i2))
545 END SUBROUTINE deriv2_sequence_weights4
546 SUBROUTINE deriv2_sequence_weights8 (self, x0, x, i0, i2, w, order)
548 integer,
optional:: order
549 real(kind=8),
dimension(:):: x, w
551 integer:: i0, i2, o, n
554 call self%interval (x0, x, i0, i2, o, order)
556 w(1:i2-i0+1) = self%deriv2_weights (x0, x(i0:i2))
557 END SUBROUTINE deriv2_sequence_weights8
562 SUBROUTINE interpolate3d (self, xi, yi, zi, fi, xo, yo, zo, fo, order)
564 real,
dimension(:):: xi, yi, zi, xo, yo, zo
565 real,
dimension(:,:,:):: fi, fo
568 real,
dimension(:,:,:),
pointer:: f1, f2
569 integer:: nix, niy, niz, nox, noy, noz, ix, iy, iz
577 allocate (f1(nix,niy,noz))
581 f1(ix,iy,iz) = self%sequence (zo(iz), zi, fi(ix,iy,:), order)
585 allocate (f2(nix,noy,noz))
589 f2(ix,iy,iz) = self%sequence (yo(iy), yi, f1(ix,:,iz), order)
597 fo(ix,iy,iz) = self%sequence (xo(ix), xi, f2(:,iy,iz), order)
602 END SUBROUTINE interpolate3d
607 SUBROUTINE trilinear3d (self, xi, yi, zi, fi, xo, yo, zo, fo, order)
609 real,
dimension(:):: xi, yi, zi, xo, yo, zo
610 real,
dimension(:,:,:):: fi, fo
613 integer:: nix, niy, niz, nox, noy, noz, ix, iy, iz, jx ,jy, jz
614 real:: px, py, pz, qx, qy, qz, dx, dy, dz
622 dx = (xi(nix)-xi(1))/(nix-1)
623 dy = (yi(niy)-yi(1))/(niy-1)
624 dz = (zi(niz)-zi(1))/(niz-1)
627 if (zi(jz+1)>zo(iz))
exit 629 pz = (zo(iz)-zi(jz))/(zi(jz+1)-zi(jz))
630 pz = max(0.0,min(1.0,pz))
633 py = (yo(iy)-yi(1))/dy + 1.
635 jy = max(1,min(niy-1,jy))
636 py = max(0.0,min(1.0,py -jy))
639 px = (xo(ix)-xi(1))/dx + 1.
641 jx = max(1,min(nix-1,jx))
642 px = max(0.0,min(1.0,px -jx))
645 qz*(qy*(qx*fi(jx ,jy , jz ) + px*fi(jx+1,jy , jz )) + &
646 py*(qx*fi(jx ,jy+1, jz ) + px*fi(jx+1,jy+1, jz ))) + &
647 pz*(qy*(qx*fi(jx ,jy , jz+1) + px*fi(jx+1,jy , jz+1)) + &
648 py*(qx*fi(jx ,jy+1, jz+1) + px*fi(jx+1,jy+1, jz+1)))
652 END SUBROUTINE trilinear3d
657 SUBROUTINE test (self)
659 integer,
parameter:: n=5
660 integer:: order, i, i0, i1
661 real(8),
dimension(n):: x, y, w
664 if (.not. io_unit%master)
return 665 print *,
'--------------- lagrange_t%test ------------------' 669 print
'(5x,a,11x,a,11x,a,4x,a)',
'x',
'y',
'y==1',
'order' 672 call self%sequence_weights (x0, x, i0, i1, w, order)
673 y0 = sum(y(i0:i1)*w(1:1+i1-i0))
674 print
'(3g12.4,i6)', x0, y0, x0**order, i1-i0
679 print
'(5x,a,11x,a,11x,a,4x,a)',
'x',
'y',
'x-4',
'order' 682 call self%sequence_weights (x0, x, i0, i1, w, order)
683 y0 = sum(y(i0:i1)*w(1:1+i1-i0))
684 print
'(3g12.4,i6)', x0, y0, x0**order, i1-i0
691 print
'(5x,a,11x,a,11x,a,4x,a)',
'x',
'y',
'x**2',
'order' 694 call self%sequence_weights (x0, x, i0, i1, w, order)
695 y0 = sum(y(i0:i1)*w(1:1+i1-i0))
696 print
'(3g12.4,i6)', x0, y0, x0**order, i1-i0
698 print *,
'--------------------------------------------------' 703 FUNCTION debug (self, verbose)
708 if (.not.self%initialized)
call self%init
709 debug = (self%verbose >= verbose)
715 SUBROUTINE init (self)
718 integer:: verbose=0, order=2
719 namelist /lagrange_params/ verbose, order
721 call trace%begin (
'lagrange_t%init')
723 if (.not.self%initialized)
then 724 self%initialized = .true.
725 rewind(io_unit%input)
726 read (io_unit%input, lagrange_params, iostat=iostat)
727 if (io_unit%master)
write (io_unit%output, lagrange_params)
728 self%verbose = verbose
738 FUNCTION interpolate1d4 (self,x0, n, x, y)
RESULT(f)
741 real(kind=4):: x0, x(n), y(n), f
745 else if (x0 > x(n))
then 748 f = self%sequence (x0, x, y)
750 END FUNCTION interpolate1d4
752 FUNCTION interpolate1d8 (self,x0, n, x, y)
RESULT(f)
755 real(kind=8):: x0, x(n), y(n), f
759 else if (x0 > x(n))
then 762 f = self%sequence (x0, x, y)
764 END FUNCTION interpolate1d8
Module for Lagrange interpolation.