Method differentiate

SRC/domain/pattern/SimpsonTimeSeriesIntegrator.cpp:128–196

TimeSeries* SimpsonTimeSeriesIntegrator::differentiate(TimeSeries *theSeries, double delta)
{
  // check for zero time step, before dividing to get number of steps
  if (delta <= 0.0)  {
    opserr << "SimpsonTimeSeriesIntegrator::differentiate() - attempting to differentiate using time step "
      << delta << "<= 0.0.\n";
    return 0;
  }
  
  // check a TimeSeries object was passed
  if (theSeries == 0)  {
    opserr << "SimpsonTimeSeriesIntegrator::differentiate() - no TimeSeries passed.\n";
    return 0;
  }
  
  // add one to get ceiling out of type cast
  long long numSteps = (long long)(theSeries->getDuration()/delta + 1.0);
  
  // create new vector for integrated values
  Vector *theDif = new Vector(numSteps);
  
  // check that the Vector was allocated properly
  if (theDif == 0 || theDif->Size() == 0)  {
    opserr << "SimpsonTimeSeriesIntegrator::differentiate() - ran out of memory allocating Vector " << endln;

    if (theDif != 0)
      delete theDif;
    
    return 0;
  }
  
  double dummyTime;
  double Fi, Fj, Fk;  //function values
  double fi, fj, fk;  //derivative values

  dummyTime = theSeries->getStartTime();

  Fi = 0.0;
  Fj = 0.0;

  fi = 0.0;
  fj = 0.0;

  for (long long i = 0; i < numSteps; i++, dummyTime += delta)  {

    Fk = theSeries->getFactor(dummyTime);
    // Apply the Simpson's rule to update the derivative
    fk = 3.0 * (Fk - Fi) / delta - fi - 4.0 * fj;

    (*theDif)[i] = fk;

    fi = fj;
    fj = fk;

    Fi = Fj;
    Fj = Fk;
  }