#!/usr/bin/env python

"""MeanderTS.py: topographic steering meandering model (S. Lancaster, 1998)
   Translated from 'meander.cpp' in Child by Sky Coyote, March 2008"""

import sys, signal
from numpy import *

# constant values
c_b7 = 1.0

def decodeData():
	lines = []
	#sys.stderr.write('MeanderTS.py.decodeData: before initial readline()\n')
	line = sys.stdin.readline()
	#sys.stderr.write('MeanderTS.py.decodeData: after initial readline()\n')
	if line == '': sys.exit(0)
	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('MeanderTS.py.decodeData: got %d lines\n' % (len(lines)))

	x = []
	y = []
	xs = []
	dels = []
	flow = []
	rerody = []
	lerody = []
	slope = []
	width = []
	depth = []
	diam = []
	lambda_ = []
	
	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 angle_(y, x):
	if x != 0 or y != 0: return arctan2(y, x)
	else: return 0.0

def initialize(stnserod, x, y):
	rho = 1000.0
	grav = 9.81
	phi  = zeros((stnserod,))
	delx = zeros((stnserod,))
	dely = zeros((stnserod,))
	
	for s in range(stnserod):
		if s < stnserod - 1:
			delx[s] = x[s + 1] - x[s]
			dely[s] = y[s + 1] - y[s]
			phi[s] = angle_(dely[s], delx[s])
		else:
			delx[s] = delx[s - 1]
			dely[s] = dely[s - 1]
			phi[s] = phi[s - 1]
			
	return (rho, grav, phi, delx, dely)

def d_sign(a, b):
	if b < 0: return -abs(a)
	else: return abs(a)

def getcurv(stnserod, delx, dely, dels):
	global c_b7
	curvature = zeros((stnserod,))
	
	# use law of cosines to find magnitude of angle, use cross
	# product to find sign of curvature
	
	for s in range(1, stnserod - 1):
		a = dels[s]
		b = dels[s - 1]
		if a == 0 or b == 0:
			sys.stderr.write("dels(s) or dels(s-1) equals zero, s = %r\n" % (s))
			sys.exit(1)
		
		c__ = sqrt((delx[s - 1] + delx[s]) * (delx[s - 1] + delx[s]) \
			+ (dely[s - 1] + dely[s]) * (dely[s - 1] + dely[s]))
		carg = (a * a + b * b - c__ * c__) / (a * 2.0 * b)
		if carg < -1: carg = -1.0
		elif carg > 1: carg = 1.0
		
		mag = (arccos(-1.0) - arccos(carg)) / ((a + b) / 2.0)
		d__1 = delx[s - 1] * dely[s] - dely[s - 1] * delx[s]
		sn = d_sign(c_b7, d__1)
		curvature[s] = mag * sn
	
	curvature[0] = 0.0
	curvature[stnserod - 1] = 0.0
	
	return curvature

