#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  7 12:41:49 2019

@author: Alexande Knebe

Sample code to read the ELG files

"""

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

#==============================================================================
# filename to be read
#==============================================================================
file = '/Users/aknebe/Office/DATA_local/Projects/UNITSIM/SAMs/SAGE/ELGs/UNITSIM1/UNITSIM1_model_z0.987_ELGs.h5'

#==============================================================================
# wrapper for reading dataset from hdf5 file
#==============================================================================
def get_h5dataset(hf, data):
    X = hf.get(data)
    X = np.array(X)
    return X

#==============================================================================
# main()
#==============================================================================
hf = h5py.File(file, 'r')
print('Available entries in file:')
print(hf.keys())
print('')

# the galaxy
Xpos             = get_h5dataset(hf, 'Xpos')                 # [Mpc/h]
Ypos             = get_h5dataset(hf, 'Ypos')                 # [Mpc/h]
Zpos             = get_h5dataset(hf, 'Zpos')                 # [Mpc/h]
Xvel             = get_h5dataset(hf, 'Xvel')                 # [km/s]
Yvel             = get_h5dataset(hf, 'Yvel')                 # [km/s]
Zvel             = get_h5dataset(hf, 'Zvel')                 # [km/s]
Mstar            = get_h5dataset(hf, 'Mstar')                # [Msun/h]
Mhalo            = get_h5dataset(hf, 'Mhalo')                # [Msun/h]
HostHaloID       = get_h5dataset(hf, 'HostHaloID')           # check Fig.A1 of Knebe et al. (2018) for definition
MainHaloID       = get_h5dataset(hf, 'MainHaloID')           # check Fig.A1 of Knebe et al. (2018) for definition
MainMhalo        = get_h5dataset(hf, 'MainMhalo')            # [Msun/h], mass of 'MainHalo' 

# the emission lines
logLHalpha       = get_h5dataset(hf, 'logLHalpha')           # [erg/s]
logLHalpha_att   = get_h5dataset(hf, 'logLHalpha_att')       # [erg/s]

logFHalpha       = get_h5dataset(hf, 'logFHalpha')           # [erg/s/cm^2]
logFHalpha_att   = get_h5dataset(hf, 'logFHalpha_att')       # [erg/s/cm^2]


#Property       = hf['Property'][:]  # this should also work to get the 'Property' out of hf[]

hf.close()



################################################################################
# check Halpha luminosity function
min = 0
max = 100

plt.figure(1)
plt.hist(logLHalpha[np.where((logLHalpha>min) & (logLHalpha<max))[0]],label='Halpha')
plt.hist(logLHalpha_att[np.where((logLHalpha_att>min) & (logLHalpha_att<max))[0]],label='Halpha_att')
plt.xlabel('log($L_{Halpha}$)')
plt.ylabel('N')
plt.legend(loc='upper left')
plt.savefig('Halpha-LF.pdf')