GammaLib

#! /usr/bin/env python

# ===========================================================================================#

# This script checks the King profile integration

# ===========================================================================================#

import math

import matplotlib.pyplot as plt

# ============================= #

# Show King profile integration #

# ============================= #

def show(sigma, gamma, color):

    """

    Show photons using matplotlib (if available).

    """

    # Set binning

    num   = 1000

    delta = 0.1

    # Create r axis

    r = []

    for i in range(num):

        r.append(i*delta)

    # Create integral curve

    integral = []

    sigma2   = sigma*sigma

    for i in range(num):

        u     = 0.5 * r[i] * r[i] / sigma2

        value = 1.0 - math.pow(1.0+u/gamma, 1.0-gamma)

        integral.append(value)

    # Compute 68% r value

    F    = 0.68

    u_68 = gamma * (math.pow(1.0-F, 1.0/(1.0-gamma)) - 1.0)

    r_68 = sigma * math.sqrt(2.0 * u_68)

    print r_68, u_68

    # Show integral

    plt.plot(r, integral, color)

    plt.plot([r_68, r_68], [0.0, 0.68], color)

    plt.plot([0.0, r_68], [0.68, 0.68], color)

    plt.xlabel("r")

    plt.ylabel("Integral")

    # Return

    return

# ======================= #

# Show slope versus gamma #

# ======================= #

def slope(F, color):

    """

    Show slope versus gamma.

    """

    # Set binning

    gmin  = 1.3

    gmax  = 2.5

    delta = 0.05

    num   = int((gmax-gmin)/delta)+2

    # Create gamma axis

    g = []

    for i in range(num):

        g.append(gmin+i*delta)

    # Compute slope

    slope = []

    for i in range(num):

        umax = g[i] * (math.pow(1.0-F, 1.0/(1.0-g[i])) - 1.0)

        s    = 1.0/math.sqrt(2.0 * umax)

        slope.append(s)

    # Show integral

    plt.plot(g, slope, color)

    plt.xlabel("Gamma")

    plt.ylabel("Slope")

    # Return

    return

#==========================#

# Main routine entry point #

#==========================#

if __name__ == '__main__':

    """

    Compute integral

    """

    plt.figure(1)

    #show(0.5, 2.0, "r-")

    #show(1.0, 1.5, "b-.")

    #show(1.0, 2.0, "b-")

    #show(1.0, 2.5, "b--")

    #show(2.0, 2.0, "g-")

    show(1.0, 1.2, "r-")

    #slope(0.68, "ro")

    plt.show()