diff --git a/source/flat.py b/source/flat.py
new file mode 100644
index 0000000..e066728
--- /dev/null
+++ b/source/flat.py
@@ -0,0 +1,738 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu July 27 10:39:40 2023
+    post-process data for cooling/heating cases
+
+@author: weibo
+"""
+
+# %% load necessary modules
+from natsort import natsorted
+import imageio as iio
+import tecplot as tp
+from tecplot.constant import *
+from tecplot.exception import *
+import variable_analysis as va
+from planar_field import PlanarField as pf
+import pytecio as pytec
+import plt2pandas as p2p
+import pandas as pd
+from scipy.interpolate import griddata
+import numpy as np
+from glob import glob
+import matplotlib.pyplot as plt
+import matplotlib
+from scipy import signal
+from line_field import LineField as lf
+from IPython import get_ipython
+
+get_ipython().run_line_magic("matplotlib", "qt")
+
+# %% set path and basic parameters
+path = "/media/weibo/Work/flat_plate_st5/"
+# path = 'E:/cases/wavy_1009/'
+p2p.create_folder(path)
+pathP = path + "probes/"
+pathF = path + "Figures/"
+pathM = path + "MeanFlow/"
+pathS = path + "SpanAve/"
+pathT = path + "TimeAve/"
+pathI = path + "Instant/"
+pathV = path + "Vortex/"
+pathSL = path + "Slice/"
+tsize = 13
+nsize = 10
+matplotlib.rcParams["xtick.direction"] = "out"
+matplotlib.rcParams["ytick.direction"] = "out"
+matplotlib.rc("font", size=tsize)
+plt.close("All")
+plt.rc("text", usetex=True)
+font = {
+    "family": "Times New Roman",  # 'color' : 'k',
+    "weight": "normal",
+    "size": "large",
+}
+cm2in = 1 / 2.54
+
+# %% mean flow contour in x-y plane
+MeanFlow = pf()
+MeanFlow.load_meanflow(path)
+MeanFlow.copy_meanval()
+# %% merge mean flow
+MeanFlow.merge_stat(pathM)
+# %% rescaled if necessary
+lh = 1.0
+MeanFlow.rescale(lh)
+x, y = np.meshgrid(np.unique(MeanFlow.x), np.unique(MeanFlow.y))
+walldist = griddata((MeanFlow.x, MeanFlow.y),
+                    MeanFlow.walldist,
+                    (x, y),
+                    method="cubic")
+corner1 = (x < 25) & (y < 0.0)
+corner2 = (x > 110) & (y < 0.0)
+corner = (walldist < 0.0)  # corner1 | corner2
+
+# %% wall buondary
+va.wall_line(MeanFlow.PlanarData, pathM, mask=None) # corner)
+# %%
+walldist = griddata((MeanFlow.x, MeanFlow.y), MeanFlow.walldist, (x, y), method="cubic")
+cover = (walldist < 0.0)  # | corner
+NewFrame = MeanFlow.PlanarData
+va.dividing_line(
+    NewFrame, pathM, show=True, mask=cover
+)  # (MeanFlow.PlanarData, pathM, loc=2.5)
+# %% Save sonic line
+va.sonic_line(MeanFlow.PlanarData, pathM,
+              option="velocity", Ma_inf=6.0, mask=corner)
+# %% enthalpy boundary layer
+va.enthalpy_boundary_edge(MeanFlow.PlanarData, pathM, Ma_inf=6.0, crit=0.986)
+# %% thermal boundary layer
+va.thermal_boundary_edge(MeanFlow.PlanarData, pathM, T_wall=3.35)
+# %% Save shock line
+va.shock_line(MeanFlow.PlanarData, pathM, var="|gradp|", val=0.05)
+# %% boundary layer 
+va.boundary_edge(MeanFlow.PlanarData, pathM, shock=False, mask=False)
+
+# %% check mesh
+"""
+    Examination of the computational mesh
+"""
+temp = MeanFlow.PlanarData  # [["x", "y"]]
+df = temp.query("x>=-100.0 & x<=90.0 & y>=0.0 & y<=30.0")
+# df = temp.query("x>=-5.0 & x<=10.0 & y>=-3.0 & y<=2.0")
+ux = np.unique(df.x)[::2]
+uy = np.unique(df.y)[::2]
+fig, ax = plt.subplots(figsize=(14*cm2in, 6*cm2in))
+matplotlib.rc("font", size=tsize)
+for i in range(np.size(ux)):
+    if i % 1 == 0:
+        df_x = df.loc[df["x"] == ux[i]]
+        ax.plot(df_x["x"], df_x["y"], color="gray",
+                linestyle="-", linewidth=0.4)
+for j in range(np.size(uy)):
+    if j % 1 == 0:  # 4
+        df_y = df.loc[df["y"] == uy[j]]
+        ax.plot(df_y["x"], df_y["y"], color="gray",
+                linestyle="-", linewidth=0.4)
+# plt.gca().set_aspect("equal", adjustable="box")
+ax.set_xlim(np.min(df.x), np.max(df.x))
+ax.set_ylim(np.min(df.y), np.max(df.y))
+ax.set_yticks(np.linspace(0.0, 30.0, 7))
+ax.tick_params(labelsize=nsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$y/l_r$", fontsize=tsize)
+plt.tight_layout(pad=0.5, w_pad=0.5, h_pad=1)
+plt.savefig(pathF + "Grid.svg", bbox_inches="tight")
+plt.show()
+
+
+# %% Plot rho contour of the mean flow field
+"""
+    mean flow field contouring by density
+"""
+var = "u"
+lh = 1.0
+rho = griddata((MeanFlow.x, MeanFlow.y), MeanFlow.rho, (x, y))
+u = griddata((MeanFlow.x, MeanFlow.y), MeanFlow.u, (x, y))
+walldist = griddata((MeanFlow.x, MeanFlow.y), MeanFlow.walldist, (x, y))
+print("rho_max=", np.max(rho))
+print("rho_min=", np.min(rho))
+cval1 = 0.3  # 1.0 #
+cval2 = 1.0  # 7.0  #
+fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
+matplotlib.rc("font", size=tsize)
+rg1 = np.linspace(cval1, cval2, 21)
+cbar = ax.contourf(x, y, rho, cmap="rainbow", levels=rg1, extend="both")
+ax.contour(x, y, u, levels=[0.0], linestyles="dotted")  # rainbow_r
+ax.set_xlim(-200, 80)
+ax.set_ylim(np.min(y), 36.0)
+ax.set_xticks(np.linspace(-200, 80.0, 8))
+ax.set_yticks(np.linspace(np.min(y), 36.0, 5))
+ax.tick_params(labelsize=nsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$y/l_r$", fontsize=tsize)
+# plt.gca().set_aspect("equal", adjustable="box")
+# Add colorbar
+rg2 = np.linspace(cval1, cval2, 3)
+cbaxes = fig.add_axes([0.17, 0.76, 0.16, 0.07])  # x, y, width, height
+cbaxes.tick_params(labelsize=nsize)
+cbar = plt.colorbar(
+    cbar, cax=cbaxes, extendrect="False", orientation="horizontal", ticks=rg2
+)
+lab = r"$\langle \rho \rangle/\rho_{\infty}$"
+cbar.set_label(lab, rotation=0, fontsize=tsize)
+Enthalpy = pd.read_csv(
+    pathM + "EnthalpyBoundaryEdge99.dat", skipinitialspace=True)
+ax.plot(Enthalpy.x / lh, Enthalpy.y / lh, "k--", linewidth=1.5)
+# Add sonic line
+sonic = pd.read_csv(pathM + "SonicLine.dat", skipinitialspace=True)
+ax.plot(sonic.x / lh, sonic.y / lh, "w--", linewidth=1.5)
+# Add dividing line(separation line)
+# dividing = pd.read_csv(pathM + "BubbleLine.dat", skipinitialspace=True)
+# ax.plot(dividing.x / lh, dividing.y / lh, "k:", linewidth=1.5)
+plt.savefig(pathF + "MeanFlow_rho.svg", bbox_inches="tight")
+plt.show()
+
+
+# %% Plot rms contour of the mean flow field
+"""
+    Root Mean Square of velocity from the statistical flow
+"""
+x, y = np.meshgrid(np.unique(MeanFlow.x), np.unique(MeanFlow.y))
+var = "<u`u`>"
+if var == "<u`u`>":
+    lab = r"$\sqrt{\langle u^\prime u^\prime \rangle}$"
+    nm = 'RMSUU'
+elif var == "<v`v`>":
+    lab = r"$\sqrt{\langle v^\prime v^\prime \rangle}$"
+    nm = 'RMSVV'
+elif var == "<w`w`>":
+    lab = r"$\sqrt{\langle w^\prime w^\prime \rangle}$"
+    nm = 'RMSWW'
+var_val = getattr(MeanFlow.PlanarData, var)
+uu = griddata((MeanFlow.x, MeanFlow.y), var_val, (x, y))
+u = griddata((MeanFlow.x, MeanFlow.y), MeanFlow.u, (x, y))
+print("uu_max=", np.max(np.sqrt(np.abs(var_val))))
+print("uu_min=", np.min(np.sqrt(np.abs(var_val))))
+
+fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
+matplotlib.rc("font", size=tsize)
+cb1 = 0.0
+cb2 = 0.01
+rg1 = np.linspace(cb1, cb2, 21)  # 21)
+cbar = ax.contourf(
+    x / lh, y / lh, np.sqrt(np.abs(uu)), cmap="coolwarm",
+    levels=rg1, extend="both"
+)  # rainbow_r # jet # Spectral_r
+ax.set_xlim(-200, 80.0)
+ax.set_ylim(np.min(y), 36.0)
+ax.set_xticks(np.linspace(-200.0, 80.0, 8))
+ax.set_yticks(np.linspace(0.0, 36.0, 5))
+ax.tick_params(labelsize=nsize)
+ax.set_xlabel(r"$x$", fontdict=font)
+ax.set_ylabel(r"$y$", fontdict=font)
+# plt.gca().set_aspect("equal", adjustable="box")
+# Add colorbar box
+rg2 = np.linspace(cb1, cb2, 3)
+cbbox = fig.add_axes([0.14, 0.53, 0.22, 0.31], alpha=0.9)
+[cbbox.spines[k].set_visible(False) for k in cbbox.spines]
+cbbox.tick_params(
+    axis="both",
+    left=False,
+    right=False,
+    top=False,
+    bottom=False,
+    labelleft=False,
+    labeltop=False,
+    labelright=False,
+    labelbottom=False,
+)
+cbbox.set_facecolor([1, 1, 1, 0.7])
+# Add colorbar
+cbaxes = fig.add_axes([0.16, 0.78, 0.18, 0.058])  # x, y, width, height
+cbaxes.tick_params(labelsize=nsize)
+cbar = plt.colorbar(
+    cbar, cax=cbaxes, orientation="horizontal", extendrect="False", ticks=rg2
+)
+
+cbar.set_label(
+    lab, rotation=0, fontsize=tsize,
+)
+# Add boundary layer
+boundary = pd.read_csv(
+    pathM + "EnthalpyBoundaryEdge99.dat", skipinitialspace=True)
+ax.plot(boundary.x / lh, boundary.y / lh, "k--", linewidth=1.0)
+# Add bubble line
+# bubble = pd.read_csv(pathM + "BubbleLine.dat", skipinitialspace=True)
+# ax.plot(bubble.x / lh, bubble.y / lh, "k:", linewidth=1.0)
+# Add sonic line
+sonic = pd.read_csv(pathM + "SonicLine.dat", skipinitialspace=True)
+ax.plot(sonic.x / lh, sonic.y / lh, "w--", linewidth=1.0)
+plt.savefig(pathF + "MeanFlow" + nm + ".svg", bbox_inches="tight")
+plt.show()
+
+
+# %% Numerical schiliren in X-Y plane
+"""
+    Numerical schiliren (gradient of density)
+"""
+flow = pf()
+file = path + "snapshots/" + "TP_2D_Z_001.szplt"
+flow.load_data(path, FileList=file, NameList='tecio')
+plane = flow.PlanarData
+# plane.to_hdf(path+'snapshots/TP_2D_Z_03_1000.h5', 'w', format='fixed')
+x, y = np.meshgrid(np.unique(plane.x), np.unique(plane.y))
+var = "|grad(rho)|"
+schlieren = griddata((plane.x, plane.y), plane[var], (x, y))
+print("rho_max=", np.max(plane[var]))
+print("rho_min=", np.min(plane[var]))
+va.sonic_line(plane, pathI, option="velocity", Ma_inf=6.0, mask=None)
+va.enthalpy_boundary_edge(plane, pathI, Ma_inf=6.0, crit=0.986)
+va.dividing_line(plane, pathI, show=True, mask=None)
+# schlieren[corner] = np.nan
+# %%
+fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
+matplotlib.rc("font", size=tsize)
+rg1 = np.linspace(0, 3.0, 21)
+cbar = ax.contourf(
+    x, y, schlieren, cmap="binary", levels=rg1, extend="both"
+)  # binary #rainbow_r# bwr
+ax.set_xlim(-200, 80.0)
+ax.set_ylim(np.min(y), 36.0)
+ax.set_xticks(np.linspace(-200, 80.0, 8))
+ax.set_yticks(np.linspace(0.0, 36.0, 5))
+ax.tick_params(labelsize=nsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$y/l_r$", fontsize=tsize)
+# plt.gca().set_aspect("equal", adjustable="box")
+# Add colorbar
+rg2 = np.linspace(0, 3.0, 3)
+cbbox = fig.add_axes([0.16, 0.50, 0.20, 0.30], alpha=0.9)
+[cbbox.spines[k].set_visible(False) for k in cbbox.spines]
+cbbox.tick_params(
+    axis="both",
+    left=False,
+    right=False,
+    top=False,
+    bottom=False,
+    labelleft=False,
+    labeltop=False,
+    labelright=False,
+    labelbottom=False,
+)
+cbbox.set_facecolor([1, 1, 1, 0.7])
+# Add colorbar
+cbaxes = fig.add_axes([0.17, 0.71, 0.18, 0.07])  # x, y, width, height
+cbaxes.tick_params(labelsize=nsize)
+cbar = plt.colorbar(
+    cbar, cax=cbaxes, orientation="horizontal", extendrect="False", ticks=rg2
+)
+cbar.set_label(
+    r"$|\nabla \rho|$", rotation=0, fontsize=tsize - 1,
+)
+# Add boundary layer
+boundary = pd.read_csv(
+    pathI + "EnthalpyBoundaryEdge.dat", skipinitialspace=True)
+# ax.plot(boundary.x, boundary.y, "g--", linewidth=1.0)
+# Add dividing line(separation line)
+dividing = pd.read_csv(pathI + "BubbleLine.dat", skipinitialspace=True)
+# ax.plot(dividing.x, dividing.y, "g:", linewidth=1.0)
+# Add sonic line
+sonic = pd.read_csv(pathI + "SonicLine.dat", skipinitialspace=True)
+# ax.plot(sonic.x, sonic.y, "b--", linewidth=1.0)
+plt.savefig(pathF + "SchlierenXY.svg", bbox_inches="tight")
+plt.show()
+# %% BL profile along streamwise
+fig, ax = plt.subplots(1, 9, figsize=(13 * cm2in, 4 * cm2in), dpi=500)
+fig.subplots_adjust(hspace=0.5, wspace=0.25)
+matplotlib.rc("font", size=nsize)
+title = [r"$(a)$", r"$(b)$", r"$(c)$", r"$(d)$", r"$(e)$"]
+matplotlib.rcParams["xtick.direction"] = "in"
+matplotlib.rcParams["ytick.direction"] = "in"
+xcoord = np.array([-120, -80, -40, -20, 0, 20, 40, 60, 80])
+for i in range(np.size(xcoord)):
+    df = MeanFlow.yprofile("x", xcoord[i])
+    y0 = df["walldist"]
+    q0 = df["u"]
+    if xcoord[i] == 0.0:
+        ind = np.where(y0 >= 0.0)[0]
+        ax[i].plot(q0[ind], y0[ind], "k-")
+        ax[i].set_ylim([0, 6])
+    else:
+        ind = np.where(y0 >= 0.0)[0]
+        ax[i].plot(q0[ind], y0[ind], "k-")
+        ax[i].set_ylim([0, 6])
+    if i != 0:
+        ax[i].set_yticklabels("")
+        ax[i].set_title(r"${}$".format(xcoord[i]), fontsize=nsize - 2)
+    ax[i].set_xticks([0, 1], minor=True)
+    ax[i].tick_params(axis="both", which="major", labelsize=nsize)
+    ax[i].grid(visible=True, which="both", linestyle=":")
+ax[0].set_title(r"$x={}$".format(xcoord[0]), fontsize=nsize - 2)
+ax[0].set_ylabel(r"$\Delta y$", fontsize=tsize)
+ax[4].set_xlabel(r"$u /u_\infty$", fontsize=tsize)
+plt.tick_params(labelsize=nsize)
+plt.savefig(pathF + "BLProfileU.svg", bbox_inches="tight", pad_inches=0.1)
+plt.show()
+
+
+# %% dataframe of the first level points
+# corner = MeanFlow.PlanarData.walldist < 0.0
+firstval = 0.03125
+ind0 = MeanFlow.PlanarData.y == 0.0
+MeanFlow.PlanarData['walldist'][ind0] = 0.0
+xy_val = va.wall_line(MeanFlow.PlanarData, path, val=firstval)  # mask=corner)
+# onelev = pd.DataFrame(data=xy_val, columns=["x", "y"])
+# onelev.drop_duplicates(subset='x', keep='first', inplace=True)
+FirstLev = pd.DataFrame(data=xy_val, columns=["x", "y"])
+# ind = FirstLev.index[FirstLev['x'] < 0]
+# FirstLev.loc[ind, 'y'] = firstval
+xx = np.arange(-200.0, 80.0, 0.25)
+FirstLev = FirstLev[FirstLev["x"].isin(xx)]
+FirstLev.to_csv(pathM + "FirstLev.dat", index=False, float_format="%9.6f")
+
+
+# %%
+firstval = 0.03125
+T_inf = 86.6
+Re = 7736
+df = MeanFlow.PlanarData
+ramp_wall = pd.read_csv(pathM + "FirstLev.dat", skipinitialspace=True)
+ramp_wall = va.add_variable(df, ramp_wall)
+# ramp1 = ramp_wall.loc[ramp_wall['x'] <= 0.0]
+# ramp2 = ramp_wall.loc[ramp_wall['x'] > 0.0]
+# for flat plate
+mu = va.viscosity(Re, ramp_wall["T"])
+xwall = ramp_wall['x']
+Cf = va.skinfriction(mu, ramp_wall['u'], firstval, factor=1)
+# %% for ramp
+angle = 15 / 180 * np.pi
+mu = va.viscosity(Re, ramp2["T"])
+ramp2_u = ramp2['u'] * np.cos(angle) + ramp2['v'] * np.sin(angle)
+xwall2 = ramp2['x']
+Cf2 = va.skinfriction(mu, ramp2_u, firstval*np.cos(angle), factor=1)
+ind = np.argmin(np.abs(Cf1[:-20]))
+xline1 = xwall1[ind]
+xline2 = np.round(np.interp(0.0, Cf2, xwall2), 2)
+yline = np.max(Cf2) * 0.92
+# xwall, Cf = va.skinfric_wavy(pathM, wavy, Re, T_inf, firstval)
+# %% skin friction
+fig2, ax2 = plt.subplots(figsize=(15*cm2in, 5*cm2in), dpi=500)
+matplotlib.rc("font", size=nsize)
+ax2.plot(xwall, Cf, "b-", linewidth=1.5)
+ax2.set_xlabel(r"$x$", fontsize=tsize)
+ax2.set_ylabel(r"$\langle C_f \rangle$", fontsize=tsize)
+ax2.set_xlim([-200.0, 80.0])
+ax2.set_xticks(np.linspace(-200.0, 80.0, 8))
+ax2.ticklabel_format(axis="y", style="sci", scilimits=(-2, 2))
+# ax2.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax2.annotate(r"$x_s={}$".format(xline1), (xline1-16, yline))
+# ax2.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+# ax2.annotate(r"$x_r={}$".format(xline2), (xline2-16, yline))
+ax2.grid(visible=True, which="both", linestyle=":")
+ax2.yaxis.offsetText.set_fontsize(nsize)
+ax2.annotate("(a)", xy=(-0.09, 1.0), xycoords="axes fraction", fontsize=nsize)
+plt.tick_params(labelsize=nsize)
+plt.savefig(pathF + "Cf.svg", bbox_inches="tight", pad_inches=0.1)
+plt.show()
+
+# %% Cp
+Ma = 6.0
+fa = 1.0  # Ma * Ma * 1.4
+b, a = signal.butter(6, Wn=1 / 12, fs=1 / 0.125)
+p_fit = signal.filtfilt(b, a, ramp_wall["p"] * fa)
+fig3, ax3 = plt.subplots(figsize=(15*cm2in, 5*cm2in), dpi=500)
+# ax3 = fig.add_subplot(212)
+ax3.plot(ramp_wall['x'], ramp_wall["p"] * fa, "b-", linewidth=1.5)
+# ax3.plot(wavy["x"], p_fit, "k", linewidth=1.5)
+ax3.set_xlabel(r"$x$", fontsize=tsize)
+# ylab = r"$\langle p_w \rangle/\rho_{\infty} u_{\infty}^2$"
+ax3.set_ylabel(r"$\langle C_p \rangle$", fontsize=tsize)
+ax3.set_xlim([-200.0, 80])
+ax3.set_xticks(np.arange(-200.0, 80.0, 40))
+ax3.ticklabel_format(axis="y", style="sci", scilimits=(-2, 2))
+# ax3.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax3.annotate(r"$x_s={}$".format(xline1), (xline1-16, yline))
+# ax3.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+# ax3.annotate(r"$x_r={}$".format(xline2), (xline2-16, yline))
+ax3.grid(visible=True, which="both", linestyle=":")
+ax3.annotate("(b)", xy=(-0.13, 1.0), xycoords="axes fraction", fontsize=nsize)
+plt.tick_params(labelsize=nsize)
+plt.savefig(pathF + "Cp.svg", bbox_inches="tight", pad_inches=0.1)
+plt.show()
+# %% Stanton number
+# firstval = 0.015625
+T_wall = 3.35
+mu = va.viscosity(Re, ramp_wall["T"], T_inf=T_inf, law="Suther")
+kt = va.thermal(mu, 0.72)
+Tt = va.stat2tot(Ma=6.0, Ts=1.0, opt="t")
+Cs = va.Stanton(ramp_wall["T"], T_wall, firstval, 
+                Re, Ma, T_inf=1.0, factor=1)
+
+# %% low-pass filter
+xwall = ramp_wall['x']
+b, a = signal.butter(6, Wn=1 / 12, fs=1 / 0.125)
+Cs_fit = signal.filtfilt(b, a, Cs)
+# ind = np.where(Cf[:] < 0.008)
+fig2, ax2 = plt.subplots(figsize=(15*cm2in, 5*cm2in), dpi=500)
+# fig = plt.figure(figsize=(8, 5), dpi=500)
+# ax2 = fig.add_subplot(211)
+matplotlib.rc("font", size=nsize)
+ax2.plot(xwall, np.abs(Cs), "b-", linewidth=1.5)
+# ax2.plot(xwall, Cs_fit, "k", linewidth=1.5)
+ax2.set_xlabel(r"$x$", fontsize=tsize)
+ax2.set_ylabel(r"$\langle C_s \rangle$", fontsize=tsize)
+ax2.set_xlim([-200, 80])
+ax2.set_xticks(np.arange(-200.0, 80.0, 40))
+ax2.ticklabel_format(axis="y", style="sci", scilimits=(-2, 2))
+# ax2.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax2.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+ax2.grid(visible=True, which="both", linestyle=":")
+ax2.yaxis.offsetText.set_fontsize(nsize)
+ax2.annotate("(c)", xy=(-0.09, 1.0), xycoords="axes fraction", fontsize=nsize)
+plt.tick_params(labelsize=nsize)
+plt.savefig(pathF + "Cs.svg", bbox_inches="tight", pad_inches=0.1)
+plt.show()
+
+# %% Time evolution of a specific variable at several streamwise locations
+"""
+    temporal evolution of signals along an axis
+"""
+probe = lf()
+lh = 1.0
+fa = 1.0
+var = "u"
+ylab = r"$u^\prime/u_\infty$"
+if var == "p":
+    fa = 6.0 * 6.0 * 1.4
+    ylab = r"$p^\prime/p_\infty$"
+timezone = np.arange(600.0, 900.0 + 0.125, 0.125)
+xloc = [-160.0, -80.0, -40.0, -20.0, 0.0]  # [-50.0, -30.0, -20.0, -10.0, 4.0]
+yloc = [0, 0, 0, 0, 0]
+zloc = 0.0
+fig, ax = plt.subplots(np.size(xloc), 1, figsize=(15*cm2in, 10*cm2in))
+fig.subplots_adjust(hspace=0.6, wspace=0.15)
+matplotlib.rc("font", size=nsize)
+matplotlib.rcParams["xtick.direction"] = "in"
+matplotlib.rcParams["ytick.direction"] = "in"
+for i in range(np.size(xloc)):
+    probe.load_probe(pathP, (xloc[i], yloc[i], zloc))
+    probe.rescale(lh)
+    probe.extract_series([timezone[0], timezone[-1]])
+    temp = (getattr(probe, var) - np.mean(getattr(probe, var))) * fa
+    ax[i].plot(probe.time, temp, "k")
+    ax[i].ticklabel_format(
+        axis="y", style="sci", useOffset=False, scilimits=(-2, 2)
+    )
+    ax[i].set_xlim([timezone[0], timezone[-1]])
+    # ax[i].set_ylim([-0.001, 0.001])
+    if i != np.size(xloc) - 1:
+        ax[i].set_xticklabels("")
+    ax[i].set_ylabel(ylab, fontsize=tsize)
+    ax[i].grid(visible=True, which="both", linestyle=":")
+    ax[i].set_title(r"$x/\delta_0={}$".format(xloc[i]), fontsize=nsize - 1)
+ax[-1].set_xlabel(r"$t u_\infty/\delta_0$", fontsize=tsize)
+
+plt.savefig(
+    pathF + var + "_TimeEvolveX" + str(zloc) + ".svg", bbox_inches="tight"
+)
+plt.show()
+
+# %%############################################################################
+
+
+# %% save and load data
+"""
+    streamwise evolution of signals along an axis
+"""
+df, SolTime = pytec.ReadINCAResults(path + 'TP_fluc_3d/',
+                                    SavePath=pathI,
+                                    SpanAve=False,
+                                    OutFile='TP_fluc_2d')
+
+# %%
+MeanFlow = pf()
+MeanFlow.merge_stat(pathI)
+# %% Streamwise evolution of a specific variable
+fa = 1.0  # 1.7 * 1.7 * 1.4
+var = "p`"
+df = pd.read_hdf(pathI + "TP_fluc_2d.h5")
+ramp_wall = pd.read_csv(pathM + "FirstLev.dat", skipinitialspace=True)
+ramp_wall = va.add_variable(df, ramp_wall)
+ramp_wall[var] = griddata((df.x, df.y), df[var],
+                          (ramp_wall.x, ramp_wall.y),
+                          method="linear")
+
+meanval = np.mean(ramp_wall[var])
+fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
+matplotlib.rc("font", size=nsize)
+ax.plot(ramp_wall.x, ramp_wall[var] - meanval, "k")
+ax.set_xlim([-200.0, 80.0])
+ax.set_xticks(np.linspace(-200.0, 80.0, 8))
+# ax.set_ylim([-0.001, 0.001])
+# ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$2p^\prime/\rho_\infty u_\infty^2$", fontsize=tsize)
+ax.ticklabel_format(axis="y", style="sci", useOffset=False, scilimits=(-2, 2))
+ax.grid(visible="both", linestyle=":")
+plt.show()
+plt.savefig(
+    pathF + var + "_streamwise.svg", bbox_inches="tight", pad_inches=0.1
+)
+
+# %% Streamwise evolution of a specific variable
+fa = 1.0  # 1.7 * 1.7 * 1.4
+var = "u`"
+ramp_wall[var] = griddata((df.x, df.y), df[var],
+                          (ramp_wall.x, ramp_wall.y),
+                          method="linear")
+meanval = np.mean(ramp_wall[var])
+fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
+matplotlib.rc("font", size=nsize)
+ax.plot(ramp_wall.x, ramp_wall[var] - meanval, "k")
+ax.set_xlim([-200.0, 80.0])
+ax.set_xticks(np.linspace(-200.0, 80.0, 8))
+# ax.set_ylim([-0.001, 0.001])
+ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+# ax.set_ylabel(r"$p^\prime/\rho_\infty u_\infty^2$", fontsize=tsize)
+# ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+ax.set_ylabel(r"$u^\prime/u_\infty$", fontsize=tsize)
+ax.ticklabel_format(axis="y", style="sci", useOffset=False, scilimits=(-2, 2))
+ax.grid(visible="both", linestyle=":")
+plt.show()
+plt.savefig(
+    pathF + var + "_streamwise.svg", bbox_inches="tight", pad_inches=0.1
+)
+
+# %% the locations with most energetic fluctuations 
+"""
+    streamwise evolution of the most energetic signals
+"""
+FlucFlow = pd.read_hdf(pathM + "MeanFlow.h5")
+varn = '<u`u`>'
+if varn == '<u`u`>':
+    savenm = "MaxRMS_u"
+    ylab = r"$\sqrt{u^{\prime 2}_\mathrm{max}}/u_{\infty}$"
+elif varn =='<p`p`>':
+    savenm = "MaxRMS_p"
+    ylab = r"$\sqrt{p^{\prime 2}_\mathrm{max}}/\rho_{\infty} u_{\infty}^2$"
+xnew = np.arange(-200.0, 80.0, 0.25)
+znew = np.zeros(np.size(xnew))
+varv = np.zeros(np.size(xnew))
+ynew = np.zeros(np.size(xnew))
+for i in range(np.size(xnew)):
+    df = va.max_pert_along_y(FlucFlow, '<u`u`>', [xnew[i], znew[i]])
+    varv[i] = df[varn]
+    ynew[i] = df['y']
+data = np.vstack((xnew, ynew, varv))
+df = pd.DataFrame(data.T, columns=['x', 'y', varn])
+# df = df.drop_duplicates(keep='last')
+df.to_csv(pathM + savenm + ".dat", sep=' ',
+          float_format='%1.8e', index=False)
+
+# %% draw maximum value curve
+fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
+matplotlib.rc('font', size=nsize)
+ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.set_ylabel(r"$y/\delta_0$", fontsize=tsize)
+ax.plot(xnew, ynew, 'k', label=r'$q^\prime_\mathrm{max}$')
+ax.set_xticks(np.linspace(-200.0, 80.0, 8))
+# ax.plot(xx, yy, 'k--', label='bubble')
+legend = ax.legend(loc='upper right', shadow=False, fontsize=nsize)
+ax.ticklabel_format(axis="y", style="sci", scilimits=(-2, 2))
+ax.grid(b=True, which="both", linestyle=":")
+plt.show()
+plt.savefig(
+    pathF + "MaxPertLoc" + savenm + ".svg", bbox_inches="tight", pad_inches=0.1
+)
+
+# %% draw maximum value curve
+ramp_max = pd.read_csv(pathM + savenm + ".dat", sep=' ', skipinitialspace=True)
+fig, ax = plt.subplots(figsize=(15*cm2in, 5.0*cm2in))
+matplotlib.rc('font', size=nsize)
+ax.set_ylabel(ylab, fontsize=tsize)
+ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.plot(ramp_max['x'], ramp_max[varn], 'k')
+ax.set_xlim([-200.0, 80])
+# ax.set_ylim([1e-7, 1e-1])
+ax.set_xticks(np.linspace(-200.0, 80.0, 8))
+ax.set_yscale('log')
+# ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+ax.grid(visible=True, which="both", linestyle=":")
+ax.tick_params(labelsize=nsize)
+plt.show()
+plt.savefig(
+    pathF + "MaxFluc" + savenm + ".svg", bbox_inches="tight", pad_inches=0.1
+)
+# %% test section
+thick = pd.read_csv(pathM + 'fit_thickness.dat', sep=' ')
+shape_fa = thick['displace'] / thick['momentum']
+xynew = pd.read_csv(pathM + 'MaxRMS_u.dat', sep=' ')
+fig, ax = plt.subplots(figsize=(6.4, 2.2))
+matplotlib.rc('font', size=nsize)
+ylab = r"$\sqrt{u^{\prime 2}_\mathrm{max}}/u_{\infty}$"
+ax.set_ylabel(ylab, fontsize=tsize)
+ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.plot(xynew['x'], xynew['<u`u`>'], 'k')
+ax.set_yscale('log')
+ax.set_xlim([-100, 20])
+# ax.ticklabel_format(axis="y", style="sci", scilimits=(-1, 1))
+ax.grid(b=True, which="both", linestyle=":")
+ax.tick_params(labelsize=nsize)
+ax.yaxis.offsetText.set_fontsize(nsize-1)
+
+ax2 = ax.twinx()
+ax2.plot(thick['x'], shape_fa, 'k--')
+ax2.set_ylabel(r"$H$", fontsize=tsize)
+
+plt.show()
+plt.savefig(
+    pathF + "MaxRMS_x.svg", bbox_inches="tight", pad_inches=0.1
+)
+
+
+# %% Draw impose mode
+
+inmode = pd.read_csv(path+"UnstableMode.inp", skiprows=5,
+                     sep=' ', index_col=False)
+fig, ax = plt.subplots(figsize=(7*cm2in, 6.5*cm2in))
+matplotlib.rc('font', size=nsize)
+xlab = r"$|q^{\prime}|/|u^\prime|_{\max}$"
+ax.set_xlabel(xlab, fontsize=tsize)
+ax.set_ylabel(r"$y/l_f$", fontsize=tsize)
+ax.plot(np.sqrt(inmode['u_r']**2+inmode['u_i']**2), inmode['y'], 'k')
+ax.plot(np.sqrt(inmode['v_r']**2+inmode['v_i']**2), inmode['y'], 'r')
+ax.plot(np.sqrt(inmode['w_r']**2+inmode['w_i']**2), inmode['y'], 'g')
+ax.plot(np.sqrt(inmode['p_r']**2+inmode['p_i']**2), inmode['y'], 'b')
+ax.plot(np.sqrt(inmode['t_r']**2+inmode['t_i']**2), inmode['y'], 'c')
+# ax.set_xlim([-100, 20])
+ax.set_ylim([0.0, 8.0])
+# ax.ticklabel_format(axis="y", style="sci", scilimits=(-1, 1))
+ax.legend(['u', 'v', 'w', 'p', 'T'])
+ax.grid(visible=True, which="both", linestyle=":")
+ax.tick_params(labelsize=nsize)
+
+plt.show()
+plt.savefig(
+    pathF + "ModeProf.svg", bbox_inches="tight", pad_inches=0.1
+)
+
+# %% animation for vortex structure
+pathF = path + "Figures/"
+pathall = glob(path + "TP_data_*/")
+tp.session.connect(port=7600)
+for i in range(np.size(pathall)):
+    dirs = glob(pathall[i] + "*.szplt")
+    # tp.session.stop()
+    print("process " + pathall[i])
+    dataset = tp.data.load_tecplot_szl(dirs, read_data_option=2)
+    soltime = np.round(dataset.solution_times[0], 2)
+    # with tp.session.suspend():
+    # tp.session.suspend_enter()
+    frame = tp.active_frame()
+    frame.width = 12.8
+    frame.height = 7.8
+    frame.position = (-1.4, 0.25)
+    # tp.macro.execute_file('/path/to/macro_file.mcr')
+    frame.load_stylesheet(path + "L2_2021.sty")
+    tp.macro.execute_command("$!Interface ZoneBoundingBoxMode = Off")
+    tp.macro.execute_command("$!FrameLayout ShowBorder = No")
+    tp.export.save_png(pathF + str(soltime) + ".png", width=2048)
+# %% create videos
+
+pathani = path + "ani/"
+prenm = "ramp"
+dirs = glob(pathani + "*.png")
+dirs = natsorted(dirs, key=lambda y: y.lower())
+flnm = pathani + prenm + "_Anima.mp4"
+recyc = np.size(dirs)
+with iio.get_writer(flnm, mode="I", fps=6, macro_block_size=None) as writer:
+    for ii in range(recyc * 6):
+        ind = ii % recyc  # mod, get reminder
+        image = iio.v3.imread(dirs[ind])
+        writer.append_data(image)
+    writer.close()
diff --git a/source/pytecio.py b/source/pytecio.py
index 91eda9c..26ec897 100644
--- a/source/pytecio.py
+++ b/source/pytecio.py
@@ -1245,14 +1245,18 @@ def flow3dto2d(files, outpath, SpanAve=False):
     #     "E:/cases/wavy_0918/TP_fluc_00100056/",
     # )
     # path = '/mnt/work/cases_new/heating2/'
-    path = '/media/weibo/VID21/ramp_st14/'
+    path = "/media/weibo/Work/flat_plate_st5/"
     pathST = path + 'TP_stat/'
-    ReadINCAResults(pathST, SpanAve=True, SavePath=path, OutFile='TP_data')
-    create_fluc_inca(path + 'TP_data_00363670/',
+    pathI = path + 'Instant/'
+    # ReadINCAResults(pathST, SpanAve=True, SavePath=path, OutFile='TP_data')
+    create_fluc_inca(path + 'TP_data_00134703/',
                      pathST,
                      path + 'TP_fluc_3d/')
-
-
+    
+    df, SolTime = ReadINCAResults(path + 'TP_fluc_3d/',
+                                  SavePath=pathI,
+                                  SpanAve=False,
+                                  OutFile='TP_fluc_2d')
 
 
     # %%
diff --git a/source/ramp.py b/source/ramp.py
index 8ac1006..5399224 100644
--- a/source/ramp.py
+++ b/source/ramp.py
@@ -120,8 +120,8 @@
 ax.set_ylim(np.min(df.y), np.max(df.y))
 ax.set_yticks(np.linspace(0.0, 30.0, 7))
 ax.tick_params(labelsize=nsize)
-ax.set_xlabel(r"$x$", fontsize=tsize)
-ax.set_ylabel(r"$y$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$y/l_r$", fontsize=tsize)
 plt.tight_layout(pad=0.5, w_pad=0.5, h_pad=1)
 plt.savefig(pathF + "Grid.svg", bbox_inches="tight")
 plt.show()
@@ -150,8 +150,8 @@
 ax.set_xticks(np.linspace(-200, 80.0, 8))
 ax.set_yticks(np.linspace(np.min(y), 36.0, 5))
 ax.tick_params(labelsize=nsize)
-ax.set_xlabel(r"$x$", fontsize=tsize)
-ax.set_ylabel(r"$y$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$y/l_r$", fontsize=tsize)
 # plt.gca().set_aspect("equal", adjustable="box")
 # Add colorbar
 rg2 = np.linspace(cval1, cval2, 3)
@@ -180,7 +180,7 @@
     Root Mean Square of velocity from the statistical flow
 """
 x, y = np.meshgrid(np.unique(MeanFlow.x), np.unique(MeanFlow.y))
