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
}