/** This program tries to help you find a Taylor series for a  particular
    function around x0.

    See:
    https://mathworld.wolfram.com/TaylorSeries.html
*/

use functionUtils.frink
use allTransforms.frink

symbolicMode[true]
showApproximations[false]

/** This function finds a Taylor series for a particular function around x0.
    x0 can be a symbol, in which case you probably want to wrap it in a
    noEval[] block.  */ 
TaylorSeries[f, x0=0, terms=6, debug=false] :=
{
   body = functionBody[f]
   vars = functionArgumentsAsSymbols[f]
   var1 = vars@0

   if debug
      println["var1 is $var1"]
   
   sum = substituteExpression[body, var1, x0]
   
   if debug
      println["sum = " + inputForm[sum]]
   
   d0 = body
   for i = 1 to terms-1
   {
      if debug
         println["Term $i:"]
      d0 = transformExpression[constructExpression["FunctionCall", ["D", d0, var1]]]
      if debug 
         println[" d0 = " + inputForm[d0]]
    
      d1 = transformExpression[substituteExpression[d0, var1, x0]]
      if debug 
         println[" d1 = " + inputForm[d1]]
      
      df = transformExpression[d1 (var1-x0)^i / i!]
      if debug
         println[" df = " + inputForm[df]]
      sum = sum + df
      if debug
         println[" sum = $sum"]
      sum = transformExpression[sum]    // This may simplify the sum.
      if debug
         println[" sum = $sum\n"]
   }

   return transformExpression[sum]
}


// This is the equation for total relativistic energy.  If you expand it with
// a Taylor series, you will get the Einsteinian mass-energy relation
//
//   E = m0 c^2
//
// where m0 is the rest mass and c is the speed of light
//
// plus the Newtonian kinetic energy
//
//   KE = 1/2 m v^2
//
// plus it also contains relativistic expansion terms which you'll see if you
// run this program through a Taylor series.  Relativity gives an infinite
// expansion of corrective terms!
//
// If you run this, you'll see that the Taylor series of this is:
// c^2 m0 + 1/2 m0 v^2 + 3/8 c^-2 m0 v^4 + 5/16 c^-4 m0 v^6 + 35/128 c^-6 m0 v^8
//
// All of this is derived through the Lorentz equation (1/sqrt[1 - v^2/c^2])

f = {|v| (m0 c^2) (1/sqrt[1 - v^2/c^2]) }


// This is the equation for relativistic kinetic energy.  It is the
// relativistic total energy (above) minus the rest-mass energy "m c^2".
// The Taylor series expansion is actually important.  The original equation
// is valid for all velocities but has numerical issues if evaluated naively
// for low velocities in the given form below.  That's why, for low velocities,
// we expand it to a Taylor series which is closer to the classical 1/2 m v^2
// (which you can actually see in the Taylor series expansion) but with
// interesting correction terms.
//
// The Taylor series expansion is:
//
// 1/2 m v^2 + 3/8 c^-2 m v^4 + 5/16 c^-4 m v^6 + 35/128 c^-6 m v^8 ...
//
// which converges quickly for low values of v, (so we choose the Taylor series
// for low values of v) but converges slowly for high values of v (so we choose
// the sqrt[...] formulation for high values of v, which has large errors for
// small values of v) which are accurate for high
// values of v which don't lose numerical precision from taking a huge
// total energy minus a huge restmass-energy.
//
// An empirical threshold is found around v = 0.13 c, depending on how many
// terms of the Taylor series are taken.

f = {|v| (m c^2) (1/sqrt[1 - v^2/c^2] - 1) }

// This is the equation for the Lorentz transformation.
f = {|v| 1 / sqrt[1 - v^2/c^2] }

// This has an elegant Taylor series expansion.
//f = {|x| 1/(1-x) }

// This is the classic example.  It's beautiful and quick to evaluate
// f = {|x| e^x }

// Arctangent has an elegant and simple Taylor series expansion and can be used
// to build lots of other inverse trigonometric functions.  It can also be used
// to calculate pi since pi = 4 arctan[1] but that would be incredibly slow.
// See https://www.johndcook.com/blog/2021/01/05/bootstrapping-math-library/
// f = {|x| arctan[x] }

println["The Taylor series for " + inputForm[functionBody[f]] + " is:"]
println[inputForm[TaylorSeries[f, 0, 15, true]]]