implement option 'xmu' to control trapezoidal integration method

default is 0.5

option xmu=0.49
  can be used to add some damping to reduce trap ringing

set xmu=0.49
  or
.option xmu=0.49
  works as well

src/frontend/com_option.c

@@ -46,6 +46,7 @@ com_option(wordlist *wl)
             printf("Unknown integration method\n");
         }
         printf("MaxOrder = %d\n", circuit->CKTmaxOrder);
+        printf("xmu = %g\n", circuit->CKTxmu);
 
         printf("\nTolerances (absolute):\n");
         printf("abstol      (current) = %g\n", circuit->CKTabstol);

src/frontend/spiceif.c

@@ -1409,6 +1409,7 @@ void com_snload(wordlist *wl)
     _t(CKTorder);
     _t(CKTmaxOrder);
     _t(CKTintegrateMethod);
+    _t(CKTxmu);
 
     _t(CKTniState);
 

src/include/ngspice/cktdefs.h

@@ -101,6 +101,7 @@ struct CKTcircuit {
     int CKTorder;               /* the integration method order */
     int CKTmaxOrder;            /* maximum integration method order */
     int CKTintegrateMethod;     /* the integration method to be used */
+    double CKTxmu;              /* for trapezoidal method */
 
 /* known integration methods */
 #define TRAPEZOIDAL 1

src/include/ngspice/optdefs.h

@@ -119,6 +119,8 @@ typedef struct {
 #define OPT_RELDV        67  /* Original: 52 (Node_Damping) */
 
 #define OPT_NOOPAC       68
+#define OPT_XMU          69
+
 
 #ifdef XSPICE
 /* gtri - begin - wbk - add new options */

src/include/ngspice/tskdefs.h

@@ -20,6 +20,7 @@ struct TSKtask {
     double TSKnomTemp;
     int TSKmaxOrder;        /* maximum integration method order */
     int TSKintegrateMethod; /* the integration method to be used */
+    double TSKxmu;          /* for trapezoidal method */
     int TSKcurrentAnalysis; /* the analysis in progress (if any) */
 
 /* defines for the value of  TSKcurrentAnalysis */

src/maths/ni/nicomcof.c

@@ -6,12 +6,11 @@ Author: 1985 Thomas L. Quarles
 #include "ngspice/ngspice.h"
 #include "ngspice/cktdefs.h"
 #include "ngspice/sperror.h"
+#include "ngspice/cpextern.h"
 
 /* xmu=0:    Backward Euler
  * xmu=0.5:  trapezoidal (standard)
- * xmu=0.48: good damping of current ringing, e.g. in R.O.s.
  */
-#define xmu 0.5
 
 
 int
@@ -41,8 +40,8 @@ NIcomCof(CKTcircuit *ckt)
             break;
 
         case 2:
-            ckt->CKTag[0] = 1.0 / ckt->CKTdelta / (1.0 - xmu);
-            ckt->CKTag[1] = xmu / (1.0 - xmu);
+            ckt->CKTag[0] = 1.0 / ckt->CKTdelta / (1.0 - ckt->CKTxmu);
+            ckt->CKTag[1] = ckt->CKTxmu / (1.0 - ckt->CKTxmu);
             break;
 
         default:

src/spicelib/analysis/cktdojob.c

@@ -52,6 +52,7 @@ CKTdoJob(CKTcircuit *ckt, int reset, TSKtask *task)
     ckt->CKTnomTemp  = task->TSKnomTemp;
     ckt->CKTmaxOrder  = task->TSKmaxOrder;
     ckt->CKTintegrateMethod  = task->TSKintegrateMethod;
+    ckt->CKTxmu  = task->TSKxmu;
     ckt->CKTbypass  = task->TSKbypass;
     ckt->CKTdcMaxIter  = task->TSKdcMaxIter;
     ckt->CKTdcTrcvMaxIter  = task->TSKdcTrcvMaxIter;

src/spicelib/analysis/cktntask.c

@@ -39,6 +39,7 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr)
         tsk->TSKnomTemp         = def->TSKnomTemp;
         tsk->TSKmaxOrder        = def->TSKmaxOrder;
         tsk->TSKintegrateMethod = def->TSKintegrateMethod;
+        tsk->TSKxmu             = def->TSKxmu;
         tsk->TSKbypass          = def->TSKbypass;
         tsk->TSKdcMaxIter       = def->TSKdcMaxIter;
         tsk->TSKdcTrcvMaxIter   = def->TSKdcTrcvMaxIter;
@@ -98,6 +99,13 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr)
         tsk->TSKdcTrcvMaxIter   = 50;
         tsk->TSKintegrateMethod = TRAPEZOIDAL;
         tsk->TSKmaxOrder        = 2;
+        /*
+         * when using trapezoidal method
+         *   xmu=0:    Backward Euler
+         *   xmu=0.5:  trapezoidal (standard)
+         *   xmu=0.49: good damping of current ringing, e.g. in R.O.s.
+         */
+        tsk->TSKxmu             = 0.5;
         tsk->TSKnumSrcSteps     = 1;
         tsk->TSKnumGminSteps    = 1;
         tsk->TSKgminFactor      = 10;

src/spicelib/analysis/cktsopt.c

@@ -114,6 +114,9 @@ CKTsetOpt(CKTcircuit *ckt, JOB *anal, int opt, IFvalue *val)
     case OPT_BYPASS:
         task->TSKbypass = val->iValue;
         break;
+    case OPT_XMU:
+        task->TSKxmu = val->rValue;
+        break;
     case OPT_MAXORD:
         task->TSKmaxOrder = val->iValue;
         /* Check options method and maxorder for consistency */
@@ -265,6 +268,7 @@ static IFparm OPTtbl[] = {
  { "lvltim", 0, IF_INTEGER,"Type of timestep control" },
  { "method", OPT_METHOD, IF_SET|IF_STRING,"Integration method" },
  { "maxord", OPT_MAXORD, IF_SET|IF_INTEGER,"Maximum integration order" },
+ { "xmu", OPT_XMU, IF_SET|IF_REAL,"Coefficient for trapezoidal method" },
  { "defm", OPT_DEFM,IF_SET|IF_REAL,"Default MOSfet Multiplier" },
  { "defl", OPT_DEFL,IF_SET|IF_REAL,"Default MOSfet length" },
  { "defw", OPT_DEFW,IF_SET|IF_REAL,"Default MOSfet width" },

src/spicelib/devices/cktinit.c

@@ -65,6 +65,7 @@ CKTinit(CKTcircuit **ckt)		/* new circuit to create */
     sckt->CKTintegrateMethod = TRAPEZOIDAL;
     sckt->CKTorder = 1;
     sckt->CKTmaxOrder = 2;
+    sckt->CKTxmu = 0.5;
     sckt->CKTpivotAbsTol = 1e-13;
     sckt->CKTpivotRelTol = 1e-3;
     sckt->CKTtemp = 300.15;