Method differentiate

SRC/domain/pattern/TrapezoidalTimeSeriesIntegrator.cpp:125–191

}

TimeSeries*
TrapezoidalTimeSeriesIntegrator::differentiate(TimeSeries *theSeries, double delta)
{	
  // Check for zero time step, before dividing to get number of steps
  if (delta <= 0.0) {
    opserr << "TrapezoidalTimeSeriesIntegrator::differentiate() Attempting to differentiate using time step" <<
      delta << "<= 0\n";
    return 0;
   }

  // check a TimeSeries object was passed
  if (theSeries == 0) {
    opserr << "TrapezoidalTimeSeriesIntegrator::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);

  Vector *theDif = new Vector (numSteps);

  // Check that the Vector was allocated properly
  if (theDif == 0 || theDif->Size() == 0) {
    opserr << "TrapezoidalTimeSeriesIntegrator::differentiate() Ran out of memory allocating Vector " << endln;

    if (theDif != 0)
      delete theDif;

    return 0;
  }

  double dummyTime;     // Dummy variable for integrating
  double Fi, Fj;        // function values
  double f;             // derivative value
      
  dummyTime = theSeries->getStartTime();

  f = 0.0;

  Fi = 0.0;

  opserr<<"differentiate()\n";
  for (long long i = 0; i < numSteps; i++, dummyTime += delta) {
    Fj = theSeries->getFactor(dummyTime);

    // Apply the trapezoidal rule to update the derivative
    f = 2.0 * (Fj - Fi) / delta - f;

    (*theDif)[i] = f;

    Fi = Fj;
    if (i < 10)
      opserr<<"data"<<f<<" "<<Fi<<" "<<Fj<<"\n";
  }

  // Set the method return value
  PathSeries *returnSeries = new PathSeries (0, *theDif, delta, 1.0, true, false, theSeries->getStartTime());
  delete theDif;