diff --git a/scripts/rys_tabulate.py b/scripts/rys_tabulate.py index 88432687..00ee7bd8 100644 --- 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 index 698de9a5..b9cd3e80 100644 --- a/scripts/sr_rys_tabulate_error.py +++ b/scripts/sr_rys_tabulate_error.py @@ -59,7 +59,7 @@ 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')) @@ -67,7 +67,7 @@ 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 index ebea454c..7ca833e3 100644 --- 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;