#!/usr/bin/python

# $Id: highwavenumber.py,v 1.1 2009/04/20 12:15:32 dintrans Exp $
# Compute Fig. 9.1 in Axel's computational paper (2003)
# http://arxiv.org/abs/astro-ph/0109497
# [only the 2nd, 4 and 6-th are plotted here]

import numpy as N
import pylab as P

k=N.linspace(0.01,N.pi,100)
pi=N.pi

#
#  first derivative
#
P.subplot(211)
P.plot(k/pi,N.sin(k)/pi,label='2nd')
P.plot(k/pi,(-N.sin(2*k)+8*N.sin(k))/(6.*pi),label='4th')
P.plot(k/pi,(N.sin(3*k)-9*N.sin(2*k)+45*N.sin(k))/(30.*pi),label='6th')
P.xlabel(r'$k/k_{\rm{Ny}}$')
P.ylabel(r'$k_{\rm{eff}}/k_{\rm{Ny}}$')
P.legend(loc='best')
P.title('First derivative')

#
#  second derivative
#
P.subplot(212)
P.plot(k/pi,(2.-2*N.cos(k))/pi**2,label='2nd')
P.plot(k/pi,(15.-16*N.cos(k)+N.cos(2*k))/(6.*pi**2),label='4th')
P.plot(k/pi,(245.-270.*N.cos(k)+27.*N.cos(2*k)-2.*N.cos(3*k))/(90.*pi**2),
label='6th')
P.xlabel(r'$k/k_{\rm{Ny}}$')
P.ylabel(r'$(k_{\rm{eff}}/k_{\rm{Ny}})^2$')
P.legend(loc='best')
P.title('Second derivative')

P.subplots_adjust(hspace=0.3,wspace=0.4,left=0.1,right=0.98,top=0.95)
P.show()