import re
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.ticker as ticker
from pyfluids import Fluid, FluidsList, Input

from pathlib import Path

plt.rcParams.update({'font.size': 30})
plt.rcParams.update({'mathtext.fontset':'cm'})
plt.rcParams.update({'font.family':'Times New Roman'})

def make_dict(paths, Dh, Pr, density, nu, Cp, k, H, A, h_area, Twall, solver):

    def find_nearest(array, value):
        array = np.asarray(array)
        idx = (np.abs(array - value)).argmin()
        return array[idx]

    plot_dict = {}

    for i, path in enumerate(paths):
        case = path.stem
        plot_dict[case] = {'periodic':{}}
        # analytical = TemperatureField(1.0, 0.5, 0.7, float(case[2::]), nu, True, 8)    
        df = pd.read_csv(path)
        x_coords = df['Points_0'].unique()
        min_x = np.min(x_coords)
        max_x = np.max(x_coords)
        plot_dict[case]['periodic']['L'] = max_x-min_x
        mid_x = find_nearest(x_coords, plot_dict[case]['periodic']['L']/2)
        
        inlet_df = df[df['Points_0']==min_x]
        outlet_df = df[df['Points_0']==max_x]
        midway_df = df[df['Points_0']==mid_x]

        plot_dict[case]['periodic']['x'] = df['Points_0'].values-min_x
        plot_dict[case]['periodic']['y'] = df['Points_1'].values   
        plot_dict[case]['periodic']['yI'] = inlet_df['Points_1'].values
        plot_dict[case]['periodic']['yM'] = midway_df['Points_1'].values
        plot_dict[case]['periodic']['yO'] = outlet_df['Points_1'].values
        
        plot_dict[case]['periodic']['Ux'] = df['UMean_0'].values
        plot_dict[case]['periodic']['Umean'] = np.mean(plot_dict[case]['periodic']['Ux'])
        plot_dict[case]['periodic']['UI'] = inlet_df['UMean_0'].values
        plot_dict[case]['periodic']['UImean'] = np.mean(plot_dict[case]['periodic']['UI'])
        plot_dict[case]['periodic']['UM'] = midway_df['UMean_0'].values
        plot_dict[case]['periodic']['UO'] = outlet_df['UMean_0'].values

        plot_dict[case]['periodic']['m_dot'] = density*A*plot_dict[case]['periodic']['UImean']
        
        plot_dict[case]['periodic']['T'] = df['TMean'].values
        plot_dict[case]['periodic']['TI'] = inlet_df['TMean'].values
        plot_dict[case]['periodic']['TM'] = midway_df['TMean'].values
        plot_dict[case]['periodic']['TO'] = outlet_df['TMean'].values

        #---------------------- h for periodic ------------------------------------------------
        periodic_TI = (np.sum(np.multiply(plot_dict[case]['periodic']['TI'], plot_dict[case]['periodic']['UI']))/
                       np.sum(plot_dict[case]['periodic']['UI']))
        periodic_TO = (np.sum(np.multiply(plot_dict[case]['periodic']['TO'], plot_dict[case]['periodic']['UO']))/
                       np.sum(plot_dict[case]['periodic']['UO']))

        temp_lm = ((Twall - periodic_TO) - (Twall - periodic_TI))/(np.log((Twall - periodic_TO)/(Twall-periodic_TI)))
        
        plot_dict[case]['periodic']['h'] = np.abs(plot_dict[case]['periodic']['m_dot']*Cp*(periodic_TI-periodic_TO)/(h_area*temp_lm))
        plot_dict[case]['periodic']['Nu'] = plot_dict[case]['periodic']['h']*2*Dh/k

        #plot_dict[case]['periodic']['Pdrop'] = np.abs(np.mean(inlet_df['pMean'].values - outlet_df['pMean'].values))*density
        plot_dict[case]['periodic']['Pdrop'] = np.abs(np.mean(solver[case]['pg'][-1000:-1]))*density

        plot_dict[case]['periodic']['f'] = ((plot_dict[case]['periodic']['Pdrop']*Dh)/
                                            (plot_dict[case]['periodic']['L']*2*density*plot_dict[case]['periodic']['Umean']**2))

        plot_dict[case]['periodic']['Re'] = (Dh*plot_dict[case]['periodic']['Umean'])/nu
    
    return plot_dict

