diff --git a/stocproc/helper.py b/stocproc/helper.py index fbade40..8d0e7d7 100644 --- a/stocproc/helper.py +++ b/stocproc/helper.py @@ -1,4 +1,6 @@ -import numpy as np +import numpy as np +from scipy.optimize import brentq +from . import class_stocproc def get_param_single_lorentz(tmax, dw_max, eta, gamma, wc, x=1e-4, verbose=0): d = gamma * np.sqrt(1/x - 1) @@ -33,3 +35,80 @@ def get_param_single_lorentz(tmax, dw_max, eta, gamma, wc, x=1e-4, verbose=0): print('-> tmax: {:.3}'.format(tmax_)) assert tmax_ > tmax, "tmax_={} > tmax={} FAILED".format(tmax_, tmax) return N, w_min, tmax + +def get_param_ohmic(t_max, spec_dens, x=1e-12, verbose=0): + fmin_spec_dens = lambda w: abs(spec_dens(w)) - spec_dens.maximum_val()*x + w_pos = spec_dens.maximum_at() + w_neg = 2*w_pos + while fmin_spec_dens(w_neg) > 0: + w_neg *= 2 + + omega_max = brentq(fmin_spec_dens, w_pos, w_neg) + if verbose > 0: + print("omega_max from threshold condition: J(w_max) = x = {:.3g} <-> w_max = {:.3g}".format(x, omega_max)) + + dw = np.pi / t_max + + if verbose > 0: + print("dw:{:.3}".format(dw)) + + ng_omega = np.ceil(omega_max / dw) # next larger integer + ng_omega = np.ceil(ng_omega / 2) * 2 - 1 # next lager odd integer + ng_t = (ng_omega + 1) / 2 # so this becomes an integer + delta_t = 2 * np.pi / (dw * ng_omega) + sp_t_max = ng_t * delta_t + if verbose > 0: + print("result ng_t: {}".format(ng_t)) + print("result sp_t_max: {:.3g}".format(sp_t_max)) + + return ng_t, sp_t_max + +def show_ohmic_sp(ng_t, sp_t_max, spec_dens, seed, ax, t_max): + try: + n = len(ng_t) + except: + n = None + + try: + m = len(sp_t_max) + except: + m = None + + if (n is None) and m is not None: + ng_t = [ng_t] * m + elif (n is not None) and m is None: + sp_t_max = [sp_t_max] * n + elif (n is not None) and (m is not None): + if n != m: + raise ValueError("len(ng_t) == len(sp_t_max) FAILED") + else: + ng_t = [ng_t] + sp_t_max = [sp_t_max] + + + for i in range(len(ng_t)): + spfft = class_stocproc.StocProc_FFT(spectral_density = spec_dens, + t_max = sp_t_max[i], + num_grid_points = ng_t[i], + seed = seed) + + spfft.new_process() + + + + t = np.linspace(0, sp_t_max[i], ng_t[i]) + t_interp = np.linspace(0, sp_t_max[i], 10*ng_t[i]) + + t = t[np.where(t < t_max)] + t_interp = t_interp[np.where(t_interp < t_max)] + + eta = spfft(t) + eta_interp = spfft(t_interp) + + p, = ax.plot(t_interp, np.real(eta_interp)) + ax.plot(t_interp, np.imag(eta_interp), color=p.get_color(), ls='--') + ax.plot(t, np.real(eta), color=p.get_color(), ls='', marker = '.') + ax.plot(t, np.imag(eta), color=p.get_color(), ls='', marker = '.') + + + \ No newline at end of file