#!/usr/bin/env python

'''Compute values for stream meandering using Johannesson-Parker 1989 method'''

# Sky Coyote, March-April 2008

import sys, signal
# must have array math package from numpy.org
from numpy import *
import Triangulate

def decodeData():
	'''Read text data input, turn into lists'''
	lines = []
	#sys.stderr.write('Meander.py.decodeData: before initial readline()\n')
	line = sys.stdin.readline()
	#sys.stderr.write('Meander.py.decodeData: after initial readline()\n')
	if line == '': sys.exit(0)  # end of input when run stand-alone
	fields = line.strip('\n').split()
	time = float(fields[0])
	stations = int(fields[1])
	stnserod = int(fields[2])
	for i in range(stnserod):
		line = sys.stdin.readline()
		lines.append(line.strip('\n'))
	#sys.stderr.write('Meander.py.decodeData: got %d lines\n' % (len(lines)))

	# initialize lists
	x = []
	y = []
	xs = []
	dels = []
	flow = []
	rerody = []
	lerody = []
	slope = []
	width = []
	depth = []
	diam = []
	lambda_ = []
	
	# append each column to appropriate list
	for i in range(len(lines)):
		fields = lines[i].split()
		
		x.append(float(fields[1]))
		y.append(float(fields[2]))
		xs.append(float(fields[3]))
		dels.append(float(fields[4]))
		flow.append(float(fields[5]))
		rerody.append(float(fields[6]))
		lerody.append(float(fields[7]))
		slope.append(float(fields[8]))
		width.append(float(fields[9]))
		depth.append(float(fields[10]))
		diam.append(float(fields[11]))
		lambda_.append(float(fields[12]))
		
	return (time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_)

def getGeometry(x, y):
	xy = column_stack((array(x), array(y)))
	tri = [None] * len(xy)
	curvature = zeros((len(xy),))
	left = zeros((len(xy), 2))
	for j in range(1, len(xy) - 1):
		pts = []
		pts.append(xy[j - 1])
		pts.append(xy[j])
		pts.append(xy[j + 1])
		pts = array(pts)
		try:
			tri[j] = Triangulate.Triangle(pts)
			curvature[j] = tri[j].orientation / tri[j].circle[2]
			left[j] = tri[j].orientation * (array((tri[j].circle[0], tri[j].circle[1])) - xy[j])
			left[j] /= sqrt(dot(left[j], left[j]))
		except Triangulate.BadTriangleException: 
			curvature[j] = 0.0
			left[j, 0] = -(xy[j + 1, 1] - xy[j, 1])
			left[j, 1] = xy[j + 1, 0] - xy[j, 0]
			left[j] /= sqrt(dot(left[j], left[j]))
	curvature[0] = curvature[1]
	#curvature[0] = 0.0
	if tri[1] != None:
		left[0] = tri[1].orientation * (array((tri[1].circle[0], tri[1].circle[1])) - xy[0])
		left[0] /= sqrt(dot(left[0], left[0]))
	else:
		left[0, 0] = -(xy[1, 1] - xy[0, 1])
		left[0, 1] = xy[1, 0] - xy[0, 0]
		left[0] /= sqrt(dot(left[0], left[0]))
	curvature[len(xy) - 1] = curvature[len(xy) - 2]
	#curvature[len(xy) - 1] = 0.0
	if tri[len(xy) - 2] != None:
		left[len(xy) - 1] = tri[len(xy) - 2].orientation \
			* (array((tri[len(xy) - 2].circle[0], tri[len(xy) - 2].circle[1])) - xy[len(xy) - 1])
		left[len(xy) - 1] /= sqrt(dot(left[len(xy) - 1], left[len(xy) - 1]))
	else:
		left[len(xy) - 1, 0] = -(xy[len(xy) - 1, 1] - xy[len(xy) - 2, 1])
		left[len(xy) - 1, 1] = xy[len(xy) - 1, 0] - xy[len(xy) - 2, 0]
		left[len(xy) - 1] /= sqrt(dot(left[len(xy) - 1], left[len(xy) - 1]))
	return (curvature, left)
	
def getK(depth, slope, diam):
	radh = depth
	if depth <= 0:
		sys.stderr.write('depth <= 0\n')
		sys.exit(1)

	shields = radh * slope / 1.65 / diam
	if diam < 1.5e-4:
		critshields = 0.08 
	elif diam < 2.5e-4:
		critshields = 0.052 
	elif diam < 3.5e-4:
		critshields = 0.039 
	elif diam < 4.5e-4:
		critshields = 0.035 
	elif diam < 5.5e-4:
		critshields = 0.034 
	elif diam < 6.5e-4:
		critshields = 0.033 
	elif diam < 7.5e-4:
		critshields = 0.034 
	elif diam < 8.5e-4:
		critshields = 0.034
	elif diam < 9.5e-4:
		critshields = 0.034 
	elif diam < 0.0015:
		critshields = 0.035 
	elif diam < 0.0025:
		critshields = 0.045 
	elif diam < 0.0035:
		critshields = 0.05 
	elif diam < 0.0045:
		critshields = 0.055 
	else:
		critshields = 0.056 

	# approximate Engelund diagram for subcritical flow:
	if shields >= 0.1 and shields < 1.0:
		# Computing 2nd power
		d__2 = log10(shields) + 1.03
		d__1 = d__2 * d__2 * 0.74 - 1.18
		grainshields = pow(10.0, d__1)
	elif shields >= 1.0 and shields < 2.0:
		# Computing 2nd power
		d__1 = shields
		grainshields = d__1 * d__1 * 0.4
	else:
		grainshields = shields

	# calculation for skin friction:
	try:
		K = (grainshields / shields) * sqrt(grainshields / critshields) * (log(grainshields * 11.0 * depth / shields / diam) * 0.5695 - 0.3606)
	except ZeroDivisionError:
		K = 0
	return K
	