def channel(stnserod, stations, flow, slope, diam, width, dels, depth, delx, dely):
	vel = zeros((stnserod,))
	transslope = zeros((stnserod,))
	acs = zeros((stnserod,))
	rightdepth = zeros((stnserod,))
	leftdepth = zeros((stnserod,))
	deln = zeros((stnserod,))
	
	# xmagtransslope is the absolute value of the transverse bedslope (transslope)
	# acs is cross-sectional area of the half-channel at the inside of the bend
	
	curvature = getcurv(stnserod, delx, dely, dels)
	
	for s in range(stnserod):
		if slope[s] <= 0:
			sys.stderr.write("neg. or zero slope: slope[s] = %r, s = %r, flow[s] = %r\n" % (slope[s], s, flow[s]))
		
		if width[s] != 0 and width[s] != depth[s] * -2:
			radh = depth[s]
			if depth[s] <= 0:
				sys.stderr.write('depth[s] = %r, s = %r\n' % (depth[s], s))
				sys.exit(1)
			
			vel[s] = flow[s] / depth[s] / width[s]
			shields = radh * slope[s] / 1.65 / diam[s]
			if diam[s] < 1.5e-4:
				critshields = .08 # s*=1.01 
			elif diam[s] < 2.5e-4:
				critshields = .052 # s*=2.84 
			elif diam[s] < 3.5e-4:
				critshields = .039 # s*=5.22 
			elif diam[s] < 4.5e-4:
				critshields = .035 # s*=8.04 
			elif diam[s] < 5.5e-4:
				critshields = .034 # s*=11.2 
			elif diam[s] < 6.5e-4:
				critshields = .033 # s*=14.8 
			elif diam[s] < 7.5e-4:
				critshields = .034 # s*=18.6 
			elif diam[s] < 8.5e-4:
				critshields = .034 # s*=22.7 
			elif diam[s] < 9.5e-4:
				critshields = .034 # s*=27.1 
			elif diam[s] < .0015:
				critshields = .035 # s*=31.8 
			elif diam[s] < .0025:
				critshields = .045 # s*=89.9 
			elif diam[s] < .0035:
				critshields = .05 # s*=165. 
			elif diam[s] < .0045:
				critshields = .055 # s*=254. 
			else:
				critshields = .056 # s*>=355 
			
			# approximate Engelund diagram for subcritical flow:
			if shields >= .1 and shields < 1.:
				# Computing 2nd power
				d__2 = log10(shields) + 1.03
				d__1 = d__2 * d__2 * .74 - 1.18
				grainshields = pow(10., d__1)
			elif shields >= 1. and shields < 2.:
				# Computing 2nd power
				d__1 = shields
				grainshields = d__1 * d__1 * .4
			else:
				grainshields = shields
			
			corner = depth[s] * 2. / width[s]
			# calculation for skin friction: */
			#print 'shields=%r critshields=%r s=%r diam[s]=%r' % (shields, critshields, s, diam[s])
			try:
				transfactor = depth[s] * (grainshields / shields) * \
					sqrt(grainshields / critshields) * (log(grainshields * 11. * \
					depth[s] / shields / diam[s]) * .5695 - .3606)
			except ZeroDivisionError:
				transfactor = 0
			rectchan = width[s] * .5 * depth[s]
			transslope[s] = transfactor * curvature[s]
			xmagtransslope = abs(transslope[s])
			if xmagtransslope <= corner:
				# Computing 2nd power
				d__1 = width[s]
				acs[s] = rectchan - d__1 * d__1 * xmagtransslope / 8.
				rightdepth[s] = depth[s] + width[s] * transslope[s] / 2.
				leftdepth[s] = depth[s] - width[s] * transslope[s] / 2.
			else:
				# Computing 2nd power
				d__1 = depth[s]
				acs[s] = d__1 * d__1 * .5 / xmagtransslope
				if curvature[s] < 0:
						rightdepth[s] = 0.
						leftdepth[s] = depth[s] * 2.
				else:
						rightdepth[s] = depth[s] * 2.
						leftdepth[s] = 0.
			
			if curvature[s] < 0.:
				deln[s] = width[s] * .5 + acs[s] * 2. / (depth[s] + rightdepth[s])
			else:
				deln[s] = width[s] * .5 + acs[s] * 2. / (depth[s] + leftdepth[s])
	
	return 	(vel, transslope, acs, rightdepth, leftdepth, deln, curvature)

