Fix rys roots

sunqm · Sep 17, 2023 · 9f93032 · 9f93032
1 parent 87635af
commit 9f93032
Show file tree

Hide file tree

Showing 3 changed files with 8 additions and 5 deletions.
diff --git a/scripts/rys_tabulate.py b/scripts/rys_tabulate.py
@@ -20,7 +20,7 @@ def clenshaw_points(n):
 ngrids = 14
 chebrt = np.array(chebyshev_roots(ngrids))
 cs = np.array(clenshaw_points(ngrids))
-TBASE = np.arange(0, 51, 2.5)
+TBASE = np.append(np.arange(0, 39, 2.5), np.arange(40, 104, 4))
 print(TBASE)
 
 def get_cheb_t_points(tbase):
@@ -93,14 +93,17 @@ def polyfit_erf(nroots, x):
 def pkl2table(prefix, pklfile):
     with open(pklfile, 'rb') as f:
         TBASE, rys_tab = pickle.load(f)
+    nt = find_tbase(81)
     TBASE = TBASE.round(6)
+    TBASE = TBASE[:nt+1]
     with open(f'{prefix}_x.dat', 'w') as fx, open(f'{prefix}_w.dat', 'w') as fw:
         fw.write(f'// DATA_TBASE[{len(TBASE)}] = ''{' + (', '.join([str(x) for x in TBASE])) + '};\n')
         fx.write(f'static double DATA_X[] = ''{\n')
         fw.write(f'static double DATA_W[] = ''{\n')
         for i, tab in enumerate(rys_tab):
             nroots = i + 6
-            for it, ttab in enumerate(tab):
+            for it in range(nt):
+                ttab = tab[it]
                 tbase = TBASE[it]
                 print(f'root {nroots}  tbase[{it}] {tbase}')
                 fx.write(f'/* root={nroots} base[{it}]={tbase} */\n')

diff --git a/scripts/sr_rys_tabulate_error.py b/scripts/sr_rys_tabulate_error.py
@@ -59,15 +59,15 @@
                 print(nroots, x, low, diff1, diff2)
 
 np.random.seed(1)
-xs = np.sort(np.append(np.random.rand(100) * 50, 3**((np.random.rand(20)-.3) * 5)))
+xs = np.sort(np.append(np.random.rand(100) * 100, 3**((np.random.rand(30)-.4) * 6)))
 xs[0] = 1e-15
 fil = 'rys_rw.pkl'
 TBASE, tab = pickle.load(open(fil, 'rb'))
 rys_tabulate.TBASE = TBASE
 
 for nroots in range(1, 15):
     for x in xs:
-        if x >= 50:
+        if x >= 100:
             continue
 
         it = np.searchsorted(rys_tabulate.TBASE, x) - 1

diff --git a/src/rys_roots.c b/src/rys_roots.c
@@ -64,7 +64,7 @@ void CINTrys_roots(int nroots, double x, double *u, double *w)
                         w[i] = POLY_SMALLX_W0[off+i] + POLY_SMALLX_W1[off+i] * x;
                 }
                 return;
-        } else if (x >= 50.) {
+        } else if (x >= 35+nroots*5) {
                 int off = nroots * (nroots - 1) / 2;
                 int i;
                 double rt;