ssa

1
package pchip
2

3
import (
4
	"math"
5
)
6

7
func Pchip2(lgl, x []float64) []float64 {
8

9
	del := make([]float64, len(lgl))
10

11
	slopes := make([]float64, len(lgl))
12
	var d float64
13
	for k := 0; k < len(del)-1; k++ {
14
		del[k] = x[k+1] - x[k]
15
	}
16

17
	for k := 0; k < len(slopes)-1; k++ {
18
		slopes[k+1] = 0.0
19
		if del[k] < 0.0 {
20
			d = del[k+1]
21
			if d <= del[k] {
22
				slopes[k+1] = del[k] / (0.5*(del[k]/d) + 0.5)
23
			} else {
24
				if d < 0.0 {
25
					slopes[k+1] = d / (0.5 + 0.5*(d/del[k]))
26
				}
27
			}
28
		} else {
29
			if del[k] > 0.0 {
30
				d = del[k+1]
31
				if d >= del[k] {
32
					slopes[k+1] = del[k] / (0.5*(del[k]/del[k+1]) + 0.5)
33
				} else {
34
					if d > 0.0 {
35
						slopes[k+1] = del[k+1] / (0.5 + 0.5*(del[k+1]/del[k]))
36
					}
37
				}
38
			}
39
		}
40
	}
41

42
	slopes[0] = exteriorSlope2(del[0], del[1], 1.0, 1.0)
43
	slopes[len(del)-1] = exteriorSlope2(del[len(del)-1], del[len(del)-2], 1.0, 1.0)
44
	return slopes
45
}
46

47
/* Function Definitions */
48
func exteriorSlope2(d1, d2, h1, h2 float64) float64 {
49
	var s float64
50
	var signd1 float64
51
	var signs float64
52
	s = ((2.0*h1+h2)*d1 - h1*d2) / (h1 + h2)
53
	signd1 = d1
54
	if d1 < 0.0 {
55
		signd1 = -1.0
56
	} else if d1 > 0.0 {
57
		signd1 = 1.0
58
	} else {
59
		if d1 == 0.0 {
60
			signd1 = 0.0
61
		}
62
	}
63

64
	signs = s
65
	if s < 0.0 {
66
		signs = -1.0
67
	} else if s > 0.0 {
68
		signs = 1.0
69
	} else {
70
		if s == 0.0 {
71
			signs = 0.0
72
		}
73
	}
74

75
	if signs != signd1 {
76
		s = 0.0
77
	} else {
78
		signs = d2
79
		if d2 < 0.0 {
80
			signs = -1.0
81
		} else if d2 > 0.0 {
82
			signs = 1.0
83
		} else {
84
			if d2 == 0.0 {
85
				signs = 0.0
86
			}
87
		}
88

89
		if (signd1 != signs) && (math.Abs(s) > math.Abs(3.0*d1)) {
90
			s = 3.0 * d1
91
		}
92
	}
93

94
	return s
95
}
96