github.com/archlabjp/eeslism-go@v0.0.0-20231109122333-4bb7bfcdf292/eeslism/spline.go (about)

     1  package eeslism
     2  
     3  // 給水温度のスプライン補間
     4  
     5  // 各月の給水温度は15日の給水温度として補間する
     6  // Nday:計算対象の通日、Tsupw:1～12月の給水温度
     7  var __Intgtsup_ic int = 0
     8  
     9  func Intgtsup(Nday int, Tsupw []float64) float64 {
    10  	var h, b, d, g, u [13]float64
    11  	var r, x, y [14]float64
    12  	var n, Mo int
    13  	var y1 float64
    14  
    15  	if __Intgtsup_ic == 0 {
    16  		n = 13
    17  		x[0], x[13] = -15.0, 380.0
    18  		for Mo = 1; Mo <= 12; Mo++ {
    19  			x[Mo] = float64(FNNday(Mo, 15))
    20  			y[Mo] = Tsupw[Mo-1]
    21  			y[0], y[13] = Tsupw[11], Tsupw[0]
    22  		}
    23  
    24  		__Intgtsup_ic = 1
    25  	}
    26  
    27  	y1 = spline(n, x[:], y[:], float64(Nday), h[:], b[:], d[:], g[:], u[:], r[:])
    28  	return y1
    29  }
    30  
    31  /**************************************/
    32  /*     ３次スプライン関数による補間   */
    33  /*          n : 区間の数              */
    34  /*          x, y : 点の座標           */
    35  /*          x1 : 補間値を求める値     */
    36  /*          h, b, d, g, u, r : 作業域 */
    37  /*          return : 補間値           */
    38  /*          coded by Y.Suganuma       */
    39  /**************************************/
    40  // http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/num/spline/C++/spline.txt
    41  func spline(n int, x, y []float64, x1 float64, h, b, d, g, u, r []float64) float64 {
    42  	i, i1, k := -1, 0, 0
    43  	var y1, qi, si, xx float64
    44  
    45  	// 区間の決定
    46  	for i1 = 1; i1 < n && i < 0; i1++ {
    47  		if x1 < x[i1] {
    48  			i = i1 - 1
    49  		}
    50  	}
    51  	if i < 0 {
    52  		i = n - 1
    53  	}
    54  
    55  	// ステップ１
    56  	for i1 = 0; i1 < n; i1++ {
    57  		h[i1] = x[i1+1] - x[i1]
    58  	}
    59  	for i1 = 1; i1 < n; i1++ {
    60  		b[i1] = 2.0 * (h[i1] + h[i1-1])
    61  		d[i1] = 3.0 * ((y[i1+1]-y[i1])/h[i1] - (y[i1]-y[i1-1])/h[i1-1])
    62  	}
    63  
    64  	// ステップ２
    65  	g[1] = h[1] / b[1]
    66  	for i1 = 2; i1 < n-1; i1++ {
    67  		g[i1] = h[i1] / (b[i1] - h[i1-1]*g[i1-1])
    68  	}
    69  	u[1] = d[1] / b[1]
    70  	for i1 = 2; i1 < n; i1++ {
    71  		u[i1] = (d[i1] - h[i1-1]*u[i1-1]) / (b[i1] - h[i1-1]*g[i1-1])
    72  	}
    73  
    74  	// ステップ３
    75  	if i > 1 {
    76  		k = i
    77  	} else {
    78  		k = 1
    79  	}
    80  	r[0] = 0.0
    81  	r[n] = 0.0
    82  	r[n-1] = u[n-1]
    83  	for i1 = n - 2; i1 >= k; i1-- {
    84  		r[i1] = u[i1] - g[i1]*r[i1+1]
    85  	}
    86  
    87  	// ステップ４
    88  	xx = x1 - x[i]
    89  	qi = (y[i+1]-y[i])/h[i] - h[i]*(r[i+1]+2.0*r[i])/3.0
    90  	si = (r[i+1] - r[i]) / (3.0 * h[i])
    91  	y1 = y[i] + xx*(qi+xx*(r[i]+si*xx))
    92  
    93  	return y1
    94  }