-var = "<w`w`>"
+var = "<u`u`>"
 if var == "<u`u`>":
     lab = r"$\sqrt{\langle u^\prime u^\prime \rangle}$"
     nm = 'RMSUU'
@@ -199,7 +199,7 @@
 fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
 matplotlib.rc("font", size=tsize)
 cb1 = 0.0
-cb2 = 0.05
+cb2 = 0.2
 rg1 = np.linspace(cb1, cb2, 21)  # 21)
 cbar = ax.contourf(
     x / lh, y / lh, np.sqrt(np.abs(uu)), cmap="coolwarm",
@@ -210,8 +210,8 @@
 ax.set_xticks(np.linspace(-200.0, 80.0, 8))
 ax.set_yticks(np.linspace(0.0, 36.0, 5))
 ax.tick_params(labelsize=nsize)
-ax.set_xlabel(r"$x$", fontdict=font)
-ax.set_ylabel(r"$y$", fontdict=font)
+ax.set_xlabel(r"$x/l_r$", fontdict=font)
+ax.set_ylabel(r"$y/l_r$", fontdict=font)
 # plt.gca().set_aspect("equal", adjustable="box")
 # Add colorbar box
 rg2 = np.linspace(cb1, cb2, 3)
@@ -313,7 +313,7 @@
 )
 # Add boundary layer
 boundary = pd.read_csv(
-    pathI + "EnthalpyBoundaryEdge99.dat", skipinitialspace=True)
+    pathI + "EnthalpyBoundaryEdge.dat", skipinitialspace=True)
 # ax.plot(boundary.x, boundary.y, "g--", linewidth=1.0)
 # Add dividing line(separation line)
 dividing = pd.read_csv(pathI + "BubbleLine.dat", skipinitialspace=True)
