#! /usr/bin/env python3.2
"""
The appropriate b value in the Sersic profile
For very low valued Sersic Index values.
Sersic Profile: 
    I=I_e*exp(-b*((r/r_e)^(1/n)-1))
Where:
      I_e: Intensity at the effective radius
      r_e: Effective radius (containing half of total light)
        n: The Sersic index
        b: A function of the Sersic Index to ensure 
           the definition of the effective radius

More explanations can be found on:
http://astr.tohoku.ac.jp/~akhlghi/Sersic_C.html


####################################################### 
#####              Mohammad Akhlaghi              #####
#####  Tohoku University Astronomical Institute   #####
#####     http://astr.tohoku.ac.jp/~akhlaghi/     #####
#####               April 27th, 2012              #####
#######################################################
"""

import numpy as np
from scipy import integrate

def b_integrand(r, re, n, C, b):
    """
    This is used in the integration to find b
    """
    return r*Sersic(r, re, n, C, b)

def Sersic(r, re, n, C, b):
    """
    Sersic Function for integration
    """
    ra=np.power((r/re), 1/n)
    return C*np.exp(-b*(ra-1))


#Limits and step of Sersic index:
low_limit=0.05
high_limit=0.06
step=0.001

#This is the value of b for n=0.36
b_test=0.426200378468
#rate of increasing b
b_drate=1.002
#rate of decreasing b
b_irate=1.001
#Effective radius:
re=2

#initiating the process
table_row=np.zeros((1,2))
out_file=open("b_of_nTable", 'w')
for n in reversed(np.arange(low_limit, high_limit, step)):
    while True:
        a=integrate.quad(b_integrand, 0, float('inf'), 
                       args=(re, n, 1, b_test))[0]
        b=integrate.quad(b_integrand, 0, re, 
                       args=(re, n, 1, b_test))[0]
        I_diff=2*np.pi*(a-(2*b))
        if I_diff>-0.00001 and I_diff<0.00001:
            break
        elif I_diff<=-0.00001:
             b_test=b_test/b_drate
        elif I_diff>=0.00001:
            b_test=b_test*b_irate
    table_row[0,0]=n
    table_row[0,1]=b_test 
    print("n={} --> b(n)={}".format(table_row[0,0], table_row[0,1]))
    print("{}\t\t{}".format(table_row[0,0], table_row[0,1]), file=out_file)
out_file.close()