#! /usr/bin/env python
#  Reads and plots euler3d xn.dat files
#  Charles O'Neill 30 April 2008
#

from numpy import *
from pylab import *
import os,fnmatch

# Load xn.dat
xn=load("xn.dat")

# Sizes
nr=(size(xn,1)-2)/3

# Open Vec
vecflag = 0
for fileName in os.listdir ( "." ):
   if fnmatch.fnmatch ( fileName, '*.vec' ):
     vec=open(fileName,'r')
     vecflag = 1
     break

# Open Con
for fileName in os.listdir ( "." ):
   if fnmatch.fnmatch ( fileName, '*.con' ):
     con=open(fileName,'r')
     break

# Find Dimensionalization Velocity in .con
mach =1.0
ainf =1.0
refdim = 1.0

a = con.readlines()
for part in a:
  part = part.split()

  if(part==[]):
    continue

  if(part[0]=="ainf"):
    part[2] = part[2].lower().replace("d","e").rstrip(',')
    ainf=float(part[2].rstrip(','))
  if(part[0].lower()=="mach"):
    part[2] = part[2].lower().replace("d","e").rstrip(',')
    mach=float(part[2].rstrip(','))
  if(part[0].lower()=="refdim"):
    part[2] = part[2].lower().replace("d","e").rstrip(',')
    refdim = float(part[2].rstrip(','))  

Vdim = mach * ainf / refdim
print Vdim

# Time
t=xn[:,1]

# Displacement
if(vecflag):
  subplot(411)
else:
  subplot(311)
ylabel('Displacement')
d_list = range(2,2+nr)
mode = 0
for d in d_list:
    plot(t,xn[:,d])
    mode = mode + 1
    text(t[-1],xn[-1,d],str(mode), fontsize=9)

# Velocity
if(vecflag):
  subplot(412)
else:
  subplot(312)
ylabel('Velocity')
v_list = range(2+nr,2+nr*2)
mode = 0
for v in v_list:
    plot(t,xn[:,v])
    mode = mode + 1
    text(t[-1],xn[-1,v],str(mode), fontsize=9)

# Force
if(vecflag):
  subplot(413)
else:
  subplot(313)

ylabel('Force')
xlabel('Non-dim Time')
f_list = range(2+nr*2,2+nr*3)
mode = 0
for f in f_list:
    plot(t,xn[:,f])
    mode = mode + 1
    text(t[-1],xn[-1,f],str(mode), fontsize=9)

# Energy
if(vecflag):
	LineVec = vec.readline()  # This reads the header
	nr = int(vec.readline())  # Number of modes

	LineVec = vec.readline()  # This reads the header
	mm =array([])
	mass_list = range(0,nr)
	for mass in mass_list:
	  LineVec = vec.readline()
	  LineVec = LineVec.lower()
	  LineVec = LineVec.replace("d","e")
	  mm = append(mm,map(float,LineVec.split()))
	mm=mm.reshape(nr,nr)

	LineVec = vec.readline()  # This reads the header
	cc =array([])
	for mass in mass_list:
	  LineVec = vec.readline()
	  LineVec = LineVec.lower()
	  LineVec = LineVec.replace("d","e")
	  cc = append(cc,map(float,LineVec.split()))
	cc=cc.reshape(nr,nr)

	LineVec = vec.readline()  # This reads the header
	kk =array([])
	for mass in mass_list:
	  LineVec = vec.readline()
	  LineVec = LineVec.lower()
	  LineVec = LineVec.replace("d","e")
	  kk = append(kk,map(float,LineVec.split()))
	kk=kk.reshape(nr,nr)
	vec.close() # Close the file

	times=range(0,size(t))
	Energy=array([])
	for time in times:
      e = 0.0
      e = dot(dot(   xn[time,2:nr+2],       kk),       xn[time,2:nr+2])  \
      + Vdim**2*dot(dot(   xn[time,2+nr:2+nr*2],  mm),  xn[time,2+nr:2+nr*2])
      Energy = append(Energy,e)


	subplot(414)
	ylabel('Energy')
	semilogy(t,Energy)

show()
V,D = linalg.eig(dot(kk,linalg.inv(mm)))
print V