root/virtexp/trunk/virtexp/tests/currents/currents.py

#!/usr/local/python272/bin/python
from glob import glob
from pprint import pformat
import cmath
import re

log_pattern = "/var/opt/gmc/gliders/ramses/logs/ramses_network_20111213*.log"
stamp_pattern = re.compile("ramses_network_(\w+)")
vx_pattern = re.compile("m_water_vx\(m/s\)=(\S+)")
vy_pattern = re.compile("m_water_vy\(m/s\)=(\S+)")

isNorth = lambda phi: (phi >= (3*cmath.pi/8))  and (phi <= (5*cmath.pi/8))
isSouth = lambda phi: (phi <= (-3*cmath.pi/8)) and (phi >= (-5*cmath.pi/8))
isNEast = lambda phi: (phi >  (cmath.pi/8))    and (phi <  (3*cmath.pi/8))
isNWest = lambda phi: (phi >  (5*cmath.pi/8))  and (phi <  (7*cmath.pi/8))
isSEast = lambda phi: (phi <  (-cmath.pi/8))   and (phi >  (-3*cmath.pi/8))
isSWest = lambda phi: (phi <  (-5*cmath.pi/8)) and (phi >  (-7*cmath.pi/8))
isEast  = lambda phi: ((phi >= 0) and (phi <=  (cmath.pi/8))) or \
                      ((phi <= 0) and (phi >=  (-cmath.pi/8)))
isWest  = lambda phi: ((phi >= 0) and (phi >=  (7*cmath.pi/8))) or \
                      ((phi <= 0) and (phi <=  (-7*cmath.pi/8)))


logs = glob(log_pattern)
logs.sort()
logs.reverse()

currents = []

for log in logs:
    stamp = stamp_pattern.search(log).groups()[0]
    handle = open(log)
    content = handle.readlines()
    handle.close()
    content = zip(content[:-1],content[1:])
    for line0,line1 in content:
        match = vx_pattern.search(line0)
        if match:
            vx = match.groups()[0]
            match = vy_pattern.search(line1)
            if match:
                vy = match.groups()[0]
                currents.append((stamp,complex(float(vx),float(vy))))
                break

ccurrents = pformat(currents)
currents = [(stamp,cmath.polar(current)) for stamp, current in currents]
pcurrents = pformat(currents)

N  = [stamp for stamp,current in currents if isNorth(current[1])]
NE = [stamp for stamp,current in currents if isNEast(current[1])]
E  = [stamp for stamp,current in currents if isEast(current[1])]
SE = [stamp for stamp,current in currents if isSEast(current[1])]
S  = [stamp for stamp,current in currents if isSouth(current[1])]
SW = [stamp for stamp,current in currents if isSWest(current[1])]
W  = [stamp for stamp,current in currents if isWest(current[1])]
NW = [stamp for stamp,current in currents if isNWest(current[1])]

print "Number of surfacings:",len(currents)
print "Current distribution:"
print "North:     %s%%" % (float(len(N))/len(currents))
print "Northeast: %s%%" % (float(len(NE))/len(currents))
print "East:      %s%%" % (float(len(E))/len(currents))
print "Southeast: %s%%" % (float(len(SE))/len(currents))
print "South:     %s%%" % (float(len(S))/len(currents))
print "Southwest: %s%%" % (float(len(SW))/len(currents))
print "West:      %s%%" % (float(len(W))/len(currents))
print "Northwest: %s%%" % (float(len(NW))/len(currents))
print
print "Cartesian currents:"
print ccurrents
print
print "Polar currents:"
print pcurrents