#!/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, "  c=[number]"
	print >> sys.stderr, "  nX=[number]"
	sys.exit(1)
	
# ----------------------------------------------------------------
def get_mean(X):
	n = len(X)
	sum = 0.0
	for i in range(0, n):
		sum += X[i]
	return sum / n

# ----------------------------------------------------------------
# X_t = (b_t + t) exp(-b_t - t/2)

T  = 2
nt = 1000
c  = 0.5
nX = 40

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

dt = T/nt
sqrtdt = sqrt(dt)
t  = 0

# ----------------------------------------------------------------
Bt = [0.0] * nX
X  = [0.0] * nX

while (t <= T):
	for k in range(0, nX):
		X[k] = (Bt[k] + t) * exp(-Bt[k] - c*t)

	print "%11.7f" % (t),
	for k in range(0, nX):
		print "%11.7f" % (X[k]),
	m = get_mean(X)
	print "%11.7f" % (m)

	for k in range(0, nX):
		dB  = random.gauss(0, sqrtdt)
		Bt[k] += dB
	t  += dt