@@ -402,7 +402,7 @@
 matplotlib.rc("font", size=nsize)
 ax2.plot(xwall1, Cf1, "b-", linewidth=1.5)
 ax2.plot(xwall2, Cf2, "b-", linewidth=1.5)
-ax2.set_xlabel(r"$x$", fontsize=tsize)
+ax2.set_xlabel(r"$x/l_r$", fontsize=tsize)
 ax2.set_ylabel(r"$\langle C_f \rangle$", fontsize=tsize)
 ax2.set_xlim([-200.0, 80.0])
 ax2.set_xticks(np.linspace(-200.0, 80.0, 8))
@@ -427,7 +427,7 @@
 # ax3 = fig.add_subplot(212)
 ax3.plot(ramp_wall['x'], ramp_wall["p"] * fa, "b-", linewidth=1.5)
 # ax3.plot(wavy["x"], p_fit, "k", linewidth=1.5)
-ax3.set_xlabel(r"$x$", fontsize=tsize)
+ax3.set_xlabel(r"$x/l_r$", fontsize=tsize)
 # ylab = r"$\langle p_w \rangle/\rho_{\infty} u_{\infty}^2$"
 ax3.set_ylabel(r"$\langle C_p \rangle$", fontsize=tsize)
 ax3.set_xlim([-200.0, 80])
