adevs
/home/rotten/adevs-2.6/include/adevs_linear_event_locator.h
00001 #ifndef _adevs_linear_event_locator_h_
00002 #define _adevs_linear_event_locator_h_
00003 #include "adevs_hybrid.h"
00004 #include <cmath>
00005 
00006 namespace adevs
00007 {
00008 
00013 template <typename X> class linear_event_locator:
00014         public event_locator<X>
00015 {
00016         public:
00023                 linear_event_locator(ode_system<X>* sys, double err_tol);
00024                 ~linear_event_locator();
00025                 bool find_events(bool* events, const double* qstart, double* qend,
00026                                 ode_solver<X>* solver, double& h);
00027         private:
00028                 double *z[2]; // State events at the start and end of [0,h]
00029                 const double err_tol; // Error tolerance
00030 };
00031 
00032 template <typename X>
00033 linear_event_locator<X>::linear_event_locator(ode_system<X>* sys,
00034                 double err_tol):
00035         event_locator<X>(sys),err_tol(err_tol)
00036 {
00037         z[0] = new double[sys->numEvents()];
00038         z[1] = new double[sys->numEvents()];
00039 }
00040 
00041 template <typename X>
00042 linear_event_locator<X>::~linear_event_locator()
00043 {
00044         delete [] z[0]; delete [] z[1];
00045 }
00046 
00047 template <typename X>
00048 bool linear_event_locator<X>::find_events(bool* events,
00049         const double* qstart, double* qend, ode_solver<X>* solver, double& h)
00050 {
00051         // Look for events at the start of the interval
00052         this->sys->state_event_func(qstart,z[0]);
00053         for (int i = 0; i < this->sys->numEvents(); i++) {
00054                 events[i] = fabs(z[0][i]) <= err_tol;
00055                 // If an event was found, the event time is zero
00056                 if (events[i]) h = 0.0;
00057         }
00058         // If an event was found at zero, put qstart in qend and return
00059         if (h == 0.0) {
00060                 for (int i = 0; i < this->sys->numVars(); i++) qend[i] = qstart[i];
00061                 return true;
00062         }
00063         // Look for events inside of the interval [0,h]
00064         for (;;) {
00065                 double tguess = h;
00066                 bool event_in_interval = false, found_event = false;
00067                 this->sys->state_event_func(qend,z[1]);
00068                 for (int i = 0; i < this->sys->numEvents(); i++) {
00069                         if ((events[i] = fabs(z[1][i]) <= err_tol)) found_event = true;
00070                         else if (z[0][i]*z[1][i] < 0.0) {
00071                                 double tcandidate = z[0][i]*h/(z[0][i]-z[1][i]);
00072                                 if (tcandidate < tguess) tguess = tcandidate;
00073                                 event_in_interval = true;
00074                         }
00075                 }
00076                 // Calculate a new solution at tguess if an event was found
00077                 if (event_in_interval) {
00078                         h = tguess;
00079                         for (int i = 0; i < this->sys->numVars(); i++)
00080                                 qend[i] = qstart[i];
00081                         solver->advance(qend,h);
00082                 }
00083                 // Stop when an event is located or is not detected in the interval
00084                 else return found_event;
00085         }
00086         // Will never reach this line
00087         return false;
00088 }
00089 
00090 } // end of namespace 
00091 
00092 #endif