# frange.py
# https://stackoverflow.com/questions/7267226/range-for-floats

from fractions import Fraction

def frange(start, stop, jump, end=False, via_str=False):
  """
  Equivalent of Python 3 range for decimal numbers.

  Notice that, because of arithmetic errors, it is safest to
  pass the arguments as strings, so they can be interpreted to exact fractions.

  >>> assert Fraction('1.1') - Fraction(11, 10) == 0.0
  >>> assert Fraction( 0.1 ) - Fraction(1, 10) == Fraction(1, 180143985094819840)

  Parameter `via_str` can be set to True to transform inputs in strings and then to fractions.
  When inputs are all non-periodic (in base 10), even if decimal, this method is safe as long
  as approximation happens beyond the decimal digits that Python uses for printing.


  For example, in the case of 0.1, this is the case:

  >>> assert str(0.1) == '0.1'
  >>> assert '%.50f' % 0.1 == '0.10000000000000000555111512312578270211815834045410'


  If you are not sure whether your decimal inputs all have this property, you are better off
  passing them as strings. String representations can be in integer, decimal, exponential or
  even fraction notation.

  >>> assert list(frange(1, 100.0, '0.1', end=True))[-1] == 100.0
  >>> assert list(frange(1.0, '100', '1/10', end=True))[-1] == 100.0
  >>> assert list(frange('1', '100.0', '.1', end=True))[-1] == 100.0
  >>> assert list(frange('1.0', 100, '1e-1', end=True))[-1] == 100.0
  >>> assert list(frange(1, 100.0, 0.1, end=True))[-1] != 100.0
  >>> assert list(frange(1, 100.0, 0.1, end=True, via_str=True))[-1] == 100.0

  """

  if via_str:
    start = str(start)
    stop = str(stop)
    jump = str(jump)

  start = Fraction(start)
  stop = Fraction(stop)
  jump = Fraction(jump)

  while start < stop:
    yield float(start)
    start += jump

  if end and start == stop:
    yield(float(start))



def linspace(start, stop, num=50, endpoint=True):
  # so I don't need numpy linspace
  vals = []
  if endpoint:
    step = (stop - start)/(num-1)
  else:
    step = (stop - start)/num
  for i in range(num):
    vals.append(start + (step * i))
  return vals



def fwdEuler(fD, xs, h, n, y0):
  ys = [0]*n
  y = y0
  for i,x in enumerate(xs):
    ys[i] = y
    y = eulerStep(x, y, fD, h)
  return ys


def eulerStep(xi, yi, fD, h):
  return yi + h*fD(xi, yi)


def rk4(fD, xs, h, n, y0):
  ys = [0]*n
  y = y0
  for i,x in enumerate(xs):
    ys[i] = y
    y = rk4Step(x, y, fD, h)
  return ys


def rk4Step(xi, yi, fD, h):
  k0 = h*fD(xi, yi)
  k1 = h*fD(xi+h/2, yi+k0/2)
  k2 = h*fD(xi+h/2, yi+k1/2)
  k3 = h*fD(xi+h, yi+k2)
  return yi + (k0 + 2*k1 + 2*k2 + k3)/6