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

import itertools

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

# Compute position of the Sun #

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

def sun(jd):

    """

    From https://en.wikipedia.org/wiki/Position_of_the_Sun

    """

    # Number of days since Greenwich noon, Terrestrial Time, on 1 January 2000

    n = jd - 2451545.0

    # Mean longitude of the Sun, corrected for the aberration of light

    L = 280.460 + 0.9856474 * n

    while L < 0.0:

        L += 360.0

    while L >= 360.0:

        L -= 360.0

    # Mean anomaly of the Sun

    g = 357.528 + 0.9856003 * n

    while g < 0.0:

        g += 360.0

    while g >= 360.0:

        g -= 360.0

    # Ecliptic longitude of the Sun

    lam = L + 1.915 * math.sin(g*gammalib.deg2rad) + \

              0.020 * math.sin(2.0*g*gammalib.deg2rad)

    # Compute Right Ascension and Declination

    ra  = math.atan(math.cos(23.43711*gammalib.deg2rad) *

                    math.tan(lam*gammalib.deg2rad))*gammalib.rad2deg

    dec = math.asin(math.sin(23.43711*gammalib.deg2rad) *

                    math.sin(lam*gammalib.deg2rad))*gammalib.rad2deg

    # Put ra in the right quadrant (same quadrant as lam)

    while lam < 0.0:

        lam += 360.0

    while lam >= 360.0:

        lam -= 360.0

    while ra < 0.0:

        ra += 360.0

    while ra >= 360.0:

        ra -= 360.0

    q_lam = int(lam / 90.0)

    q_ra  = int(ra / 90.0)

    ra   += (q_lam-q_ra)*90.0

    print(q_lam,q_ra,q_lam-q_ra,(q_lam-q_ra)*90.0)

    # Return

    return (ra, dec)

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

# Compute Sun constraint map #

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

def sun_constraint(lim=100.0):

    """

    Compute Sun constraint map

    """

    # Initialise sky map

    map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180)

    # Set Julian days for one "orbit" of the Sun

    nsteps = 1000

    jds    = [float(i)*365.25636/float(nsteps) for i in range(nsteps)]

    norm   = 1.0/float(nsteps)

    # Loop over JDs

    for jd in jds:

        # Get position of sun

        ra, dec = sun(jd)

        print(jd,ra,dec)

        sundir  = gammalib.GSkyDir()

        sundir.radec_deg(ra, dec)

        # Compute for all pixels distance to Sun

        for inx in range(map.npix()):

            dir = map.inx2dir(inx)

            if sundir.dist_deg(dir) > lim:

                map[inx] += norm

    # Save map

    map.save('sun.fits', True)

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

# Plot zenith angle #

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

def plot_zenith_angle(lat=-24.4, lim=60.0):

    """

    Plot zenith angle.

    The zenith angle of a position (ra,dec) depends on the declination and

    the hour angle h and is given by

    zenith(ra,dec) = arccos( sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(h) )

    The hour angle (or local hour angle, LHA) is defined as the difference

    between local siderial time (LST) and the Right Ascension

    h = LST - ra

    LST is the difference between the Greenwhich siderial time (GST) and the

    geographic longitude of the observer

    LST = GST - lon

    where lon is counted positive west from the meridian. Hence

    h = GST - lon - ra

    and

    zenith(ra,dec) = arccos( sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(GST-lon-ra) )

    The zenith angle of the Sun depends of the hour angle in the local solar

    time h_sol, the current declination of the Sun dec_sol and the geographic

    latitude

    zenith_sol = arccos( sin(lat) * sin(dec_sol) + cos(lat) * cos(dec_sol) * cos(h_sol) )

    """

    # Initialise sky map

    map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180,60)

    # Set hour angle array in radians

    hours = [float(i)*0.1 for i in range(3600)]

    norm  = 24.0 / float(len(hours)) * 365.0    # hours per year per time step

    # Create figure

    plt.figure(1)

    plt.title('Zenith angles')

    # Loop over declinations

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

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

    for k, dec in enumerate(decs):

        # Compute zenith angles for a given declination as function of hour

        # angle

        cos_dec = math.cos(dec*gammalib.deg2rad)

        sin_dec = math.sin(dec*gammalib.deg2rad)

        cos_phi = math.cos(lat*gammalib.deg2rad)

        sin_phi = math.sin(lat*gammalib.deg2rad)

        zeniths = [math.acos(sin_phi*sin_dec +

                             cos_phi*cos_dec*math.cos(h*gammalib.deg2rad))*

                   gammalib.rad2deg for h in hours]

        # Select only zenith angles < lim

        x = []

        y = []

        for i in range(len(hours)):

            if zeniths[i] <= lim:

                x.append(hours[i])

                y.append(zeniths[i])

        # Update visibility

        visible[k] = float(len(x)) * norm

        # Plot zenith angle

        if len(x) > 0:

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

        # Fill visibility cube for all hour angles. For each hour angle compute

        # the appropriate zenith angle index iz, and set all Right Ascensions.

        # This implies that we loop uniformly over all local siderial times.

        #

        # In reality, however, the Sun constraint (Sun under horizon) translates

        # into a set of time intervals during which observations are performed.

        for i in range(len(hours)):

            iz = int(zeniths[i]/1.0)

            if iz < 60:

                for ira in range(360):

                    index = ira + 360 * k

                    map[index,iz] += norm

    # Set axes

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

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

    # Create figure

    plt.figure(2)

    plt.title('Visibility')

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

    plt.xlabel('Declination (degrees)')

    plt.ylabel('Visibility (hours)')

    # Save sky map

    map.save('viscube_north.fits', True)

    # Show plot

    plt.show()

    # Return

    return

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

