adevs
|
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