#########################################################################
#
# Script: coords.py
#
# Author: Igor Volobouev, March 2005
#
# Purpose: contains code relevant to coordinate systems used
#          in cosmic ray generation
#
#########################################################################

from basic_kinematics import *

# Tags for basic coordinate systems
LOCAL = 0
CELESTIAL = 1
GALACTIC = 2
ZENITH = 3
GEOGRAPHIC = 4

# Rotation to transform celectial coordinates into galactic.
# Transformation constants are taken from
# http://astro.estec.esa.nl/Hipparcos/CATALOGUE_VOL1/sect1_05.pdf
__r1 = Rotation(V3(0, 0, 1), 32.93192/180.0*math.pi)
__r2 = Rotation(__r1(V3(1, 0, 0)), (27.12825 - 90.0)/180.0*math.pi)
__r3 = Rotation((__r2*__r1)(V3(0, 0, 1)), -(192.85948 + 90.0)/180.0*math.pi)
__celectial_to_galactic = __r3 * __r2 * __r1

# Rotation to transform galactic coordinates into celestial
__galactic_to_celestial = ~__celectial_to_galactic

def rotation(fromsys, tosys, orbit=None, t=None, recursion_level=0):
    "Rotation needed to convert vectors between two coordinate systems."
    if (fromsys == tosys):
        return Rotation()
    elif (fromsys == LOCAL and tosys == CELESTIAL):
        return orbit.orientation(t)
    elif (fromsys == GALACTIC and tosys == CELESTIAL):
        return __galactic_to_celestial
    elif (fromsys == CELESTIAL and tosys == GALACTIC):
        return __celectial_to_galactic
    elif (tosys == GALACTIC):
        return __celectial_to_galactic * rotation(fromsys, CELESTIAL,
                                                  orbit, t, recursion_level)
    elif (fromsys == GALACTIC):
        return rotation(CELESTIAL, tosys, orbit, t, recursion_level) * \
               __galactic_to_celestial
    elif (fromsys == GEOGRAPHIC and tosys == CELESTIAL):
        zenith, alt = orbit.zenith_and_altitude(t)
        # For now, we are going to neglect all polar motion
        # (precession, nutation, etc.). This is good enough
        # for the purpose of simulating the rigidity cutoffs
        # which are done to a few degrees anyway. However,
        # the precision of this conversion is probably too
        # low for any other purpose since the earth's pole
        # direction changes by about 50 arcseconds per year.
        pole_direction = V3(0.0, 0.0, 1.0)
        eastward = vprod3(pole_direction, zenith)
        return Rotation(eastward.ra(), eastward.dec(),
                        zenith.ra(), zenith.dec())
    elif (recursion_level == 0):
        return ~rotation(tosys, fromsys, orbit, t, 1)
    else:
        raise NotImplementedError, "coords.rotation"