def forcelag(stnserod, stations, rho, vel, depth, width, curvature, acs, dels, deln):
	global c_b7
	latforce = zeros((stnserod,))
	lag = zeros((stnserod,))
	
	# calculate lateral force and downstream lag:
	if stations != stnserod:
		for s in range(stations - 1):
			if width[s] != 0. and width[s] > depth[s] * 2.: 
				# Computing 2nd power
				d__1 = vel[s]
				forcefactor = rho * (d__1 * d__1) / depth[s]
				latforce[s] = 0.
				lag[s] = 0.
				xmagcurvep1 = abs(curvature[s + 1])
				xmagcurve = abs(curvature[s])
				if (curvature[s + 1] * curvature[s] > 0. and xmagcurvep1 > \
					xmagcurve and xmagcurvep1 * width[s] <= 2.) or \
					(curvature[s] == 0. and xmagcurvep1 * width[s] <= 2.):
					delacs = acs[s + 1] - acs[s];
				elif curvature[s + 1] * curvature[s] < 0. and \
					xmagcurvep1 * width[s] <= 2.:
					delacs = acs[s + 1] - depth[s] * width[s] / 2.
				else:
					delacs = 0.
				
				d__1 = -curvature[s + 1]
				latforce[s] = forcefactor * delacs * delacs / dels[s] * \
					d_sign(c_b7, d__1);
				latvel = vel[s] * -1. * delacs / depth[s] / dels[s];
				if (latvel > 0.):
					lag[s] = vel[s] * (deln[s + 1] + deln[s]) / 2. / latvel;
				else:
					latforce[s] = 0.
					lag[s] = 0.
			
			else:
				latforce[s] = 0.
				lag[s] = 0.
	
	else:
		for s in range(stations - 1):
			if width[s] != 0.: 
				# Computing 2nd power
				d__1 = vel[s]
				forcefactor = rho * (d__1 * d__1) / depth[s]
				latforce[s] = 0.
				lag[s] = 0.
				xmagcurvep1 = abs(curvature[s + 1])
				xmagcurve = abs(curvature[s])
				if (curvature[s + 1] * curvature[s] > 0. and xmagcurvep1 > \
					xmagcurve and xmagcurvep1 * width[s] <= 2.) or \
					(curvature[s] == 0. and xmagcurvep1 * width[s] <= 2.):
					delacs = acs[s + 1] - acs[s];
				elif curvature[s + 1] * curvature[s] < 0. and \
					xmagcurvep1 * width[s] <= 2.:
					delacs = acs[s + 1] - depth[s] * width[s] / 2.
				else:
					delacs = 0.
				
				d__1 = -curvature[s + 1]
				latforce[s] = forcefactor * delacs * delacs / dels[s] * \
					d_sign(c_b7, d__1);
				latvel = vel[s] * -1. * delacs / depth[s] / dels[s];
				if (latvel > 0.):
					lag[s] = vel[s] * (deln[s + 1] + deln[s]) / 2. / latvel;
				else:
					latforce[s] = 0.
					lag[s] = 0.
			
			else:
				latforce[s] = 0.
				lag[s] = 0.
				
	return (latforce, lag)

def forcedist(stnserod, stations, lambda_, lag, latforce, phi, curvature, depth, rightdepth, leftdepth, xs):
	tauwall = zeros((stnserod,))
	spreaddelta_x__ = zeros((stnserod,))
	spreaddelta_y__ = zeros((stnserod,))
	
	for s in range(stnserod):
		tauwall[s] = 0.
		spreaddelta_x__[s] = 0.
		spreaddelta_y__[s] = 0.
	
	# For each latforce(s), scan sp = s++ until the
	# distance downstream is greater than xs(s) + lag(s) 
	# + 2*lambda; for each sp, add to tauwall(sp): 
	# latforce(s) * gaussian(xs(s) + lag(s) - xs(sp)) 
	# find deeper bank; calc. 
	# normalized (unit vector divided by downstream 
	# increment and bank depth) lateral direction vectors: 
	for s in range(stations):
		tenlambda = lambda_[s] * 10.
		if tenlambda > xs[stations - 1]:
			tenlambda = xs[stations - 1]
		
		if lag[s] <= tenlambda:
			xdest = xs[s] + lag[s]
			xtrmnt = xdest + lambda_[s] * 2.
			xstrt = xdest - lambda_[s] * 2.
			if xstrt < xs[s]:
				xstrt = xs[s]
			
			sp = s
			while sp < stnserod and xs[sp] <= xtrmnt:
				if xs[sp] >= xstrt:
					xdel = abs(xdest - xs[sp])
					if lambda_[s] != 0.:
						# Computing 2nd power
						d__1 = xdel
						# Computing 2nd power
						d__2 = lambda_[s]
						gaussian = exp(d__1 * d__1 * -1. / 2. / (d__2 * d__2)) \
							 / sqrt(2. * 3.1416) / lambda_[s]
						gaussfactor = gaussian * latforce[s]
					else:
						gaussfactor = 0.
					
					tauwall[sp] += gaussfactor
				sp += 1
	
	for s in range(stnserod):
		if depth[s] != 0.:
			tauwall[s] /= depth[s]
			spreaddelta_x__[s] = tauwall[s] * -1. * sin(phi[s])
			spreaddelta_y__[s] = tauwall[s] * cos(phi[s])
		else:
			tauwall[s] = 0.
			spreaddelta_x__[s] = 0.
			spreaddelta_y__[s] = 0.
			
	return (spreaddelta_x__, spreaddelta_y__)

