#!/usr/bin/env python
import numpy as np
import sord, time

# parameters
path = 'infile'
dtype = '<'
dix_x = 1
dit_x = 25
dix_t = 2
dit_t = 2

# Graves' broadband timeslice format.
fh = open( path, 'br' )
start = tuple( np.fromfile( fh, dtype + 'i', 4 ) )
shape = tuple( np.fromfile( fh, dtype + 'i', 4 ) )
delta = tuple( np.fromfile( fh, dtype + 'f', 4 ) )
rot = np.fromfile( fh, dtype + 'f', 1 )
origin = tuple( np.fromfile( fh, dtype + 'f', 2 ) )
meta = dict( start=start, shape=shape, delta=delta, rot=rot, origin=origin )
print meta
sord.util.save( path + '-meta.py' )

# decimate arrays and compute PGV, PGD
nn = shape[:-1]
nt = shape[-1]
n  = np.prod( nn )
u1 = np.zeros( nn )
u2 = np.zeros( nn )
u3 = np.zeros( nn)
pgv = np.zeros( nn )
pgd = np.zeros( nn )
x1 = open( 'v1_x', 'wb' )
x2 = open( 'v2_x', 'wb' )
x3 = open( 'v3_x', 'wb' )
t1 = open( 'v1_t', 'wb' )
t2 = open( 'v2_t', 'wb' )
t3 = open( 'v3_t', 'wb' )
t0 = time.time()
for i in xrange( nt ):
    v1 = np.fromfile( fh, dtype, n ).reshape( nn[::-1] ).T
    v2 = np.fromfile( fh, dtype, n ).reshape( nn[::-1] ).T
    v3 = np.fromfile( fh, dtype, n ).reshape( nn[::-1] ).T
    u1 = u1 + dt * v1
    u2 = u2 + dt * v2
    u3 = u3 + dt * v3
    pgv = np.maximum( pgv, v1*v1 + v2*v2 + v3*v3 )
    pgd = np.maximum( pgd, u1*u1 + u2*u2 + u3*u3 )
    if np.mod( i, dit_x ) == 0:
        v1[::dix_x,::dix_x].T.tofile( x1 )
        v2[::dix_x,::dix_x].T.tofile( x2 )
        v3[::dix_x,::dix_x].T.tofile( x3 )
    if np.mod( i, dit_t ) == 0:
        v1[::dix_t,::dix_t].T.tofile( t1 )
        v2[::dix_t,::dix_t].T.tofile( t2 )
        v3[::dix_t,::dix_t].T.tofile( t3 )
    sord.util.progress( i+1, nt, time.time()-t0 )
x1.close()
x2.close()
x3.close()
t1.close()
t2.close()
t3.close()
np.array( np.sqrt( pgv ), dtype ).T.tofile( 'pgv' )
np.array( np.sqrt( pgd ), dtype ).T.tofile( 'pgd' )
del( v1, v2, v3, u1, u2, u3, pgv, pgd )

# transpose t-arrays
nn = t_shape[0], t_shape[1] * t_shape[2] * t_shape[3]
n = nn[0] * nn[1]
np.fromfile( 'v1_t', dtype, n ).reshape( nn ).T.tofile( 'v1_t' )
np.fromfile( 'v2_t', dtype, n ).reshape( nn ).T.tofile( 'v2_t' )
np.fromfile( 'v3_t', dtype, n ).reshape( nn ).T.tofile( 'v3_t' )