forked from gsn/predictor
53 lines
1.6 KiB
Go
53 lines
1.6 KiB
Go
package grib
|
||
|
||
import "math"
|
||
|
||
func lerp(a, b, t float64) float64 { return a + t*(b-a) }
|
||
|
||
// Interpolate 16‑point (time, p, lat, lon)
|
||
func (d *dataset) uv(lat, lon, alt float64, tHours float64) (float64, float64) {
|
||
if lon < 0 {
|
||
lon += 360
|
||
}
|
||
iy := (lat + 90) * 2
|
||
y0 := int(math.Floor(iy))
|
||
y1 := y0 + 1
|
||
wy := iy - float64(y0)
|
||
ix := lon * 2
|
||
x0 := int(math.Floor(ix)) % d.cube.lon
|
||
x1 := (x0 + 1) % d.cube.lon
|
||
wx := ix - float64(x0)
|
||
it0 := int(math.Floor(tHours / 3.0))
|
||
wt := (tHours - float64(it0*3)) / 3.0
|
||
p := pressureFromAlt(alt)
|
||
ip0 := 0
|
||
for ip0+1 < len(pressureLevels) && pressureLevels[ip0+1] > p {
|
||
ip0++
|
||
}
|
||
ip1 := ip0 + 1
|
||
wp := (pressureLevels[ip0] - p) / (pressureLevels[ip0] - pressureLevels[ip1])
|
||
fetch := func(ti, pi int) (float64, float64) {
|
||
u00 := d.cube.val(0, ti, pi, y0, x0)
|
||
u10 := d.cube.val(0, ti, pi, y0, x1)
|
||
u01 := d.cube.val(0, ti, pi, y1, x0)
|
||
u11 := d.cube.val(0, ti, pi, y1, x1)
|
||
v00 := d.cube.val(1, ti, pi, y0, x0)
|
||
v10 := d.cube.val(1, ti, pi, y0, x1)
|
||
v01 := d.cube.val(1, ti, pi, y1, x0)
|
||
v11 := d.cube.val(1, 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
|
||
}
|