def loop_plots_T(plots, H, colormap, Twall, dx, dy):
    fig_x, fig_y = 12, 10
    figs = {}
    axs = {}

    H2 = H/2

    for idx, case in enumerate(plots):
        figs[case], axs[case] = plt.subplots(figsize=(fig_x, fig_y), sharex=True)
        plt.subplots_adjust(hspace=0.1)
        figs[case].tight_layout()
        figs[case].subplots_adjust(top=0.945)
        figs[case].suptitle('Re ' + str(np.rint(plot_dict[case]['periodic']['Re'])))
        
        periodic_T = (plots[case]['periodic']['T']-Twall)/np.mean(plots[case]['periodic']['TI']-Twall)
        
        cf3 = axs[case].tricontour(plots[case]['periodic']['x']/plots[case]['periodic']['L'], 
                                   plots[case]['periodic']['y']/H2, periodic_T, levels=20, cmap=colormap)
        
        axs[case].plot(dx, dy/H, 'k', linewidth=3)
        axs[case].plot(dx, -dy/H, 'k', linewidth=3)        
        
        axs[case].clabel(cf3, inline=True, fontsize=22)

        axs[case].set_xlabel('x/L')
        axs[case].set_ylabel('$\mathrm{y/H_{avg}}$')
    
    return figs, axs

def loop_plots_Ux(plots, H, colormap, dx, dy):
    fig_x, fig_y = 12, 10
    figs = {}
    axs = {}

    H2 = H/2

    for idx, case in enumerate(plots):
        figs[case], axs[case] = plt.subplots(figsize=(fig_x, fig_y), sharex=True)
        plt.subplots_adjust(hspace=0.1)
        figs[case].tight_layout()
        figs[case].subplots_adjust(top=0.945)
        figs[case].suptitle('Re ' + str(np.rint(plot_dict[case]['periodic']['Re'])))
        
        cf3 = axs[case].tricontour(plots[case]['periodic']['x']/plots[case]['periodic']['L'], plots[case]['periodic']['y']/H2, 
                                    plots[case]['periodic']['Ux']/plots[case]['periodic']['UImean'], levels=20, cmap=colormap)
        
        axs[case].plot(dx, dy/H, 'k', linewidth=3)
        axs[case].plot(dx, -dy/H, 'k', linewidth=3)
        
        axs[case].clabel(cf3, inline=True, fontsize=22)

        axs[case].set_xlabel('x/L')
        axs[case].set_ylabel('y/H')
    
    return figs, axs

