#include <PadeTableODE.h>

Public Member Functions
void	calculate (double tau, double inputCurrent, double dIdt, double d2Id2t, const double x, unsigned lenX, unsigned firstNode, double derivative) const

unsigned	getPadeColumn () const

unsigned	getPadeRow () const

unsigned	nParameters () const

	PadeTableODE (unsigned padeRow, unsigned padeColumn)

void	setParameters (const double *pars, unsigned nPars)

Private Attributes
unsigned	col_

unsigned	row_

Detailed Description

Definition at line 11 of file PadeTableODE.h.

Constructor & Destructor Documentation

PadeTableODE::PadeTableODE	(	unsigned	padeRow,
		unsigned	padeColumn
	)

Definition at line 6 of file PadeTableODE.cc.

References col_, Exception, and row_.

                                                                             : row_(padeRow), col_(padeColumn) {
   if (row_ > 2U)
     throw cms::Exception("In PadeTableODE constructor: Pade table row number out of range");
   if (col_ > 3U)
     throw cms::Exception("In PadeTableODE constructor: Pade table column number out of range");
 }

Member Function Documentation

void PadeTableODE::calculate	(	double	tau,
		double	inputCurrent,
		double	dIdt,
		double	d2Id2t,
		const double *	x,
		unsigned	lenX,
		unsigned	firstNode,
		double *	derivative
	)		const

Definition at line 13 of file PadeTableODE.cc.

References cms::cuda::assert(), col_, Exception, row_, and metsig::tau.

                                                        {
   // Check input sanity
   if (lenX < firstNode + col_)
     throw cms::Exception("In PadeTableODE::calculate: insufficient number of variables");
   if (tau <= 0.0)
     throw cms::Exception("In PadeTableODE::calculate: delay time is not positive");
   if (col_)
     assert(x);
   assert(derivative);
 
   switch (col_) {
     case 0U:
       // Special case: no ODE to solve
       derivative[firstNode] = 0.0;
       switch (row_) {
         case 2U:
           derivative[firstNode] += 0.5 * tau * tau * d2Id2t;
           [[fallthrough]];
         case 1U:
           derivative[firstNode] -= tau * dIdt;
           [[fallthrough]];
         case 0U:
           derivative[firstNode] += currentIn;
           break;
 
         default:
           assert(0);
       }
       break;
 
     case 1U:
       // First order ODE to solve
       switch (row_) {
         case 0U:
           derivative[firstNode] = (currentIn - x[firstNode]) / tau;
           break;
 
         case 1U:
           derivative[firstNode] = 2.0 * (currentIn - x[firstNode]) / tau - dIdt;
           break;
 
         case 2U:
           derivative[firstNode] = 3.0 * (currentIn - x[firstNode]) / tau - 2.0 * dIdt + 0.5 * tau * d2Id2t;
           break;
 
         default:
           assert(0);
       }
       break;
 
     case 2U:
       // Second order ODE to solve
       derivative[firstNode] = x[firstNode + 1];
       switch (row_) {
         case 0U:
           derivative[firstNode + 1] = 2.0 * (currentIn - x[firstNode] - tau * x[firstNode + 1]) / tau / tau;
           break;
 
         case 1U:
           derivative[firstNode + 1] =
               (6.0 * (currentIn - x[firstNode]) - 2.0 * tau * dIdt - 4.0 * tau * x[firstNode + 1]) / tau / tau;
           break;
 
         case 2U:
           derivative[firstNode + 1] =
               12.0 * (currentIn - x[firstNode]) / tau / tau - 6.0 * (x[firstNode + 1] + dIdt) / tau + d2Id2t;
           break;
 
         default:
           assert(0);
       }
       break;
 
     case 3U:
       // Third order ODE to solve
       derivative[firstNode] = x[firstNode + 1];
       derivative[firstNode + 1] = x[firstNode + 2];
       switch (row_) {
         case 0U:
           derivative[firstNode + 2] =
               6.0 * (currentIn - x[firstNode] - tau * x[firstNode + 1] - 0.5 * tau * tau * x[firstNode + 2]) / tau /
               tau / tau;
           break;
 
         case 1U:
           derivative[firstNode + 2] = 24.0 / tau / tau / tau *
                                       (currentIn - x[firstNode] - 0.25 * tau * dIdt - 0.75 * tau * x[firstNode + 1] -
                                        0.25 * tau * tau * x[firstNode + 2]);
           break;
 
         case 2U:
           derivative[firstNode + 2] = 60.0 / tau / tau / tau *
                                       (currentIn - x[firstNode] - 0.4 * tau * dIdt + 0.05 * tau * tau * d2Id2t -
                                        0.6 * tau * x[firstNode + 1] - 0.15 * tau * tau * x[firstNode + 2]);
           break;
 
         default:
           assert(0);
       }
       break;
 
     default:
       //
       // In principle, it is possible to proceed a bit further, but
       // we will soon encounter difficulties. For example, row 0 and
       // column 4 is going to generate a 4th order differential
       // equation for which all roots of the characteristic equation
       // still have negative real parts. The most "inconvenient" pair
       // of roots there is (-0.270556 +- 2.50478 I) which leads
       // to oscillations with damping. The characteristic equation
       // of 5th and higher order ODEs are going to have roots with
       // positive real parts. Unless additional damping is
       // purposefully introduced into the system, numerical
       // solutions of such equations will just blow up.
       //
       assert(0);
   }
 }