# plots petrosian functions for sersic profiles 
# as per Fig 6 of
# Graham and Driver (astro-ph/0503176)
#
import numpy, os, math, ppgplot, sys

def get_b_n(m):
#   find sersic b_n coefficient
#   this is more accurate than the usual b_n = 1.9992*n_sersic - 0.3271
    if m == 0.0: return -0.3271
    b_n = 2.0*m - 1.0/3.0 + 4.0/m/405.0 + 46.0/m/m/25515.0 \
                + 131.0/m/m/m/1148175.0  \
                - 2194697.0/m/m/m/m/30690717750.0
    return b_n

def gammln(xx):
#  Logarithm of the gamma function.
        gammln_cof = numpy.array([76.18009173, -86.50532033, 24.01409822,
	                          -1.231739516e0, 0.120858003e-2, -0.536382e-5])
        gammln_stp = 2.50662827465
	x = xx - 1.
	tmp = x + 5.5
	tmp = (x + 0.5)*math.log(tmp) - tmp
	ser = 1.
	for j in range(6):
		x = x + 1.
		ser = ser + gammln_cof[j]/x
	return tmp + math.log(gammln_stp*ser)

def gser(a, x, itmax=700, eps=3.e-7):
	"""Series approx'n to the incomplete gamma function."""
	gln = gammln(a)
	if (x < 0.):
		raise bad_arg, x
	if (x == 0.):
		return(0.)
	ap = a
	sum = 1. / a
	delta = sum
	n = 1
	while n <= itmax:
		ap = ap + 1.
		delta = delta * x / ap
		sum = sum + delta
		if (abs(delta) < abs(sum)*eps):
			return (sum * math.exp(-x + a*math.log(x) - gln), gln)
		n = n + 1
	raise max_iters, str((abs(delta), abs(sum)*eps))


def gcf(a, x, itmax=200, eps=3.e-7):
#"""Continued fraction approx'n of the incomplete gamma function."""
    gln = gammln(a)
    gold = 0.
    a0 = 1.
    a1 = x
    b0 = 0.
    b1 = 1.
    fac = 1.
    n = 1
    while n <= itmax:
        an = n
	ana = an - a
	a0 = (a1 + a0*ana)*fac
	b0 = (b1 + b0*ana)*fac
	anf = an*fac
	a1 = x*a0 + anf*a1
	b1 = x*b0 + anf*b1
	if (a1 != 0.):
           fac = 1. / a1
	   g = b1*fac
	   if (abs((g-gold)/g) < eps):
              return (g*math.exp(-x+a*math.log(x)-gln), gln)
	   gold = g
	   n = n + 1
    raise max_iters, str(abs((g-gold)/g))

def gammp(a, x):
	"""Incomplete gamma function."""
	if (x < 0. or a <= 0.):
		raise ValueError, (a, x)
	if (x < a+1.):
#                print gser(a,x)
		return gser(a,x)[0]
	else:
		return 1.-gcf(a,x)[0]

def plot_curve(n):
    one_over_n = 1./n
    b_n = get_b_n(n)
    two_n = n*2.0
    for ix in range(1,401):
        xx = float(ix)/100.
        xxx = b_n*xx**one_over_n
        incgam =  gammp(two_n, xxx)*math.exp(gammln(two_n))
        eta = two_n*incgam/math.exp(-1.*xxx)/(xxx**two_n)
        yy = 2.5*math.log10(eta)
#    print ix, xx, incgam,  eta, yy
        if ix == 1: ppgplot.pgmove(xx,yy)                                
        ppgplot.pgdraw(xx,yy)  

os.environ["PGPLOT_FOREGROUND"] = "black"
os.environ["PGPLOT_BACKGROUND"] = "white"
os.environ["PGPLOT_FONT"] = "/usr/local/pgplot/grfont.dat"

ppgplot.pgopen('petro.ps/ps')                         
ppgplot.pgvsiz(1.0,6.0,1.0,4.0)                        
ppgplot.pgslw(2)                                       
ppgplot.pgswin(0.0,4.0,0.0,3.0)                     
ppgplot.pgbox('bcnts',0.0,0,'bcnvts',0.0,0)            
ppgplot.pglab('R/R\de\u','2.5 log \gy [mag arcsec \u-2\d]',' ')  
ppgplot.pgslw(1) 

plot_curve(0.5)
plot_curve(1.0)
plot_curve(2.0)
plot_curve(3.0)
plot_curve(4.0)
plot_curve(5.0)
plot_curve(6.0)
plot_curve(7.0)
plot_curve(8.0)
plot_curve(9.0)
plot_curve(10.0)

ppgplot.pgclos()