#!/usr/bin/env python

# ----------------------------------------------------------------
# John Kerl
# kerl.john.r@gmail.com
# 2006-05-19
#
# Contour integral around a (closed) curve C of f(z)/(z-z0)^n dz.
#
# Example:
#
# for f in 1 sin cos exp log; do
#   for p in 0 1 2 3 4; do
#     echo -n "$f $p: "
#     numintc1 -f $f -p $p
#   done
#   echo ""
# done
# ----------------------------------------------------------------
# This software is released under the terms of the GNU GPL.
# Please see LICENSE.txt in the same directory as this file.
# ================================================================

from __future__ import division # 1/2 = 0.5, not 0.
# Import the cmath module rather than the math module so that exp et al. may
# take complex arguments.
from cmath import *
import sys

# ----------------------------------------------------------------
# Path
def pathz(t, z0):
	R = 0.01
	return z0 + R * exp(1j * t)

# ----------------------------------------------------------------
# Non-library-function integrands.  
def one_func(z):
	return 1.0

# ----------------------------------------------------------------
def usage():
	print >> sys.stderr, "Usage: %s [options]" % (sys.argv[0])
	print >> sys.stderr, "Integrates w(z) = f(z)/(z-z0)^n around a small circle centered at z0."
	print >> sys.stderr, "Options:"
	print >> sys.stderr, "  -t:         Tabulate t, z, dz, w, and accumulated sum over the mash."
	print >> sys.stderr, "  -nt {n}:    Specify the number of mesh points."
	print >> sys.stderr, "  -z0 {z0}:   Specify z0 in 1/(z-z0)^n in the integrand."
	print >> sys.stderr, "  -p  {n}:    Specify exponent n in 1/(z-z0)^n in the integrand."
	print >> sys.stderr, "  -f  {name}: Specify f(z) of integrand: 1, sin, cos, exp, or log."
	sys.exit(1)

# ----------------------------------------------------------------
# Main script

# Integrand parameters
z0 = 0.0
power = 1
f = one_func

# Mesh over path
tlo = 0.0
thi = 2.0 * pi
nt  = 10000
tabulate = 0

# Command-line parsing
argc = len(sys.argv)
argi = 1
while (argi < argc):
	if (sys.argv[argi] == '-t'):
		tabulate = 1
	elif (sys.argv[argi] == '-nt'):
		if ((argc-argi) < 2):
			usage()
		argi += 1
		nt = int(sys.argv[argi])
	elif (sys.argv[argi] == '-p'):
		if ((argc-argi) < 2):
			usage()
		argi += 1
		power = int(sys.argv[argi])
	elif (sys.argv[argi] == '-z0'):
		if ((argc-argi) < 2):
			usage()
		argi += 1
		z0 = complex(sys.argv[argi])
	elif (sys.argv[argi] == '-f'):
		if ((argc-argi) < 2):
			usage()
		argi += 1
		if   (sys.argv[argi] == '1'):
			f = one_func
		elif (sys.argv[argi] == 'one'):
			f = one_func
		elif (sys.argv[argi] == 'sin'):
			f = sin
		elif (sys.argv[argi] == 'cos'):
			f = cos
		elif (sys.argv[argi] == 'exp'):
			f = exp
		elif (sys.argv[argi] == 'log'):
			f = log
		else:
			usage()
	else:
		usage()
	argi += 1

# Start of numerical approximation to the integral
sum = 0.0
t   = tlo
dt  = (thi - tlo) / nt
for it in range(0, nt):
	t1 = t
	t2 = t + dt

	z  = pathz(t1, z0)
	dz = pathz(t2, z0) - z

	w = f(z) / ((z - z0) ** power)
	sum += w * dz

	if tabulate:
		print ">> t = %8.4f/2pi z = %8.4f +%8.4fi dz = %8.4f +%8.4fi w = %8.4f +%8.4fi sum = %8.4f +%8.4fi" % \
			(t/2/pi, z.real, z.imag, dz.real, dz.imag, w.real, w.imag, sum.real, sum.imag)

	t += dt

print "%11.7f +%11.7fi" % (sum.real, sum.imag)