main
1package harmonica
2
3// This file defines a simplified damped harmonic oscillator, colloquially
4// known as a spring. This is ported from Ryan Juckett’s simple damped harmonic
5// motion, originally written in C++.
6//
7// Example usage:
8//
9// // Run once to initialize.
10// spring := NewSpring(FPS(60), 6.0, 0.2)
11//
12// // Update on every frame.
13// pos := 0.0
14// velocity := 0.0
15// targetPos := 100.0
16// someUpdateLoop(func() {
17// pos, velocity = spring.Update(pos, velocity, targetPos)
18// })
19//
20// For background on the algorithm see:
21// https://www.ryanjuckett.com/damped-springs/
22
23/******************************************************************************
24
25 Copyright (c) 2008-2012 Ryan Juckett
26 http://www.ryanjuckett.com/
27
28 This software is provided 'as-is', without any express or implied
29 warranty. In no event will the authors be held liable for any damages
30 arising from the use of this software.
31
32 Permission is granted to anyone to use this software for any purpose,
33 including commercial applications, and to alter it and redistribute it
34 freely, subject to the following restrictions:
35
36 1. The origin of this software must not be misrepresented; you must not
37 claim that you wrote the original software. If you use this software
38 in a product, an acknowledgment in the product documentation would be
39 appreciated but is not required.
40
41 2. Altered source versions must be plainly marked as such, and must not be
42 misrepresented as being the original software.
43
44 3. This notice may not be removed or altered from any source
45 distribution.
46
47*******************************************************************************
48
49 Ported to Go by Charmbracelet, Inc. in 2021.
50
51******************************************************************************/
52
53import (
54 "math"
55 "time"
56)
57
58// FPS returns a time delta for a given number of frames per second. This
59// value can be used as the time delta when initializing a Spring. Note that
60// game engines often provide the time delta as well, which you should use
61// instead of this function, if possible.
62//
63// Example:
64//
65// spring := NewSpring(FPS(60), 5.0, 0.2)
66//
67func FPS(n int) float64 {
68 return (time.Second / time.Duration(n)).Seconds()
69}
70
71// In calculus ε is, in vague terms, an arbitrarily small positive number. In
72// the original C++ source ε is represented as such:
73//
74// const float epsilon = 0.0001
75//
76// Some Go programmers use:
77//
78// const epsilon float64 = 0.00000001
79//
80// We can, however, calculate the machine’s epsilon value, with the drawback
81// that it must be a variable versus a constant.
82var epsilon = math.Nextafter(1, 2) - 1
83
84// Spring contains a cached set of motion parameters that can be used to
85// efficiently update multiple springs using the same time step, angular
86// frequency and damping ratio.
87//
88// To use a Spring call New with the time delta (that's animation frame
89// length), frequency, and damping parameters, cache the result, then call
90// Update to update position and velocity values for each spring that neeeds
91// updating.
92//
93// Example:
94//
95// // First precompute spring coefficients based on your settings:
96// var x, xVel, y, yVel float64
97// deltaTime := FPS(60)
98// s := NewSpring(deltaTime, 5.0, 0.2)
99//
100// // Then, in your update loop:
101// x, xVel = s.Update(x, xVel, 10) // update the X position
102// y, yVel = s.Update(y, yVel, 20) // update the Y position
103//
104type Spring struct {
105 posPosCoef, posVelCoef float64
106 velPosCoef, velVelCoef float64
107}
108
109// NewSpring initializes a new Spring, computing the parameters needed to
110// simulate a damped spring over a given period of time.
111//
112// The delta time is the time step to advance; essentially the framerate.
113//
114// The angular frequency is the angular frequency of motion, which affects the
115// speed.
116//
117// The damping ratio is the damping ratio of motion, which determines the
118// oscillation, or lack thereof. There are three categories of damping ratios:
119//
120// Damping ratio > 1: over-damped.
121// Damping ratio = 1: critlcally-damped.
122// Damping ratio < 1: under-damped.
123//
124// An over-damped spring will never oscillate, but reaches equilibrium at
125// a slower rate than a critically damped spring.
126//
127// A critically damped spring will reach equilibrium as fast as possible
128// without oscillating.
129//
130// An under-damped spring will reach equilibrium the fastest, but also
131// overshoots it and continues to oscillate as its amplitude decays over time.
132func NewSpring(deltaTime, angularFrequency, dampingRatio float64) (s Spring) {
133 // Keep values in a legal range.
134 angularFrequency = math.Max(0.0, angularFrequency)
135 dampingRatio = math.Max(0.0, dampingRatio)
136
137 // If there is no angular frequency, the spring will not move and we can
138 // return identity.
139 if angularFrequency < epsilon {
140 s.posPosCoef = 1.0
141 s.posVelCoef = 0.0
142 s.velPosCoef = 0.0
143 s.velVelCoef = 1.0
144 return s
145 }
146
147 if dampingRatio > 1.0+epsilon {
148 // Over-damped.
149 var (
150 za = -angularFrequency * dampingRatio
151 zb = angularFrequency * math.Sqrt(dampingRatio*dampingRatio-1.0)
152 z1 = za - zb
153 z2 = za + zb
154
155 e1 = math.Exp(z1 * deltaTime)
156 e2 = math.Exp(z2 * deltaTime)
157
158 invTwoZb = 1.0 / (2.0 * zb) // = 1 / (z2 - z1)
159
160 e1_Over_TwoZb = e1 * invTwoZb
161 e2_Over_TwoZb = e2 * invTwoZb
162
163 z1e1_Over_TwoZb = z1 * e1_Over_TwoZb
164 z2e2_Over_TwoZb = z2 * e2_Over_TwoZb
165 )
166
167 s.posPosCoef = e1_Over_TwoZb*z2 - z2e2_Over_TwoZb + e2
168 s.posVelCoef = -e1_Over_TwoZb + e2_Over_TwoZb
169
170 s.velPosCoef = (z1e1_Over_TwoZb - z2e2_Over_TwoZb + e2) * z2
171 s.velVelCoef = -z1e1_Over_TwoZb + z2e2_Over_TwoZb
172
173 } else if dampingRatio < 1.0-epsilon {
174 // Under-damped.
175 var (
176 omegaZeta = angularFrequency * dampingRatio
177 alpha = angularFrequency * math.Sqrt(1.0-dampingRatio*dampingRatio)
178
179 expTerm = math.Exp(-omegaZeta * deltaTime)
180 cosTerm = math.Cos(alpha * deltaTime)
181 sinTerm = math.Sin(alpha * deltaTime)
182
183 invAlpha = 1.0 / alpha
184
185 expSin = expTerm * sinTerm
186 expCos = expTerm * cosTerm
187 expOmegaZetaSin_Over_Alpha = expTerm * omegaZeta * sinTerm * invAlpha
188 )
189
190 s.posPosCoef = expCos + expOmegaZetaSin_Over_Alpha
191 s.posVelCoef = expSin * invAlpha
192
193 s.velPosCoef = -expSin*alpha - omegaZeta*expOmegaZetaSin_Over_Alpha
194 s.velVelCoef = expCos - expOmegaZetaSin_Over_Alpha
195
196 } else {
197 // Critically damped.
198 var (
199 expTerm = math.Exp(-angularFrequency * deltaTime)
200 timeExp = deltaTime * expTerm
201 timeExpFreq = timeExp * angularFrequency
202 )
203
204 s.posPosCoef = timeExpFreq + expTerm
205 s.posVelCoef = timeExp
206
207 s.velPosCoef = -angularFrequency * timeExpFreq
208 s.velVelCoef = -timeExpFreq + expTerm
209 }
210
211 return s
212}
213
214// Update updates position and velocity values against a given target value.
215// Call this after calling NewSpring to update values.
216func (s Spring) Update(pos, vel float64, equilibriumPos float64) (newPos, newVel float64) {
217 oldPos := pos - equilibriumPos // update in equilibrium relative space
218 oldVel := vel
219
220 newPos = oldPos*s.posPosCoef + oldVel*s.posVelCoef + equilibriumPos
221 newVel = oldPos*s.velPosCoef + oldVel*s.velVelCoef
222
223 return newPos, newVel
224}