trunk/libf/regr1_lint.f90

module regr1_lint_m

  ! Author: Lionel GUEZ

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