#!/usr/bin/env python

#Copyright 2007-2008 John C. Vernaleo
#Unfortunately, I am not comfortable putting unobscured email addresses
#on the web, but these shouldn't be too hard to figure out.
#
#		(my_first_name)@netpurgatory.com
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#02/26/2007
#The only thing you want to use from this is
#get_hdf which takes a string with the filename and path
#If you have anything (other than 0) as the 2nd arguement
#it reads the (not so often plotted) velocity variables also.
#hdf=get_hdf(filename)

import commands
import re
import sys
import math
import os

try:
    from Numeric import *
except ImportError:
    try:
        from numarray import *
    except ImportError:
        print "No array package present"
        sys.exit()


class hdf:
    def __init__(self):
        self.x3=[]
        self.x2=[]
        self.x1=[]
        self.v1=[]
        self.v2=[]
        self.v3=[]
        self.b1=[]
        self.b2=[]
        self.b3=[]
        self.d=[]
        self.e=[]
        self.gp=[]
        self.gamma=5.0/3.0
    def info(self):
        return(self.x3,self.x2,self,x3,self.v1,self.v2,self,v3,self.b1,self.b2,self.b3,self,d,self.e,self.gp,self.gamma)
    def calc_ent(self):
        """Calculate specific Entropy
	Cannot use both entropy and log entropy!
        """
        self.a=((self.gamma-1.0)*self.e)/(self.d**self.gamma)
    def calc_logent(self):
        """Calculate log of specific Entropy
	Cannot use both entropy and log entropy!
        """
        self.a=log10(((self.gamma-1.0)*self.e)/(self.d**self.gamma))
    def calc_tkev(self):
        """Calculate Temperature in KeV
        """
        self.tkev=(self.e/self.d)*9.5293
    def calc_gp(self):
        """Calculate gravitational potential for our
        cluster simulations
        """
        cs=1
        #need gp to have x1 in it!
        self.gp=self.d
        for k in range(0,size(self.x3)):
            for j in range(0,size(self.x2)):
                for i in range(0,size(self.x1)):
                    self.gp[k,j,i]=self.x1[i]
        self.gp=((cs**2)/self.gamma)*log(1.0/(1+(self.gp/2.0)**2)**(3.0/4.0))
    def calc_volume_2d(self):
	"""Calculate per cell volume for 2D simulations in spherical
	"""
	x1size=len(self.x1)
	x2size=len(self.x2)
	x3size=len(self.x3)

	self.volume=zeros((x3size,x2size,x1size),Float)
	for j in range(0,x2size):
	    if j == 0:
		dt=self.x2[0]+((self.x2[1]-self.x2[0])/2.0)-0.001
	    elif j == (x2size-1):
		dt = math.pi - self.x2[x2size-2]+((self.x2[x2size-1]-self.x2[x2size-2])/2.0)
	    else:
		dt=(self.x2[j+1]-self.x2[j-1])/2.0
	    for i in range(0,x1size):
		if i == 0:
		    dr = self.x1[0]+((self.x1[1]-self.x1[0])/2.0)-0.1
		elif i == (x1size-1):
		    dr = 30 - self.x1[x1size-2]+((self.x1[x1size-1]-self.x1[x1size-2])/2.0)
		else:
		    dr=(self.x1[i+1]-self.x1[i-1])/2.0
                #Make sure no negative volumes slip in.
                test(dt,'dt',j)
                test(dr,'dr',i)
		self.volume[0,j,i]=2.0*dr*dt*math.pi*(self.x1[i]**2)*math.sin(self.x2[j])

    def calc_volume_3d(self):
	"""Calculate per cell volume for 3D simulations in spherical
	"""
	x1size=len(self.x1)
	x2size=len(self.x2)
	x3size=len(self.x3)

	self.volume=zeros((x3size,x2size,x1size),Float)
	for k in range(0,x3size):
            if k == 0:
                dp=self.x3[0]+((self.x3[1]-self.x3[0])/2.0)-0.003
            elif k == (x3size-1):
                dp=2.0*math.pi - self.x3[x3size-2]+((self.x3[x3size-1]-self.x3[x3size-2])/2.0)
            else:
                dp=(self.x3[k+1]-self.x3[k-1])/2.0
            for j in range(0,x2size):
                if j == 0:
                    dt=self.x2[0]+((self.x2[1]-self.x2[0])/2.0)-0.003
                elif j == (x2size-1):
                    dt=(math.pi/2.0) - self.x2[x2size-2]+((self.x2[x2size-1]-self.x2[x2size-2])/2.0)
                else:
                    dt=(self.x2[j+1]-self.x2[j-1])/2.0
                for i in range(0,x1size):
                    if i == 0:
                        dr=self.x1[0]+((self.x1[1]-self.x2[0])/2.0)-0.2
                    elif i == (x1size-1):
                        dr=20.0 - self.x1[x1size-2]+((self.x1[x1size-1]-self.x1[x1size-2])/2.0)
                    else:
                        dr=(self.x1[i+1]-self.x1[i-1])/2.0
                    #Make sure no negative volumes slip in.
                    test(dp,'dp',k)
                    test(dt,'dt',j)
                    test(dr,'dr',i)
                    self.volume[k,j,i]=(self.x1[i]**2)*math.sin(self.x2[j])*dt*dt*dp