def loop_plots_all(plots, H, Tw, colormap, dx, dy, bahaidarah_velocity, reference_dict):
    fig_x, fig_y = 20, 12
    figs = {}
    axs = {}

    H2 = H/2

    for idx, case in enumerate(plots):
        figs[case], axs[case] = plt.subplots(1,3, figsize=(fig_x, fig_y), width_ratios=[2,1,1], sharey=True)
        plt.subplots_adjust(hspace=0.1)
        figs[case].tight_layout()
        figs[case].subplots_adjust(top=0.945)
        figs[case].suptitle('Re ' + str(np.rint(plot_dict[case]['periodic']['Re'])))
        
        periodic_T = plots[case]['periodic']['T']-Tw
        periodic_TI = plots[case]['periodic']['TI']-Tw
        periodic_TO = plots[case]['periodic']['TO']-Tw
        periodic_TM = plots[case]['periodic']['TM']-Tw
        
        cf2 = axs[case][0].tricontour(plots[case]['periodic']['x']/plots[case]['periodic']['L'], plots[case]['periodic']['y']/H2, 
                                     (periodic_T/np.mean(periodic_TI)), levels=20, cmap=colormap)
        
        axs[case][0].plot(dx, dy/H, 'k', linewidth=3)
        axs[case][0].plot(dx, -dy/H, 'k', linewidth=3)

        periodic_UI = plots[case]['periodic']['UI']
        periodic_UO = plots[case]['periodic']['UO']
        periodic_UM = plots[case]['periodic']['UM']

        step = 3
        #---------------- Velocity Profile ------------------------------------
        axs[case][1].plot((periodic_UI[::step]/plots[case]['periodic']['UImean']), plots[case]['periodic']['yI'][::step]/H2, 'o', color='C0', label = 'x=0')
        axs[case][1].plot((periodic_UM[::step]/plots[case]['periodic']['UImean']), plots[case]['periodic']['yM'][::step]/H2, '^', color='C1', label='x=L/2')
        
        if np.rint(plot_dict[case]['periodic']['Re']) == 25.0:
            if flow_type == 'arc':
                b_y = bahaidarah_velocity['y_25'].values
                b_u = bahaidarah_velocity['u_25'].values
                axs[case][1].plot(b_u, (b_y*3)-1.5, color='C4', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']')
            else:
                b_y = bahaidarah_velocity['y_25'].values
                b_u = bahaidarah_velocity['u_25'].values
                axs[case][1].plot(b_u, (b_y*3)-1.5, color='C4', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']')
        elif np.rint(plot_dict[case]['periodic']['Re']) == 100.0:
            if flow_type == 'arc':
                b_y = bahaidarah_velocity['y_100'].values
                b_u = bahaidarah_velocity['u_100'].values
                axs[case][1].plot(b_u, (b_y*3)-1.5, color='C4', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']')
            else:
                b_y = bahaidarah_velocity['y_100'].values
                b_u = bahaidarah_velocity['u_100'].values
                axs[case][1].plot(b_u, (b_y*3)-1.5, color='C4', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']')
        elif np.rint(plot_dict[case]['periodic']['Re']) == 400.0:
            if flow_type == 'arc':
                b_y = bahaidarah_velocity['y_400'].values
                b_u = bahaidarah_velocity['u_400'].values
                axs[case][1].plot(b_u, (b_y*3)-1.5, color='C4', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']')
            else:
                b_y = bahaidarah_velocity['y_400'].values
                b_u = bahaidarah_velocity['u_400'].values
                axs[case][1].plot(b_u, (b_y*3)-1.5, color='C4', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']')

        axs[case][1].plot((periodic_UO[::step]/plots[case]['periodic']['UImean']), plots[case]['periodic']['yI'][::step]/H2, 'd', color='C2', label='x=L')

        #------------------ Temperature Profile -----------------------------------------
        axs[case][2].plot((periodic_TI[::step]/np.mean(periodic_TI)), plots[case]['periodic']['yI'][::step]/H2, 'o', color='C0', label='x=0')
        axs[case][2].plot((periodic_TM[::step]/np.mean(periodic_TI)), plots[case]['periodic']['yM'][::step]/H2, '^', color='C1', label='x=L/2')
        axs[case][2].plot((periodic_TO[::step]/np.mean(periodic_TI)), plots[case]['periodic']['yI'][::step]/H2, 'd', color='C2', label='x=L')

        axs[case][0].clabel(cf2, inline=True, fontsize=22)

        axs[case][0].set_xlabel('x/L')
        axs[case][0].set_ylabel('$y/H_{avg}$')

        axs[case][1].set_xlabel('$U/U_0$')

        x_ticks = np.arange(0, 2.5, 0.5)
        # Rewrite the y labels
        axs[case][1].set_xticks(x_ticks)
        axs[case][1].xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))

        axs[case][2].set_xlabel('$T/T_0$')
        # Rewrite the y labels
        axs[case][2].set_xticks(x_ticks)
        axs[case][2].xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))

        axs[case][1].legend(loc=1)
        axs[case][2].legend(loc=1)
    
    return figs, axs

def search_str(file_path):
    time = []
    pg = []
    with open(file_path, 'r') as f:
        for line in f:
            if re.match(r'Time = ', line):
                time.append(float(line.split()[2]))
            elif re.match(r'Pressure gradient', line):
                pg.append(float(line.split()[-1]))
    
    return np.array(time), np.array(pg)

flow_types = ['arc', 'sine']
solver_type = ['steady', 'unsteady']
fluid_types = ['air','water']

