00001 #ifndef _adevs_event_locators_h_
00002 #define _adevs_event_locators_h_
00003 #include "adevs_hybrid.h"
00004 #include <cmath>
00005 #include <iostream>
00006
00007 namespace adevs
00008 {
00009
00014 template <typename X> class event_locator_impl:
00015 public event_locator<X>
00016 {
00017 public:
00024 enum Mode { INTERPOLATE, BISECTION };
00025 event_locator_impl(ode_system<X>* sys, double err_tol, Mode mode);
00026 ~event_locator_impl();
00027 bool find_events(bool* events, const double* qstart, double* qend,
00028 ode_solver<X>* solver, double& h);
00029 private:
00030 double *z[2];
00031 const double err_tol;
00032 Mode mode;
00033
00034 int sign(double x) const
00035 {
00036 if (x < 0.0) return -1;
00037 else if (x > 0.0) return 1;
00038 else return 0;
00039 }
00040 };
00041
00042 template <typename X>
00043 event_locator_impl<X>::event_locator_impl(ode_system<X>* sys,
00044 double err_tol, Mode mode):
00045 event_locator<X>(sys),
00046 err_tol(err_tol),
00047 mode(mode)
00048 {
00049 z[0] = new double[sys->numEvents()];
00050 z[1] = new double[sys->numEvents()];
00051 }
00052
00053 template <typename X>
00054 event_locator_impl<X>::~event_locator_impl()
00055 {
00056 delete [] z[0]; delete [] z[1];
00057 }
00058
00059 template <typename X>
00060 bool event_locator_impl<X>::find_events(bool* events,
00061 const double* qstart, double* qend, ode_solver<X>* solver, double& h)
00062 {
00063
00064
00065 this->sys->state_event_func(qstart,z[0]);
00066
00067 for (;;)
00068 {
00069 double tguess = h;
00070 bool event_in_interval = false, found_event = false;
00071 this->sys->state_event_func(qend,z[1]);
00072
00073 for (int i = 0; i < this->sys->numEvents(); i++)
00074 {
00075 events[i] = false;
00076 if (sign(z[1][i]) != sign(z[0][i]))
00077 {
00078
00079 if (fabs(z[1][i]) <= err_tol)
00080 {
00081 events[i] = found_event = true;
00082 }
00083
00084 else
00085 {
00086 if (mode == INTERPOLATE)
00087 {
00088 double tcandidate = z[0][i]*h/(z[0][i]-z[1][i]);
00089
00090 if (tcandidate < h/4.0) tcandidate = h/4.0;
00091 if (tcandidate < tguess) tguess = tcandidate;
00092 }
00093 event_in_interval = true;
00094 }
00095 }
00096 }
00097
00098 if (event_in_interval)
00099 {
00100 if (mode == BISECTION) h /= 2.0;
00101 else h = tguess;
00102 for (int i = 0; i < this->sys->numVars(); i++)
00103 qend[i] = qstart[i];
00104 solver->advance(qend,h);
00105 }
00106 else return found_event;
00107 }
00108
00109 return false;
00110 }
00111
00115 template <typename X>
00116 class bisection_event_locator:
00117 public event_locator_impl<X>
00118 {
00119 public:
00120 bisection_event_locator(ode_system<X>* sys, double err_tol):
00121 event_locator_impl<X>(
00122 sys,err_tol,event_locator_impl<X>::BISECTION){}
00123 };
00124
00128 template <typename X>
00129 class linear_event_locator:
00130 public event_locator_impl<X>
00131 {
00132 public:
00133 linear_event_locator(ode_system<X>* sys, double err_tol):
00134 event_locator_impl<X>(
00135 sys,err_tol,event_locator_impl<X>::INTERPOLATE){}
00136 };
00137
00138 }
00139
00140 #endif