#!/usr/bin/env python

#
# $Id: p_power.py,v 1.2 2009/04/29 09:29:49 dintrans Exp $
# Plot the power spectra
#

import numpy as N
import pylab as P
import pencil as pc

P.figure(figsize=(8, 8))
# kinetic power spectra
t, data = pc.read_power('data/power_kin.dat')
nt=len(t)
P.subplot(221)
for i in range(nt):
  P.semilogy(data[i, :])
P.semilogy(data[0, :], 'o', label='t=%.2f'%t[0])
P.semilogy(data[nt-1, :], 'o', label='t=%.2f'%t[nt-1])
P.legend(loc='best')
P.title('Kinetic energy spectra')
P.xlabel('k')

# magnetic spectra
t, data = pc.read_power('data/power_mag.dat')
P.subplot(222)
for i in range(nt):
  P.semilogy(data[i, :])
P.semilogy(data[0, :], 'o', label='t=%.2f'%t[0])
P.semilogy(data[nt-1, :], 'o', label='t=%.2f'%t[nt-1])
P.legend(loc='best')
P.title('Magnetic energy spectra')
P.xlabel('k')

# kinetic helicity spectra
t, data = pc.read_power('data/powerhel_kin.dat')
P.subplot(223)
for i in range(nt):
  P.semilogy(data[i, :])
P.semilogy(data[0, :], 'o', label='t=%.2f'%t[0])
P.semilogy(data[nt-1, :], 'o', label='t=%.2f'%t[nt-1])
P.legend(loc='best')
P.title('Kinetic helicity spectra')
P.xlabel('k')

# magnetic helicity spectra
t, data = pc.read_power('data/powerhel_mag.dat')
P.subplot(224)
for i in range(nt):
  P.semilogy(data[i, :])
P.semilogy(data[0, :], 'o', label='t=%.2f'%t[0])
P.semilogy(data[nt-1, :], 'o', label='t=%.2f'%t[nt-1])
P.legend(loc='best')
P.title('Magnetic helicity spectra')
P.xlabel('k')

P.subplots_adjust(hspace=0.25, bottom=0.05, top=0.95, left=0.1, right=0.95)
P.show()