for flow_type in flow_types:
    for subcase in solver_type:
        for fluid in fluid_types:

            data_directory = Path(r'./').joinpath(flow_type, 'Re_runs_'+ subcase, fluid)
            periodic_paths = list(data_directory.glob('*.csv'))
            solver_files = sorted(list(data_directory.rglob('log.cyclic*Foam')))

            save_path = data_directory

            # grab the pressure gradient
            solver = {}
            for file in solver_files:
                case = file.parent.stem
                solver[case] = {}
                solver[case]['time'], solver[case]['pg'] = search_str(file)
            
            if fluid == 'air':
                working_fluid = Fluid(FluidsList.Air).with_state(
                    Input.pressure(101325), Input.temperature(35)
                )
            else:
                working_fluid = Fluid(FluidsList.Water).with_state(
                    Input.pressure(101325), Input.temperature(35)
                )

            Twall = 293.15

            plot_dict = {}

            L = 1.0
            if flow_type == 'arc':
                Hmax = 0.7138
                Hmin = 0.214286
                Havg = Hmax/2 + Hmin/2
                Dh = Havg
                h_area = 0.0455905
                A = 0.00428571
            else:
                Hmax = 0.7134
                Hmin = 0.214286
                Havg = Hmax/2 + Hmin/2
                Dh = Havg
                h_area = 0.0455769
                A = 0.00428571

            Pr = working_fluid.prandtl
            density = working_fluid.density
            mu = working_fluid.dynamic_viscosity
            nu = working_fluid.dynamic_viscosity/working_fluid.density 
            Cp = working_fluid.specific_heat
            alpha = nu/Pr
            k = Cp*mu/Pr

            plot_dict = make_dict(periodic_paths, Dh, Pr, density, nu, Cp, k, 1.0, A, h_area, Twall, solver)

            validation_data_path = Path(r'./wavy_channel_data.xlsx')

            bahaidarah = pd.read_excel(validation_data_path, sheet_name='Bahaidarah')
            bahaidarah_velocity_arc = pd.read_excel(validation_data_path, sheet_name='Bahaidarah_velocity_arc')
            bahaidarah_velocity_sine = pd.read_excel(validation_data_path, sheet_name='Bahaidarah_velocity_sine')
            niceno = pd.read_excel(validation_data_path, sheet_name='Niceno')
            wang_vanaka = pd.read_excel(validation_data_path, sheet_name='Wang_Vanaka')
            ramgadia = pd.read_excel(validation_data_path, sheet_name='Ramgadia')
            harikrishnan = pd.read_excel(validation_data_path, sheet_name='Harikrishnan_Tiwari')
            nishimura = pd.read_excel(validation_data_path, sheet_name='Nishimura')
            
            reference_dict = {'bahaidarah':21, 'niceno':24,'wang_vanaka':23,'ramgadia':22,'harikrishnan':3,'nishimura':25}

            from matplotlib.lines import Line2D

            markers = ['o', '*', '^', 'd', 'v', '>', '<']

            fig1, ax1 = plt.subplots(figsize=(12, 10), facecolor='w')

            if fluid == 'air':
                if flow_type == 'arc':
                        b_Re = bahaidarah['Re'].values
                        b_nu = (bahaidarah['arc_Nu'].values*k)/Dh
                        n_Re = niceno['Re_Nu_arc'].values
                        n_nu = (niceno['Nu_arc'].values*k)/Dh
                else:
                        b_Re = bahaidarah['Re'].values
                        b_nu = (bahaidarah['sine_Nu'].values*k)/Dh
                        n_Re = niceno['Re_Nu_sine'].values
                        n_nu = (niceno['Nu_sine'].values*k)/Dh
                        w_Re = wang_vanaka['Re_Nu'].values
                        w_nu = (wang_vanaka['Nu'].values*k)/Dh
                        r_Re = ramgadia['Re'].values
                        r_nu = (ramgadia['Nu'].values*k)/Dh
                        h_Re = harikrishnan['Re'].values
                        h_nu = (harikrishnan['Nu'].values*k)/Dh

            sim_Re = np.zeros(len(plot_dict))
            sim_h = np.zeros(len(plot_dict))

            for i, case in enumerate(plot_dict):
                sim_Re[i] = plot_dict[case]['periodic']['Re']
                sim_h[i] = np.abs(plot_dict[case]['periodic']['h'])

            inds = sim_Re.argsort()

            ax1.plot(sim_Re[inds], sim_h[inds], marker=markers[0], markersize=10, color='C0')

            if fluid == 'air':
                ax1.scatter(b_Re, b_nu, marker=markers[1], color='C1')
                ax1.scatter(n_Re, n_nu, marker=markers[2], color='C2')

                legend_elements = [Line2D([0], [0], marker=markers[0], color='w', label='Present Work',
                                            markerfacecolor='C0', markersize=12),
                                        Line2D([0], [0], marker=markers[1], color='w', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']',
                                            markerfacecolor='C1', markersize=12),
                                        Line2D([0], [0], marker=markers[2], color='w', label='Ničeno and Nobile ' + '[' + str(reference_dict['niceno']) + ']',
                                            markerfacecolor='C2', markersize=12)]

                if flow_type == 'sine':
                    ax1.scatter(w_Re, w_nu, marker=markers[3], color='C3')
                    ax1.scatter(r_Re, r_nu, marker=markers[4], color='C4')
                    ax1.scatter(h_Re, h_nu, marker=markers[5], color='C5')
                    
                    legend_elements = [Line2D([0], [0], marker=markers[0], color='w', label='Present Work',
                                        markerfacecolor='C0', markersize=12),
                                    Line2D([0], [0], marker=markers[1], color='w', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']',
                                        markerfacecolor='C1', markersize=12),
                                    Line2D([0], [0], marker=markers[2], color='w', label='Ničeno and Nobile ' + '[' + str(reference_dict['niceno']) + ']',
                                        markerfacecolor='C2', markersize=12),
                                    Line2D([0], [0], marker=markers[3], color='w', label='Wang and Vanka ' + '[' + str(reference_dict['wang_vanaka']) + ']',
                                        markerfacecolor='C3', markersize=12),
                                    Line2D([0], [0], marker=markers[4], color='w', label='Ramgadia and Saha ' + '[' + str(reference_dict['ramgadia']) + ']',
                                        markerfacecolor='C4', markersize=12),
                                    Line2D([0], [0], marker=markers[5], color='w', label='Harikrishnan and Tiwari ' + '[' + str(reference_dict['harikrishnan']) + ']',
                                        markerfacecolor='C5', markersize=12)]

            # plt.xticks(rotation=45)
            ax1.set_xlabel('Re')
            ax1.set_ylabel('Local Heat Transfer Coefficient [W/mK]')
            ax1.set_yscale('log')
            ax1.set_xscale('log')
            if fluid == 'air':
                if flow_type == 'sine':
                        ax1.legend(handles=legend_elements, loc=1, ncol=2, bbox_to_anchor=(1.01, 1.20))
                else:
                        ax1.legend(handles=legend_elements, loc=1, ncol=2, bbox_to_anchor=(1.01, 1.15))
            #fig1.savefig(save_path.joinpath(flow_type + '_' + subcase + '_' + fluid + '_ht_coefficient_steady.pdf'), facecolor='w', dpi=300, bbox_inches='tight')
            plt.close()

            fig2, ax2 = plt.subplots(figsize=(12, 10), facecolor='w')

            if flow_type == 'arc':
                b_Re = bahaidarah['Re'].values
                b_nu = bahaidarah['arc_Nu'].values
                n_Re = niceno['Re_Nu_arc'].values
                n_nu = niceno['Nu_arc'].values
            else:
                b_Re = bahaidarah['Re'].values
                b_nu = bahaidarah['sine_Nu'].values
                n_Re = niceno['Re_Nu_sine'].values
                n_nu = niceno['Nu_sine'].values
                w_Re = wang_vanaka['Re_Nu'].values
                w_nu = wang_vanaka['Nu'].values
                r_Re = ramgadia['Re'].values
                r_nu = ramgadia['Nu'].values
                h_Re = harikrishnan['Re'].values
                h_nu = harikrishnan['Nu'].values


            sim_Re = np.zeros(len(plot_dict))
            sim_nu = np.zeros(len(plot_dict))

            for i, case in enumerate(plot_dict):
                sim_Re[i] = plot_dict[case]['periodic']['Re']
                sim_nu[i] = np.abs(plot_dict[case]['periodic']['Nu'])

            inds = sim_Re.argsort()

            ax2.plot(sim_Re[inds], sim_nu[inds], marker=markers[0], markersize=10, color='C0')

            if fluid == 'air':
                ax2.scatter(b_Re, b_nu, marker=markers[1], color='C1')
                ax2.scatter(n_Re, n_nu, marker=markers[2], color='C2')


                legend_elements = [Line2D([0], [0], marker=markers[0], color='w', label='Present Work',
                                            markerfacecolor='C0', markersize=12),
                                        Line2D([0], [0], marker=markers[1], color='w', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']',
                                            markerfacecolor='C1', markersize=12),
                                        Line2D([0], [0], marker=markers[2], color='w', label='Ničeno and Nobile ' + '[' + str(reference_dict['niceno']) + ']',
                                            markerfacecolor='C2', markersize=12)]

                if flow_type == 'sine':
                    ax2.scatter(w_Re, w_nu, marker=markers[3], color='C3')
                    ax2.scatter(r_Re, r_nu, marker=markers[4], color='C4')
                    ax2.scatter(h_Re, h_nu, marker=markers[5], color='C5')
                    
                    legend_elements = [Line2D([0], [0], marker=markers[0], color='w', label='Present Work',
                                        markerfacecolor='C0', markersize=12),
                                    Line2D([0], [0], marker=markers[1], color='w', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']',
                                        markerfacecolor='C1', markersize=12),
                                    Line2D([0], [0], marker=markers[2], color='w', label='Ničeno and Nobile ' + '[' + str(reference_dict['niceno']) + ']',
                                        markerfacecolor='C2', markersize=12),
                                    Line2D([0], [0], marker=markers[3], color='w', label='Wang and Vanka ' + '[' + str(reference_dict['wang_vanaka']) + ']',
                                        markerfacecolor='C3', markersize=12),
                                    Line2D([0], [0], marker=markers[4], color='w', label='Ramgadia and Saha ' + '[' + str(reference_dict['ramgadia']) + ']',
                                        markerfacecolor='C4', markersize=12),
                                    Line2D([0], [0], marker=markers[5], color='w', label='Harikrishnan and Tiwari ' + '[' + str(reference_dict['harikrishnan']) + ']',
                                        markerfacecolor='C5', markersize=12)]

            dx = np.linspace(25, 400, 10)
            ax2.plot(dx, np.ones(len(dx))*7.54, 'k', linewidth=2, label='Plane Channel')
            
            # plt.xticks(rotation=45)
            ax2.set_xlabel('Re')
            ax2.set_ylabel('Nu')
            ax2.set_xscale('log')
            ax2.set_yscale('log')
            if fluid == 'air':
                if flow_type == 'sine':
                        ax2.legend(handles=legend_elements, loc=1, ncol=2, bbox_to_anchor=(1.01, 1.2))
                else:
                        ax2.legend(handles=legend_elements, loc=1, ncol=2, bbox_to_anchor=(1.01, 1.15))
            fig2.savefig(save_path.joinpath(flow_type + '_' + subcase + '_' + fluid + '_Nusselt.pdf'), facecolor='w', dpi=300, bbox_inches='tight')
            plt.close()

            fig3, ax3 = plt.subplots(figsize=(12, 10), facecolor='w')

            if flow_type == 'arc':
                b_Re = bahaidarah['Re'].values
                b_f = bahaidarah['arc_f'].values
                n_Re = niceno['Re_f_arc'].values
                n_f = niceno['f_arc'].values
            else:
                b_Re = bahaidarah['Re'].values
                b_f = bahaidarah['sine_f'].values
                n_Re = niceno['Re_f_sine'].values
                n_f = niceno['f_sine'].values
                w_Re = wang_vanaka['Re_f'].values
                w_f = wang_vanaka['f'].values
                r_Re = ramgadia['Re'].values
                r_f = ramgadia['f'].values
                h_Re = harikrishnan['Re'].values
                h_f = harikrishnan['f'].values
                nish_Re = nishimura['Re'].values
                nish_f = nishimura['f'].values

            sim_Re = np.zeros(len(plot_dict))
            sim_f = np.zeros(len(plot_dict))

            for i, case in enumerate(plot_dict):
                sim_Re[i] = plot_dict[case]['periodic']['Re']
                sim_f[i] = np.abs(plot_dict[case]['periodic']['f'])

            inds = sim_Re.argsort()

            ax3.plot(sim_Re[inds], sim_f[inds], marker=markers[0], markersize=10, color='C0')

            if fluid == 'air':
                ax3.scatter(b_Re, b_f, marker=markers[1], color='C1')
                ax3.scatter(n_Re, n_f, marker=markers[2], color='C2')

                legend_elements = [Line2D([0], [0], marker=markers[0], color='w', label='Present Work',
                                            markerfacecolor='C0', markersize=12),
                                        Line2D([0], [0], marker=markers[1], color='w', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']',
                                            markerfacecolor='C1', markersize=12),
                                        Line2D([0], [0], marker=markers[2], color='w', label='Ničeno and Nobile ' + '[' + str(reference_dict['niceno']) + ']',
                                            markerfacecolor='C2', markersize=12)]

                if flow_type == 'sine':
                        ax3.scatter(w_Re, w_f, marker=markers[3], color='C3')
                        ax3.scatter(r_Re, r_f, marker=markers[4], color='C4')
                        ax3.scatter(h_Re, h_f, marker=markers[5], color='C5')
                        ax3.scatter(nish_Re, nish_f, marker=markers[6], color='C6')
                        
                        legend_elements = [Line2D([0], [0], marker=markers[0], color='w', label='Present Work',
                                            markerfacecolor='C0', markersize=12),
                                        Line2D([0], [0], marker=markers[1], color='w', label='Bahaidarah et al. ' + '[' + str(reference_dict['bahaidarah']) + ']',
                                            markerfacecolor='C1', markersize=12),
                                        Line2D([0], [0], marker=markers[2], color='w', label='Ničeno and Nobile ' + '[' + str(reference_dict['niceno']) + ']',
                                            markerfacecolor='C2', markersize=12),
                                        Line2D([0], [0], marker=markers[3], color='w', label='Wang and Vanka ' + '[' + str(reference_dict['wang_vanaka']) + ']',
                                            markerfacecolor='C3', markersize=12),
                                        Line2D([0], [0], marker=markers[4], color='w', label='Ramgadia and Saha ' + '[' + str(reference_dict['ramgadia']) + ']',
                                            markerfacecolor='C4', markersize=12),
                                        Line2D([0], [0], marker=markers[5], color='w', label='Harikrishnan and Tiwari ' + '[' + str(reference_dict['harikrishnan']) + ']',
                                            markerfacecolor='C5', markersize=12),
                                        Line2D([0], [0], marker=markers[6], color='w', label='Nishimura et al. ' + '[' + str(reference_dict['nishimura']) + ']',
                                            markerfacecolor='C6', markersize=12)]
            
            dx = np.linspace(10, 400, 100)
            ax3.plot(dx, 24/dx, color='k', linewidth=2)
            
            # plt.xticks(rotation=45)
            ax3.set_xlabel('Re')
            ax3.set_ylabel('f')
            ax3.set_xscale('log')
            ax3.set_yscale('log')
            if fluid == 'air':
                if flow_type == 'sine':
                        ax3.legend(handles=legend_elements, loc=1, ncol=2, bbox_to_anchor=(1.01, 1.255))
                else:
                        ax3.legend(handles=legend_elements, loc=1, ncol=2, bbox_to_anchor=(1.01, 1.15))
            fig3.savefig(save_path.joinpath(flow_type + '_' + subcase + '_' + fluid + '_f.pdf'), facecolor='w', dpi=300, bbox_inches='tight')
            plt.close()

            # make outline of geometry
            Lx = plot_dict[case]['periodic']['L'] 
            x_top = np.linspace(0, Lx, 100)
            a = (Hmax - Hmin)/2

            if flow_type == 'sine':
                y_top = Hmin + 2*a*np.sin(np.pi*x_top/Lx)**2
            else:
                y_top = Hmin + 2*a*np.sin(np.pi*x_top/Lx)

            # Plot the temperature field
            # figs, axs = loop_plots_T(plot_dict, Havg, 'twilight', Twall, x_top, y_top)
            # for case in figs:
            #     figs[case].savefig(save_path.joinpath(flow_type + '_' + case + '_' + subcase + '_' + fluid + '_T.pdf'), facecolor='w', dpi=300, bbox_inches='tight')
            #     plt.close()

            # Plot the velocity field
            # figs, axs = loop_plots_Ux(plot_dict, Havg, 'twilight', x_top, y_top)
            # for case in figs:
            #     figs[case].savefig(save_path.joinpath(flow_type + '_' + case + '_' + subcase + '_' + fluid + '_Ux.pdf'), facecolor='w', dpi=300, bbox_inches='tight')
            #     plt.close()
            
            # Plot the three fields combined
            if flow_type == 'arc':
                figs, axs = loop_plots_all(plot_dict, Havg, Twall,'twilight', x_top, y_top, bahaidarah_velocity_arc, reference_dict)
            else:
                figs, axs = loop_plots_all(plot_dict, Havg, Twall,'twilight', x_top, y_top, bahaidarah_velocity_sine, reference_dict)
            for case in figs:
                figs[case].savefig(save_path.joinpath(flow_type + '_' + case + '_' + subcase + '_' + fluid + '_combined.pdf'), facecolor='w', dpi=300, bbox_inches='tight')
                plt.close()
            
