lib/e_elemnt.cc

   1 /*                             -*- C++ -*-
   2  * Copyright (C) 2001 Albert Davis
   3  * Author: Albert Davis <aldavis@gnu.org>
   4  *
   5  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20  * 02110-1301, USA.
  21  *------------------------------------------------------------------
  22  * Base class for elements of a circuit
  23  */
  24 #include "m_divdiff.h"
  25 #include "u_xprobe.h"
  26 #include "e_aux.h"
  27 #include "e_elemnt.h"
  28 /*--------------------------------------------------------------------------*/
  29 ELEMENT::ELEMENT()
  30   :COMPONENT(),
  31    _input_order(1),
  32    _trsteporder(-1),
  33    _loaditer(0),
  34    _m0(),
  35    _m1(),
  36    _loss0(0.),
  37    _loss1(0.),
  38    _acg(0.),
  39    _ev(0.),
  40    _dt(0.),
  41    _discont(disNONE)
  42 {
  43   _n = _nodes;
  44   assert(_y[0].x == 0. && _y[0].f0 == 0. && _y[0].f1 == 0.);
  45   assert(_y1 == _y[0]);
  46
  47   std::fill_n(_time, int(OPT::_keep_time_steps), 0.);
  48 }
  49 /*--------------------------------------------------------------------------*/
  50 ELEMENT::ELEMENT(const ELEMENT& p)
  51   :COMPONENT(p),
  52    _input_order(p._input_order),
  53    _trsteporder(-1),
  54    _loaditer(0),
  55    _m0(),
  56    _m1(),
  57    _loss0(p._loss0),
  58    _loss1(p._loss1),
  59    _acg(0.),
  60    _ev(0.),
  61    _dt(0.),
  62    _discont(disNONE)
  63 {
  64   trace0(long_label().c_str());
  65   _n = _nodes;
  66   if (p._n == p._nodes) {
  67     for (int ii = 0;  ii < NODES_PER_BRANCH;  ++ii) {
  68       _n[ii] = p._n[ii];
  69     }
  70   }else{
  71     assert(p._nodes);
  72     // the constructor for a derived class will take care of it
  73   }
  74   assert(_y[0].x == 0. && _y[0].f0 == 0. && _y[0].f1 == 0.);
  75   assert(_y1 == _y[0]);
  76
  77   notstd::copy_n(p._time, int(OPT::_keep_time_steps), _time);
  78 }
  79 /*--------------------------------------------------------------------------*/
  80 bool ELEMENT::skip_dev_type(CS& cmd)
  81 {
  82   return cmd.umatch(dev_type() + ' ');
  83 }
  84 /*--------------------------------------------------------------------------*/
  85 void ELEMENT::expand()
  86 {
  87   trace3("ELEMENT::expand", long_label(), input_order(), hp(common()));
  88   assert(!subckt() || subckt()->size() < 2);
  89   COMPONENT::expand();
  90
  91 #if 0
  92   if(input_order()>1){
  93     if (!subckt()) new_subckt();
  94     trace3("DEV_VCVS::expand branch", long_label(), hp(common()), input_order());
  95     ELEMENT* more = prechecked_cast<ELEMENT*>(clone());
  96     node_t nodes[input_order()*2];
  97     nodes[0] = _n[0];
  98     nodes[1] = _n[1];
  99     notstd::copy_n(_n+4, input_order()*2-2, nodes+2);
 100     assert(mutable_common()==common());
 101     more->set_parameters("branch", this, mutable_common(), 0., 0, NULL, input_order()*2, nodes);
 102     more->set_input_order(input_order()-1);
 103     subckt()->push_front(more);
 104     subckt()->set_slave();
 105
 106     assert(more->common() == common());
 107     subckt()->expand();
 108     trace2("DEV_VCVS::expand branch", hp(more->common()), hp(common()));
 109     assert(more->common() == common());
 110   }
 111 #else
 112
 113   if (input_order()>1) {
 114     if (!subckt()) new_subckt();
 115     trace3("ELEMENT::expand branch", long_label(), hp(common()), input_order());
 116     ELEMENT* more = prechecked_cast<ELEMENT*>(clone());
 117     node_t* nodes = new node_t[input_order()*2];
 118     nodes[0] = _n[0];
 119     nodes[1] = _n[1];
 120     notstd::copy_n(_n+4, input_order()*2-2, nodes+2);
 121     delete[] nodes;
 122     assert(mutable_common()==common());
 123     more->set_input_order(input_order()-1);
 124     more->set_parameters("branch", this, mutable_common(), 0., 0, NULL, 2+(input_order()-1)*ports_per_input(), nodes);
 125
 126     for (unsigned i=1; i<num_current_ports(); i++) {
 127       string a = current_port_value(i);
 128       more->set_port_by_index(i+1, a);
 129     }
 130
 131 //    for(unsigned i=4; i<net_nodes(); i++){untested();
 132 //      string a = port_value(i);
 133 //      trace2("ELEMENT::expand ", i, a);
 134 //      more->set_port_by_index(i-2, a);
 135 //    }
 136
 137     subckt()->push_front(more);
 138     subckt()->set_slave();
 139
 140     assert(more->common() == common());
 141     trace2("DEV_CCSRC::expand branch", hp(more->common()), hp(common()));
 142     subckt()->expand();
 143     // assert(more->common() == common()); maybe not.
 144   }else{
 145   }
 146 #endif
 147 }
 148 /*--------------------------------------------------------------------------*/
 149 void ELEMENT::precalc_last()
 150 {
 151   COMPONENT::precalc_last();
 152
 153   //BUG// This is needed for AC analysis without doing op (or dc or tran ...) first.
 154   // Something like it should be moved to ac_begin.
 155   if (_sim->has_op() == s_NONE) {
 156     _y[0].x  = 0.;
 157     _y[0].f0 = LINEAR;
 158     _y[0].f1 = value();
 159   }else{
 160   }
 161 }
 162 /*--------------------------------------------------------------------------*/
 163 void ELEMENT::tr_begin()
 164 {
 165   _trsteporder = -1;
 166   trace4("ELEMENT::tr_begin", long_label(), value(), has_common(), tr_needs_eval());
 167   if (_sim->_time0) { untested();
 168   }
 169   _time[0] = _sim->_time0;
 170   _y[0].x  = 0.;
 171   _y[0].f0 = LINEAR;
 172   _y[0].f1 = value();
 173   _y1 = _y[0];
 174   for (int i=1; i<OPT::_keep_time_steps; ++i) {
 175     _time[i] = _sim->_time0;
 176     _y[i] = FPOLY1(0., 0., 0.);
 177   }
 178   _dt = NOT_VALID;
 179 }
 180 /*--------------------------------------------------------------------------*/
 181 void ELEMENT::tr_restore()
 182 {
 183   if (_time[0] > _sim->_time0) {
 184     // _freezetime
 185     incomplete();
 186     //BUG// wrong values in _time[]
 187     for (int i=0  ; i<OPT::_keep_time_steps-1; ++i) {
 188       _time[i] = _time[i+1];
 189       _y[i] = _y[i+1];
 190     }
 191     _time[OPT::_keep_time_steps-1] = 0.;
 192     _y[OPT::_keep_time_steps-1]    = FPOLY1(0., 0., 0.);
 193   }else if (_time[0] == _sim->_time0) {
 194     // the usual continue where the last one left off
 195   }else{unreachable();
 196     // skipping ahead, not implemented
 197   }
 198
 199   //assert(_time[0] == _sim->_time0);
 200   if (is_constant()){
 201
 202   }else if (_time[0] != _sim->_time0) {
 203     error(bDANGER, "//BUG// %s restore time mismatch.  t0=%.12f, s->t=%.12f\n", long_label().c_str(),
 204         _time[0], _sim->_time0);
 205
 206     _trsteporder = -1;
 207     for (int i=OPT::_keep_time_steps-1; i>=0; --i) {
 208       _time[i] = _sim->_time0;
 209     }
 210     for (int i=1  ; i<OPT::_keep_time_steps-1; ++i) {
 211       _y[i] = _y[0];
 212     }
 213     //BUG// happens when continuing after a ^c,
 214     // when the last step was not printed
 215     //BUG// also happens with _freezetime
 216     // _time[0] is the non-printed time.  _sim->_time0 is the printed time.
 217   }else{
 218   }
 219
 220   for (int i=OPT::_keep_time_steps-1; i>0; --i) {
 221     assert(_time[i] < _time[i-1] || _time[i] == 0. || _time[i] == _sim->_time0);
 222   }
 223 }
 224 /*--------------------------------------------------------------------------*/
 225 void ELEMENT::dc_advance()
 226 {
 227   assert(_sim->_time0 == 0.); // DC
 228
 229   for (int i=OPT::_keep_time_steps-1; i>=0; --i) {
 230     if(_time[i] == _sim->_time0){
 231     }else{
 232       // in _dc_cont case allow to skip dc analysis....
 233       _time[i] = _sim->_time0;
 234     }
 235   }
 236
 237   _dt = NOT_VALID;
 238 }
 239 /*--------------------------------------------------------------------------*/
 240 void ELEMENT::tr_advance()
 241 {
 242   assert(_time[0] < _sim->_time0); // moving forward
 243
 244   for (int i=OPT::_keep_time_steps-1; i>0; --i) {
 245     assert(_time[i] < _time[i-1] || _time[i] == 0. || _time[i] < _sim->_time0);
 246     _time[i] = _time[i-1];
 247     _y[i] = _y[i-1];
 248   }
 249   _time[0] = _sim->_time0;
 250
 251   _dt = _time[0] - _time[1];
 252
 253   DISCONT d = disNONE;
 254   for (unsigned i=0;i < net_nodes(); ++i) {
 255     d |= _n[i]->discont();
 256     assert(_n[i].n_() != &ground_node || !_n[i]->discont());
 257   }
 258
 259   if (d & 1) {
 260     _trsteporder = -1;
 261   }else if (d & 2) {
 262     _trsteporder = 0;
 263   }else if (_trsteporder < (int)OPT::trsteporder) {
 264     ++_trsteporder;
 265   }
 266 }
 267 /*--------------------------------------------------------------------------*/
 268 void ELEMENT::tr_regress()
 269 {
 270   trace4("ELEMENT::tr_regress", long_label(), _sim->_time0, _time[0], _time[1]);
 271   assert(is_constant() || _time[0] >= _sim->_time0); // moving backwards
 272   assert(_time[1] <= _sim->_time0); // but not too far backwards
 273
 274   for (int i=OPT::_keep_time_steps-1; i>0; --i) {
 275     trace1("i", i);
 276     assert(_time[i] < _time[i-1] || _time[i] == 0. || _time[i] == _sim->_time0);
 277   }
 278   _time[0] = _sim->_time0;
 279
 280   _dt = _time[0] - _time[1];
 281 }
 282 /*--------------------------------------------------------------------------*/
 283 /*--------------------------------------------------------------------------*/
 284 /* some regress hack??
 285 void ELEMENT::tr_regress()
 286 {
 287   trace3("ELEMENT::tr_regress", _time[0], _time[1], _sim->_time0);
 288   assert(_time[0] >= _sim->_time0); // moving backwards
 289 //  assert(_time[1] <= _sim->_time0); // but not too far backwards
 290   if (  _time[1] <  _sim->_time0){
 291     error(bDANGER, "regress %s //BUG// some time mismatch.  _time[1]=%.17f, _sim->_time0=%.17f\n",
 292         (short_label()).c_str(),
 293     _time[1], _sim->_time0);
 294   }
 295
 296   for (int i=OPT::_keep_time_steps-1; i>0; --i) {
 297     assert(_time[i] < _time[i-1] || _time[i] == 0.);
 298   }
 299   _time[0] = _sim->_time0;
 300
 301   _dt = _time[0] - _time[1];
 302 }
 303 */
 304 /*--------------------------------------------------------------------------*/
 305 TIME_PAIR ELEMENT::tr_review()
 306 {
 307   trace2("ELEMENT::tr_review", long_label(), order());
 308   COMPONENT::tr_review();
 309   if (order() > 0 && _y[0].f0 != LINEAR) {
 310     double timestep = tr_review_trunc_error(_y);
 311     double newtime = tr_review_check_and_convert(timestep);
 312     _time_by.min_error_estimate(newtime);
 313   }else{
 314   }
 315
 316   return _time_by;
 317 }
 318 /*--------------------------------------------------------------------------*/
 319 void ELEMENT::tr_accept()
 320 {
 321   COMPONENT::tr_accept();
 322   if (_discont) {
 323     for (unsigned int i=0; i < matrix_nodes(); ++i) {
 324       _n[i]->discont(_discont);
 325     }
 326   } else {
 327   }
 328 }
 329 /*--------------------------------------------------------------------------*/
 330 void ELEMENT::tr_iwant_matrix_passive()
 331 {
 332   assert(matrix_nodes() == 2);
 333   assert(is_device());
 334   //assert(!subckt()); ok for subckt to exist for logic
 335   trace3("ELEMENT::tr_iwant_matrix_passive", long_label().c_str(), _n[OUT1].m_(), _n[OUT2].m_());
 336
 337   if (_n[OUT1].m_() == INVALID_NODE) {
 338     std::cerr << "ELEMENT::tr_iwant_matrix_passive " << long_label() << "\n";
 339     exit(4);
 340   }
 341   assert(_n[OUT2].m_() != INVALID_NODE);
 342
 343   _sim->_aa.iwant(_n[OUT1].m_(),_n[OUT2].m_());
 344   _sim->_lu.iwant(_n[OUT1].m_(),_n[OUT2].m_());
 345 }
 346 /*--------------------------------------------------------------------------*/
 347 void ELEMENT::tr_iwant_matrix_active()
 348 {
 349   trace3("ELEMENT::tr_iwant_matrix_active", long_label(), _n[OUT1].m_(), _n[OUT2].m_());
 350   assert(matrix_nodes() >= 4);
 351   assert(is_device());
 352   // assert(!subckt());
 353
 354   assert(_n[OUT1].m_() != INVALID_NODE);
 355   assert(_n[OUT2].m_() != INVALID_NODE);
 356   assert(_n[IN1].m_() != INVALID_NODE);
 357   assert(_n[IN2].m_() != INVALID_NODE);
 358
 359   //_sim->_aa.iwant(_n[OUT1].m_(),_n[OUT2].m_());
 360   _sim->_aa.iwant(_n[OUT1].m_(),_n[IN1].m_());
 361   _sim->_aa.iwant(_n[OUT1].m_(),_n[IN2].m_());
 362   _sim->_aa.iwant(_n[OUT2].m_(),_n[IN1].m_());
 363   _sim->_aa.iwant(_n[OUT2].m_(),_n[IN2].m_());
 364   //_sim->_aa.iwant(_n[IN1].m_(),_n[IN2].m_());
 365
 366   //_sim->_lu.iwant(_n[OUT1].m_(),_n[OUT2].m_());
 367   _sim->_lu.iwant(_n[OUT1].m_(),_n[IN1].m_());
 368   _sim->_lu.iwant(_n[OUT1].m_(),_n[IN2].m_());
 369   _sim->_lu.iwant(_n[OUT2].m_(),_n[IN1].m_());
 370   _sim->_lu.iwant(_n[OUT2].m_(),_n[IN2].m_());
 371   //_sim->_lu.iwant(_n[IN1].m_(),_n[IN2].m_());
 372 }
 373 /*--------------------------------------------------------------------------*/
 374 void ELEMENT::tr_iwant_matrix_extended()
 375 {
 376   trace4("ELEMENT::tr_iwant_matrix_extended", long_label(), dev_type(),  ext_nodes(), int_nodes());
 377   assert(is_device());
 378   assert(ext_nodes() + int_nodes() == matrix_nodes());
 379
 380   for (uint_t ii = 0;  ii < matrix_nodes();  ++ii) {
 381       trace2("ELEMENT::tr_iwant_matrix_extended", ii, _n[ii].m_() );
 382   }
 383   for (uint_t ii = 0;  ii < matrix_nodes();  ++ii) {
 384     if (_n[ii].m_()  != INVALID_NODE) {
 385       for (uint_t jj = 0;  jj < ii ;  ++jj) {
 386         _sim->_aa.iwant(_n[ii].m_(),_n[jj].m_());
 387         _sim->_lu.iwant(_n[ii].m_(),_n[jj].m_());
 388       }
 389     }else{itested();
 390       // node 1 is grounded or invalid
 391     }
 392   }
 393 }
 394 /*--------------------------------------------------------------------------*/
 395 void ELEMENT::ac_iwant_matrix_passive()
 396 {
 397   trace2(long_label().c_str(), _n[OUT1].m_(), _n[OUT2].m_());
 398   _sim->_acx.iwant(_n[OUT1].m_(),_n[OUT2].m_());
 399 }
 400 /*--------------------------------------------------------------------------*/
 401 void ELEMENT::ac_iwant_matrix_active()
 402 {
 403   //_sim->_acx.iwant(_n[OUT1].m_(),_n[OUT2].m_());
 404   _sim->_acx.iwant(_n[OUT1].m_(),_n[IN1].m_());
 405   _sim->_acx.iwant(_n[OUT1].m_(),_n[IN2].m_());
 406   _sim->_acx.iwant(_n[OUT2].m_(),_n[IN1].m_());
 407   _sim->_acx.iwant(_n[OUT2].m_(),_n[IN2].m_());
 408   //_sim->_acx.iwant(_n[IN1].m_(),_n[IN2].m_());
 409 }
 410 /*--------------------------------------------------------------------------*/
 411 void ELEMENT::ac_iwant_matrix_extended()
 412 {
 413   assert(is_device());
 414   assert(ext_nodes() + int_nodes() == matrix_nodes());
 415
 416   for (uint_t ii = 0;  ii < matrix_nodes();  ++ii) {
 417     if (_n[ii].m_() !=INVALID_NODE) {
 418       for (uint_t jj = 0;  jj < ii ;  ++jj) {
 419         _sim->_acx.iwant(_n[ii].m_(),_n[jj].m_());
 420       }
 421     }else{itested();
 422       // node 1 is grounded or invalid
 423     }
 424   }
 425 }
 426 /*--------------------------------------------------------------------------*/
 427 hp_float_t ELEMENT::tr_amps()const
 428 {
 429   hp_float_t root = fixzero((_loss0 * tr_outvolts() + _m0.c1 * tr_involts() + _m0.c0),
 430       _m0.c0);
 431
 432   assert (root==root);
 433   return (root);
 434 }
 435 /*--------------------------------------------------------------------------*/
 436 void ELEMENT::tr_save_amps(int n){
 437
 438   trace2( "ELEMENT::tr_save_amps ",  n, _sim->_tt_accepted);
 439   hp_float_t tramps = tr_amps();
 440
 441   hp_float_t _tr_amps_diff_cur;
 442   hp_float_t _tr_amps_sum_cur;
 443
 444   if( _amps_new == 0 ){
 445     incomplete(); // somebody forgot to initilize
 446     return;
 447   }
 448
 449   if(_amps!=0) {
 450     trace3( "saving _amps[ ",  n ,  net_nodes(), _amps[n] );
 451     _tr_amps_sum_cur = fabs( _amps[n] ) + fabs( tramps );
 452     _tr_amps_diff_cur = fabs( _amps[n] - tramps );
 453     //          std::cerr << "ELEMENT::tr_save_amps: not first. diff= " << _tr_amps_diff_cur << "\n";
 454   } else {
 455     trace2( "dummy _amps[ ",  n ,  net_nodes() );
 456     _tr_amps_diff_cur = 0;
 457     _tr_amps_sum_cur = 1;
 458   }
 459
 460   //    std::cerr << short_label() << ": saving _amps[ " << n << " ]" << _amps << " \n";
 461   trace2( (short_label() + " have ").c_str(), tramps, n);
 462   _amps_new[ n ]= tramps;
 463
 464
 465   trace1("ELEMENT::tr_save_amps", _tr_amps_sum_cur);
 466   tr_behaviour_del = _tr_amps_diff_cur;
 467   tr_behaviour_rel = _tr_amps_diff_cur / _tr_amps_sum_cur;
 468
 469
 470   //    tr_behaviour_del = 0;
 471   //    tr_behaviour_rel = 0;
 472
 473   COMPONENT::tr_behaviour();
 474
 475 }
 476 /*--------------------------------------------------------------------------*/
 477 double ELEMENT::tr_probe_num(const std::string& x)const
 478 {
 479   if (Umatch(x, "v{out} ")) {
 480     return tr_outvolts();
 481   }else if (Umatch(x, "amps_diff ")) {
 482     return tr_amps_diff();
 483   }else if (Umatch(x, "vi{n} |v0 ")) {
 484     return tr_involts();
 485 //  }else if (Umatch(x, "dv ")) { untested();
 486 //      return ( _n[OUT2].vt1() - _n[OUT1].vt1() );
 487   }else if (Umatch(x, "v1 ")) {
 488     if(input_order()>1){
 489       const ELEMENT* more = prechecked_cast<const ELEMENT*>(*(subckt()->begin()));
 490       return more->tr_involts();
 491     }
 492     return NOT_VALID;
 493   }else if (Umatch(x, "i ")) {
 494     return tr_amps();
 495   }else if (Umatch(x, "i1 ")) {
 496     const ELEMENT* more = prechecked_cast<const ELEMENT*>(*(subckt()->begin()));
 497     if(input_order()>1){
 498       assert(more);
 499       return more->tr_amps();
 500     }
 501     return NOT_VALID;
 502   }else if (Umatch(x, "p ")) {
 503     return tr_amps() * tr_outvolts();
 504   }else if (Umatch(x, "pd ")) {
 505     double p = tr_amps() * tr_outvolts();
 506     return (p > 0.) ? p : 0.;
 507   }else if (Umatch(x, "ps ")) {
 508     double p = tr_amps() * tr_outvolts();
 509     return (p < 0.) ? -p : 0.;
 510   }else if (Umatch(x, "in{put} ")) {
 511     return _y[0].x;
 512   }else if (Umatch(x, "f ")) {
 513     if (_y[0].f0 == LINEAR) {itested();
 514       return _y[0].x * _y[0].f1;
 515     }else{
 516       return _y[0].f0;
 517     }
 518   }else if (Umatch(x, "ev |df ")) {
 519     return _y[0].f1;
 520   }else if (Umatch(x, "dx ")) {
 521     return _y1.x-_y[0].x;
 522   }else if (Umatch(x, "res ")) {
 523     return _y1.f0-_y[0].f0;
 524   }else if (Umatch(x, "conv{erged} ")) {
 525     // assert(converged()==conv_check()); no.
 526     return converged() + 10*conv_check();
 527   }else if (Umatch(x, "nv ")) {
 528     return value();
 529   }else if (Umatch(x, "eiv ")) {untested();
 530     return _m0.x;
 531   }else if (Umatch(x, "y |m0c1 ")) {
 532     return _m0.c1;
 533   }else if (Umatch(x, "is{tamp} ")) {
 534     return _m0.f0();
 535   }else if (Umatch(x, "iof{fset} ")) {
 536     return _m0.c0;
 537   }else if (Umatch(x, "ip{assive} ")) {itested();
 538     return _m0.c1 * tr_involts();
 539   }else if (Umatch(x, "il{oss} ")) {untested();
 540     return _loss0 * tr_outvolts();
 541   }else if (Umatch(x, "dt ")) {
 542     return _dt;
 543   }else if (Umatch(x, "dtr{equired} ")) {
 544     return ((_time_by._error_estimate - _time[0]) > 0)
 545       ? _time_by._error_estimate - _time[0]
 546       : _time_by._error_estimate - _time[1];
 547   }else if (Umatch(x, "time ")) {untested();
 548     return _time[0];
 549   }else if (Umatch(x, "timeo{ld} ")) {untested();
 550     return _time[1];
 551   //}else if (Umatch(x, "didt ")) {untested();
 552     //double i0  = (_m0.c1  * _m0.x  + _m0.c0);
 553     //double it1 = (mt1.f1 * mt1.x + mt1.c0);
 554     //return  (i0 - it1) / (_sim->_time0 - _time1);
 555   }else if (Umatch(x, "r ")) {
 556     return (_m0.c1!=0.) ? 1/_m0.c1 : MAXDBL;
 557   }else if (Umatch(x, "ic ")) {
 558     if (!has_ic()){ untested();
 559       return NOT_VALID;
 560     } else if (const EVAL_BM_ACTION_BASE* x = dynamic_cast<const EVAL_BM_ACTION_BASE*>(common())) {
 561       return x->_ic;
 562     } else { unreachable();
 563       return NOT_VALID;
 564     }
 565   }else if (Umatch(x, "z ")) {
 566     return port_impedance(_n[OUT1], _n[OUT2], _sim->_lu, (double) (mfactor()*(_m0.c1+_loss0)));
 567 #ifndef NDEBUG
 568   }else if (Umatch(x, "mp ")) {
 569     return _n[OUT1].m_();
 570   }else if (Umatch(x, "mn ")) {
 571     return _n[OUT2].m_();
 572 #endif
 573   }else if (Umatch(x, "common ")) {
 574     return has_common();
 575   }else if (Umatch(x, "treval ")) {
 576     return has_tr_eval();
 577   }else if (Umatch(x, "utreval ")) {
 578     return using_tr_eval();
 579   }else if (Umatch(x, "zraw ")) {
 580     return port_impedance(_n[OUT1], _n[OUT2], _sim->_lu, 0.);
 581   }else if (Umatch(x, "ord{er} ")) {
 582     return order();
 583   }else{
 584     return COMPONENT::tr_probe_num(x);
 585   }
 586 }
 587 /*--------------------------------------------------------------------------*/
 588 COMPLEX ELEMENT::ac_amps()const
 589 {
 590   assert(!is_source());
 591   return (ac_involts() * _acg + ac_outvolts() * (double)_loss0);
 592 }
 593 /*--------------------------------------------------------------------------*/
 594 XPROBE ELEMENT::ac_probe_ext(const std::string& x)const
 595 {
 596   COMPLEX admittance = (is_source()) ? _loss0 : _acg+ (double)_loss0;
 597
 598   if (Umatch(x, "v{out} ")) {                   /* volts (out) */
 599     return XPROBE(ac_outvolts());
 600   }else if (Umatch(x, "vin ")) {                /* volts (in) */
 601     return XPROBE(ac_involts());
 602   }else if (Umatch(x, "i ")) {                  /* amps */
 603     return XPROBE(ac_amps());
 604   }else if (Umatch(x, "p ")) {                  /* complex "power" */
 605     return XPROBE(ac_outvolts() * conj(ac_amps()), mtREAL, 10.);
 606   }else if (Umatch(x, "nv ")) {untested();      /* nominal value */
 607     return XPROBE(value());
 608   }else if (Umatch(x, "ev ")) {                 /* effective value */
 609     return XPROBE((COMPLEX)_ev);
 610   }else if (Umatch(x, "y ")) {
 611     /* admittance */
 612     return XPROBE(admittance, mtREAL);
 613   }else if (Umatch(x, "r ")) {                  /* complex "resistance" */
 614     if (admittance == 0.) {untested();
 615       return XPROBE(MAXDBL);
 616     }else{
 617       return XPROBE(1. / admittance);
 618     }
 619   }else if (Umatch(x, "z ")) {                  /* port impedance */
 620     return XPROBE(port_impedance(_n[OUT1], _n[OUT2], _sim->_acx, mfactor()*admittance));
 621   }else if (Umatch(x, "zraw ")) {               /* port impedance, raw */
 622     return XPROBE(port_impedance(_n[OUT1], _n[OUT2], _sim->_acx, COMPLEX(0.)));
 623   }else if (Umatch(x, "n ")) {                  /* port noise, raw */
 624     return XPROBE(port_noise(_n[OUT1], _n[OUT2]));
 625   }else{                                        /* bad parameter */
 626     return COMPONENT::ac_probe_ext(x);
 627   }
 628 }
 629 /*--------------------------------------------------------------------------*/
 630 double ELEMENT::tr_review_trunc_error(const FPOLY1* q)
 631 {
 632   int error_deriv;
 633   if (order() >= OPT::_keep_time_steps - 2) {
 634     error_deriv = OPT::_keep_time_steps - 1;
 635   }else if (order() < 0) {itested();
 636     error_deriv = 1;
 637   }else{
 638     error_deriv = order()+1;
 639   }
 640
 641   double timestep;
 642   trace1("ELEMENT::tr_review_trunc_error", error_deriv);
 643   trace2("ELEMENT::tr_review_trunc_error", _time[0], error_deriv);
 644   // if (_time[0] <= _sim->_time0)
 645   if (_sim->analysis_is_tran_restore()) {
 646     timestep = NEVER;
 647   }else if (_time[0] == 0.) {
 648     // DC, I know nothing
 649     timestep = NEVER;
 650   }else if (error_deriv - 1 - OPT::initsc < 0 || _time[error_deriv - 1 - OPT::initsc] <= 0 ) {
 651     // first few steps, I still know nothing
 652     // repeat whatever step was used the first time
 653     timestep = _dt;
 654   }else{
 655     for (int i=error_deriv-1; i>0; --i) {
 656       assert(_time[i] < _time[i-1]); // || _time[i] == 0.);
 657     }
 658
 659     double c[OPT::_keep_time_steps];
 660     for (int i=0; i<OPT::_keep_time_steps; ++i) {
 661       c[i] = q[i].f0;
 662     }
 663     assert(error_deriv < OPT::_keep_time_steps);
 664
 665     // better use divdiff.
 666     // better, only compute up to error_div
 667     if (error_deriv && 0. == _time[error_deriv-1]) {
 668       // vile hack, probably not a good idea
 669       _time[error_deriv] = -_time[error_deriv-2];
 670       derivatives(c, OPT::_keep_time_steps, _time);
 671       _time[error_deriv] = 0;
 672     } else {
 673       derivatives(c, OPT::_keep_time_steps, _time);
 674     }
 675     // now c[i] is i'th derivative
 676
 677     assert(OPT::_keep_time_steps >= 5);
 678     trace0(("ts" + long_label()).c_str());
 679     trace5("time", _time[0], _time[1], _time[2], _time[3], _time[4]);
 680     trace5("charge", q[0].f0, q[1].f0, q[2].f0, q[3].f0, q[4].f0);
 681     trace5("deriv", c[0], c[1], c[2], c[3], c[4]);
 682
 683     if (c[error_deriv] == 0) {
 684       timestep = NEVER;
 685     }else{
 686       double chargetol = std::max(OPT::chgtol,
 687                                 OPT::reltol * std::max((double)std::abs(q[0].f0), (double) std::abs(q[1].f0)));
 688       double tol = OPT::trtol * chargetol;
 689       double denom = error_factor() * std::abs(c[error_deriv]);
 690       assert(tol > 0.);
 691       assert(denom > 0.);
 692       switch (error_deriv) { // pow is slow.
 693       case 1:   timestep = tol / denom;         break;
 694       case 2:   timestep = sqrt(tol / denom);   break;
 695       case 3:   timestep = cbrt(tol / denom);   break;
 696       default:  timestep = pow((tol / denom), 1./(error_deriv)); break;
 697       }
 698       trace4("", chargetol, tol, denom, timestep);
 699     }
 700   }
 701   assert(timestep > 0.);
 702   return timestep;
 703 }
 704 /*--------------------------------------------------------------------------*/
 705 double ELEMENT::tr_review_check_and_convert(double timestep)
 706 {
 707   double time_future;
 708   if (timestep == NEVER) {
 709     time_future = NEVER;
 710   }else{
 711     if (timestep < _sim->_dtmin) {
 712       timestep = _sim->_dtmin;
 713     }else{
 714     }
 715
 716     if (timestep < _dt * OPT::trreject) {
 717       if (_time[order()] == 0) {
 718         error(bTRACE, "initial step rejected:" + long_label() + '\n');
 719         error(bTRACE, "new=%g  old=%g  required=%g\n",
 720               timestep, _dt, _dt * OPT::trreject);
 721       }else{
 722         error(bTRACE, "step rejected:" + long_label() + '\n');
 723         error(bTRACE, "new=%g  old=%g  required=%g\n",
 724               timestep, _dt, _dt * OPT::trreject);
 725       }
 726       time_future = _time[1] + timestep;
 727       trace3("reject", timestep, _dt, time_future);
 728     }else{
 729       time_future = _time[0] + timestep;
 730       trace3("accept", timestep, _dt, time_future);
 731     }
 732   }
 733   assert(time_future > 0.);
 734   assert(time_future > _time[1]);
 735   return time_future;
 736 }
 737 /*--------------------------------------------------------------------------*/
 738 void ELEMENT::set_ic(double x)
 739 {
 740   // tr_unload();
 741   trace1("ELEMENT::set_ic " + long_label(), x);
 742   *(_value.pointer_hack()) = x;
 743   // do_tr();
 744   // tr_load();
 745   q_eval();
 746   // q_load();
 747 }
 748 /*--------------------------------------------------------------------------*/
 749 void ELEMENT::keep_ic()
 750 {
 751   unreachable();
 752 }
 753 /*--------------------------------------------------------------------------*/
 754 void ELEMENT::tt_advance()
 755 {
 756   trace3("ELEMENT::tt_advance", short_label(), _sim->_time0, _sim->_dt0);
 757   if (_time[0] > 0.) {
 758     for (int i=0; i<OPT::_keep_time_steps-1; ++i) {
 759       _time[i] = 0.;
 760       _y[i] = _y[i+1];
 761     }
 762     _time[OPT::_keep_time_steps-1] = 0.;
 763     _y[OPT::_keep_time_steps-1]    = FPOLY1(0., 0., 0.);
 764   }else if (_time[0] == _sim->_time0) {
 765   }else{
 766   }
 767
 768   for (int i=OPT::_keep_time_steps-1; i>=0; --i) {
 769     assert(!_time[i]);
 770   }
 771 }
 772 /*--------------------------------------------------------------------------*/
 773 // hmmmm like tr_begin? start from scratch?
 774 void ELEMENT::tt_regress()
 775 {
 776   assert(_sim->_time0 == 0.);
 777   trace3("ELEMENT::tt_regress", short_label(), _sim->_time0, _sim->_dt0);
 778   if (_time[0] > _sim->_time0) {
 779     for (int i=0; i<OPT::_keep_time_steps-1; ++i) {
 780       _time[i] = 0.;
 781       _y[i] = _y[i+1];
 782     }
 783     _time[OPT::_keep_time_steps-1] = 0.;
 784     _y[OPT::_keep_time_steps-1]    = FPOLY1(0., 0., 0.);
 785   }else if (_time[0] == _sim->_time0) {
 786
 787   }else{
 788
 789   }
 790
 791   for (int i=OPT::_keep_time_steps-1; i>=0; --i) {
 792     assert(!_time[i]);
 793   }
 794 }
 795 /*--------------------------------------------------------------------------*/
 796 // probably ineffective?
 797 void ELEMENT::set_param_by_name(string Name, string Value)
 798 {
 799   trace3("ELEMENT::set_param_by_name", Name, Value, value_name());
 800   bool isvalue = false;
 801   if(value_name() == "") {
 802   }else if(Umatch(Name,value_name())) {
 803     set_value(Value);
 804     isvalue = true;
 805   } else {
 806   }
 807
 808   if (has_common()) {
 809     COMMON_COMPONENT* c = mutable_common()->clone();
 810     try{
 811       c->set_param_by_name(Name,Value);
 812     }catch(Exception_No_Match){
 813       if(!isvalue){
 814         delete c;
 815         throw;
 816       }
 817     }
 818     if(isvalue){ itested();
 819       c->set_modelname(Value);
 820     }
 821     attach_common(c);
 822   }else if (Umatch (Name,"type")) { untested();
 823     COMMON_COMPONENT* c = bm_dispatcher[Value];
 824     if(!c) throw Exception_No_Match(Value);
 825     attach_common(c->clone());
 826   } else {
 827     COMPONENT::set_param_by_name(Name,Value);
 828   }
 829 }
 830 /*--------------------------------------------------------------------------*/
 831 string ELEMENT::dev_type()const
 832 {
 833   if (common()) {
 834     return common()->modelname();
 835   }else{
 836     return element_type();
 837   }
 838 }
 839 /*--------------------------------------------------------------------------*/
 840 /*--------------------------------------------------------------------------*/
 841 /*--------------------------------------------------------------------------*/
 842 // vim:ts=8:sw=2:noet: