zeus_python_old/util.src at master · Higginbottom/zeus_python_old · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
	module interpolation

	contains

		function binarysearch(length,array,value,delta)
        ! Given an array and a value, returns the index of the element that
        ! is closest to, but less than, the given value.
        ! Uses a binary search algorithm.
        ! "delta" is the tolerance used to determine if two values are equal
        ! if ( abs(x1 - x2) <= delta) then
        !    assume x1 = x2
        ! endif

        implicit none
        integer, intent(in) :: length
        real*8, dimension(501), intent(in) :: array
        !f2py depend(length) array
        real*8, intent(in) :: value
        real*8, intent(in), optional :: delta

        integer :: binarysearch

        integer :: left, middle, right,n
        real :: d


c		print *,"In binary search length=",length

c		do 30 n=1,length
c			print *,value,n,array(n)
c 30	continue


        if (present(delta) .eqv. .true.) then
            d = delta
        else
            d = 1e-9
        endif

        left = 1
        right = length
        do
            if (left > right) then
                exit
            endif
            middle = nint((left+right) / 2.0)
            if ( abs(array(middle) - value) <= d) then
                binarySearch = middle
                return
            else if (array(middle) > value) then
                right = middle - 1
            else
                left = middle + 1
            end if
        end do
        binarysearch = right

	end function binarysearch

	real function interpolate(x_len,x_array,y_len,y_array,f,x,y,delta)
        ! This function uses bilinear interpolation to estimate the value
        ! of a function f at point (x,y)
        ! f is assumed to be sampled on a regular grid, with the grid x values specified
        ! by x_array and the grid y values specified by y_array
        ! Reference: http://en.wikipedia.org/wiki/Bilinear_interpolation
        implicit none
        integer, intent(in) :: x_len, y_len
        real*8, dimension(501), intent(in) :: x_array
        real*8, dimension(501), intent(in) :: y_array
        real*8, dimension(501,501), intent(in) :: f
        real*8, intent(in) :: x,y
        real*8, intent(in), optional :: delta
        !f2py depend(x_len) x_array, f
        !f2py depend(y_len) y_array, f

        real :: denom, x1, x2, y1, y2
        integer :: i,j,n

c		do 30 n=1,x_lenc
c			print *,'init',n,x_array(n)
c 30		continue

c		print *,"Going to binarysearch"

        i = binarysearch(x_len, x_array, x)
        j = binarysearch(y_len, y_array, y)

c		print *,"Back from binarysearch"


        x1 = x_array(i)
        x2 = x_array(i+1)

        y1 = y_array(j)
        y2 = y_array(j+1)

c		print *,x,y,i,j,f(i,j),f(i+1,j),f(i,j+1),f(i+1,j+1)

        denom = (x2 - x1)*(y2 - y1)

        interpolate = (f(i,j)*(x2-x)*(y2-y) + f(i+1,j)*(x-x1)*(y2-y) +
     +   f(i,j+1)*(x2-x)*(y-y1) + f(i+1, j+1)*(x-x1)*(y-y1))/denom

c		print *,x1,x,x2,y1,y,y2,denom,interpolate


		end function interpolate

	end module interpolation