#!/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, "  A=[number]"
	print >> sys.stderr, "  B=[number]"
	sys.exit(1)

# ----------------------------------------------------------------
# Brownian motion on an ellipse.
# SDE:
#   dX_t = -1/2 X dt - A/B Y_t db_t
#   dY_t = -1/2 Y dt + B/A X_t db_t
# Solution:
#   X_t  = A cos(b_t)
#   Y_t  = B sin(b_t).

T  = 2.0
nt = 20000
A  = 2.0
B  = 2.0

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'^A=',  sys.argv[argi]): A  = float(sys.argv[argi][2:])
	elif re.match(r'^B=',  sys.argv[argi]): B  = float(sys.argv[argi][2:])
	else: usage()

dt = T/nt
t  = 0.0
Bt = 0.0
sqrtdt = sqrt(dt)

Xapprox = A * cos(Bt)
Yapprox = B * sin(Bt)

while (t <= T):
	Xexact = A * cos(Bt)
	Yexact = B * sin(Bt)
	print "%11.7f %11.7f %11.7f %11.7f %11.7f %9.3e %9.3e" % \
		(t, Xexact, Xapprox, Xexact-Xapprox, Yexact, Yapprox, Yexact-Yapprox)

	dB  = random.gauss(0, sqrtdt)

	Xapprox += -0.5 * Xapprox * dt - A/B * Yapprox * dB
	Yapprox += -0.5 * Yapprox * dt + B/A * Xapprox * dB

	Bt += dB
	t  += dt