adevs
/home/rotten/adevs-2.6/include/adevs_bisection_event_locator.h
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