1 files changed, 136 insertions, 0 deletions
diff --git a/script/gen_phi0 b/script/gen_phi0
new file mode 100755
index 0000000..af01a31
--- /dev/null
+++ b/script/gen_phi0
@@ -0,0 +1,136 @@
+#!/usr/bin/python3
+
+import math
+
+# Save the points in separate lists 
+points_t1 = []  # A list of tuples (x, r)
+points_t2 = []  # A list of tuples (x, r)
+points_t3 = []  # A list of tuples (x, r)
+
+######################################################
+# Generate
+######################################################
+
+adj = 0.50  # Offset ratio for step of 1, estimate
+
+def loop_even(upper, lower, step, ary):
+  result = []
+  for i in range((int)(upper * step), (int)(lower * step), -1):
+    x = i / step
+    z = math.exp(x + (adj / step))
+    r = math.log((z+1)/(z-1))
+    ary.append( (x, r) )
+
+#####
+# Create initial list which will be sorted desending
+
+# Tbl-1 9:5 in linear steps of 0.5
+loop_even(9.999, 4.999, 2, points_t1)
+
+# Tbl-2 5:1 in linear steps of 1/16
+loop_even(4.999, 0.999, 16, points_t2)
+
+# Tbl-3 log steps:  1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
+for i in range(1, 28):
+    e = (2 ** (i/2))
+    x = 1.0 / e
+    z = math.exp(x+ (.19 / e))
+    r = math.log((z+1)/(z-1))
+    points_t3.append( (x, r) )
+
+
+######################################################
+# Output
+######################################################
+
+print("""
+// phi0.c
+//
+// An approximation of the function
+//
+// This file is generated by the gen_phi0 scripts
+// Any changes should be made to that file, not this one
+
+#include <stdint.h>
+
+#define SI16(f) ((int32_t)(f * (1<<16)))
+
+float phi0( float xf ) {
+
+    int32_t x = SI16(xf);
+
+    if (x >= SI16(10.0f)) return(0.0f);
+""", end="")
+
+
+##################################
+# Tbl-1 9:5 in linear steps of 0.5  (9.5, 9.0, 8.5, 8.0, .., 5.0)
+print("""    else {
+      if (x >= SI16(5.0f)) {
+        int i = 19 - (x >> 15);
+        switch (i) {
+""", end="")
+for i in range(len(points_t1)):
+    #assert(points_t1[i][0] == ((18 - i) / 2))
+    print("{}case {}: return({:.9f}f); // ({:.1f})".format(
+        (10*" "), i, points_t1[i][1], points_t1[i][0]))
+print("{}}}".format(( 8*" ")))
+print("{}}}".format(( 6*" ")))
+
+
+##################################
+# Tbl-2 5:1 in linear steps of 1/16  (4-15/16, 4-7/7, ... 1)
+print("""      else {
+        if (x >= SI16(1.0f)) {
+          int i = 79 - (x >> 12);
+          switch (i) {
+""", end="")
+for i in range(len(points_t2)):
+    #assert(points_t2[i][0] == ((18 - i) / 2))
+    print("{}case {}: return({:.9f}f); // ({:.4f})".format(
+        (12*" "), i, points_t2[i][1], points_t2[i][0]))
+print("{}}}".format((10*" ")))
+print("{}}}".format(( 8*" ")))
+
+
+##################################
+# Tbl-3 log steps:  1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
+# Output as a balanced search
+
+def prnt_cmp(x, ind):
+    print("{}if (x > SI16({:.6f}f)) {{".format((" "*ind), x))
+
+def prnt_rtn(r, ind):
+    print("{}return({:.9f}f);".format((" "*ind), r))
+
+def one_level(pts, ind, dft):
+    #print("# One_Level({})".format(pts))
+    mid = (int)(len(pts)/2)
+    lft = pts[:mid]
+    rgt = pts[mid+1:]
+    x = pts[mid][0]
+    r = pts[mid][1]
+    #print("## {}, {}, {}".format(lft, x, rgt))
+    prnt_cmp(x, ind)
+    if (len(lft)): 
+        one_level(lft, ind+2, r)
+    else:
+        prnt_rtn(r, (ind+2))
+    print("{}}} else {{".format((" "*(ind))))
+    if (len(rgt)): 
+        one_level(rgt, ind+2, dft)
+    else:
+        prnt_rtn(dft, (ind+2))
+    print("{}}}".format((" "*(ind))))
+    
+
+# Start recursive process
+print("{}else {{".format((8*" ")))
+indent = 10
+one_level(points_t3, indent, 10)
+print("{}}}".format((8*" ")))  # End of tb1_1
+print("{}}}".format((6*" ")))  # End of tb1_2
+print("{}}}".format((4*" ")))  # End of tb1_3
+print("{}return(10.0f);".format((4*" ")))
+print("}")
+