#!/usr/bin/env python
# script to convert from the velocity profile in vm.out to
# one that is actually usable in GMT.
#
# Modified to handle several velocity profiles at different
# x offsets, and output to an xyz file.
#
# Chris LeBlanc
# Dalhousie University, 2003
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU program 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
# program General Public License for more details.
#
# You should have received a copy of the GNU program 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

import sys
import fileinput
from string import *

if len(sys.argv) != 3:
    print """
    
    Error: usage is 'vmodel2xyz.py [input] [output]'.
    
    Used to convert the velocity profile(s) in vm.out (rayinvr) 
    to a simple 'x offset(km), depth(km), velocity(km/s)' format.
    Will also handle 'x offset(km), 2wayTime(s), velocity(km/s)'.
    Rayinvr must be set up to output a velocity profile (ie: set 
    ivp, izort, and xvp in vm.in and run 'vmodel'.  See Rayinvr 
    docs).
    
    """
    #print len(sys.argv)
    sys.exit()

inFile, outFile = sys.argv[1:3]
inputData = []
resultData = []

# Checking of input data, should be 5 or 7 (arrrg, zelt!) columns
# of data.
xoffset = -99999
for line in fileinput.input(inFile):
    fields = line.split()
    
    # Obtaining the x offset
    try:
        if str(fields[0]) == "x" and str(fields[1]) == "position:":
            xoffset = float(fields[2])
            #print xoffset
    except:
        continue
    
    # Only taking lines with 5 or 7 columns
    if len(fields) < 5:
        continue
    if len(fields) > 7:
        continue
    if len(fields) == 6:
        continue

    # Eliminating all lines that dont have a number in the
    # second column
    try:
        float(fields[1])
    except ValueError:
        continue
    
    # Deleting the "layer# x" fields
    if str(fields[0]) == "layer#":
        del fields[0:2]

    inputData.append((xoffset, fields[1], fields[3]))

outFileObj = open(outFile, "w")

outFileObj.write("Rayinvr 1D velocity profile \n\
[x offset(km), depth(km) or 2wayTime, velocity(km/s)]\n")

for idx in range(1, len(inputData), 2):
    zOffset, z1, z2 = inputData[idx - 1]
    vOffset, vp1, vp2 = inputData[idx]
    #print (z1, vp1, z2, vp2)
    outFileObj.write("%06.1f %s %s\n%06.1f %s %s\n" % \
                    (zOffset, z1, vp1, vOffset, z2, vp2))

outFileObj.close()