#!/usr/bin/env python
#
# ----------------------------------------------------------------------
#
#                           Brad T. Aagaard
#                        U.S. Geological Survey
#
# ----------------------------------------------------------------------
#

cell = "tet4"
dx = 200

inputRoot = "output/%s_%3dm-fault" % (cell,dx)
outdir = "scecfiles/%s_%3dm/" % (cell,dx)

# ----------------------------------------------------------------------
import h5py
import numpy
import time

# ----------------------------------------------------------------------
threshold = 0.001 # threshold for detecting slip has started
maxTime = 1.0e+10 # Large value for default rupture time

h5 = h5py.File("%s.h5" % inputRoot, 'r', driver="sec2")
vertices = h5['geometry/vertices'][:]
slip = h5['vertex_fields/slip'][:]
slipRate = h5['vertex_fields/slip_rate'][:]
timeStamps =  h5['time'][:].ravel()
dt = timeStamps[1] - timeStamps[0]

h5.close()

nsteps = timeStamps.shape[0]
npts = vertices.shape[0]

rupTime = maxTime * numpy.ones( (npts,), dtype=numpy.float64)

# Create buffer for current rupture time
tmpTime = numpy.zeros( (npts,), dtype=numpy.float64)

itime = 0
for timestamp in timeStamps:
    t = timeStamps[itime]

    # Compute magnitude of slip rate
    slipRateMag = (slipRate[itime,:,0]**2 + slipRate[itime,:,1]**2)**0.5

    # Set rupture time at locations where threshold is exceeded
    mask = slipRateMag > threshold
    tmpTime[:] = maxTime
    tmpTime[mask] = t

    indices = numpy.where(tmpTime < rupTime)[0]
    rupTime[indices] = t

    itime += 1

# Write data
headerA = \
    "# problem = TPV102\n" + \
    "# author = Brad Aagaard\n" + \
    "# date = %s\n" % (time.asctime()) + \
    "# code = PyLith\n" + \
    "# code_version = 1.8.1a (scecdynrup branch)\n" + \
    "# element_size = %s m\n" % dx + \
    "# Contour data in 3 columns of E14.6:\n" + \
    "# Column #1 = Distance along strike from hypocenter (m)\n" + \
    "# Column #2 = Distance down-dip from surface (m)\n" + \
    "# Column #3 = Rupture time (s)\n" + \
    "#\n" + \
    "# Data fields\n" + \
    "j k t\n" + \
    "#\n"


distDip = -vertices[:,2]
distStrike = vertices[:,1]

filename = "%sruptime.dat" % (outdir)
fout = open(filename, "w")
fout.write(headerA)
data = numpy.transpose((distStrike, distDip, rupTime))
numpy.savetxt(fout, data, fmt='%14.6e')
fout.close()
