00001
00031 #ifndef _adevs_event_locators_h_
00032 #define _adevs_event_locators_h_
00033 #include "adevs_hybrid.h"
00034 #include <cmath>
00035 #include <iostream>
00036
00037 namespace adevs
00038 {
00039
00044 template <typename X> class event_locator_impl:
00045 public event_locator<X>
00046 {
00047 public:
00054 enum Mode { INTERPOLATE, BISECTION, DISCONTINUOUS };
00055 event_locator_impl(ode_system<X>* sys, double err_tol, Mode mode);
00056 ~event_locator_impl();
00057 bool find_events(bool* events, const double* qstart, double* qend,
00058 ode_solver<X>* solver, double& h);
00059 private:
00060 double *z[2];
00061 const double err_tol;
00062 Mode mode;
00063
00064 int sign(double x) const
00065 {
00066 if (x < 0.0) return -1;
00067 else if (x > 0.0) return 1;
00068 else return 0;
00069 }
00070 };
00071
00072 template <typename X>
00073 event_locator_impl<X>::event_locator_impl(ode_system<X>* sys,
00074 double err_tol, Mode mode):
00075 event_locator<X>(sys),
00076 err_tol(err_tol),
00077 mode(mode)
00078 {
00079 z[0] = new double[sys->numEvents()];
00080 z[1] = new double[sys->numEvents()];
00081 }
00082
00083 template <typename X>
00084 event_locator_impl<X>::~event_locator_impl()
00085 {
00086 delete [] z[0]; delete [] z[1];
00087 }
00088
00089 template <typename X>
00090 bool event_locator_impl<X>::find_events(bool* events,
00091 const double* qstart, double* qend, ode_solver<X>* solver, double& h)
00092 {
00093
00094
00095 this->sys->state_event_func(qstart,z[0]);
00096
00097 for (;;)
00098 {
00099 double tguess = h;
00100 bool event_in_interval = false, found_event = false;
00101 this->sys->state_event_func(qend,z[1]);
00102
00103 for (int i = 0; i < this->sys->numEvents(); i++)
00104 {
00105 events[i] = false;
00106 if (sign(z[1][i]) != sign(z[0][i]))
00107 {
00108
00109 if (
00110
00111 ((mode != DISCONTINUOUS) && (fabs(z[1][i]) <= err_tol)) ||
00112
00113 ((mode == DISCONTINUOUS) && (h <= err_tol))
00114 )
00115 {
00116 events[i] = found_event = true;
00117 }
00118
00119 else
00120 {
00121 if (mode == INTERPOLATE)
00122 {
00123 double tcandidate = z[0][i]*h/(z[0][i]-z[1][i]);
00124
00125 if (tcandidate < h/4.0) tcandidate = h/4.0;
00126 if (tcandidate < tguess) tguess = tcandidate;
00127 }
00128 event_in_interval = true;
00129 }
00130 }
00131 }
00132
00133 if (event_in_interval)
00134 {
00135 if (mode == BISECTION || mode == DISCONTINUOUS) h /= 2.0;
00136 else h = tguess;
00137 for (int i = 0; i < this->sys->numVars(); i++)
00138 qend[i] = qstart[i];
00139 solver->advance(qend,h);
00140 }
00141 else return found_event;
00142 }
00143
00144 return false;
00145 }
00146
00150 template <typename X>
00151 class bisection_event_locator:
00152 public event_locator_impl<X>
00153 {
00154 public:
00155 bisection_event_locator(ode_system<X>* sys, double err_tol):
00156 event_locator_impl<X>(
00157 sys,err_tol,event_locator_impl<X>::BISECTION){}
00158 };
00159
00163 template <typename X>
00164 class linear_event_locator:
00165 public event_locator_impl<X>
00166 {
00167 public:
00168 linear_event_locator(ode_system<X>* sys, double err_tol):
00169 event_locator_impl<X>(
00170 sys,err_tol,event_locator_impl<X>::INTERPOLATE){}
00171 };
00172
00176 template <typename X>
00177 class discontinuous_event_locator:
00178 public event_locator_impl<X>
00179 {
00180 public:
00181 discontinuous_event_locator(ode_system<X>* sys, double err_tol):
00182 event_locator_impl<X>(
00183 sys,err_tol,event_locator_impl<X>::DISCONTINUOUS){}
00184 };
00185
00186 }
00187
00188 #endif