MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
interpolation.f90
Go to the documentation of this file.
1 
7 module interpolation
8  use constants, only: dp
9  implicit none
10  private
11  public pwl_interp_1d
12 contains
13 !!********************************BEGINNING*************************************
14 !!minimal working example (you have to supply xy values in "refl.dat" and change value of xi)
15 !subroutine test_pwl_interp_1d
16 ! integer :: fi=154 !file index
17 ! integer :: i !counter
18 ! integer, parameter :: nd=101 !number of data points
19 ! real(dp), dimension(nd) :: xd,yd !data for interpolation
20 ! integer :: ni=1 !number of points, where we want to interpolate
21 ! real(dp) :: xi(1) !x-values of points, where we want to interpolate
22 ! real(dp) :: yi(1) !interpolated y-values
23 ! open(fi,file="refl.dat")
24 ! do i=1,nd
25 ! read(fi,*) xd(i),yd(i)
26 ! enddo
27 ! close(fi)
28 ! xi=2e-6_dp
29 ! call pwl_interp_1d ( nd, xd, yd, ni, xi, yi )
30 ! write(*,*) xi,yi
31 !end subroutine test_pwl_interp_1d
32 !!***********************************END****************************************
33 
34 
35 subroutine pwl_interp_1d ( nd, xd, yd, ni, xi, yi )
36 
37 !*****************************************************************************80
38 !
39 !! PWL_INTERP_1D evaluates the piecewise linear interpolant.
40 !
41 ! Discussion:
42 !
43 ! The piecewise linear interpolant L(ND,XD,YD)(X) is the piecewise
44 ! linear function which interpolates the data (XD(I),YD(I)) for I = 1
45 ! to ND.
46 !
47 ! Licensing:
48 !
49 ! This code is distributed under the GNU LGPL license.
50 !
51 ! Modified:
52 !
53 ! 22 September 2012
54 !
55 ! Author:
56 !
57 ! John Burkardt
58 !
59 ! Parameters:
60 !
61 ! Input, integer ( kind = 4 ) ND, the number of data points.
62 ! ND must be at least 1.
63 !
64 ! Input, real ( kind = 8 ) XD(ND), the data points.
65 !
66 ! Input, real ( kind = 8 ) YD(ND), the data values.
67 !
68 ! Input, integer ( kind = 4 ) NI, the number of interpolation points.
69 !
70 ! Input, real ( kind = 8 ) XI(NI), the interpolation points.
71 !
72 ! Output, real ( kind = 8 ) YI(NI), the interpolated values.
73 !
74  implicit none
75 
76  integer ( kind = 4 ) nd
77  integer ( kind = 4 ) ni
78 
79  integer ( kind = 4 ) i
80  integer ( kind = 4 ) k
81  real ( kind = dp ) t
82  real ( kind = dp ) xd(nd)
83  real ( kind = dp ) yd(nd)
84  real ( kind = dp ) xi(ni)
85  real ( kind = dp ) yi(ni)
86 
87  yi(1:ni) = 0.0e+00_dp
88 
89  if ( nd == 1 ) then
90  yi(1:ni) = yd(1)
91  return
92  end if
93 
94  do i = 1, ni
95 
96  if ( xi(i) <= xd(1) ) then
97 
98  t = ( xi(i) - xd(1) ) / ( xd(2) - xd(1) )
99  yi(i) = ( 1.0e+00_dp - t ) * yd(1) + t * yd(2)
100 
101  else if ( xd(nd) <= xi(i) ) then
102 
103  t = ( xi(i) - xd(nd-1) ) / ( xd(nd) - xd(nd-1) )
104  yi(i) = ( 1.0e+00_dp - t ) * yd(nd-1) + t * yd(nd)
105 
106  else
107 
108  do k = 2, nd
109 
110  if ( xd(k-1) <= xi(i) .and. xi(i) <= xd(k) ) then
111 
112  t = ( xi(i) - xd(k-1) ) / ( xd(k) - xd(k-1) )
113  yi(i) = ( 1.0e+00_dp - t ) * yd(k-1) + t * yd(k)
114  exit
115 
116  end if
117 
118  end do
119 
120  end if
121 
122  end do
123 
124  return
125 end subroutine pwl_interp_1d
126 end module interpolation