import sys import numpy as np from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, sinh, cosh, tanh, sqrt, abs from scipy import signal from random import random as rand import pandas as pd import matplotlib.pyplot as plt from tklib.tkutils import replace_path, getarg, getintarg, getfloatarg, pint, pfloat from tklib.tksci.tksci import kB, e from tklib.tkvariousdata import tkVariousData from tklib.tkapplication import tkApplication from tklib.tkparams import tkParams from tklib.tksci.tkintegration import gauss_legendre, simpson from tklib.tksci.tksci import asin, acos, atan, degcos, degsin, degtan, degacos, degasin, degatan, arcsin, arccos, arctan from tklib.tksci.tksci import factorial, log10, gamma, combination, eVTonm, nmToeV, Bn from tklib.tksci.tksci import h, h_bar, hbar, e, kB, NA, c, pi, pi2, torad, todeg, basee from tklib.tksci.tksci import me, mp, mn from tklib.tksci.tksci import u0, e0, a0, R, F, g, G from tklib.tksci.tksci import HartreeToeV, RyToeV, KToeV, eVToK, JToeV, eVToJ, Debye app = tkApplication() cparams = tkParams() cparams.func_str = 'exp(-x) + exp(-2.0 * x)' cparams.mode = 'forward' cparams.xmin = 0.0 cparams.xmax = 10.0 cparams.nx = 100 cparams.smin = 0.0 cparams.smax = 10.0 cparams.ns = 100 cparams.sigma = 0.3 cparams.nsmooth = 7 cparams.nsmoothorder = 5 cparams.figsize = [12, 8] cparams.fontsize = 12 cparams.legend_fontsize = 10 logfile = 'fft_func-out.txt' #app.replace_path(cparams.infile) print(f"Open logfile [{logfile}]") app.redirect(targets = ["stdout", logfile], mode = 'w') def usage(app): print(f"Usage: python {sys.argv[0]} func_str x0 x1 nx") def update_vars(app): cparams.mode = getarg(1, cparams.mode) cparams.func_str = getarg(2, cparams.func_str) cparams.xmin = getfloatarg(3, cparams.xmin) cparams.xmax = getfloatarg(4, cparams.xmax) cparams.nx = getintarg (5, cparams.nx) cparams.smin = getfloatarg(6, cparams.smin) cparams.smax = getfloatarg(7, cparams.smax) cparams.ns = getintarg (8, cparams.ns) cparams.sigma = getfloatarg(9, cparams.sigma) cparams.nsmooth = getintarg(10, cparams.nsmooth) cparams.nsmoothorder = getintarg(11, cparams.nsmoothorder) def func(x): y = eval(cparams.func_str, globals(), {"x": x}) return y def integ_func(x, s): f = exp(complex(0.0, cparams.sigma) * x) * exp(-s * x) * func(x) return f def common_func(): print("") print("========================================") print(f"mode: {cparams.mode}") print(f"func: {cparams.func_str}") print(f"smearing sigma: {cparams.sigma}") print(f"x range : {cparams.xmin} - {cparams.xmax} for {cparams.nx} points") print(f"s range : {cparams.smin} - {cparams.smax} for {cparams.ns} points") print(f"Smoothing:") print(f" Number of points: {cparams.nsmooth}") print("a=", cparams.nsmooth % 2) if cparams.nsmooth % 2 == 0: cparams.nsmooth += 1 print(f" Number of points must be odd. Changed to {cparams.nsmooth}") print(f" Order of polynomial: {cparams.nsmoothorder}") if '***' in cparams.func_str: app.terminate("\nError: choose function\n", usage = usage, pause = True) print("") print(f"Build function [{cparams.func_str}]") cparams.xstep = (cparams.xmax - cparams.xmin) / (cparams.nx - 1) x = [cparams.xmin + i * cparams.xstep for i in range(cparams.nx)] y = [func(_x) for _x in x] return x, y def forward(): x, y = common_func() xmin = min(x) xmax = max(x) sstep = (cparams.smax - cparams.smin) / (cparams.ns - 1) s = [cparams.smin + i * sstep for i in range(cparams.ns)] Fr = [] Fi = [] Fabs = [] invF = [] for i in range(cparams.ns): S = simpson(lambda x1: integ_func(x1, s[i]), xmin, xmax, cparams.nx) Fr.append(S.real) Fi.append(S.imag) Sabs = abs(S) Fabs.append(Sabs) invF.append(1.0 / Sabs) print("") print(f"{'s':10} {'y(LT)(real)':10} {'y(LT)(imag)':10} {'|y(LT)|':10}") for i in range(cparams.ns): print(f"{s[i]:10.4g} {Fr[i]:10.4g} {Fi[i]:10.4g} {Fabs[i]:10.4g}") invsmooth = signal.savgol_filter(invF, cparams.nsmooth, cparams.nsmoothorder, deriv = 0) ds = s[1] - s[0] invdiff = [(invsmooth[1] - invsmooth[0]) / ds] for i in range(1, cparams.ns - 1): # d = (invsmooth[i] - invsmooth[i - 1]) / ds d = (invsmooth[i + 1] - invsmooth[i - 1]) / 2.0 / ds invdiff.append(d) invdiff.append((invsmooth[cparams.ns - 1] - invsmooth[cparams.ns - 2]) / ds) print("") print("First differential") for i in range(cparams.ns): print(f"{s[i]:10.4g} {invdiff[i]:10.4g}") xintercept = [] for i in range(cparams.ns): if abs(invdiff[i]) > 1.0e-5: # y = y0 + (s - s0) * slope = 0: y = 1/L yi = invsmooth[i] xi = s[i] - yi / invdiff[i] else: xi = None xintercept.append(xi) print("") print("s intercept") for i in range(cparams.ns): print(f"{s[i]:10.4g} {xintercept[i]}") outfile = 'fft_func-Laplace.xlsx' #app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-Laplace.xlsx"]) df = pd.DataFrame(np.array([x, y]).T, columns = ['x', 'y']) print("") print(f"Save to [{outfile}]") df.to_excel(outfile, index = False) outfile = 'fft_func-LaplaceTransformed.xlsx' #app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-LaplaceTransformed.xlsx"]) df = pd.DataFrame(np.array([s, Fr, Fi, Fabs]).T, columns = ['s', 'y(LT)(real)', 'y(LT)(imag)', '|y(LT)|']) print("") print(f"Save to [{outfile}]") df.to_excel(outfile, index = False) print("") print("Plot") fig, axes = plt.subplots(2, 2, figsize = cparams.figsize) axes = axes.flatten() ax2 = axes[2].twinx() axes[0].tick_params(labelsize = cparams.fontsize) axes[1].tick_params(labelsize = cparams.fontsize) axes[2].tick_params(labelsize = cparams.fontsize) axes[3].tick_params(labelsize = cparams.fontsize) axes[0].plot(x, y, label = 'input data', color = 'black', linewidth = 0.5) axes[0].set_xlabel('x', fontsize = cparams.fontsize) axes[0].set_ylabel(f'y = {cparams.func_str}', fontsize = cparams.fontsize) axes[0].legend(fontsize = cparams.legend_fontsize) axes[1].plot(s, Fr, label = 'real', color = 'black', linewidth = 0.5) axes[1].plot(s, Fi, label = 'imag', color = 'red', linewidth = 0.5) axes[1].plot(s, Fabs, label = 'abs', color = 'blue', linewidth = 0.5) axes[1].plot(axes[1].get_xlim(), [0.0, 0.0], linestyle = 'dashed', color = 'red', linewidth = 0.5) axes[1].set_xlabel('s', fontsize = cparams.fontsize) axes[1].set_ylabel('Laplace transformation', fontsize = cparams.fontsize) axes[1].legend(fontsize = cparams.legend_fontsize) axes[1].set_ylim([-100.0, 100.0]) ins1 = axes[2].plot(s, invF, label = '1/LT', linestyle = '', marker = 'o', markersize = 0.5, markerfacecolor = 'black', markeredgecolor = 'black') ins2 = axes[2].plot(s, invsmooth, label = '1/LT(smoothed)', color = 'blue', linewidth = 0.5) axes[2].plot(axes[2].get_xlim(), [0.0, 0.0], linestyle = 'dashed', color = 'red', linewidth = 0.5) ins3 = ax2.plot(s, invdiff, label = '$d(1/LT)/ds$', color = 'red', linewidth = 0.5) axes[2].set_xlabel('s', fontsize = cparams.fontsize) axes[2].set_ylabel('Inverse of Laplace transformation', fontsize = cparams.fontsize) ax2.set_ylabel('$d(1/LT)/ds$', fontsize = cparams.fontsize) axes[2].legend(fontsize = cparams.legend_fontsize) ins = ins1 + ins2 + ins3 axes[2].legend(ins, [l.get_label() for l in ins], fontsize = cparams.legend_fontsize) axes[3].plot(s, xintercept, label = 'real', color = 'black', linewidth = 0.5) axes[3].set_xlabel('s', fontsize = cparams.fontsize) axes[3].set_ylabel('intercept in 1/LT', fontsize = cparams.fontsize) plt.tight_layout() plt.pause(0.001) app.terminate("\n", usage = usage, pause = True) def main(): print("") print( "#######################################") print(f"# {sys.argv[0]}") print( "#######################################") update_vars(app) # cparams.print_parameters() if cparams.mode == 'forward': forward() else: app.terminate(f"\nmode [{cparams.mode}] is not implemented.\n", usage = usage, pause = True) if __name__ == "__main__": main()