star-wf

commit 66220095b441c72c7eb24d5792bbe03018c2fb11
parent 7479559cea358a9f29145a31b078a2ffec8a5ef3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 12 Mar 2024 09:19:35 +0100

Delete the calculation of the 3rd derivative of H

It was not used anywhere. Only the first and second derivatives are used
for quadratic extrapolation.

Diffstat:
M
src/swf_H.c
 | 
29
+++++++++--------------------
M
src/swf_tabulation.h
 | 
3
+--

2 files changed, 10 insertions(+), 22 deletions(-)
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -48,16 +48,14 @@ check_swf_H_tabulate_args(const struct swf_H_tabulate_args* args)
   return RES_OK;
 }
 
-/* H(x)    = 1 + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x)]
- * H'(x)   =   - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^2]
- * H''(x)  =   + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^4]
- * H'''(x) =   - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^6] */
+/* H(x)  = 1 + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x)]
+ * H'(x) =   - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^2]
+ * H"(x) =   + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^4] */
 static INLINE double
 H3d
   (const double x,
    double* dHx, /* Derivative. NULL <=> do not compute it */
-   double* d2Hx, /* Second derivative. NULL <=> do not compute it */
-   double* d3Hx) /* Third derivative. NULL <=> do not compute it */
+   double* d2Hx) /* Second derivative. NULL <=> do not compute it */
 {
   double Hx = 0; /* Value of the function */
 
@@ -65,7 +63,6 @@ H3d
   double sum = 0; /* Sum of the terms */
   double sumd = 0; /* Sum of the derivative terms */
   double sumd2 = 0; /* Sum of the second derivative terms */
-  double sumd3 = 0; /* Sum of the third derivative terms */
 
   double sign = -1; /* Sign of the serie term */
   unsigned k = 1; /* Index of the serie term */
@@ -74,7 +71,6 @@ H3d
   int error = 1;
   int errord = dHx != NULL;
   int errord2 = d2Hx != NULL;
-  int errord3 = d3Hx != NULL;
 
   ASSERT(x >= 0); /* Check pre-condition */
 
@@ -106,19 +102,12 @@ H3d
       sumd2 = sumd2_next;
     }
 
-    if(d3Hx) { /* Third derivative */
-      const double termd3 = term * t2*t2*t2;
-      const double sumd3_next = sumd3 + termd3;
-      errord3 = sumd3_next != sumd3;
-      sumd3 = sumd3_next;
-    }
-  } while((error || errord || errord2 || errord3) && ++k < UINT_MAX);
+  } while((error || errord || errord2) && ++k < UINT_MAX);
 
   Hx = 1.0 + 2.0 * sum;
 
   if(dHx)  *dHx  = -2.0 * sumd;  /* Derivative */
   if(d2Hx) *d2Hx =  2.0 * sumd2; /* Second derivative */
-  if(d3Hx) *d3Hx = -2.0 * sumd3; /* Third derivative */
   return Hx;
 }
 
@@ -128,7 +117,7 @@ H3d
 double
 swf_H3d_eval(const double x)
 {
-  return H3d(x, NULL, NULL, NULL);
+  return H3d(x, NULL, NULL);
 }
 
 double
@@ -171,13 +160,13 @@ swf_H3d_tabulate
 
   /* Setup the tabulation */
   i = 0;
-  for(x = args->x_min; x < args->x_max; x += x*args->step) {
+  for(x = args->x_min; x <= args->x_max; x += x*args->step) {
     const struct item* pitem = darray_item_cdata_get(&tab->items) + i - 1;
     struct item item = ITEM_NULL;
 
-    /* Evaluate H(x) and its derivatives up to order 3 */
+    /* Evaluate H(x) and its derivatives up to order 2 */
     item.x = MMIN(x, args->x_max);
-    item.fx = H3d(item.x, &item.dfx, &item.d2fx, &item.d3fx);
+    item.fx = H3d(item.x, &item.dfx, &item.d2fx);
 
     if(i > 0) {
      /* Detect oscillations due to numerical problems.
diff --git a/src/swf_tabulation.h b/src/swf_tabulation.h
@@ -24,9 +24,8 @@ struct item {
   double fx; /* Value of the function */
   double dfx; /* Value of the derivative of the function */
   double d2fx; /* Value of the second derivative of the function */
-  double d3fx; /* Value of the third derivative of the function */
 };
-#define ITEM_NULL__ {0, 0, 0, 0, 0}
+#define ITEM_NULL__ {0, 0, 0, 0}
 static const struct item ITEM_NULL = ITEM_NULL__;
 
 /* Declare the dynamic array of items */