# This subroutine moves the meander nodes. 
# xcp is the cross-product of the flow direction and the direction of 
# bank erosion (perpendicular to the flow direction); the sign tells 
# you whether to use right or left bank erodibility. 

def changeposition(stnserod, lerody, rerody, spreaddelta_x__, spreaddelta_y__, delx, dely):
	delta_x__ = zeros((stnserod,))
	delta_y__ = zeros((stnserod,))
	
	for s in range(stnserod):
		xcp = delx[s] * spreaddelta_y__[s] - spreaddelta_x__[s] * dely[s]
		if xcp < 0.:
			delta_x__[s] = rerody[s] * spreaddelta_x__[s]
			delta_y__[s] = rerody[s] * spreaddelta_y__[s]
		else:
			delta_x__[s] = lerody[s] * spreaddelta_x__[s]
			delta_y__[s] = lerody[s] * spreaddelta_y__[s]
			
	return (delta_x__, delta_y__)

def meander(time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_):
	
	# declare parameter values; get initial channel config.:
	rho, grav, phi, delx, dely = initialize(stnserod, x, y)
	
	# define channel characteristics:
	vel, transslope, acs, rightdepth, leftdepth, deln, curvature = \
		channel(stnserod, stations, flow, slope, diam, width, dels, depth, delx, dely)
	
	# calculate lateral force and downstream lag:
	latforce, lag = forcelag(stnserod, stations, rho, vel, depth, width, curvature, acs, dels, deln)
	
	# propagate and smooth lateral force:
	spreaddelta_x__, spreaddelta_y__ = \
		forcedist(stnserod, stations, lambda_, lag, latforce, phi, curvature, depth, rightdepth, leftdepth, xs)
	
	# change channel position:
	delta_x__, delta_y__ = changeposition(stnserod, lerody, rerody, spreaddelta_x__, spreaddelta_y__, delx, dely)
	
	# adjust for fortran indexing
	delta_x__ = roll(delta_x__, 1)
	delta_y__ = roll(delta_y__, 1)
	rightdepth = roll(rightdepth, 1)
	leftdepth = roll(leftdepth, 1)
	delta_x__[0] = 0.0
	delta_y__[0] = 0.0
	rightdepth[0] = 0.0
	leftdepth[0] = 0.0
	
	return (delta_x__, delta_y__, rightdepth, leftdepth)

def encodeData(delx, dely, rightdepth, leftdepth):
	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('MeanderTS.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:
		time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_ = decodeData()
		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])
		delx, dely, rightdepth, leftdepth = \
			meander(time, stations, stnserod, x, y, xs, dels, flow, rerody, lerody, slope, width, depth, diam, lambda_)
		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])
		encodeData(delx, dely, rightdepth, leftdepth)