Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module bipm_db
21 :
22 : use const_def, only: dp
23 :
24 : implicit none
25 :
26 : contains
27 :
28 7 : subroutine do_mkbipm_db(x, nx, y, ny, f1, nf2, ierr)
29 : use interp_1d_def
30 : use interp_1d_lib
31 : integer, intent(in) :: nx ! length of x vector
32 : integer, intent(in) :: ny ! length of y vector
33 : real(dp), intent(in), pointer :: x(:) ! (nx) ! x vector, strict ascending
34 : real(dp), intent(in), pointer :: y(:) ! (ny) ! y vector, strict ascending
35 : integer, intent(in) :: nf2 ! 2nd dimension of f, nf2.ge.nx
36 : real(dp), intent(inout), pointer :: f1(:) ! (4,nf2,ny) ! data & interpolant coefficients
37 : integer, intent(out) :: ierr ! =0 on exit if there is no error.
38 : integer, parameter :: nwork = pm_work_size
39 : integer :: i
40 94 : real(dp), target :: work_ary(nx*nwork)
41 1 : real(dp), pointer :: work(:), f(:)
42 1 : work => work_ary
43 1 : ierr = 0
44 36 : do i = 1, ny
45 35 : f(1:4*nf2) => f1(1 + (i - 1)*4*nf2:i*4*nf2)
46 35 : call interp_pm(x, nx, f, nwork, work, 'do_mkbipm_db', ierr)
47 36 : if (ierr /= 0) exit
48 : end do
49 2 : end subroutine do_mkbipm_db
50 :
51 6 : subroutine do_evbipm_db(xget, yget, x, nx, y, ny, fin1, nf2, f, ierr)
52 1 : use num_lib, only: binary_search
53 : use interp_1d_def
54 : use interp_1d_lib
55 : integer, intent(in) :: nx, ny
56 : real(dp), intent(in) :: xget, yget ! target of this interpolation
57 : real(dp), intent(in), pointer :: x(:) ! (nx) ! ordered x grid
58 : real(dp), intent(in), pointer :: y(:) ! (ny) ! ordered y grid
59 : integer, intent(in) :: nf2
60 : real(dp), intent(in), pointer :: fin1(:) ! fin(4,nf2,ny) ! function data
61 : real(dp), intent(out) :: f
62 : integer, intent(out) :: ierr ! error code =0 ==> no error
63 :
64 : integer, parameter :: nwork = pm_work_size
65 84 : real(dp) :: x0, x1, dx, y0, y1, dy, alfa, beta
66 48 : real(dp) :: ys(4), ynew(1), val(1)
67 : integer :: j, jlo, i, ix, jy, ii
68 180 : real(dp), target :: work_ary(4*nwork), ff_ary(4*4)
69 6 : real(dp), pointer :: work(:), fin(:, :, :), ff(:, :), ff1(:)
70 6 : work => work_ary
71 6 : ff1 => ff_ary
72 6 : ff(1:4, 1:4) => ff_ary(1:4*4)
73 6 : fin(1:4, 1:nf2, 1:ny) => fin1(1:4*nf2*ny)
74 :
75 6 : ierr = 0
76 :
77 6 : ix = binary_search(nx, x, 0, xget) ! x(ix) <= xget < x(ix+1)
78 6 : if (ix < 1 .or. ix >= nx) then
79 0 : ierr = -1
80 0 : return
81 : end if
82 6 : jy = binary_search(ny, y, 0, yget) ! y(jy) <= yget < y(jy+1)
83 6 : if (jy < 1 .or. jy >= ny) then
84 0 : ierr = -1
85 0 : return
86 : end if
87 :
88 6 : x0 = x(ix); x1 = x(ix + 1)
89 6 : y0 = y(jy); y1 = y(jy + 1)
90 6 : dx = xget - x0
91 6 : dy = yget - y0
92 6 : beta = dy/(y1 - y0) ! fraction of y1 result
93 6 : alfa = 1 - beta ! fraction of y0 result
94 :
95 6 : ynew(1) = yget
96 6 : if (jy == 1) then
97 6 : jlo = 1; ; ii = 1
98 6 : else if (jy >= ny - 1) then
99 0 : jlo = jy - 2; ii = 3
100 : else
101 6 : jlo = jy - 1; ii = 2
102 : end if
103 :
104 30 : do i = 1, 4
105 24 : j = jlo + i - 1
106 24 : ys(i) = y(j)
107 30 : ff(1, i) = fin(1, ix, j) + dx*(fin(2, ix, j) + dx*(fin(3, ix, j) + dx*fin(4, ix, j)))
108 : end do
109 :
110 6 : call interp_pm(ys, 4, ff1, nwork, work, 'do_evbipm_db', ierr)
111 6 : if (ierr /= 0) return
112 :
113 6 : call interp_values(ys, 4, ff1, 1, ynew, val, ierr)
114 6 : if (ierr /= 0) return
115 :
116 6 : f = val(1)
117 :
118 12 : end subroutine do_evbipm_db
119 :
120 : end module bipm_db
|
