ctools

#! /usr/bin/env python

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

# This script computes a visibility cube

#

# Copyright (C) 2016 Juergen Knoedlseder

#

# This program is free software: you can redistribute it and/or modify

# it under the terms of the GNU General Public License as published by

# the Free Software Foundation, either version 3 of the License, or

# (at your option) any later version.

#

# This program is distributed in the hope that it will be useful,

# but WITHOUT ANY WARRANTY; without even the implied warranty of

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

# GNU General Public License for more details.

#

# You should have received a copy of the GNU General Public License

# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#

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

import math

import gammalib

import matplotlib.pyplot as plt

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

# Plot visibility #

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

def plot_visibility(ra):

    """

    Plot visibility

    """

    # Load visibility cube

    cube = gammalib.GSkyMap(file)

    # Create figure

    plt.figure(1)

    plt.title(title)

    # Plot visibility

    #decs = [float(i*10)-90.0 for i in range(18)]

    decs = [-64.4, -54.4, -44.4, -34.4, -24.4, -14.4, -4.4, 5.6, 15.6]

    for dec in decs:

        # Get RA,DEC index

        index = int(ra) + int((dec+90.0))*360

        # Extract vectors

        zeniths    = [float(i) for i in range(cube.nmaps())]

        visibility = [cube[index, i] for i in range(cube.nmaps())]

        print(dec, sum(visibility))

        # Plot visibility

        plt.plot(zeniths, visibility, 'r-')

    # Set axes

    plt.xlabel('Zenith angle (degrees)')

    plt.ylabel('Visibility (hours)')

    # Return

    return

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

# Plot zenith angle #

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

def plot_zenith_angle(ra):

    """

    Plot zenith angle

    """

    # Load visibility cube

    cube = gammalib.GSkyMap(file)

    # Create figure

    plt.figure(2)

    plt.title(title)

    # Plot visibility

    zeniths    = [float(i) for i in range(60)]

    visibility = [0.0 for i in range(60)]

    for izenith, zenith in enumerate(zeniths):

        # Extract vectors

        visibility[izenith] = sum([cube[i,izenith] for i in range(cube.npix())])

    # Plot visibility

    plt.plot(zeniths, visibility, 'r-')

    # Set axes

    plt.xlabel('Zenith (degrees)')

    plt.ylabel('Visibility * Solid angle (hours deg$^2$)')

    plt.xlim([0,60])

    # Return

    return

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

# Plot declination #

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

def plot_declination(ra):

    """

    """

    # Load visibility cube

    cube = gammalib.GSkyMap(file)

    # Create figure

    plt.figure(3)

    plt.title(title)

    # Plot visibility

    decs       = [float(i)-89.5 for i in range(180)]

    visibility = [0.0 for i in range(180)]

    for idec, dec in enumerate(decs):

        # Get RA,DEC index

        index = int(ra) + idec*360

        # Extract vectors

        visibility[idec] = sum([cube[index, i] for i in range(cube.nmaps())])

    # Plot visibility

    plt.plot(decs, visibility, 'r-')

    # Set axes

    plt.xlabel('Declination (degrees)')

    plt.ylabel('Visibility (hours)')

    plt.xlim([-90,90])

    # Return

    return

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

# Plot minimum zenith angle #

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

def plot_min_zenith(ra):

    """

    """

    # Load visibility cube

    cube = gammalib.GSkyMap(file)

    # Create figure

    plt.figure(4)

    plt.title(title)

    # Determine minimum zenith angle for each declination

    decs = [float(i)-89.5 for i in range(180)]

    x    = []

    y    = []

    for idec, dec in enumerate(decs):

        # Get RA,DEC index

        index = int(ra) + idec*360

        # Search first non-zero zenith angle

        for iz in range(cube.nmaps()):

            if cube[index, iz] > 0.0:

                x.append(dec)

                y.append(float(iz))

                break

    # Plot visibility

    plt.plot(x, y, 'r-')

    # Set axes

    plt.xlabel('Declination (degrees)')

    plt.ylabel('Minimum zenith angle (deg)')

    plt.xlim([-90,90])

    # Return

    return

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Set parameters

    #file  = 'viscube_north.fits'

    #title = 'CTA North'

    file  = 'viscube_south.fits'

    title = 'CTA South'

    ra    = 220.0

    # Plot zenith angle

    plot_visibility(ra)

    plot_zenith_angle(ra)

    plot_declination(ra)

    plot_min_zenith(ra)

    # Show plot

    plt.show()