GammaLib

#! /usr/bin/env python

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

# Plot difference between GSkyDir.sun() results and PyEphem results

#

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

import ephem

import matplotlib.pyplot as plt

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

# Compute Sun position #

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

def compute_difference(time):

    """

    """

    # Compute ephem position

    sun_ephem     = ephem.Sun(time.datetime())

    sun_ephem_ra  = sun_ephem.a_ra*gammalib.rad2deg

    sun_ephem_dec = sun_ephem.a_dec*gammalib.rad2deg

    # Compute ctools position

    sun_ctools=gammalib.GSkyDir()

    sun_ctools.sun(time)

    sun_ctools_ra  = sun_ctools.ra_deg()

    sun_ctools_dec = sun_ctools.dec_deg()

    if sun_ctools_ra < 0.0:

        sun_ctools_ra += 360.0

    # Compute difference

    ra_diff  = sun_ephem_ra  - sun_ctools_ra

    dec_diff = sun_ephem_dec - sun_ctools_dec

    if ra_diff > 180.0:

        ra_diff -= 360.0

    if ra_diff < -180.0:

        ra_diff += 360.0

    # Print results

    #print(sun_ephem_ra, sun_ephem_dec)

    #print(sun_ctools_ra, sun_ctools_dec)

    #print(ra_diff, dec_diff)

    # Return

    return ra_diff, dec_diff

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

# Get difference #

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

def get_difference(dt=50000.0, nt=31560):

    """

    """

    # Initialise lists

    list_mjd  = []

    list_dra  = []

    list_ddec = []

    # Set start time

    time = gammalib.GTime()

    time.utc('2000-01-01T00:00:00')

    print('Start time .: %s' % time.utc())

    # Loop over times

    for i in range(nt):

        # Compute difference

        ra_diff, dec_diff = compute_difference(time)

        # Print difference

        #print('%s (%f): %f %f' % (time.utc(), time.mjd(), ra_diff, dec_diff))

        # Append result to lists

        list_mjd.append(time.mjd())

        list_dra.append(ra_diff)

        list_ddec.append(dec_diff)

        # Increment time

        time += dt

    # Print stop time

    print('Stop time ..: %s' % time.utc())

    # Return

    return list_mjd, list_dra, list_ddec

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

# Plot difference #

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

def plot_difference(list_mjd, list_dra, list_ddec):

    """

    """

    # Create figure

    fig = plt.figure('Sun position difference', (12, 4))

    fig.subplots_adjust(left=0.12, right=0.98, top=0.94, bottom=0.14)

    ax = plt.subplot()

    # Plot the differeces

    ax.plot(list_mjd, list_dra,  'b-', label='Right Ascension')

    ax.plot(list_mjd, list_ddec, 'r-', label='Declination')

    # Set labels

    ax.set_xlabel('MJD (days)', fontsize=14)

    ax.set_ylabel('Difference (deg)', fontsize=14)

    # Add legend

    ax.legend()

    # Show figure

    plt.show()

    # Return

    return

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Get difference

    list_mjd, list_dra, list_ddec = get_difference()

    # Plot difference

    plot_difference(list_mjd, list_dra, list_ddec)