trunk/libf/regr1_lint.f90

module regr1_lint_m

  implicit none

  interface regr1_lint
     ! Each procedure regrids by linear interpolation.
     ! The regridding operation is done on the first dimension of the
     ! input array.
     ! The difference betwwen the procedures is the rank of the first argument.
     module procedure regr11_lint, regr12_lint
  end interface

  private
  public regr1_lint

contains

  function regr11_lint(vs, xs, xt) result(vt)

    ! "vs" has rank 1.

    use nrutil, only: assert_eq
    use interpolation, only: hunt, polint

    real, intent(in):: vs(:)
    ! (values of the function at source points "xs")

    real, intent(in):: xs(:)
    ! (abscissas of points in source grid, in strictly monotonic order)

    real, intent(in):: xt(:)
    ! (abscissas of points in target grid)

    real vt(size(xt)) ! values of the function on the target grid

    ! Variables local to the procedure:
    integer is, it, ns
    integer is_b ! "is" bound between 1 and "ns - 1"

    !--------------------------------------

    ns = assert_eq(size(vs), size(xs), "regr11_lint ns")

    is = -1 ! go immediately to bisection on first call to "hunt"

    do it = 1, size(xt)
       call hunt(xs, xt(it), is)
       is_b = min(max(is, 1), ns - 1)
       call polint(xs(is_b:is_b+1), vs(is_b:is_b+1), xt(it), vt(it))
!!$       vt(it) = ((xs(is_b+1) - xt(it)) * vs(is_b) &
!!$            + (xt(it) - xs(is_b)) * vs(is_b+1)) / (xs(is_b+1) - xs(is_b))
    end do

  end function regr11_lint

  !*********************************************************

  function regr12_lint(vs, xs, xt) result(vt)

    ! "vs" has rank 2.

    use nrutil, only: assert_eq
    use interpolation, only: hunt

    real, intent(in):: vs(:, :)
    ! (values of the function at source points "xs")

    real, intent(in):: xs(:)
    ! (abscissas of points in source grid, in strictly monotonic order)

    real, intent(in):: xt(:)
    ! (abscissas of points in target grid)

    real vt(size(xt), size(vs, 2)) ! values of the function on the target grid

    ! Variables local to the procedure:
    integer is, it, ns
    integer is_b ! "is" bound between 1 and "ns - 1"

    !--------------------------------------

    ns = assert_eq(size(vs, 1), size(xs), "regr12_lint ns")

    is = -1 ! go immediately to bisection on first call to "hunt"

    do it = 1, size(xt)
       call hunt(xs, xt(it), is)
       is_b = min(max(is, 1), ns - 1)
       vt(it, :) = ((xs(is_b+1) - xt(it)) * vs(is_b, :) &
            + (xt(it) - xs(is_b)) * vs(is_b+1, :)) / (xs(is_b+1) - xs(is_b))
    end do

  end function regr12_lint

end module regr1_lint_m
1	guez	9	module regr1_lint_m
2
3			implicit none
4
5			interface regr1_lint
6			! Each procedure regrids by linear interpolation.
7			! The regridding operation is done on the first dimension of the
8			! input array.
9			! The difference betwwen the procedures is the rank of the first argument.
10			module procedure regr11_lint, regr12_lint
11			end interface
12
13			private
14			public regr1_lint
15
16			contains
17
18			function regr11_lint(vs, xs, xt) result(vt)
19
20			! "vs" has rank 1.
21
22			use nrutil, only: assert_eq
23			use interpolation, only: hunt, polint
24
25			real, intent(in):: vs(:)
26			! (values of the function at source points "xs")
27
28			real, intent(in):: xs(:)
29			! (abscissas of points in source grid, in strictly monotonic order)
30
31			real, intent(in):: xt(:)
32			! (abscissas of points in target grid)
33
34			real vt(size(xt)) ! values of the function on the target grid
35
36			! Variables local to the procedure:
37			integer is, it, ns
38			integer is_b ! "is" bound between 1 and "ns - 1"
39
40			!--------------------------------------
41
42			ns = assert_eq(size(vs), size(xs), "regr11_lint ns")
43
44			is = -1 ! go immediately to bisection on first call to "hunt"
45
46			do it = 1, size(xt)
47			call hunt(xs, xt(it), is)
48			is_b = min(max(is, 1), ns - 1)
49			call polint(xs(is_b:is_b+1), vs(is_b:is_b+1), xt(it), vt(it))
50			!!$ vt(it) = ((xs(is_b+1) - xt(it)) * vs(is_b) &
51			!!$ + (xt(it) - xs(is_b)) * vs(is_b+1)) / (xs(is_b+1) - xs(is_b))
52			end do
53
54			end function regr11_lint
55
56			!*********************************************************
57
58			function regr12_lint(vs, xs, xt) result(vt)
59
60			! "vs" has rank 2.
61
62			use nrutil, only: assert_eq
63			use interpolation, only: hunt
64
65			real, intent(in):: vs(:, :)
66			! (values of the function at source points "xs")
67
68			real, intent(in):: xs(:)
69			! (abscissas of points in source grid, in strictly monotonic order)
70
71			real, intent(in):: xt(:)
72			! (abscissas of points in target grid)
73
74			real vt(size(xt), size(vs, 2)) ! values of the function on the target grid
75
76			! Variables local to the procedure:
77			integer is, it, ns
78			integer is_b ! "is" bound between 1 and "ns - 1"
79
80			!--------------------------------------
81
82			ns = assert_eq(size(vs, 1), size(xs), "regr12_lint ns")
83
84			is = -1 ! go immediately to bisection on first call to "hunt"
85
86			do it = 1, size(xt)
87			call hunt(xs, xt(it), is)
88			is_b = min(max(is, 1), ns - 1)
89			vt(it, :) = ((xs(is_b+1) - xt(it)) * vs(is_b, :) &
90			+ (xt(it) - xs(is_b)) * vs(is_b+1, :)) / (xs(is_b+1) - xs(is_b))
91			end do
92
93			end function regr12_lint
94
95			end module regr1_lint_m