-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathS3M_Module_Tools_Interp.f90
104 lines (80 loc) · 3.46 KB
/
S3M_Module_Tools_Interp.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
!------------------------------------------------------------------------------------------
! File: S3M_Module_Tools_Interp.f90
! Author(s): Fabio Delogu.
!
! Created on April 28, 2017, 4:12 PM
!
! Module to define interpolation tools
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Module header
module S3M_Module_Tools_Interp
!------------------------------------------------------------------------------------------
! External module(s) and implicit none
use S3M_Module_Tools_Debug
implicit none
!------------------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------------------
! Subroutine to space linearly an array to two numbers
subroutine linspace(xmin,xmax,x)
real(kind = 4), intent(in) :: xmin,xmax
real(kind = 4), intent(out) :: x(:)
integer(kind = 4) :: i,n
n = size(x)
if (n == 1) then
if(xmin /= xmax) then
call mprintf(.true., iERROR, 'Cannot call linspace with n=1 and xmin /= xmax')
stop
else
x = xmin
end if
else
do i=1,n
x(i) = (xmax-xmin) * real(i-1,kind = 4) / real(n-1,kind = 4) + xmin
end do
end if
end subroutine linspace
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Subroutine to evaluate the nearest neighbor interpolant
subroutine nearest_interp_1d ( nd, xd, yd, ni, xi, yi, ki)
! The nearest neighbor interpolant L(ND,XD,YD)(X) is the piecewise
! constant function which interpolates the data (XD(I),YD(I)) for I = 1
! to ND.
! Input, integer ND, the number of data points.
! ND must be at least 1.
! Input, double precision XD(ND), the data points.
! Input, double precision YD(ND), the data values.
! Input, integer NI, the number of interpolation points.
! Input, double precision XI(NI), the interpolation points.
! Output, double precision YI(NI), the interpolated values.
integer(kind = 4) :: nd
integer(kind = 4) :: ni
real(kind = 4) :: d
real(kind = 4) :: d2
integer(kind = 4) :: i
integer(kind = 4) :: j
integer(kind = 4) :: k
real(kind = 4) :: xd(nd)
real(kind = 4) :: xi(ni)
real(kind = 4) :: yd(nd)
real(kind = 4) :: yi(ni)
integer(kind = 4) :: ki(ni)
do i = 1, ni
k = 1
d = abs ( xi(i) - xd(k) )
do j = 2, nd
d2 = abs ( xi(i) - xd(j) )
if ( d2 .lt. d ) then
k = j
d = d2
end if
end do
yi(i) = yd(k)
ki(i) = k
end do
end subroutine nearest_interp_1d
!------------------------------------------------------------------------------------
end module S3M_Module_Tools_Interp
!------------------------------------------------------------------------------------------