#! /usr/bin/env python
#
# Python script to add calculate the "offset" for a horizon.
# This is designed to assist in creating a 2d Rayinvr velocity
# model, based on a 1d velocity profile.  Very narrow scope script.
#
# Todo: possibly include decimation.  Maybe use getopt or optparse?
#
#
# 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


if len(sys.argv) != 4:
	print """
	Error: usage is 'horizonOffset.py [originTrace] [input] [output]'.
	
	Simple script to calculate the offset for a horizon, in relation to
	a given trace.  Input horizon must be in Geoquest card image 7 format.
	"""
	#print len(sys.argv)
	sys.exit()

originTrace, inFile, outFile = sys.argv[1:4]
inputData = []
resultData = []
originX, originY = 0.0, 0.0

# Looping through the whole file to get the origin
for line in file(inFile, "r").readlines():
	# splitting the line into columns
	fields = line.split()
	
	# taking only lines of 8 columns
	if len(fields) != 8:
		continue
	
	# Associating the command line origin with a trace in the file.
	try:
		originTrace = int(originTrace)
	except:
		print "Please choose an integer value for the originTrace."
	
	try:
		if originTrace == int(fields[4]):
			originX = float(fields[0])
			originY = float(fields[1])
	except:
		continue

# Checking to see that the origins have been read from the selected trace
# Would prefer to use an exception
if originX == 0.0 or originY == 0.0:
	print "\n  Error, trace number does not exist in input file."
	print ("%s %s\n" % ("  Trace number:", originTrace))
	sys.exit()

# Now that we have the origin, we can calculate the offsets for each trace.
outFileObj = open(outFile, "w")
outFileObj.write("Converted horizon\ntrace, travelTime(s), offset(km)\n")

for line in file(inFile, "r").readlines():	
	fields = line.split()
	
	# taking only lines of 8 columns
	if len(fields) != 8:
		continue
	
	# Assigning names to columns, throwing out those that dont work.
	try:
		traceX = float(fields[0])
		traceY = float(fields[1])
		time   = float(fields[2]) / 1000
		trace  = int(fields[4])
	except:
		continue
	
	# Using pythagoreous to calculate offset:
	offset = ((traceX - originX) **2 + (traceY - originY) **2) **0.5 / 1000
	
	# Simple way of calculating offset sign.  Make sure it matches the OBS SEGY.
	if originX - traceX > 0:
		offset = -offset
	
	outFileObj.write("%4d %06.3f %07.3f\n" % (trace, time, offset))

outFileObj.close()