def meander(time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_, printDetails=False, migration_rate=1.0):
	delxList = [] 
	delyList = [] 
	rightdepthList = [] 
	leftdepthList = [] 
	
	#print 'MeanderJP.meander'
	g = 9.8
	curvature, left = getGeometry(x, y)
	curvature_e = curvature
	integral = 0.0
	integral_e = 0.0
	integrands = []
	integrands_e = []
	dt = 1.0
	
	seterr(all='raise')
	
	if printDetails:
		print
		print 'migration_rate = %f' % (migration_rate)
		print
		print "    i              s              U            C_f            A_s              K      curvature      integrand       integral           u_b1            acc"
		print "----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------"

	for i in range(len(x)):
		H = depth[i]
		I = slope[i]
		U = flow[i] / (width[i] * depth[i])
		C_f = (g * H * I) / U**2
		
		chi_1 = 0.077 / sqrt(C_f)
		chi = chi_1 - (1.0/3.0)
		chi_20 = (1.0 / chi_1**3) * (chi**3 + chi**2 + (2.0/5.0) * chi + (2.0/35.0))
		
		b = width[i]
		A_s = (181.0 / chi_1) * (2.0 * H / b)**2 * (2.0 * chi**2 + (4.0/5.0) * chi + (1.0/15.0))
		
		K = getK(depth[i], slope[i], diam[i])
		
		s = xs[i]
		C = abs(curvature[i])
		C_e = abs(curvature_e[i])
		
		integrand = C * exp(2.0 * C_f * s / H)
		integrands.append(integrand)
		integrand_e = C_e * exp(2.0 * C_f * s / H)
		integrands_e.append(integrand_e)
		if i > 0:
			integral += ((integrands[i] + integrands[i - 1]) / 2.0) * (xs[i] - xs[i - 1])
			integral_e += ((integrands_e[i] + integrands_e[i - 1]) / 2.0) * (xs[i] - xs[i - 1])
			
		u_b1 = chi_20 * U * b * C \
			+ (C_f * U * b**2 / H) * (chi_20 * (U**2 / (g * H) + 2.0) - 1.0) * exp(-2.0 * C_f * s / H) * integral \
			+ (C_f * U * b / H) * (K + A_s) * exp(-2.0 * C_f * s / H) * integral_e
			
		if curvature[i] > 0:
			delx, dely = -left[i] * u_b1 * rerody[i] * migration_rate * dt
		else:
			delx, dely = left[i] * u_b1 * lerody[i] * migration_rate * dt
	
		acc = U**2 * curvature_e[i]
		del_depth = (acc / g) * (b / 2.0)
		rightdepth = depth[i] + del_depth
		leftdepth = depth[i] - del_depth

		delxList.append(delx) 
		delyList.append(dely) 
		rightdepthList.append(rightdepth) 
		leftdepthList.append(leftdepth) 
		
		if printDetails:
			print "%5d %14.6f %14.6f %14.6f %14.6f %14.6f %14.9f %14.6e %14.6e %14.6f %14.6f" \
				% (i, s, U, C_f, A_s, K, curvature[i], integrand, integral, u_b1, acc)

	return (delxList, delyList, rightdepthList, leftdepthList)

def encodeData(delx, dely, rightdepth, leftdepth):
	'''Write outputs to text'''
	for i in range(len(delx)):
		print '%5d %14.6f %14.6f %14.6f %14.6f' % (i, delx[i], dely[i], rightdepth[i], leftdepth[i])
	sys.stdout.flush()
	#sys.stderr.write('Meander.py.encodeData: wrote %d lines\n' % (len(delx)))

def sigintHandler(a, b):
	'''Control-C handler'''
	sys.exit(0)

if __name__ == '__main__':
	# install ctrl-c handler
	signal.signal(signal.SIGINT, sigintHandler)
	
	while True:
		# get data from standard input
		time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_ = decodeData()
		# print inputs
		if False:
			print '%g %d %d' % (time, stations, stnserod)
			for i in range(len(x)):
				print "$$%5d %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f" \
					% (i, x[i], y[i], xs[i], dels[i], flow[i], rerody[i], lerody[i], slope[i], width[i], depth[i], diam[i], lambda_[i])
		# calculate meandering
		delx, dely, rightdepth, leftdepth = \
			meander(time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_)
		# print outputs
		if False:
			for i in range(len(delx)):
				print "$$%5d %14.6f %14.6f %14.6f %14.6f" % (i, delx[i], dely[i], rightdepth[i], leftdepth[i])
		# send results to standard output
		encodeData(delx, dely, rightdepth, leftdepth)