#!/usr/bin/env python

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

# Sky Coyote, April 2008

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

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] = Triangle.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 Triangle.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 sigma_s_deriv(phi, sigma_s, delta, sigma, r):
	return delta * (sigma - sigma_s) / r
	
def u_1b_deriv(phi, u_1b, r, chi_20, F, A, A_s, sigma, sigma_s, d_sigma_d_phi):
	return (-r * chi_20 * d_sigma_d_phi + (F**2 * chi_20 - 1.0) * sigma + (A + A_s) * sigma_s - 2.0 * u_1b) / r
	
def meander(time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_, \
	mode=0, printDetails=False, mrate=1.0, mexp=1.0, ldist=0.0, sdist=0.0, caller=None):

	seterr(all='raise')
	
	delxList = [] 
	delyList = [] 
	rightdepthList = [] 
	leftdepthList = [] 
	
	curvature, left = getGeometry(x, y)
	g = 9.8
	b = width[0] / 2.0
	H = depth[0]
	I = slope[0]
	U = flow[0] / (width[0] * depth[0])
	C_f = (g * H * I) / U**2
	r = H / C_f
	F = U / sqrt(g * H)
	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))
	delta = chi_1**2 * (chi + 0.25) / ((1.0/12.0) * chi**2 + (11.0/360.0) * chi + (1.0/504.0))
	A_s = (181.0 / chi_1) * (H / b)**2 * (2.0 * chi**2 + (4.0/5.0) * chi + (1.0/15.0))
	if diam[0] > 0:
		A = getK(depth[0], slope[0], diam[0])
	else:
		A = 0.0
	
	if printDetails:
		print
		print 'MeanderJP5: mode=%d  mrate=%f  mexp=%f  ldist=%f  sdist=%f' % (mode, mrate, mexp, ldist, sdist)
		print '             b              H              I              U            C_f              r              F'
		print '-------------- -------------- -------------- -------------- -------------- -------------- --------------'
		print '%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f' \
			% (b, H, I, U, C_f, r, F)	
		print
		print '         chi_1            chi         chi_20          delta            A_s              A'
		print '-------------- -------------- -------------- -------------- -------------- --------------'
		print '%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f' \
			% (chi_1, chi, chi_20, delta, A_s, A)	
		print

	for i in range(0, len(x)):
		if caller != None and caller.runState == 2: return (delxList, delyList, rightdepthList, leftdepthList)

		if i == 0:
			phi = 0.0
			sigma = b * abs(curvature)
			sigma_s = sigma[0]
			u_1b = sigma[0]
			bedslope = A * H * abs(curvature[i])

			if printDetails:
				print '    i              x              y            phi      curvature       bedslope          sigma        sigma_s           u_1b'
				print '----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------'
				print '%5d %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f' \
					% (i, x[i], y[i], phi, curvature[i], bedslope, sigma[i], sigma_s, u_1b)	
			
		else:
			d_phi = sqrt((x[i] - x[i - 1])**2 + (y[i] - y[i - 1])**2)
			d_phi_2 = d_phi / 2.0

			if i == 1:
				d_sigma_d_phi = (sigma[i] - sigma[i - 1]) / d_phi
			else:
				d_sigma_d_phi = (sigma[i] - sigma[i - 2]) / (2.0 * d_phi)

			u_1b_k1 = u_1b_deriv(phi, u_1b, r, chi_20, F, A, A_s, sigma[i - 1], sigma_s, d_sigma_d_phi)
			sigma_s_k1 = sigma_s_deriv(phi, sigma_s, delta, sigma[i - 1], r)

			u_1b_k2 = u_1b_deriv(phi + d_phi_2, u_1b + d_phi_2 * u_1b_k1, r, chi_20, F, A, A_s, \
				(sigma[i] + sigma[i - 1]) / 2.0, sigma_s + d_phi_2 * sigma_s_k1, (sigma[i] - sigma[i - 1]) / d_phi)
			sigma_s_k2 = sigma_s_deriv(phi + d_phi_2, sigma_s + d_phi_2 * sigma_s_k1, delta, (sigma[i] + sigma[i - 1]) / 2.0, r)

			u_1b_k3 = u_1b_deriv(phi + d_phi_2, u_1b + d_phi_2 * u_1b_k2, r, chi_20, F, A, A_s, \
				(sigma[i] + sigma[i - 1]) / 2.0, sigma_s + d_phi_2 * sigma_s_k2, (sigma[i] - sigma[i - 1]) / d_phi)
			sigma_s_k3 = sigma_s_deriv(phi + d_phi_2, sigma_s + d_phi_2 * sigma_s_k2, delta, (sigma[i] + sigma[i - 1]) / 2.0, r)

			if i == len(x) - 1:
				d_sigma_d_phi = (sigma[i] - sigma[i - 1]) / d_phi
			else:
				d_sigma_d_phi = (sigma[i] - sigma[i - 2]) / (2.0 * d_phi)

			u_1b_k4 = u_1b_deriv(phi + d_phi, u_1b + d_phi * u_1b_k3, r, chi_20, F, A, A_s, \
				sigma[i], sigma_s + d_phi * sigma_s_k3, d_sigma_d_phi)
			sigma_s_k4 = sigma_s_deriv(phi + d_phi, sigma_s + d_phi * sigma_s_k3, delta, sigma[i], r)

			u_1b += (u_1b_k1 + 2.0 * u_1b_k2 + 2.0 * u_1b_k3 + u_1b_k4) * (d_phi / 6.0)
			sigma_s += (sigma_s_k1 + 2.0 * sigma_s_k2 + 2.0 * sigma_s_k3 + sigma_s_k4) * (d_phi / 6.0)
			
			phi += d_phi

			if printDetails:
				print '%5d %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f' \
					% (i, x[i], y[i], phi, curvature[i], bedslope, sigma[i], sigma_s, u_1b)	
			
		if curvature[i] > 0:
			delx, dely = -left[i] * u_1b * U * rerody[i] * mrate
			rightdepth = depth[i] + bedslope * b
			leftdepth = depth[i] - bedslope * b
		else:
			delx, dely = left[i] * u_1b * U * lerody[i] * mrate
			rightdepth = depth[i] - bedslope * b
			leftdepth = depth[i] + bedslope * b
			
		if rightdepth < 0:
			rightdepth = 0.0
		elif rightdepth > depth[i] * 2.0:
			rightdepth = depth[i] * 2.0 
		if leftdepth < 0:
			leftdepth = 0.0
		elif leftdepth > depth[i] * 2.0:
			leftdepth = depth[i] * 2.0 
			
		delxList.append(delx) 
		delyList.append(dely) 
		rightdepthList.append(rightdepth) 
		leftdepthList.append(leftdepth) 

	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)