adevs
|
00001 #ifndef _adevs_bisection_event_locator_h_ 00002 #define _adevs_bisection_event_locator_h_ 00003 #include "adevs_hybrid.h" 00004 #include <cmath> 00005 00006 namespace adevs 00007 { 00015 template <typename X> class bisection_event_locator: 00016 public event_locator<X> 00017 { 00018 public: 00025 bisection_event_locator(ode_system<X>* sys, double err_tol); 00026 ~bisection_event_locator(); 00027 bool find_events(bool* events, const double* qstart, double* qend, 00028 ode_solver<X>* solver, double& h); 00029 private: 00030 double *z[2]; // State events at the start and end of [0,h] 00031 const double err_tol; // Error tolerance 00032 }; 00033 00034 template <typename X> 00035 bisection_event_locator<X>::bisection_event_locator(ode_system<X>* sys, 00036 double err_tol): 00037 event_locator<X>(sys),err_tol(err_tol) 00038 { 00039 z[0] = new double[sys->numEvents()]; 00040 z[1] = new double[sys->numEvents()]; 00041 } 00042 00043 template <typename X> 00044 bisection_event_locator<X>::~bisection_event_locator() 00045 { 00046 delete [] z[0]; delete [] z[1]; 00047 } 00048 00049 template <typename X> 00050 bool bisection_event_locator<X>::find_events(bool* events, 00051 const double* qstart, double* qend, ode_solver<X>* solver, double& h) 00052 { 00053 // Calculate the state event functions at the start of the interval 00054 this->sys->state_event_func(qstart,z[0]); 00055 // Look for events inside of the interval [0,h] 00056 for (;;) 00057 { 00058 bool event_in_interval = false, found_event = false; 00059 // Get the state event functions at the end of the interval 00060 this->sys->state_event_func(qend,z[1]); 00061 // Do any of the z functions change sign? Have we found an event 00062 for (int i = 0; i < this->sys->numEvents(); i++) 00063 { 00064 events[i] = false; 00065 if (z[1][i]*z[0][i] <= 0.0) 00066 { 00067 // There is an event at h 00068 if (fabs(z[1][i]) <= err_tol) 00069 events[i] = found_event = true; 00070 // There is an event prior to h 00071 else event_in_interval = true; 00072 } 00073 } 00074 // Guess at a new h and calculate qend for that time 00075 if (event_in_interval) 00076 { 00077 h /= 2.0; 00078 for (int i = 0; i < this->sys->numVars(); i++) 00079 qend[i] = qstart[i]; 00080 solver->advance(qend,h); 00081 } 00082 else return found_event; 00083 } 00084 // Will never reach this line 00085 return false; 00086 } 00087 00088 } // end of namespace 00089 00090 #endif