@@ -462,7 +462,7 @@
 matplotlib.rc("font", size=nsize)
 ax2.plot(xwall, np.abs(Cs), "b-", linewidth=1.5)
 # ax2.plot(xwall, Cs_fit, "k", linewidth=1.5)
-ax2.set_xlabel(r"$x$", fontsize=tsize)
+ax2.set_xlabel(r"$x/l_r$", fontsize=tsize)
 ax2.set_ylabel(r"$\langle C_s \rangle$", fontsize=tsize)
 ax2.set_xlim([-200, 80])
 ax2.set_xticks(np.arange(-200.0, 80.0, 40))
@@ -512,7 +512,7 @@
         ax[i].set_xticklabels("")
     ax[i].set_ylabel(ylab, fontsize=tsize)
     ax[i].grid(visible=True, which="both", linestyle=":")
-    ax[i].set_title(r"$x/\delta_0={}$".format(xloc[i]), fontsize=nsize - 1)
+    ax[i].set_title(r"$x/l_r={}$".format(xloc[i]), fontsize=nsize - 1)
 ax[-1].set_xlabel(r"$t u_\infty/\delta_0$", fontsize=tsize)
 
 plt.savefig(
@@ -531,11 +531,13 @@
                                     SavePath=pathI,
                                     SpanAve=False,
                                     OutFile='TP_fluc_2d')
