#!/usr/bin/python

# $Id: presult.py,v 1.4 2009/04/23 08:37:12 dintrans Exp $
# Effective wavenumbers for first and second derivatives (Section 9.3.1 in 
# Brandenburg (2003), Computational aspects of astrophysical MHD and turbulence 
# http://arxiv.org/abs/astro-ph/0109497

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

dim=pc.read_dim()
# Box size Lx = 2pi by default --> kNy=pi/dx=pi/(Lx/Nx)=Nx/2
kNy=dim.nx/2

dat=N.loadtxt('result.dat',comments='-')

k=dat[:,0]
abm=dat[:,1]
ajm=dat[:,2]
brms=dat[:,3]
arms=dat[:,4]

keff=abm/arms**2
keff2=ajm/arms**2

k=k/kNy
keff=keff/kNy
keff2=keff2/kNy**2

P.figure(figsize=(10,5))
P.subplot(121)
P.plot(k,k,label='theory')
P.plot(k,k,'o')
P.plot(k,keff,label='numerical')
P.plot(k,keff,'o')
P.xlabel(r'$k/k_{\rm{Ny}}$')
P.ylabel(r'$k_{\rm{eff}}/k_{\rm{Ny}}$')
P.title('First derivative')
P.legend(loc='best')

P.subplot(122)
P.plot(k,k**2,label='theory')
P.plot(k,k**2,'o')
P.plot(k,keff2,label='numerical')
P.plot(k,keff2,'o')
P.xlabel(r'$k/k_{\rm{Ny}}$')
P.ylabel(r'$(k_{\rm{eff}}/k_{\rm{Ny}})^2$')
P.title('Second derivative')
P.legend(loc='best')

P.subplots_adjust(wspace=0.3, left=0.1, right=0.98, top=0.9)
P.show()