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_sg
21 :
22 : implicit none
23 :
24 : contains
25 :
26 1 : subroutine do_mkbipm_sg(x, nx, y, ny, f1, nf2, ierr)
27 : use interp_1d_def
28 : use interp_1d_lib_sg
29 : integer, intent(in) :: nx ! length of x vector
30 : integer, intent(in) :: ny ! length of y vector
31 : real, intent(in) :: x(:) ! (nx) ! x vector, strict ascending
32 : real, intent(in) :: y(:) ! (ny) ! y vector, strict ascending
33 : integer, intent(in) :: nf2 ! 2nd dimension of f, nf2.ge.nx
34 : real, intent(inout), pointer :: f1(:) ! (4,nf2,ny) ! data & interpolant coefficients
35 : integer, intent(out) :: ierr ! =0 on exit if there is no error.
36 : integer, parameter :: nwork = pm_work_size
37 : integer :: i
38 94 : real, target :: work_ary(nx*nwork)
39 1 : real, pointer :: work(:), f(:)
40 1 : work => work_ary
41 1 : ierr = 0
42 36 : do i = 1, ny
43 35 : f(1:4*nf2) => f1(1 + (i - 1)*4*nf2:i*4*nf2)
44 35 : call interp_pm_sg(x, nx, f, nwork, work, 'do_mkbipm_sg', ierr)
45 36 : if (ierr /= 0) exit
46 : end do
47 1 : end subroutine do_mkbipm_sg
48 :
49 6 : subroutine do_evbipm_sg(xget, yget, x, nx, y, ny, fin1, nf2, f, ierr)
50 : use num_lib, only: binary_search_sg
51 : use interp_1d_def
52 : use interp_1d_lib_sg
53 : integer, intent(in) :: nx, ny
54 : real, intent(in) :: xget, yget ! target of this interpolation
55 : real, intent(in), pointer :: x(:) ! (nx) ! ordered x grid
56 : real, intent(in), pointer :: y(:) ! (ny) ! ordered y grid
57 : integer, intent(in) :: nf2
58 : real, intent(in), pointer :: fin1(:) ! fin(4,nf2,ny) ! function data
59 : real, intent(out) :: f
60 : integer, intent(out) :: ierr ! error code =0 ==> no error
61 :
62 : integer, parameter :: nwork = pm_work_size
63 84 : real :: x0, x1, dx, y0, y1, dy, alfa, beta
64 48 : real :: ys(4), ynew(1), val(1)
65 : integer :: j, jlo, i, ix, jy, ii
66 180 : real, target :: work_ary(4*nwork), ff_ary(4*4)
67 6 : real, pointer :: work(:), fin(:, :, :), ff(:, :), ff1(:)
68 6 : work => work_ary
69 6 : ff1 => ff_ary
70 6 : ff(1:4, 1:4) => ff_ary(1:4*4)
71 6 : fin(1:4, 1:nf2, 1:ny) => fin1(1:4*nf2*ny)
72 :
73 6 : ierr = 0
74 :
75 6 : ix = binary_search_sg(nx, x, 0, xget) ! x(ix) <= xget < x(ix+1)
76 6 : if (ix < 1 .or. ix >= nx) then
77 0 : ierr = -1
78 0 : return
79 : end if
80 6 : jy = binary_search_sg(ny, y, 0, yget) ! y(jy) <= yget < y(jy+1)
81 6 : if (jy < 1 .or. jy >= ny) then
82 0 : ierr = -1
83 0 : return
84 : end if
85 :
86 6 : x0 = x(ix); x1 = x(ix + 1)
87 6 : y0 = y(jy); y1 = y(jy + 1)
88 6 : dx = xget - x0
89 6 : dy = yget - y0
90 6 : beta = dy/(y1 - y0) ! fraction of y1 result
91 6 : alfa = 1 - beta ! fraction of y0 result
92 :
93 6 : ynew(1) = yget
94 6 : if (jy == 1) then
95 6 : jlo = 1; ; ii = 1
96 6 : else if (jy >= ny - 1) then
97 0 : jlo = jy - 2; ii = 3
98 : else
99 6 : jlo = jy - 1; ii = 2
100 : end if
101 :
102 30 : do i = 1, 4
103 24 : j = jlo + i - 1
104 24 : ys(i) = y(j)
105 30 : ff(1, i) = fin(1, ix, j) + dx*(fin(2, ix, j) + dx*(fin(3, ix, j) + dx*fin(4, ix, j)))
106 : end do
107 :
108 6 : call interp_pm_sg(ys, 4, ff1, nwork, work, 'do_evbipm_sg', ierr)
109 6 : if (ierr /= 0) return
110 :
111 6 : call interp_values_sg(ys, 4, ff1, 1, ynew, val, ierr)
112 6 : if (ierr /= 0) return
113 :
114 6 : f = val(1)
115 :
116 12 : end subroutine do_evbipm_sg
117 :
118 : end module bipm_sg
|
