accretion_disk_spectrum

def main():

    import sympy
    import numpy
    import pylab
    import math

    f = sympy.Symbol('f')
    c1 = 1.
    c2 = 1e3
    f_range = numpy.logspace(-6,2,100)
    xi = sympy.Symbol('xi')
    beta = 3./4.
    integrand = lambda x: 0 if x>50 else (2./beta)*x**(2./beta-1)/(math.exp(x)-1)
    F_list = [f**(3.-2./beta)*sympy.mpmath.quad(integrand,[f*c1,f*c2]) for f in f_range]
    F_list = numpy.array([float(x) for x in F_list])
    f_mid = 0.5*(f_range[1:]+f_range[:-1])
    dlnFdlnf_list = numpy.diff(numpy.log(F_list))/numpy.diff(numpy.log(f_range))
    pylab.subplot(211)
    pylab.loglog(f_range, F_list)
    pylab.ylim((0.5*F_list[0],2*numpy.max(F_list)))
    pylab.title('Spectrum of an accretion disk')
    pylab.ylabel(r'$F_f$')
    pylab.subplot(212)
    pylab.semilogx(f_mid, dlnFdlnf_list)
    pylab.xlabel('f')
    pylab.ylabel(r'$d \ln F_f / d \ln f$')
    pylab.ylim((0,2.1))
    pylab.show()

if __name__ == '__main__':

    main()