# plots surface brightness profiles as per Fig 1 (top panel) of
# Graham and Driver (astro-ph/0503176) plus a random profile (dashed line)
#
import numpy, os, math, ppgplot

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 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 plot_it(r_sersic,sb_e,n_sersic):
  b_n = get_b_n(n_sersic)
  n2 = n_sersic*2.0
  f_n = n_sersic*math.exp(b_n)*math.exp(gammln(n2))/(b_n**n2)
  con2 = math.log(10.)
  for i in range(801):
      r = float(i)/100.0
      sb = sb_e + 2.5*b_n*((r/r_sersic)**(1./n_sersic) - 1.)/con2
      if i == 0: ppgplot.pgmove(r,sb)                                
      ppgplot.pgdraw(r,sb)                                

ppgplot.pgopen('prof.ps/ps')                         
ppgplot.pgvsiz(1.0,6.0,1.0,4.0)                        
ppgplot.pgslw(2)                                       
ppgplot.pgswin(0.0,8.0,28.0,12.0)                     
ppgplot.pgbox('bcnts',0.0,0,'bcnvts',0.0,0)            
ppgplot.pglab('R/R\de\u','\gm [mag arcsec\u-2\d]',' ')  
ppgplot.pgslw(1) 

sb_e = 20.0
r_sersic = 1.0
n_sersic = 10.0
plot_it(r_sersic,sb_e,n_sersic)
n_sersic = 4.0
plot_it(r_sersic,sb_e,n_sersic)
n_sersic = 2.0
plot_it(r_sersic,sb_e,n_sersic)
n_sersic = 1.0
plot_it(r_sersic,sb_e,n_sersic)
n_sersic = 0.5
plot_it(r_sersic,sb_e,n_sersic)


# for a galaxy with
tot_mag  =  16.0
r_sersic = 10.0
n_sersic = 4.0

const = 2.5*math.log10(2.*math.pi)
sb_av = tot_mag + const + 5.0*math.log10(r_sersic) 
b_n = get_b_n(n_sersic)
n2 = n_sersic*2.0
f_n = n_sersic*math.exp(b_n)*math.exp(gammln(n2))/(b_n**n2)
sb_e = sb_av + 2.5*math.log10(f_n)
ppgplot.pgsls(2) 
plot_it(r_sersic,sb_e,n_sersic)


ppgplot.pgclos()