CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
PadeTableODE.cc
Go to the documentation of this file.
1 #include <cassert>
3 
5 
6 PadeTableODE::PadeTableODE(const unsigned padeRow, const unsigned padeColumn) : row_(padeRow), col_(padeColumn) {
7  if (row_ > 2U)
8  throw cms::Exception("In PadeTableODE constructor: Pade table row number out of range");
9  if (col_ > 3U)
10  throw cms::Exception("In PadeTableODE constructor: Pade table column number out of range");
11 }
12 
13 void PadeTableODE::calculate(const double tau,
14  const double currentIn,
15  const double dIdt,
16  const double d2Id2t,
17  const double* x,
18  const unsigned lenX,
19  const unsigned firstNode,
20  double* derivative) const {
21  // Check input sanity
22  if (lenX < firstNode + col_)
23  throw cms::Exception("In PadeTableODE::calculate: insufficient number of variables");
24  if (tau <= 0.0)
25  throw cms::Exception("In PadeTableODE::calculate: delay time is not positive");
26  if (col_)
27  assert(x);
29 
30  switch (col_) {
31  case 0U:
32  // Special case: no ODE to solve
33  derivative[firstNode] = 0.0;
34  switch (row_) {
35  case 2U:
36  derivative[firstNode] += 0.5 * tau * tau * d2Id2t;
37  [[fallthrough]];
38  case 1U:
39  derivative[firstNode] -= tau * dIdt;
40  [[fallthrough]];
41  case 0U:
42  derivative[firstNode] += currentIn;
43  break;
44 
45  default:
46  assert(0);
47  }
48  break;
49 
50  case 1U:
51  // First order ODE to solve
52  switch (row_) {
53  case 0U:
54  derivative[firstNode] = (currentIn - x[firstNode]) / tau;
55  break;
56 
57  case 1U:
58  derivative[firstNode] = 2.0 * (currentIn - x[firstNode]) / tau - dIdt;
59  break;
60 
61  case 2U:
62  derivative[firstNode] = 3.0 * (currentIn - x[firstNode]) / tau - 2.0 * dIdt + 0.5 * tau * d2Id2t;
63  break;
64 
65  default:
66  assert(0);
67  }
68  break;
69 
70  case 2U:
71  // Second order ODE to solve
72  derivative[firstNode] = x[firstNode + 1];
73  switch (row_) {
74  case 0U:
75  derivative[firstNode + 1] = 2.0 * (currentIn - x[firstNode] - tau * x[firstNode + 1]) / tau / tau;
76  break;
77 
78  case 1U:
79  derivative[firstNode + 1] =
80  (6.0 * (currentIn - x[firstNode]) - 2.0 * tau * dIdt - 4.0 * tau * x[firstNode + 1]) / tau / tau;
81  break;
82 
83  case 2U:
84  derivative[firstNode + 1] =
85  12.0 * (currentIn - x[firstNode]) / tau / tau - 6.0 * (x[firstNode + 1] + dIdt) / tau + d2Id2t;
86  break;
87 
88  default:
89  assert(0);
90  }
91  break;
92 
93  case 3U:
94  // Third order ODE to solve
95  derivative[firstNode] = x[firstNode + 1];
96  derivative[firstNode + 1] = x[firstNode + 2];
97  switch (row_) {
98  case 0U:
99  derivative[firstNode + 2] =
100  6.0 * (currentIn - x[firstNode] - tau * x[firstNode + 1] - 0.5 * tau * tau * x[firstNode + 2]) / tau /
101  tau / tau;
102  break;
103 
104  case 1U:
105  derivative[firstNode + 2] = 24.0 / tau / tau / tau *
106  (currentIn - x[firstNode] - 0.25 * tau * dIdt - 0.75 * tau * x[firstNode + 1] -
107  0.25 * tau * tau * x[firstNode + 2]);
108  break;
109 
110  case 2U:
111  derivative[firstNode + 2] = 60.0 / tau / tau / tau *
112  (currentIn - x[firstNode] - 0.4 * tau * dIdt + 0.05 * tau * tau * d2Id2t -
113  0.6 * tau * x[firstNode + 1] - 0.15 * tau * tau * x[firstNode + 2]);
114  break;
115 
116  default:
117  assert(0);
118  }
119  break;
120 
121  default:
122  //
123  // In principle, it is possible to proceed a bit further, but
124  // we will soon encounter difficulties. For example, row 0 and
125  // column 4 is going to generate a 4th order differential
126  // equation for which all roots of the characteristic equation
127  // still have negative real parts. The most "inconvenient" pair
128  // of roots there is (-0.270556 +- 2.50478 I) which leads
129  // to oscillations with damping. The characteristic equation
130  // of 5th and higher order ODEs are going to have roots with
131  // positive real parts. Unless additional damping is
132  // purposefully introduced into the system, numerical
133  // solutions of such equations will just blow up.
134  //
135  assert(0);
136  }
137 }
138 
139 void PadeTableODE::setParameters(const double* /* pars */, const unsigned nPars) { assert(nPars == 0U); }
void calculate(double tau, double inputCurrent, double dIdt, double d2Id2t, const double *x, unsigned lenX, unsigned firstNode, double *derivative) const
Definition: PadeTableODE.cc:13
Derivative< X, A >::type derivative(const A &_)
Definition: Derivative.h:18
void setParameters(const double *pars, unsigned nPars)
assert(be >=bs)
PadeTableODE(unsigned padeRow, unsigned padeColumn)
Definition: PadeTableODE.cc:6
unsigned col_
Definition: PadeTableODE.h:31
float x
unsigned row_
Definition: PadeTableODE.h:30