import numpy as np
from sys import argv
import matplotlib.pyplot as plt

z_dipoles = []
with open("ch2o_pbe_rtp-output-moments.dat") as f:
#with open(argv[1]) as f:
    lines = f.readlines()
    
    max_value = 0
    for line in lines:
        line = line.split(' ')
        line = filter(None, line)
        flip = 0
        for string in line:
            if string == 'Z=':
                flip = 1
                continue
            if flip == 1:
                z_dipoles.append(float(string))
                break

signal = np.array(z_dipoles, dtype=float)
fourier = np.fft.fft(signal)
n = signal.size
timestep = np.float64(0.025 * 10E-15) #2.419E-17
freq = np.fft.fftfreq(n, d=timestep) * np.float64(6.626E-34 * 6.24E18)


print(len(z_dipoles))
plt.plot(z_dipoles)
plt.show()
#y = range(len(z_dipoles))
#x = numpy.array(y, z_dipoles)
x = np.fft.rfft(z_dipoles)
plt.plot(freq[0:len(x[1:])],abs(x[1:]))
plt.show()