#!/usr/bin/python -Wall

# ================================================================
# John Kerl
# kerl at math dot arizona dot edu
# 2008-05-12
# ================================================================

from   __future__ import division # 1/2 = 0.5, not 0.
import random
from   math import *
import copy
import sys
import re

# ----------------------------------------------------------------
def usage():
	print >> sys.stderr, "Usage: %s [options]" % (sys.argv[0])
	print >> sys.stderr, "Options:"
	print >> sys.stderr, "  T=[number]"
	print >> sys.stderr, "  nt=[number]"
	print >> sys.stderr, "  X0=[number]"
	print >> sys.stderr, "  Y0=[number]"
	print >> sys.stderr, "  alpha=[number]"
	print >> sys.stderr, "  beta=[number]"
	sys.exit(1)

# ----------------------------------------------------------------
# Stochastically forced vibrating string equation:
#   dX =  Y dt + alpha da
#   dY = -X dt + beta  db
# where a and b are independent Brownian motions.

T  = 100
nt = 10000
X0 =  1
Y0 =  0
alpha = 0.200
beta  = 0.200

argc = len(sys.argv)
for argi in range(1, argc):
	if   re.match(r'^T=',     sys.argv[argi]): T     = float(sys.argv[argi][2:])
	elif re.match(r'^nt=',    sys.argv[argi]): nt    = int  (sys.argv[argi][3:])
	elif re.match(r'^X0=',    sys.argv[argi]): X0    = float(sys.argv[argi][3:])
	elif re.match(r'^Y0=',    sys.argv[argi]): Y0    = float(sys.argv[argi][3:])
	elif re.match(r'^alpha=', sys.argv[argi]): alpha = float(sys.argv[argi][6:])
	elif re.match(r'^beta=',  sys.argv[argi]): beta  = float(sys.argv[argi][5:])
	else: usage()

t  =  0
dt =  T/nt
At = 0.0
Bt = 0.0
X  = X0
Y  = Y0
sqrtdt = sqrt(dt)

while (t <= T):
	print "%11.7f %11.7f %11.7f" % (t, X, Y)
	dA  = random.gauss(0, sqrtdt)
	dB  = random.gauss(0, sqrtdt)
	X += Y * dt + alpha * dA
	Y -= X * dt + beta  * dB
	At += dA
	Bt += dB
	t  += dt