|
| 1 | +package main |
| 2 | + |
| 3 | +import ( |
| 4 | + "fmt" |
| 5 | + |
| 6 | + "gonum.org/v1/gonum/mat" |
| 7 | +) |
| 8 | + |
| 9 | +// Given: |
| 10 | +// (1) p¹ = p₀¹ + v¹⋅t |
| 11 | +// (2) p² = p₀² + v²⋅t |
| 12 | +// (3) p³ = p₀³ + v³⋅t |
| 13 | +// ... |
| 14 | +// |
| 15 | +// Find: |
| 16 | +// p = p₀ + v⋅t |
| 17 | +// |
| 18 | +// such as: |
| 19 | +// p¹=p when t¹=t |
| 20 | +// p²=p when t²=t |
| 21 | +// p³=p when t³=t |
| 22 | +// ... |
| 23 | +// |
| 24 | +// where |
| 25 | +// pⁿ = [pₓⁿ, pᵧⁿ, pₛⁿ] |
| 26 | +// p₀ⁿ = [pₓ₀ⁿ, pᵧⁿ, pₛ₀ⁿ] |
| 27 | +// vⁿ = [vₓⁿ, vᵧⁿ, vₛⁿ] |
| 28 | + |
| 29 | +// p¹ + v¹⋅t = p + v⋅t |
| 30 | +// p² + v²⋅t = p + v⋅t |
| 31 | +// p³ + v³⋅t = p + v⋅t |
| 32 | +// pⁿ + vⁿ⋅t = p + v⋅t |
| 33 | +// |
| 34 | +// (p - p¹) = (v¹ - v)t¹ |
| 35 | +// (p - p²) = (v² - v)t² |
| 36 | +// (p - p³) = (v³ - v)t³ |
| 37 | +// (p - pⁿ) = (vⁿ - v)tⁿ |
| 38 | +// and so: |
| 39 | +// (p - pⁿ) ⋅ (v-vⁿ) = 0 |
| 40 | +func main() { |
| 41 | + hail := parse() |
| 42 | + |
| 43 | + // Ax = b |
| 44 | + // |
| 45 | + // A = [ 0, -v⁰ₛ +v¹ₛ, v⁰ᵧ - v¹ᵧ, 0, -p¹ₛ +p⁰ₛ, p¹ᵧ -p⁰ᵧ |
| 46 | + // v⁰ₛ -v¹ₛ, 0, -v⁰ₓ + v¹ₓ, p¹ₛ -p⁰ₛ, 0, -p¹ₓ +p⁰ₓ |
| 47 | + // -v⁰ᵧ +v¹ᵧ, v⁰ₓ -v¹ₓ, 0, -p¹ᵧ +p⁰ᵧ, p¹ₓ -p⁰ₓ, 0 |
| 48 | + // 0, -v⁰ₛ +v²ₛ, v⁰ᵧ - v²ᵧ, 0, -p²ₛ +p⁰ₛ, p²ᵧ -p⁰ᵧ |
| 49 | + // v⁰ₛ -v²ₛ, 0, -v⁰ₓ + v²ₓ, p²ₛ -p⁰ₛ, 0, -p²ₓ +p⁰ₓ |
| 50 | + // -v⁰ᵧ +v²ᵧ, v⁰ₓ -v²ₓ, 0, -p²ᵧ +p⁰ᵧ, p²ₓ -p⁰ₓ, 0 ] |
| 51 | + // |
| 52 | + // x = [ pₓ |
| 53 | + // pᵧ |
| 54 | + // pₛ |
| 55 | + // vₓ |
| 56 | + // vᵧ |
| 57 | + // vₛ ] |
| 58 | + // |
| 59 | + // b = [ p¹ᵧ⋅v¹ₛ -p¹ₛ⋅v¹ᵧ -p⁰ᵧ⋅v⁰ₛ +p⁰ₛ⋅v⁰ᵧ, |
| 60 | + // p¹ₛ⋅v¹ₓ -p¹ₓ⋅v¹ₛ -p⁰ₛ⋅v⁰ₓ +p⁰ₓ⋅v⁰ₛ, |
| 61 | + // p¹ₓ⋅v¹ᵧ -p¹ᵧ⋅v¹ₓ -p⁰ₓ⋅v⁰ᵧ +p⁰ᵧ⋅v⁰ₓ, |
| 62 | + // p²ᵧ⋅v²ₛ -p²ₛ⋅v²ᵧ -p⁰ᵧ⋅v⁰ₛ +p⁰ₛ⋅v⁰ᵧ, |
| 63 | + // p²ₛ⋅v²ₓ -p²ₓ⋅v²ₛ -p⁰ₛ⋅v⁰ₓ +p⁰ₓ⋅v⁰ₛ, |
| 64 | + // p²ₓ⋅v²ᵧ -p²ᵧ⋅v²ₓ -p⁰ₓ⋅v⁰ᵧ +p⁰ᵧ⋅v⁰ₓ ] |
| 65 | + // |
| 66 | + |
| 67 | + // hail[0] with hail[1] |
| 68 | + // (p0 - p[1]) x (v0 - v[1]) == 0 |
| 69 | + A00 := diff(crossMatrix(hail[0].Vel), crossMatrix(hail[1].Vel)) |
| 70 | + A03 := diff(crossMatrix(hail[1].Pos), crossMatrix(hail[0].Pos)) |
| 71 | + |
| 72 | + // hail[0] with hail[2] |
| 73 | + // (p0 - p[2]) x (v0 - v[2]) == 0 |
| 74 | + A30 := diff(crossMatrix(hail[0].Vel), crossMatrix(hail[2].Vel)) |
| 75 | + A33 := diff(crossMatrix(hail[2].Pos), crossMatrix(hail[0].Pos)) |
| 76 | + |
| 77 | + A := mat.NewDense(6, 6, []float64{ |
| 78 | + A00[0], A00[1], A00[2], A03[0], A03[1], A03[2], |
| 79 | + A00[3], A00[4], A00[5], A03[3], A03[4], A03[5], |
| 80 | + A00[6], A00[7], A00[8], A03[6], A03[7], A03[8], |
| 81 | + A30[0], A30[1], A30[2], A33[0], A33[1], A33[2], |
| 82 | + A30[3], A30[4], A30[5], A33[3], A33[4], A33[5], |
| 83 | + A30[6], A30[7], A30[8], A33[6], A33[7], A33[8], |
| 84 | + }) |
| 85 | + |
| 86 | + b0 := diff(hail[1].Pos.cross(hail[1].Vel).toF(), hail[0].Pos.cross(hail[0].Vel).toF()) |
| 87 | + b3 := diff(hail[2].Pos.cross(hail[2].Vel).toF(), hail[0].Pos.cross(hail[0].Vel).toF()) |
| 88 | + |
| 89 | + b := mat.NewVecDense(6, []float64{b0[0], b0[1], b0[2], b3[0], b3[1], b3[2]}) |
| 90 | + |
| 91 | + var x mat.VecDense |
| 92 | + _ = x.SolveVec(A, b) |
| 93 | + |
| 94 | + rock := HailStone{ |
| 95 | + Pos: P{x.At(0, 0), x.At(1, 0), x.At(2, 0)}, |
| 96 | + Vel: P{x.At(3, 0), x.At(4, 0), x.At(5, 0)}, |
| 97 | + } |
| 98 | + |
| 99 | + fmt.Printf("%.0f\n", rock.Pos.x+rock.Pos.y+rock.Pos.z) |
| 100 | +} |
| 101 | + |
| 102 | +func crossMatrix(p P) []float64 { |
| 103 | + return []float64{ |
| 104 | + 0, -p.z, p.y, |
| 105 | + p.z, 0, -p.x, |
| 106 | + -p.y, p.x, 0, |
| 107 | + } |
| 108 | +} |
| 109 | + |
| 110 | +func diff(a, b []float64) []float64 { |
| 111 | + res := make([]float64, len(a)) |
| 112 | + for i := 0; i < len(a); i++ { |
| 113 | + res[i] = a[i] - b[i] |
| 114 | + } |
| 115 | + return res |
| 116 | +} |
| 117 | + |
| 118 | +// https://wikimedia.org/api/rest_v1/media/math/render/svg/3242bd71d63c393d02302c7dbe517cd0ec352d31 |
| 119 | +// https://en.wikipedia.org/wiki/Cross_product#Coordinate_notation |
| 120 | +func (p P) cross(p2 P) P { |
| 121 | + return P{ |
| 122 | + p.y*p2.z - p.z*p2.y, |
| 123 | + p.z*p2.x - p.x*p2.z, |
| 124 | + p.x*p2.y - p.y*p2.x, |
| 125 | + } |
| 126 | +} |
| 127 | + |
| 128 | +func (p P) toF() []float64 { |
| 129 | + return []float64{p.x, p.y, p.z} |
| 130 | +} |
0 commit comments