# Compute zenith angle map #

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

def zenith_angle_map(lat=-24.4):

    """

    Compute zenith angle map

    The zenith angle of a position (ra,dec) depends on the declination and

    the hour angle h and is given by

    zenith(h,dec) = arccos( sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(h) )

    The hour angle h (or local hour angle, LHA) is defined as the difference

    between local siderial time (LST) and the Right Ascension

    h = LST - ra

    LST is the difference between the Greenwich siderial time (GST) and the

    geographic longitude of the observer

    LST = GST - lon

    where lon is counted positive west from the meridian. Hence

    h = GST - lon - ra

    The map is computed for h=-ra which is equivalent to GST=lon (or LST=0).

    In other words, the map corresponds to the time when ra=0 passes the

    local meridian.

    """

    # Initialise zenith angle map

    map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180)

    # Set hour angle and declination vectors

    hours = [float(i) for i in range(360)]

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

    # Precompute latitude terms

    cos_lat = math.cos(lat*gammalib.deg2rad)

    sin_lat = math.sin(lat*gammalib.deg2rad)

    # Loop over all declination and hour angles and compute the zenith

    # angle. Store the zenith angle in the map

    index = 0

    for dec in decs:

        cos_dec = math.cos(dec*gammalib.deg2rad)

        sin_dec = math.sin(dec*gammalib.deg2rad)

        for h in hours:

            cos_h = math.cos(h*gammalib.deg2rad)

            zenith = math.acos(sin_lat*sin_dec +

                               cos_lat*cos_dec*cos_h)*gammalib.rad2deg

            map[index] = zenith

            index     += 1

    # Save zenith angle map

    map.save('zenith_angles.fits', True)

    # Return

    return

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

# Compute visibility cube #

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

def visibility_cube():

    """

    Compute visibility cube by displacing the zenith angle map for all

    hour angles. The visibility cube contains the number of hours a given

    celestial position is visible under a given zenith angle interval.

    Summing over all zenith angle intervals specifies for how long a given

    celestial position will be visible.

    """

    # Set zenith angle range and bin size

    dz = 1.0

    nz = int(60.0/dz)

    # Initialise visibility cube

    cube = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180,nz)

    # Load zenith angle map

    map = gammalib.GSkyMap('zenith_angles.fits')

    # Setup hour angle weights. They specify for how many hours a given

    # hour angle will be observed during one year.

    hours = hour_angle_weight()

    # Compute normalisation factor

    norm = 365.0/float(len(hours))

    # Loop over all hour angle weights

    for h, weight in enumerate(hours):

        # Compute temporal shift

        shift = 360.0/float(len(hours)) * float(h)

        # Initialise shifted zenith angle map

        map_shift = gammalib.GSkyMap('CAR','CEL',shift,0.0,-1.0,1.0,360,180)

        # Merge zenith angle map into shifted map

        map_shift += map

        map_shift.save('zenith_angles_shifted.fits', True)

        # Loop over all pixels of shifted map

        for i, zenith in enumerate(map_shift):

            iz     = int(zenith/dz)

            if iz < nz:

                cube[i,iz] += weight * norm

    # Save visibility cube

    cube.save('viscube.fits', True)

    # Return

    return

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

# Compute visibility cube #

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

def hour_angle_weight():

    """

    Compute during which time the array was observing due to hour angle

    constraints. We take here a very simple model where the night lasts

    between 6 and 10 hours (or 90-150 degrees), with a sinusoidal variation

    along the year. Every day the hour angle of the Sun set is advancing by 

    4 minutes.

    """

    # Set number of hour angle bins

    nsteps = 1000

    d_min  =  6.0/24.0*nsteps

    d_max  = 10.0/24.0*nsteps

    # Initialise hour angle array and weight

    hours  = [0.0 for i in range(nsteps)]

    weight = 24.0/float(nsteps)  # Weight per hour angle step

    # Loop over all days in one year

    for day in range(nsteps):

        h_set  = day

        h_rise = day + int(d_min + (d_max-d_min)*math.sin(day/float(nsteps)*gammalib.pi))

        for i in range(h_set,h_rise):

            if i >= nsteps:

                hours[i-nsteps] += weight

            else:

                hours[i] += weight

    # Return

    return hours

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Plot zenith angle

    #zenith_angle_map(lat=-24.4)        # South

    zenith_angle_map(lat=+28.7569)     # North

    #visibility_cube()

    visibility_cube()