GammaLib

#! /usr/bin/env python

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

# Check BVC file

#

# Copyright (C) 2021 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 gammalib

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

# Check BVC file #

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

def check_bvc(filename):

    """

    Check BVC file

    """

    # Open BVC file

    fits = gammalib.GFits(filename)

    # Get BVC table

    bvc = fits.table('COMPTEL_BVC')

    # Extract information

    for irow in range(bvc.nrows()):

        # Read information

        tjd    = bvc['TJD'][irow]

        tics   = bvc['TICS'][irow]

        ssb_x  = bvc['SSB_X'][irow]

        ssb_y  = bvc['SSB_Y'][irow]

        ssb_z  = bvc['SSB_Z'][irow]

        tdelta = bvc['TDELTA'][irow]

        #print(tjd,tics,ssb_x,ssb_y,ssb_z,tdelta)

        # Convert into sky direction

        time = gammalib.com_time(tjd,tics)            # Get time of row

        ssb  = gammalib.GVector(ssb_x, ssb_y, ssb_z)  # Setup vector pointing from SSB -> CGRO

        ssb /= ssb.norm()                             # Normalise vector to unity

        ssb  = -ssb                                   # Invert vector to point from CGRO -> SSB

        dir  = gammalib.GSkyDir()                     #

        dir.celvector(ssb)                            # Interpret vector in celestial coordinates

        sun  = gammalib.GSkyDir()

        sun.sun(time)                                 # Get coordinates of Sun

        # Show result

        print('%4.4d:%9.9d: (%13.2f,%13.2f,%13.2f) %10.7f  CGRO->SSB(%9.4f,%9.4f)  Sun(%9.4f,%9.4f)' %

              (tjd,tics,ssb_x,ssb_y,ssb_z,tdelta,dir.ra_deg(),dir.dec_deg(),sun.ra_deg(),sun.dec_deg()))

    # Return

    return

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Set filename

    filename = '/project-data/comptel/data/phase01/vp0001_0/s10150_bvc.fits'

    # Check BVC file

    check_bvc(filename)

    # Create figure

    #fig = plt.figure(figsize=(12,8))

    # Show figure

    #plt.show()