import "dart:math";
import "package:breaker/lambert_w_function.dart";
///
/// ExponentialRate models the position and velocity over time of an object
/// that experiences an acceleration proportional to its velocity.
///
/// The system of differential equations that govern our objects movement is:
///
/// dV(t)/dt = -lambda V(t)
/// dX(t)/dt = V(t)
///
class ExponentialRate
{
final double initial_value;
final double initial_rate;
final double total_time;
final double lambda;
///
/// We know how far the object should travel, how long it should take to
/// travel that far, and how fast it should be moving at the end of its
/// travels. Given these factors, the initial speed, and the
/// proportionality coefficient (lambda), are completely defined. But
/// computing the initial speed is tricky.
///
ExponentialRate(double initial_value,
double final_value,
double final_rate,
double total_time) :
this._internal(initial_value,
_initial_rate(final_value - initial_value,
final_rate,
total_time),
final_rate,
total_time);
///
/// We cache the value of lambda instead of recomputing it each time it is
/// needed.
///
ExponentialRate._internal(this.initial_value,
double initial_rate,
double final_rate,
double total_time) :
initial_rate = initial_rate,
total_time = total_time,
lambda = _lambda(final_rate, initial_rate, total_time);
/// start function rate
///
/// Return the rate of change of the value at a given time.
///
/// Time can be negative, positive, or zero, it should also be able to be
/// positive or negative infinity, but I haven't verified that Dart's exp
/// function correctly handles those cases.
///
double rate(double time) => initial_rate * exp(-lambda * time);
/// end function rate
/// start function value
///
/// Return the value at a given time.
///
/// Time can be negative, positive, or zero, it should also be able to be
/// positive or negative infinity, but I haven't verified that Dart's exp
/// function correctly handles those cases.
///
double value(double time)
{
//
// This version of hex maintains about 14 digits of precision, not full
// machine precision (~16 digits for doubles).
//
double hex(double x)
{
//
// A threshold of 10^-2 seems like a good compromise. The smaller
// the threshold is, the fewer terms we need to use to accurately
// evaluate hex. But smaller thresholds increase the loss of
// significance when using the exponential form below. The
// truncation error introduced by a threshold of 10^-2 is O(x^6),
// or on the order of 1e-12. In fact, the first term that we lose
// from the sum is x^6/6!, or 10^-12/720, or 1.38e-15, so our
// truncation error is closer to 1e-15.
//
if (x.abs() < 1e-2)
return 1 + (1 + (1 + (1 + x / 5) * x / 4) * x / 3) * x / 2;
//
// The worst case loss of significance from the subtraction below
// will be when x = 10^-2.
//
// log2(1-(1/e^(10^-2))) ~ -6.6
//
// So we can expect to loose 7 bits of precision, from our doubles
// 53 bit mantissa. Or ~2 digits from our ~16 digits of precision.
//
return (exp(x) - 1) / x;
}
return initial_value + initial_rate * time * hex(-lambda * time);
}
/// end function value
/// start function _initial_rate
///
/// Return the initial rate that will result in the given change in value
/// and rate, over the specified total time.
///
/// delta_value in [ 0 .. infinity)
/// final_rate in (-infinity .. infinity)
/// total_time in [ 0 .. infinity)
///
static double _initial_rate(double delta_value,
double final_rate,
double total_time)
{
///
/// If delta_value is zero, then final_rate must also be zero. It is
/// not possible to have zero change in value with a non-zero
/// final_rate because the exponential is monotonic.
///
assert(delta_value != 0.0 || final_rate == 0.0);
///
/// Since we're not expected to change the value at all, and final_rate
/// is zero, _initial_rate should also be zero.
///
if (delta_value == 0.0)
return 0.0;
///
/// If total_time is zero, then delta_value must also be zero. If
/// delta_value were not zero we would have to somehow give two
/// different values for the same instant in time.
///
assert(total_time != 0.0 || delta_value == 0.0);
///
/// If total_time is zero then _initial_rate will match final_rate.
///
if (total_time == 0.0)
return final_rate;
final double A = -total_time * final_rate / delta_value;
if (A > -1.0)
return W_negative_one(A * exp(A)) * final_rate / A;
if (A < -1.0)
return W_zero(A * exp(A)) * final_rate / A;
return final_rate;
}
/// end function _initial_rate
/// start function _lambda
///
/// Compute the proportionality coefficient lambda given the start and end
/// rates as well as total time taken.
///
static double _lambda(double final_rate,
double initial_rate,
double total_time)
{
//
// If final_rate, or initial_rate is zero then the other must also be
// zero. This is because the only way to actually get to zero with an
// exponential is to already be there. And if you're there, there is
// no way to leave. In this case lambda should be zero, because our
// rate is not changing.
//
assert((final_rate == 0.0) == (initial_rate == 0.0));
if (final_rate == 0.0)
return 0.0;
//
// Lambda is not a complex number, we are only interested in solving
// for real values. So final_rate and initial_rate must also have the
// same sign.
//
assert((final_rate > 0.0) == (initial_rate > 0.0));
return -log(final_rate / initial_rate) / total_time;
}
/// end function _lambda
}