From 6c1be283a757d829f1a3799a73893d0ad9bb58e3 Mon Sep 17 00:00:00 2001
From: Holger Vogt <holger.vogt@uni-due.de>
Date: Tue, 29 Nov 2022 15:36:31 +0100
Subject: [PATCH] Add a function ddt (derivative versus time) to the B-source
 funtion parser.

---
 src/include/ngspice/inpptree.h |  1 +
 src/spicelib/parser/inpptree.c | 31 ++++++++++++++++++
 src/spicelib/parser/inpxx.h    |  1 +
 src/spicelib/parser/ptfuncs.c  | 57 ++++++++++++++++++++++++++++++++--
 4 files changed, 87 insertions(+), 3 deletions(-)

diff --git a/src/include/ngspice/inpptree.h b/src/include/ngspice/inpptree.h
index 81686a492..7fc00cae6 100644
--- a/src/include/ngspice/inpptree.h
+++ b/src/include/ngspice/inpptree.h
@@ -129,6 +129,7 @@ void INPptPrint(char *str, IFparseTree * ptree);
 #define PTF_CEIL	34
 #define PTF_FLOOR	35
 #define PTF_NINT	36
+#define PTF_DDT	37
 
 /* The following things are used by the parser -- these are the token types the
  * lexer returns.
diff --git a/src/spicelib/parser/inpptree.c b/src/spicelib/parser/inpptree.c
index 4d725077b..1f97cf082 100644
--- a/src/spicelib/parser/inpptree.c
+++ b/src/spicelib/parser/inpptree.c
@@ -171,6 +171,7 @@ static struct func {
     { "pwr",    PTF_PWR,    (void(*)(void)) PTpwr},
     { "min",    PTF_MIN,    (void(*)(void)) PTmin},
     { "max",    PTF_MAX,    (void(*)(void)) PTmax},
+    { "ddt",    PTF_DDT,    (void(*)(void)) PTddt},
 } ;
 
 #define NUM_FUNCS (int)NUMELEMS(funcs)
@@ -566,6 +567,11 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum)
             arg1 = mkcon(0.0);
             break;
 
+        case PTF_DDT:
+            arg1 = mkcon(0.0);
+            arg1->data = p->data;
+            break;
+
         case PTF_MIN:
         case PTF_MAX:
         /* min(a,b) -->   (a<b)       ? a : b
@@ -1080,6 +1086,20 @@ static INPparseNode *prepare_PTF_PWL(INPparseNode *p)
     return (p);
 }
 
+static INPparseNode* prepare_PTF_DDT(INPparseNode* p)
+{
+    struct ddtdata { int n; double* vals; } *data;
+    int i, ii;
+    /* store 3 recent times and 3 recent values in pairs t0, v0, t1, v1, t2, v2  */
+    i = 0;
+    data = TMALLOC(struct ddtdata, 1);
+    data->vals = TMALLOC(double, 7);
+    for (ii = 0; ii < 7; ii++) {
+        data->vals[ii] = 0;
+    }
+    p->data = (void*)data;
+    return (p);
+}
 
 INPparseNode *PT_mkfnode(const char *fname, INPparseNode * arg)
 {
@@ -1140,6 +1160,9 @@ INPparseNode *PT_mkfnode(const char *fname, INPparseNode * arg)
     if(p->funcnum == PTF_PWL)
         p = prepare_PTF_PWL(p);
 
+    if (p->funcnum == PTF_DDT)
+        p = prepare_PTF_DDT(p);
+
     return (p);
 }
 
@@ -1573,6 +1596,14 @@ void free_tree(INPparseNode *pt)
         }
     }
 
+    if (pt->type == PT_FUNCTION && (pt->funcnum == PTF_DDT)) {
+        struct ddtdata { int n; double* vals; } *data = (struct ddtdata*)(pt->data);
+        if (data) {
+            txfree(data->vals);
+            txfree(data);
+        }
+    }
+
     txfree(pt);
 }
 
diff --git a/src/spicelib/parser/inpxx.h b/src/spicelib/parser/inpxx.h
index 77e4a2143..81e9355db 100644
--- a/src/spicelib/parser/inpxx.h
+++ b/src/spicelib/parser/inpxx.h
@@ -78,5 +78,6 @@ double PTle0(double arg);
 double PTceil(double arg);
 double PTfloor(double arg);
 double PTnint(double arg);
+double PTddt(double arg, void* data);
 
 #endif
diff --git a/src/spicelib/parser/ptfuncs.c b/src/spicelib/parser/ptfuncs.c
index e0211e74c..2ea1d8ba4 100644
--- a/src/spicelib/parser/ptfuncs.c
+++ b/src/spicelib/parser/ptfuncs.c
@@ -13,12 +13,10 @@ Author: 1987 Wayne A. Christopher, U. C. Berkeley CAD Group
 #include "ngspice/fteext.h"
 #include "ngspice/ifsim.h"
 #include "ngspice/inpptree.h"
+#include "ngspice/cktdefs.h"
 #include "inpxx.h"
 #include "ngspice/compatmode.h"
 
-/* XXX These should be in math.h */
-
-
 double PTfudge_factor;
 
 #define MODULUS(NUM,LIMIT) ((NUM) - ((int) ((NUM) / (LIMIT))) * (LIMIT))
@@ -414,3 +412,56 @@ PTnint(double arg1)
      */
     return nearbyint(arg1);
 }
+
+double
+PTddt(double arg, void* data)
+{
+    struct ddtdata { int n; double* vals; } *thing = (struct ddtdata*)data;
+    double y, time;
+
+    double debarr[7];
+    int i;
+
+    CKTcircuit* ckt = ft_curckt->ci_ckt;
+
+    time = ckt->CKTtime;
+
+    if (time == 0) {
+        thing->vals[3] = arg;
+        return 0;
+    }
+
+    if (!(ckt->CKTmode & MODETRAN) /*&& !(ckt->CKTmode & MODEINITTRAN)*/)
+        return 0;
+
+
+    if (time > thing->vals[0]) {
+        thing->vals[4] = thing->vals[2];
+        thing->vals[5] = thing->vals[3];
+        thing->vals[2] = thing->vals[0];
+        thing->vals[3] = thing->vals[1];
+        thing->vals[0] = time;
+        thing->vals[1] = arg;
+
+        for (i = 0; i < 7; i++)
+            debarr[i] = thing->vals[i];
+
+        if (thing->n > 1) {
+            //       if (thing->vals[2] > 0) {
+           //            debarr[6] = thing->vals[6] = 0.5 * ((arg - thing->vals[3]) / (time - thing->vals[2]) + thing->vals[6]);
+            debarr[6] = thing->vals[6] = (thing->vals[1] - thing->vals[3]) / (thing->vals[2] - thing->vals[4]);
+        }
+        else {
+            debarr[6] = thing->vals[6] = 0;
+            debarr[3] = thing->vals[3] = arg;
+        }
+        thing->n += 1;
+    }
+
+    for (i = 0; i < 7; i++)
+        debarr[i] = thing->vals[i];
+
+    y = thing->vals[6];
+
+    return y;
+}