forked from gsn/predictor
126 lines
3.2 KiB
Go
126 lines
3.2 KiB
Go
package grib
|
|
|
|
import "math"
|
|
|
|
func lerp(a, b, t float64) float64 { return a + t*(b-a) }
|
|
|
|
// ghInterp returns interpolated geopotential height at given time/pressure/lat/lon
|
|
func (d *dataset) ghInterp(ti, pi int, y0, y1, x0, x1 int, wy, wx float64) float64 {
|
|
g00 := d.cube.val(0, ti, pi, y0, x0)
|
|
g10 := d.cube.val(0, ti, pi, y0, x1)
|
|
g01 := d.cube.val(0, ti, pi, y1, x0)
|
|
g11 := d.cube.val(0, ti, pi, y1, x1)
|
|
return (1-wy)*((1-wx)*float64(g00)+wx*float64(g10)) + wy*((1-wx)*float64(g01)+wx*float64(g11))
|
|
}
|
|
|
|
// searchAltLevel uses geopotential height to find pressure level bracket for target altitude.
|
|
func (d *dataset) searchAltLevel(alt float64, ti, y0, y1, x0, x1 int, wy, wx float64) (int, float64) {
|
|
levels := d.ds.Levels
|
|
nLevels := len(levels)
|
|
|
|
lo, hi := 0, nLevels-1
|
|
for lo < hi-1 {
|
|
mid := (lo + hi) / 2
|
|
ghMid := d.ghInterp(ti, mid, y0, y1, x0, x1, wy, wx)
|
|
if ghMid < alt {
|
|
lo = mid
|
|
} else {
|
|
hi = mid
|
|
}
|
|
}
|
|
|
|
ghLo := d.ghInterp(ti, lo, y0, y1, x0, x1, wy, wx)
|
|
ghHi := d.ghInterp(ti, hi, y0, y1, x0, x1, wy, wx)
|
|
|
|
wp := 0.0
|
|
if ghHi != ghLo {
|
|
wp = (alt - ghLo) / (ghHi - ghLo)
|
|
}
|
|
if wp < 0 {
|
|
wp = 0
|
|
}
|
|
if wp > 1 {
|
|
wp = 1
|
|
}
|
|
|
|
return lo, wp
|
|
}
|
|
|
|
// uv выполняет интерполяцию ветра по 4 измерениям (time, pressure, lat, lon).
|
|
func (d *dataset) uv(lat, lon, alt float64, tHours float64) (float64, float64) {
|
|
if lon < 0 {
|
|
lon += 360
|
|
}
|
|
|
|
inv := d.ds.InvResolution()
|
|
|
|
// GRIB scan north→south: index 0 = 90°N
|
|
iy := (90 - lat) * inv
|
|
y0 := int(math.Floor(iy))
|
|
if y0 < 0 {
|
|
y0 = 0
|
|
}
|
|
if y0 >= d.cube.lat-1 {
|
|
y0 = d.cube.lat - 2
|
|
}
|
|
y1 := y0 + 1
|
|
wy := iy - float64(y0)
|
|
|
|
ix := lon * inv
|
|
x0 := int(math.Floor(ix)) % d.cube.lon
|
|
x1 := (x0 + 1) % d.cube.lon
|
|
wx := ix - float64(x0)
|
|
|
|
// Время: tHours делим на шаг, чтобы получить индекс в кубе
|
|
tIdx := tHours / float64(d.ds.TimeStep)
|
|
it0 := int(math.Floor(tIdx))
|
|
if it0 < 0 {
|
|
it0 = 0
|
|
}
|
|
if it0 >= d.cube.t-1 {
|
|
it0 = d.cube.t - 2
|
|
}
|
|
wt := tIdx - float64(it0)
|
|
|
|
// ISA: высота → давление → индекс уровня
|
|
levels := d.ds.Levels
|
|
p := pressureFromAlt(alt)
|
|
ip0 := 0
|
|
for ip0+1 < len(levels) && levels[ip0+1] > p {
|
|
ip0++
|
|
}
|
|
ip1 := ip0 + 1
|
|
if ip1 >= len(levels) {
|
|
ip1 = len(levels) - 1
|
|
}
|
|
wp := 0.0
|
|
if levels[ip0] != levels[ip1] {
|
|
wp = (levels[ip0] - p) / (levels[ip0] - levels[ip1])
|
|
}
|
|
|
|
fetch := func(ti, pi int) (float64, float64) {
|
|
u00 := d.cube.val(1, ti, pi, y0, x0)
|
|
u10 := d.cube.val(1, ti, pi, y0, x1)
|
|
u01 := d.cube.val(1, ti, pi, y1, x0)
|
|
u11 := d.cube.val(1, ti, pi, y1, x1)
|
|
v00 := d.cube.val(2, ti, pi, y0, x0)
|
|
v10 := d.cube.val(2, ti, pi, y0, x1)
|
|
v01 := d.cube.val(2, ti, pi, y1, x0)
|
|
v11 := d.cube.val(2, ti, pi, y1, x1)
|
|
uxy := (1-wy)*((1-wx)*float64(u00)+wx*float64(u10)) + wy*((1-wx)*float64(u01)+wx*float64(u11))
|
|
vxy := (1-wy)*((1-wx)*float64(v00)+wx*float64(v10)) + wy*((1-wx)*float64(v01)+wx*float64(v11))
|
|
return uxy, vxy
|
|
}
|
|
|
|
u0p0, v0p0 := fetch(it0, ip0)
|
|
u0p1, v0p1 := fetch(it0, ip1)
|
|
u1p0, v1p0 := fetch(it0+1, ip0)
|
|
u1p1, v1p1 := fetch(it0+1, ip1)
|
|
uLow := lerp(u0p0, u0p1, wp)
|
|
vLow := lerp(v0p0, v0p1, wp)
|
|
uHig := lerp(u1p0, u1p1, wp)
|
|
vHig := lerp(v1p0, v1p1, wp)
|
|
u := lerp(uLow, uHig, wt)
|
|
v := lerp(vLow, vHig, wt)
|
|
return u, v
|
|
}
|