//Integração Numérica
/*  
    Adaptive Simpson's Quadrature
	(Desenvolvido por Lucio Tavernini -www.tavernini.com  -  Adaptado por Alfredo J G A Borba - WebCalc)
*/

var integralError = false
var integralBound = 0
var ru = roundingUnit()

function Integral(f, a, b, errorBound) {
  var left, right, max, mstart, j
  var jend, step, m, m1, fa, fb, v1, v2
  var error, bound, h, h6, value, result
  var alerta = false;
  
  if (a == b) return 0
  integralBound = errorBound
  integralError = false 
  if (integralError) return Math.NaN

  //  Limites trocados?
  if (a > b)  return -Integral(f, b, a, errorBound)

  //  Integrar de -infinito a +infinito? 
  if (a == Number.NEGATIVE_INFINITY && 
    b == Number.POSITIVE_INFINITY) 
    return Integral(f, 0, b, 0.5*errorBound) 
             + Integral(f, a, 0, 0.5*errorBound)

  //  Integrar de a até infinito? 
  if (b == Number.POSITIVE_INFINITY) {
    h = 5
    left = a
    right = a + h
    result = 0
    do {
      value = Integral(f, left, right, errorBound/h)
      result += value
      h = 2*h
      left = right
      right = left + h
    }
    while (!isNaN(value) && 
      Math.abs(value) >= 0.5*errorBound/h && left < right)
    if (integralError) result = Number.NaN
    return result
  }

  //  Integrar de -infinito até b? 
  if (a == Number.NEGATIVE_INFINITY) {
    h = 5
    left = b - h
    right = b
    result = 0
    do {
      value = Integral(f, left, right, errorBound/h)
      result += value
      h = 2*h
      right = left
      left = right - h
    }
    while (!isNaN(value) && 
      Math.abs(value) >= 0.5*errorBound/h && left < right)
    if (integralError) result = Number.NaN
    return result
  }

  //  Integrar de a até b?  Inicializar: 
  if (integralError) return Number.NaN
  max = 1024 
  var x = new Array(max)
  var f1 = new Array(max)
  var f2 = new Array(max)
  var f3 = new Array(max)
  var v = new Array(max)                   
  step = 1
  m = 1
  bound = errorBound
  value = 0
  h = b - a
  x[0] = a
  f1[0] = f(a)
  f2[0] = f(0.5*(a + b))
  f3[0] = f(b)
  v[0] = h*(f1[0] + 4*f2[0] + f3[0])/6
  do {
    //  Estamos indo para frente ou para trás? 
    if (step == -1) {
      //  Para frente: j = m,...,max 
      step = 1
      j = m + 1
      jend = max
      m = 0
      mstart = 0 
    }
    else {
      //  Para trás: j = m,...,1
      step = -1
      j = m - 1
      jend = -1 
      m = max - 1
      mstart = max - 1
    }
    h = 0.5*h
    h6 = h/6
    bound = 0.5*bound
    do {
      left = x[j]
      right = x[j] + 0.5*h
      //  Perda de sentido?
      if (left >= right) {
        if (alerta) alert('integral: Erro 1')
        return value
      }
      fa = f(x[j] + 0.5*h)
      fb = f(x[j] + 1.5*h)
      v1 = h6*(f1[j] + 4*fa + f2[j])
      v2 = h6*(f2[j] + 4*fb + f3[j])
      error = (v[j] - v1 - v2)/15
      if (Math.abs(error) <= bound || 
        Math.abs(v1 + v2) < Math.abs(value)*ru) {
        value = ((v1 + v2) + value) - error
      }
      else {
        if (integralError) return Number.NaN
        //  Estamos sem memória?
        if (m == j) {
          left = x[j]
          right = x[j] + 0.5*h
          //  Perda de sentido?
          if (left >= right) {
            if (alerta) alert('integral: Erro 2')
            return value
          }
          value += Integral(f, left, x[j] + 2*h, bound)
        }
        else {
          //  Não
          left = x[j]
          right = x[j] + 0.125*h
          if (left >= right) {
            msg = ' O limite de erro especificado (' + 
                  integralBound + ') é muito pequeno.'
            if (alerta) alert(msg);
            integralError = true
            return Math.NaN
          }
          m1 = m + step
          x[m] = x[j]
          x[m1] = x[j] + h
          v[m] = v1
          v[m1] = v2
          f1[m] = f1[j]
          f2[m] = fa
          f3[m] = f2[j]
          f1[m1] = f2[j]
          f2[m1] = fb
          f3[m1] = f3[j]
          m += 2*step
        }
      }
      j += step
    }
    while (j != jend)
  }
  while (m != mstart)
  result = value
  if (integralError) result = Math.NaN
  return result
}

function roundingUnit() {
  var ru = 1
  do {
    ru = 0.5*ru
    var value = 1 + ru
  }
  while (value != 1)
  return 2*ru
}
//Fim de Integração Numérica