import argparse import configparser import logging import pathlib import matplotlib.pyplot as plt import numpy as np import scipy.signal as sgl from .reflex3s import reflex3s from .pa.reflex3S import reflex3S from .pa.PUV import PUV_dam parser = argparse.ArgumentParser(description="Post-process swash output") parser.add_argument("-v", "--verbose", action="count", default=0) parser.add_argument("-c", "--config", default="config.ini") args = parser.parse_args() logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) log = logging.getLogger("post") log.info("Starting post-processing") config = configparser.ConfigParser() config.read(args.config) inp = pathlib.Path(config.get("post", "inp")) root = pathlib.Path(config.get("swash", "out")) log.info(f"Reading data from '{inp}'") x = np.load(inp.joinpath("xp.npy")) t = np.load(inp.joinpath("tsec.npy")) botl = np.load(inp.joinpath("botl.npy")) watl = np.load(inp.joinpath("watl.npy")) vel = np.load(inp.joinpath("vel.npy")) x0 = config.getint("post", "x0") arg_x0 = np.abs(x - x0).argmin() t0 = config.getfloat("post", "t0") arg_t0 = np.abs(t - t0).argmin() dt = np.diff(t).mean() f = 1 / dt idx = [arg_x0 - 5, arg_x0, arg_x0 + 7] wl = watl[arg_t0:, idx] # Reflex3S res = reflex3S(*wl.T, *x[idx], botl[idx].mean(), f, 0.02, 0.2) f_, ai_, ar_ = reflex3s(wl.T, x[idx], botl[idx].mean(), f) ai, ar, _, _, _, _, fre = res # Custom PUV eta = watl[t > t0, arg_x0] u = vel[t > t0, 0, arg_x0] phi_eta = sgl.welch(eta, f) phi_u = sgl.welch(u, f) phi_eta_u = sgl.csd(eta, u, f) R = np.sqrt( (np.abs(phi_eta[1]) + np.abs(phi_u[1]) - 2 * phi_eta_u[1].real) / (np.abs(phi_eta[1]) + np.abs(phi_u[1]) + 2 * phi_eta_u[1].real) ) out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}") log.info(f"Saving plots in '{out}'") out.mkdir(parents=True, exist_ok=True) plt.plot(fre, np.asarray(ar) / np.asarray(ai)) plt.xlim((0, 0.3)) plt.ylim((0, 1)) plt.savefig(out.joinpath("reflex3s_pa.pdf")) plt.clf() plt.plot(f_, ar_ / ai_) plt.plot(phi_eta[0], R) plt.xlim((0, 0.3)) plt.ylim((0, 1)) plt.savefig(out.joinpath("reflex3s.pdf")) # pressk = np.load(inp.joinpath("pressk.npy")) # press = pressk[:, 1, :] # res = PUV_dam(vel[arg_t0:, 0, 500], press[arg_t0:, 500], botl[500], f, botl[500]/2, botl[500]/2) # print(res)