import matplotlib.pyplot as plt
%pylab inline
%matplotlib inline
Populating the interactive namespace from numpy and matplotlib
xtmp = np.arange(-4.,4.,0.05)
plt.semilogy(xtmp,N1d(xtmp,0.,1.))
print intT1d(xtmp,N1d(xtmp,0.,1.)), intT1d(xtmp,N1d(xtmp,0.,0.2**-2)), intT1d(xtmp,N1d(xtmp,0.,2.**-2))
mdev10={0.00139:0.00087, \
0.00941:0.00296, \
0.04441:0.00792, \
0.16162:0.01902, \
0.48121:0.04289, \
1.20357:0.09351, \
2.54182:0.20168, \
4.46441:0.44126, \
6.22820:1.01833, \
6.15393:2.74555}
ss = 0
for k,v in mdev10.iteritems():
ss += k
print ss, k
print ss
0.999929007421 1.0 0.95310433342 4.46441 4.46441 7.00623 2.54182 8.2098 1.20357 8.21921 0.00941 8.38083 0.16162 8.38222 0.00139 14.61042 6.2282 14.65483 0.04441 20.80876 6.15393 21.28997 0.48121 21.28997
def N1d(x,m,sig2inv):
"""
weird input: sig2inv to follow the 2d Gaussian definition.
sig2inv = 1./sig^2 where sig2 is the dispersion.
"""
exparg = (x-m)**2*sig2inv
return (0.5*sig2inv/np.pi)**0.5*np.exp(-0.5*exparg)
def deVauc(r,Ie,re,aopt=1):
"""
aopt = 0: straight-up deVauc formula.
aopt = 1: SDSS approximation ('luv') Eqn 17-18 of Hogg and Lang.
aopt = 2: Hogg and Lang approximation to luv (sum over 10 Gaussians).
"""
xi = np.fabs(r/re)
if aopt == 0:
return Ie*np.exp(-7.66924944*((xi)**(0.25)-1))
if aopt == 1:
result = np.zeros(len(r))
result = Ie*np.exp(-7.66925*(((xi)**2 + 0.0004)**(0.125)-1))
xx = np.where(xi > 8.)[0]
if len(xx) > 0:
result[xx] = 0.
xx = np.where((xi >= 7.) & (xi < 8.))[0]
if len(xx) > 0.:
result[xx] = result[xx] * (1.-(xi[xx]-7)**2)**2
return result
if aopt == 2:
## hard code Hogg/Lang coefficients.
mdev10={0.00139:0.00087, \
0.00941:0.00296, \
0.04441:0.00792, \
0.16162:0.01902, \
0.48121:0.04289, \
1.20357:0.09351, \
2.54182:0.20168, \
4.46441:0.44126, \
6.22820:1.01833, \
6.15393:2.74555}
result = np.zeros(len(r))
xx = np.where(xi > 8.)[0]
result += N1d(r,)
if len(xx) > 0:
result[xx] = 0.
def intT1d(r,pofr):
"""
Trapezoid rule for integration.
"""
return 0.5*((r[1:]-r[:-1])*(pofr[1:]+pofr[:-1])).sum()
## Plot of the deVauc and LUV profiles. Can you add the plot of the Gaussian approximation?
r=10**np.arange(-3,2,0.001)
print len(r)
# plt.loglog(r,deVauc(r,1.,1.,0),'b')
plt.loglog(r,deVauc(r,1.,1.,0),'b')
plt.loglog(r,deVauc(r,1.,1.,1),'g')
# plt.loglog(r,deVauc(r,1.,1.,2),'r')
5000
[<matplotlib.lines.Line2D at 0x10898b510>]
x=r
y=deVauc(r,1.,1.,0)*2*np.pi*x
y2=deVauc(r,1.,1.,1)*2*np.pi*x
ihalf = np.where(r > 1.)[0].min()
print 'half light radius correct?', x[ihalf], intT1d(x[:ihalf],y[:ihalf])/intT1d(x,y)
print 'half light approx correct for LUV profile?',intT1d(x[:ihalf],y2[:ihalf])/intT1d(x,y2)
half light radius correct? 1.00230523808 0.499975031558 half light approx correct for LUV profile? 0.52884049292