8
func exteriorSlope(d1, d2, h1, h2 float64) float64 {
12
s = ((2.0*h1+h2)*d1 - h1*d2) / (h1 + h2)
46
if (signd1 != signs) && (math.Abs(s) > math.Abs(3.0*d1)) {
53
type PchipCoefs struct {
57
func Pchip(x, y, new_x []float64, x_len, new_x_len int) ([]float64, []float64, PchipCoefs) {
59
new_y := make([]float64, new_x_len)
62
del := make([]float64, x_len-1)
63
slopes := make([]float64, x_len)
64
h := make([]float64, x_len-1)
69
pp_coefs := make([]float64, (x_len-1)+(3*(x_len-1)))
74
for low_ip1 := 0; low_ip1 < x_len-1; low_ip1++ {
75
hs = x[low_ip1+1] - x[low_ip1]
76
del[low_ip1] = (y[low_ip1+1] - y[low_ip1]) / hs
80
for low_ip1 := 0; low_ip1 < x_len-2; low_ip1++ {
81
hs = h[low_ip1] + h[low_ip1+1]
83
w1 = (h[low_ip1] + hs) / hs3
84
hs = (h[low_ip1+1] + hs) / hs3
86
if del[low_ip1] < 0.0 {
87
if del[low_ip1+1] <= del[low_ip1] {
88
hs3 = del[low_ip1] / (w1*(del[low_ip1]/del[low_ip1+1]) + hs)
90
if del[low_ip1+1] < 0.0 {
91
hs3 = del[low_ip1+1] / (w1 + hs*(del[low_ip1+1]/del[low_ip1]))
95
if del[low_ip1] > 0.0 {
96
if del[low_ip1+1] >= del[low_ip1] {
97
hs3 = del[low_ip1] / (w1*(del[low_ip1]/del[low_ip1+1]) + hs)
99
if del[low_ip1+1] > 0.0 {
100
hs3 = del[low_ip1+1] / (w1 + hs*(del[low_ip1+1]/del[low_ip1]))
106
slopes[low_ip1+1] = hs3
109
return new_y, pp_coefs, PchipCoefs{}
111
fmt.Println("x_len", x_len)
113
slopes[0] = exteriorSlope(del[0], del[1], h[0], h[1])
114
slopes[x_len-1] = exteriorSlope(del[x_len-2], del[x_len-3], h[x_len-2], h[x_len-3])
115
for low_ip1 := 0; low_ip1 < x_len-1; low_ip1++ {
116
hs = (del[low_ip1] - slopes[low_ip1]) / h[low_ip1]
117
hs3 = (slopes[low_ip1+1] - del[low_ip1]) / h[low_ip1]
118
pp_coefs[low_ip1] = (hs3 - hs) / h[low_ip1]
119
pp_coefs[low_ip1+x_len-1] = 2.0*hs - hs3
120
pp_coefs[low_ip1+(2*(x_len-1))] = slopes[low_ip1]
121
pp_coefs[low_ip1+(3*(x_len-1))] = y[low_ip1]
124
for ix := 0; ix < new_x_len; ix++ {
128
for high_i > low_ip1 {
129
mid_i = ((low_i + high_i) + 1) >> 1
130
if new_x[ix] >= x[mid_i-1] {
138
hs = new_x[ix] - x[low_i]
139
hs3 = pp_coefs[low_i]
140
for low_ip1 := 0; low_ip1 < 3; low_ip1++ {
141
hs3 = hs*hs3 + pp_coefs[low_i+(low_ip1+1)*(x_len-1)]
155
lenCoefs := len(pp_coefs) / 4
156
coefs.A = make([]float64, lenCoefs)
157
coefs.B = make([]float64, lenCoefs)
158
coefs.C = make([]float64, lenCoefs)
159
coefs.D = make([]float64, lenCoefs)
160
for i := 0; i < lenCoefs; i++ {
161
coefs.A[i] = pp_coefs[i]
162
coefs.B[i] = pp_coefs[1*(len(pp_coefs)/4)+i]
163
coefs.C[i] = pp_coefs[2*(len(pp_coefs)/4)+i]
164
coefs.D[i] = pp_coefs[3*(len(pp_coefs)/4)+i]
167
return new_y, pp_coefs, coefs