#!/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 hemisphericPower.py'
"""

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

def distance(p0, p1):
  """ Calculate the distance between p0 & p1--two points in R**3 """
  return( numpy.sqrt( (p0[0]-p1[0])**2 + 
                      (p0[1]-p1[1])**2 + 
                      (p0[2]-p1[2])**2 ) )

def calcFaceAreas(x,y,z):
  """ 
  Calculate the area of each face in a quad mesh.
 
  >>> calcFaceAreas(numpy.array([[ 0., 1.],[ 1., 0.]]), numpy.array([[ 0.,1.],[ 1.,0.]]), numpy.array([[ 0., 0.],[ 0.,0.]]))
  array([[ 2.]])
  """
  (nLonP1, nLatP1) = x.shape
  (nLon, nLat) = (nLonP1-1, nLatP1-1)

  area = numpy.zeros((nLon, nLat))

  for i in range(nLon):
    for j in range(nLat):
      left  = distance( (x[i,j],   y[i,j],   z[i,j]),   (x[i,j+1],   y[i,j+1],   z[i,j+1]) )
      right = distance( (x[i+1,j], y[i+1,j], z[i+1,j]), (x[i+1,j+1], y[i+1,j+1], z[i+1,j+1]) )
      top =   distance( (x[i,j+1], y[i,j+1], z[i,j+1]), (x[i+1,j+1], y[i+1,j+1], z[i+1,j+1]) )
      bot =   distance( (x[i,j],   y[i,j],   z[i,j]),   (x[i+1,j],   y[i+1,j],   z[i+1,j]) )
      
      area[i,j] = 0.5*(left+right) * 0.5*(top+bot)

  return area

if __name__ == '__main__':
    # Read the data
    filename='AGU_CMITs_mix_2006-12-14T08-00-00Z.hdf'
    hdffile = SD(filename,SDC.READ)
    x=hdffile.select('Grid X').get()
    y=hdffile.select('Grid Y').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()
    
    # Calculate area of MIX grid:
    x[:,0] = 0.0
    y[:,0] = 0.0
    z = numpy.sqrt(1.0-x**2-y**2)
    ri = 6370.0e3  # Radius of ionosphere
    areaMixGrid = calcFaceAreas(x,y,z)*ri*ri
    
    # --- Hemispheric Power
    hp_n = areaMixGrid*energy_n[:-1,:-1] * flux_n[:-1,:-1]
    hp_s = areaMixGrid*energy_s[:-1,:-1] * flux_s[:-1,:-1]

    # KeV/cm^2s to mW/m^2 to GW
    hp_n = hp_n.sum() * 1.6e-21    
    hp_s = hp_s.sum() * 1.6e-21

    print 'File='+filename
    print 'Hemispheric power in North='+str(hp_n)
    print 'Hemispheric power in South='+str(hp_s)