+
+# %%
 MeanFlow = pf()
-MeanFlow.merge_stat(pathM)
+MeanFlow.merge_stat(pathI)
 # %% Streamwise evolution of a specific variable
 fa = 1.0  # 1.7 * 1.7 * 1.4
-var = "p`"
+var = "u`"
 df = pd.read_hdf(pathI + "TP_fluc_2d.h5")
 ramp_wall = pd.read_csv(pathM + "FirstLev.dat", skipinitialspace=True)
 ramp_wall = va.add_variable(df, ramp_wall)
@@ -552,7 +554,7 @@
 # ax.set_ylim([-0.001, 0.001])
 ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
 ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
-ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
 ax.set_ylabel(r"$2p^\prime/\rho_\infty u_\infty^2$", fontsize=tsize)
 ax.ticklabel_format(axis="y", style="sci", useOffset=False, scilimits=(-2, 2))
 ax.grid(visible="both", linestyle=":")
@@ -563,7 +565,7 @@
 
 # %% Streamwise evolution of a specific variable
 fa = 1.0  # 1.7 * 1.7 * 1.4
-var = "u`"
+var = "p`"
 ramp_wall[var] = griddata((df.x, df.y), df[var],
                           (ramp_wall.x, ramp_wall.y),
                           method="linear")
@@ -574,7 +576,7 @@
 ax.set_xlim([-200.0, 80.0])
 ax.set_xticks(np.linspace(-200.0, 80.0, 8))
 # ax.set_ylim([-0.001, 0.001])
-ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
 # ax.set_ylabel(r"$p^\prime/\rho_\infty u_\infty^2$", fontsize=tsize)
 ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
 ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
@@ -615,8 +617,8 @@
 # %% draw maximum value curve
 fig, ax = plt.subplots(figsize=(15*cm2in, 5*cm2in))
 matplotlib.rc('font', size=nsize)
-ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
-ax.set_ylabel(r"$y/\delta_0$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
+ax.set_ylabel(r"$y/l_r$", fontsize=tsize)
 ax.plot(xnew, ynew, 'k', label=r'$q^\prime_\mathrm{max}$')
 ax.set_xticks(np.linspace(-200.0, 80.0, 8))
 # ax.plot(xx, yy, 'k--', label='bubble')
@@ -633,14 +635,14 @@
 fig, ax = plt.subplots(figsize=(15*cm2in, 5.0*cm2in))
 matplotlib.rc('font', size=nsize)
 ax.set_ylabel(ylab, fontsize=tsize)
-ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
 ax.plot(ramp_max['x'], ramp_max[varn], 'k')
 ax.set_xlim([-200.0, 80])
 # ax.set_ylim([1e-7, 1e-1])
 ax.set_xticks(np.linspace(-200.0, 80.0, 8))
 ax.set_yscale('log')
-ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
-ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
+# ax.axvline(x=xline1, color="gray", linestyle=":", linewidth=1.5)
+# ax.axvline(x=xline2, color="gray", linestyle=":", linewidth=1.5)
 ax.grid(visible=True, which="both", linestyle=":")
 ax.tick_params(labelsize=nsize)
 plt.show()
@@ -655,7 +657,7 @@
 matplotlib.rc('font', size=nsize)
 ylab = r"$\sqrt{u^{\prime 2}_\mathrm{max}}/u_{\infty}$"
 ax.set_ylabel(ylab, fontsize=tsize)
-ax.set_xlabel(r"$x/\delta_0$", fontsize=tsize)
+ax.set_xlabel(r"$x/l_r$", fontsize=tsize)
 ax.plot(xynew['x'], xynew['<u`u`>'], 'k')
 ax.set_yscale('log')
 ax.set_xlim([-100, 20])
@@ -699,29 +701,6 @@
 plt.savefig(
     pathF + "ModeProf.svg", bbox_inches="tight", pad_inches=0.1
 )
-# %%
-ampl = 0.001
-x = 0.0
-z = 0.0
-t = 804
-omeg_r = 0.443259
-alpha = 0.633513
-beta = 1.147607
-theta = alpha*x+beta*z-omeg_r*t
-q_perturb = ampl*(qr*np.cos(theta)-qi*np.sin(theta))
-fig, ax = plt.subplots()
-ax.plot(q_perturb, y, 'k', linewidth=1.5)
-# ax.set_xlabel (r'$p^{\prime}/(\rho_{\infty} u_{\infty}^{2})$', fontdict = font3)
-ax.set_xlabel(r'$u^{\prime}/u_{\infty}$', fontdict=font3)
-ax.set_ylabel(r'$y/\delta_0$', fontdict=font3)
-# ax.set_xlim ([-0.00002, 0.00018])
-ax.set_ylim([0, 5])
-# ax.xaxis.set_major_locator (xmajorLocator)
-# ax.xaxis.set_minor_locator (xminorLocator)
-ax.grid(b=True, which='both', linestyle=':')
-plt.tight_layout(pad=0.5, w_pad=0.2, h_pad=0.2)
-plt.savefig(path + 'PertUProfile.pdf', dpi=300)
-plt.show()
 
 # %% animation for vortex structure
 pathF = path + "Figures/"