#! /usr/bin/env python3
#-------------------------------*- python -*-----------------------------------
# ==  == ====== ====   ====    |
#                   \\     ||  | Multiphase Code Repository by HZDR
# ======   //   ||  || ===//   | Website: https://doi.org/10.14278/rodare.767
# ||  ||  //    ||  // || \\   | License: GPL-3.0-or-later
# ==  == ====== ====   ==  ==  |
#------------------------------------------------------------------------------
# License
#     This file is part of the Multiphase Code Repository by HZDR.
#
#     Copyright (C) 2025 by Helmholtz-Zentrum Dresden-Rossendorf e.V. (HZDR),
#     Website: https://hzdr.de
#
#     Multiphase Code Repository by HZDR is based on the free software for
#     computational fluid dynamics (CFD) from the OpenFOAM Foundation.
#     Copyright (C) 2025 by OpenFOAM Foundation, Website: https://openfoam.org
#
#     If you are interested in which files are original OpenFOAM Foundation
#     files, which OpenFOAM Foundation files were modified, and which files were
#     newly created, see FILES.md.
#
#     Multiphase Code Repository by HZDR is free software: you can redistribute
#     it and/or modify it under the terms of the GNU General Public License as
#     published by the Free Software Foundation, either version 3 of the
#     License, or (at your option) any later version.
#
#     Multiphase Code Repository by HZDR is distributed in the hope that it will
#     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
#     Public License for more details.
#
#     You should have received a copy of the GNU General Public License along
#     with the Multiphase Code Repository by HZDR. If not, see
#     <http://www.gnu.org/licenses/>.
#
#------------------------------------------------------------------------------

""" Plot OpenFOAM patch surface data from cylinder as 2D image via heatmap """

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import NearestNDInterpolator

# Load data
film_thickness = np.loadtxt('../postProcessing/filmThickness/0/patch.xy')
x_values = film_thickness[:,0]
y_values = film_thickness[:,1]
z_values = film_thickness[:,2]
film_thickness_values = film_thickness[:,3]

# Calculate angle
angle_values = np.arctan2(-y_values, x_values)/np.pi + 2.5

# Closed-loop condition
angle_values = np.append(angle_values, angle_values-2)
z_values = np.append(z_values, z_values)
film_thickness_values = np.append(film_thickness_values, film_thickness_values)

# Map to 2D grid
X = np.linspace(min(z_values), max(z_values), num = 100)
Y = np.linspace(min(angle_values), max(angle_values), num = 100)
X, Y = np.meshgrid(X, Y) # 2D grid for interpolation
interp = NearestNDInterpolator(list(zip(z_values, angle_values)),
                              film_thickness_values)
Z = interp(X, Y)*1000 #[mm]

# Plot
fig, ax = plt.subplots()
plt.pcolormesh(X, Y, Z, shading='auto', cmap = 'viridis', vmin=0)
cb1 = plt.colorbar(orientation='horizontal')
cb1.set_label(r'$h_\mathrm{f}$ [mm]')
plt.title('Film thickness distribution', fontsize=12)
plt.xlabel('z [m]')
plt.ylabel(r'$\varphi$ [$\pi$]')
plt.ylim(0, 2)
plt.savefig('filmThickness')
