#################################################################################
#
# Tutorial file how to read in hdf5 data from SP-Wind Windfarm simulations,
# generated in the context of the TotalControl project.
#
# A thorough description of the simulations and data can be found at in the deliverable report
#
# author: Ishaan Sood
# date: 24 Feb 2020
# email: sood.ishaan26@gmail.com
#
#################################################################################

'''
Read and plot timeseries snapshots of velocity and temperature field data from hdf5 file at xy and xz planes between
time 905s and 4495s at intervals of 10 seconds

'''


import matplotlib.animation as animation
import numpy as np
import matplotlib.pyplot as plt
import h5py as hp

#======================================================================================================================#

#Simulation grid parameters for plotting

Nx = 1200
Ny = 1200
Nz = 225

Lx = 16000
Ly = 16000
Lz = 1500

xvec = np.linspace(0.0,Lx,Nx+1)
yvec = np.linspace(0.0,Ly,Ny+1)
zvec = np.linspace(0.0,Lz,Nz+1)

plt.close("all")

# open the file

ftser = hp.File('cross_sections_timeseries.h5', 'r')

###########################################################################################################################
#
# FILE: cross_sections_timeseries.h5
#
#   group: xz_plane
#       description:    Cross section velocity (and temperature) plane at top turbine row
#       subgroups: u, v, w, (th)
#           datasets:  t_905.0 (e.g. at t = 905 sec)
#
#   group: xy_plane
#       description:    Cross section velocity (and temperature) plane at hub height
#       subgroups: u, v, w, (th)
#           datasets:  t_905.0 (e.g. at t = 905 sec)
###########################################################################################################################

plane = 'xy_plane'
component = 'u'
t = 4005

inflow_group = ftser[plane]

setname = "t_{:04.1f}".format(t)

data = inflow_group[component][setname]

plt.figure()
if plane == 'xz_plane':
    plt.imshow(np.flipud(data),extent=[xvec.min(), 0.9*xvec.max(), zvec.min(), zvec.max()])
else: 
    plt.imshow(np.flipud(data),extent=[xvec.min(), 0.9*xvec.max(), yvec.min(), yvec.max()])
plt.savefig('xy_u_plane.png')
plt.show()