from pycrtools import *
import matplotlib.pyplot as plt

plt.figure()

filename_lofar = LOFARSOFT + "/data/lofar/VHECR_example.h5"

datafile = open(filename_lofar)

fftdata = datafile.empty("FFT_DATA")

nblocks = datafile["MAXIMUM_READ_LENGTH"] / datafile["BLOCKSIZE"]

avspectrum = hArray(float, dimensions=fftdata, name="Average spectrum")

# Calculate spectral power
for block in range(nblocks):
    datafile["BLOCK"] = block
    fftdata.read(datafile, "FFT_DATA")
    hSpectralPower(avspectrum[...], fftdata[...])

frequencies = datafile["FREQUENCY_DATA"].setUnit("M", "")

# Plotting
avspectrum.par.xvalues = frequencies
avspectrum.par.title = "Average spectrum"
avspectrum[0].plot(logplot="y", clf=True)