def test(val,name,it):
    if(val<=0):
        print "Error, negative volume element."
        print name,"at",it,val
        sys.exit()

def readaxis(size,label,hdffile):
    cmd="hdp dumpsds -d -n "+label+" "+hdffile
    temp=commands.getoutput(cmd)
    temp=temp.split()
    out=zeros(size,Float)
    i=0
    for x in temp:
        out[i]=float(x)
        i+=1
    return out

def readvar(size1,size2,size3,label,hdffile):
    cmd="hdp dumpsds -d -n "+label+" "+hdffile
    temp=commands.getoutput(cmd)
    temp=temp.split()
    size=size1*size2*size3
    out=zeros((size),Float)
    i=0
    for x in temp:
        out[i]=float(x)
        i+=1
    out=resize(out,(size1,size2,size3))
    return out

def get_hdf(file,opts=""):
    data=hdf()
    #Velocities
    if opts.find("v") != -1:
        vel=1
    else:
        vel=0
    #Magnetic fields
    if opts.find("b") != -1:
        b=1
    else:
        b=0
    #Potential
    if opts.find("p") != -1:
        p=1
    else:
        p=0
    status=os.path.exists(file)
    if not status:
	print "Cannot find "+file
	sys.exit()
    cmd="hdp dumpsds -h "+file
    header=commands.getoutput(cmd)
    
    finddims=re.compile('Dimension Variable Name')
    findvars=re.compile('Data-Set-')

    header=header.splitlines()

    i=0
    indices=[]
    for x in header:
        if finddims.search(x):
            indices.append(i)
            
        if findvars.search(x):
            indices.append(i)
        i+=1
    findsize=re.compile('Size')
    getnum=re.compile('\d+')
    findname=re.compile('Value')
    axissize=[]
                
    for y in range(0,3):
        for x in range(indices[y],indices[y+1]):
            if findsize.search(header[x]):
                num=getnum.findall(header[x])
                num=int(num[0])
                axissize.append(num)

    data.x3=readaxis(axissize[0],"fakeDim0",file)
    data.x2=readaxis(axissize[1],"fakeDim1",file)
    data.x1=readaxis(axissize[2],"fakeDim2",file)
    if vel:
        data.v1=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-2",file)
        data.v2=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-3",file)
        data.v3=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-4",file)
    if b:
        data.b1=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-5",file)
        data.b2=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-6",file)
        data.b3=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-7",file)
        data.d=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-8",file)
        data.e=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-9",file)
    else:
        data.d=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-5",file)
        data.e=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-6",file)
    if p:
        data.gp=readvar(axissize[0],axissize[1],axissize[2],"Data-Set-7",file)

    return data