| 55 | |
| 56 | |
| 57 | TimeSeries* SimpsonTimeSeriesIntegrator::integrate(TimeSeries *theSeries, double delta) |
| 58 | { |
| 59 | // check for zero time step, before dividing to get number of steps |
| 60 | if (delta <= 0.0) { |
| 61 | opserr << "SimpsonTimeSeriesIntegrator::integrate() - attempting to integrate using time step " |
| 62 | << delta << "<= 0.0.\n"; |
| 63 | return 0; |
| 64 | } |
| 65 | |
| 66 | // check a TimeSeries object was passed |
| 67 | if (theSeries == 0) { |
| 68 | opserr << "SimpsonTimeSeriesIntegrator::integrate() - no TimeSeries passed.\n"; |
| 69 | return 0; |
| 70 | } |
| 71 | |
| 72 | // add one to get ceiling out of type cast |
| 73 | long long numSteps = (long long)(theSeries->getDuration()/delta + 1.0); |
| 74 | |
| 75 | // create new vector for integrated values |
| 76 | Vector *theInt = new Vector(numSteps); |
| 77 | |
| 78 | // check that the Vector was allocated properly |
| 79 | if (theInt == 0 || theInt->Size() == 0) { |
| 80 | opserr << "SimpsonTimeSeriesIntegrator::integrate() - ran out of memory allocating Vector " << endln; |
| 81 | |
| 82 | if (theInt != 0) |
| 83 | delete theInt; |
| 84 | |
| 85 | return 0; |
| 86 | } |
| 87 | |
| 88 | double dummyTime; |
| 89 | double fi, fj, fk; //function values |
| 90 | double Fi, Fj, Fk; //intergral values |
| 91 | |
| 92 | dummyTime = theSeries->getStartTime(); |
| 93 | |
| 94 | Fi = 0.0; |
| 95 | Fj = 0.0; |
| 96 | |
| 97 | fi = 0.0; |
| 98 | fj = 0.0; |
| 99 | |
| 100 | for (long long i = 0; i < numSteps; i++, dummyTime += delta) { |
| 101 | |
| 102 | fk = theSeries->getFactor(dummyTime); |
| 103 | // Apply the Simpson's rule to update the integral |
| 104 | Fk = Fi + delta / 3.0 * (fi + 4.0 * fj + fk); |
| 105 | |
| 106 | (*theInt)[i] = Fk; |
| 107 | |
| 108 | fi = fj; |
| 109 | fj = fk; |
| 110 | |
| 111 | Fi = Fj; |
| 112 | Fj = Fk; |
| 113 | } |
| 114 |
nothing calls this directly
no test coverage detected