#!/bin/env python
"""
Instructions:

1.  Configure your environment for Python.  In bash on the HAO
    workstations, you would want to add this to your ~/.bashrc file:

    export PROCESSOR=`uname -p`
    if [ "$PROCESSOR" == "x86_64" ]; then
    export PYTHONPATH=$PYTHONPATH:/home/schmitt/opt/64/lib64/python2.4/site-packages
    export PYTHONPATH=$PYTHONPATH:/home/schmitt/opt/64/lib/python2.4/site-packages
    else
    export PYTHONPATH=$PYTHONPATH:/home/schmitt/opt/32/lib/python2.5/site-packages
    fi

2.  Make sure you have a data file.  This script assumes you have the
    file 'AGU_CMITs_mix_2006-12-14T08-00-00Z.hdf' in the same
    directory as this script.

3.  Execute  'python plot.py'
"""

from pyhdf.SD import SD, SDC
import pylab
import numpy

# Read the data
hdffile = SD('/hao/aim/schmitt/data/LTR-2_1_4/Agu2006/CMITs/AGU_CMITs_mix_2006-12-14T22-05-00Z.hdf',SDC.READ)
x=hdffile.select('Grid X').get()
y=hdffile.select('Grid Y').get()
psi_n=hdffile.select('Potential North [V]').get()
psi_s=hdffile.select('Potential South [V]').get()
fac_n=hdffile.select('FAC North [A/m^2]').get()
fac_s=hdffile.select('FAC South [A/m^2]').get()
sigmap_n=hdffile.select('Pedersen conductance North [S]').get()
sigmap_s=hdffile.select('Pedersen conductance South [S]').get()
sigmah_n=hdffile.select('Hall conductance North [S]').get()
sigmah_s=hdffile.select('Hall conductance South [S]').get()
energy_n=hdffile.select('Average energy North [keV]').get()
energy_s=hdffile.select('Average energy South [keV]').get()
flux_n=hdffile.select('Number flux North [1/cm^2 s]').get()
flux_s=hdffile.select('Number flux South [1/cm^2 s]').get()
x=hdffile.select('Grid X').get()
y=hdffile.select('Grid Y').get()
psi_n=hdffile.select('Potential North [V]').get()/1000.
psi_s=hdffile.select('Potential South [V]').get()/1000.
fac_n=hdffile.select('FAC North [A/m^2]').get()*1.e6
fac_s=hdffile.select('FAC South [A/m^2]').get()*1.e6
sigmap_n=hdffile.select('Pedersen conductance North [S]').get()
sigmap_s=hdffile.select('Pedersen conductance South [S]').get()
sigmah_n=hdffile.select('Hall conductance North [S]').get()
sigmah_s=hdffile.select('Hall conductance South [S]').get()
energy_n=hdffile.select('Average energy North [keV]').get()
energy_s=hdffile.select('Average energy South [keV]').get()
flux_n=hdffile.select('Number flux North [1/cm^2 s]').get()
flux_s=hdffile.select('Number flux South [1/cm^2 s]').get()

# Go from (x,y) to cylindrical (r, theta).  You could also do
# spherical coordinates with r=1... but you'll get the same results
# when projecting onto a 2d plane.
theta=numpy.arctan2(y,x)
theta[theta<0]=theta[theta<0]+2*numpy.pi
r=numpy.sqrt(x**2+y**2)

# Create a polar plot.
pylab.figure()
pylab.axes(polar=True)
contours = numpy.linspace(psi_n.min(),psi_n.max(),100)
pylab.contour(theta,r,psi_n,contours,extend='both')
pylab.contourf(theta,r,psi_n,contours,extend='both')
pylab.colorbar()

# Create a cartesian plot
pylab.figure()
pylab.contour(theta,r,psi_n,contours,extend='both')
pylab.contourf(theta,r,psi_n,contours,extend='both')
pylab.colorbar()

# Obtain all values at a fixed UT.
print 'psi_n at 0 UT is:'
print psi_n[0,:]
print 'psi_n at 6 UT is:'
print psi_n[135,:]
print 'psi_n at 12 UT is:'
print psi_n[90,:]
print 'psi_n at 18 UT is:'
print psi_n[45,:]

# You can also plot a range of UT's by slicing into the arrays.  Try this on for size:
pylab.figure()
pylab.axes(polar=True)
contours = numpy.linspace(psi_n.min(),psi_n.max(),100)
pylab.contour(theta[0:10,:],r[0:10,:],psi_n[0:10,:],contours,extend='both')
pylab.contourf(theta[0:10,:],r[0:10,:],psi_n[0:10,:],contours,extend='both')
pylab.colorbar()

pylab.show()