2次元スプライン関数を使用して、緯度と経度のペアに対応する目的地までの時間を計算

2024年4月7日

/*
	c:\users\ebata\tomioka3b\src\others\main41.go
	2次元スプライン関数を使用して、緯度と経度のペアに対応する目的地までの時間を計算します。
*/

package main

import (
	"fmt"
	"sort"
)

// スプライン補間用の関数
func splineInterpolation(dataSet [][]float64) func(float64, float64) float64 {
	n := len(dataSet)

	// xの値を昇順にソートする
	sort.Slice(dataSet, func(i, j int) bool {
		return dataSet[i][0] < dataSet[j][0]
	})

	// トリディアゴナル行列を作成
	h := make([]float64, n-1)
	for i := 0; i < n-1; i++ {
		h[i] = dataSet[i+1][0] - dataSet[i][0]
	}

	// 2階微分の値を計算
	delta := make([]float64, n)
	for i := 1; i < n-1; i++ {
		delta[i] = (dataSet[i+1][1]-dataSet[i][1])/h[i] - (dataSet[i][1]-dataSet[i-1][1])/h[i-1]
	}

	// トリディアゴナル行列を解く
	m := make([]float64, n)
	l := make([]float64, n)
	zArray := make([]float64, n)
	l[0] = 1
	for i := 1; i < n-1; i++ {
		l[i] = 2*(dataSet[i+1][0]-dataSet[i-1][0]) - h[i-1]*m[i-1]
		m[i] = h[i] / l[i]
		zArray[i] = (delta[i] - h[i-1]*zArray[i-1]) / l[i]
	}
	l[n-1] = 1
	zArray[n-1] = 0

	c := make([]float64, n)
	b := make([]float64, n)
	d := make([]float64, n)
	for j := n - 2; j >= 0; j-- {
		c[j] = zArray[j] - m[j]*c[j+1]
		b[j] = (dataSet[j+1][1]-dataSet[j][1])/h[j] - h[j]*(c[j+1]+2*c[j])/3
		d[j] = (c[j+1] - c[j]) / (3 * h[j])
	}

	// 補間関数を返す
	return func(xVal, yVal float64) float64 {
		// xの範囲を確認
		if xVal < dataSet[0][0] || xVal > dataSet[n-1][0] {
			panic("x value is out of range")
		}

		// 対応するiを探す
		i := 0
		for i < n-1 && dataSet[i+1][0] <= xVal {
			i++
		}

		// スプライン補間を計算
		dx := xVal - dataSet[i][0]
		return dataSet[i][2] + b[i]*dx + c[i]*dx*dx + d[i]*dx*dx*dx
	}
}

func main() {
	// データ点の定義
	dataSet := [][]float64{
		{35.7281578, 139.7680665, 73},
		{35.7214573, 139.7754384, 63},
		{35.7141672, 139.7748342, 57},
		{35.7075171, 139.7564025, 67},
		{35.698383, 139.7704968, 65},
		{35.6997175, 139.7618655, 61},
		{35.7020484, 139.7509267, 66},
		{35.6918216, 139.7683569, 63},
		{35.6812362, 139.7645499, 57},
		{35.6750133, 139.7604455, 59},
		{35.666379, 139.7557649, 55},
		{35.6553809, 139.754554, 54},
		{35.6457361, 139.7449875, 54},
		{35.6284713, 139.7361848, 43},
		{35.6208763, 139.7405387, 50},
		{35.6147448, 139.7412068, 49},
		{35.6086671, 139.7426731, 45},
		{35.6052508, 139.7421627, 50},
		{35.5983959, 139.7391046, 50},
		{35.58754, 139.7355111, 51},
		{35.5788485, 139.7348961, 40},
		{35.5721351, 139.7318494, 53},
		{35.5666392, 139.7281759, 53},
		{35.5613144, 139.7215761, 33},
		{35.5446458, 139.7493338, 46},
		{35.5329877, 139.6982966, 30},
		{35.5348941, 139.7449763, 42},
		{35.5229857, 139.6889874, 41},
		{35.5179115, 139.6811867, 40},
		{35.5071467, 139.677526, 34},
		{35.4998507, 139.6713032, 31},
		{35.4923606, 139.6622968, 30},
		{35.4852628, 139.6544734, 29},
		{35.4846744, 139.6452731, 32},
		{35.4808759, 139.6394986, 26},
		{35.4763238, 139.631979, 30},
		{35.4688634, 139.6268306, 32},
		{35.4659811, 139.6194871, 20},
		{35.4571602, 139.6199226, 27},
		{35.4464751, 139.6235656, 25},
		{35.436673, 139.6217708, 24},
		{35.4334892, 139.6206713, 18},
		{35.4314711, 139.6030542, 16},
		{35.424238, 139.5927802, 17},
		{35.4201765, 139.586678, 14},
		{35.413768, 139.582819, 10},
		{35.3819518, 139.6165374, 3},
		{35.4091204, 139.5781752, 7},
		{35.3966138, 139.6051573, 4},
		{35.3645904, 139.6267988, 1},
		{35.3428573, 139.6190711, 3},
		{35.3314185, 139.6176347, 6},
		{35.330416, 139.6295796, 17},
		{35.3337822, 139.630869, 18},
		{35.3407534, 139.6332478, 18},
		{35.3400425, 139.6390968, 20},
		{35.323276, 139.6151587, 13},
		{35.2973885, 139.5758056, 31},
		{35.3158184, 139.6222558, 15},
		{35.312573, 139.6250891, 16},
		{35.308405, 139.628248, 17},
		{35.2825803, 139.6405151, 19},
		{35.2916167, 139.6283632, 17},
		{35.277848, 139.6379121, 21},
		{35.2802815, 139.6599103, 19},
		{35.2786985, 139.6674653, 20},
		{35.2616259, 139.6679659, 26},
		{35.2635438, 139.6841348, 23},
		{35.2482969, 139.6775595, 25},
		{35.2509382, 139.712402, 30},
		{35.2315701, 139.699661, 29},
		{35.2120725, 139.6824547, 33},
		{35.2055645, 139.671602, 44},
		{35.1986888, 139.663162, 37},
		{35.1880928, 139.6507295, 39},
		{35.1775554, 139.6306392, 42},
	}

	// スプライン補間関数を作成
	interpolatedFunction := splineInterpolation(dataSet)

	// 特定の座標で補間された値を計算
	xValue := 35.4688634
	yValue := 139.6268306

	xValue, yValue = 35.7214573, 139.7754384

	interpolatedValue := interpolatedFunction(xValue, yValue)

	// 結果の出力
	fmt.Printf("補間された値: %.2f\n", interpolatedValue)
}