# Source code for cosmo_utils.utils.geometry

#! /usr/bin/env python
# -*- coding: utf-8 -*-

# Victor Calderon
# Created      : 2018-05-03
from __future__ import absolute_import, division, print_function
__author__     = ['Victor Calderon']
__email__      = ['victor.calderon@vanderbilt.edu']
__maintainer__ = ['Victor Calderon']
__all__        = [  "flip_angles",
"Ang_Distance",
"Coord_Transformation"]
"""
Set of geometrical definitions for translations, coordinate tranformations,
etc.
"""

## Import modules
import numpy as np
import pandas as pd
from   cosmo_utils.utils import file_utils as fd
from   cosmo_utils.custom_exceptions import LSSUtils_Error

## Functions

# Restricting angles to be between 0 and 360
[docs]def flip_angles(ang, unit='deg'):
"""
Ensures that an angle is always between 0 and 360 degrees.

Parameters
----------
ang : float, int, numpy.ndarray
Angle in units of degrees.

unit : {'deg', 'rad'} str, optional
Unit of the angle. This variable is set to 'deg' by default.

Returns
----------
ang_out : float
Convertions of ang, ranging from 0 to 360 degrees.

Raises
----------
LSSUtils_Error : Exception
This error gets raised when ang is not a digit or a number.

Examples
----------
>>> flip_angles(-50, unit='deg')
310.0

>>> flip_angles(110, unit='deg')
110.0
"""
file_msg = fd.Program_Msg(__file__)
# Checking type of ang
valid_types = (float, int, list, np.ndarray)
if not (isinstance(ang, valid_types)):
msg = '{0} ang ({1}) is not a number! Exiting!'.format(
file_msg, ang)
raise LSSUtils_Error(msg)
# Checking unit
if not (unit.lower() in ['deg', 'rad']):
msg = '{0} unit ({1}) is not a valid option! Exiting!'.format(
file_msg, unit)
raise LSSUtils_Error(msg)
##
## Converting angle
if isinstance(ang, float) or isinstance(ang, int):
ang = float(ang)
# Converting from radians to degrees, if applicable
ang_converted = np.degrees(ang)
elif unit == 'deg':
ang_converted = ang
# Checking the range of ang
if ang_converted < 0:
ang_converted += 360
# Converting back to radians, if applicable
elif unit == 'deg':
ang_final = float(ang_converted)
else:
try:
ang = np.asarray(ang)
# Converting to degrees, if applicable
ang_converted = np.degrees(ang)
elif unit == 'deg':
ang_converted = ang
# Checking the range of ang
ang_converted = np.asarray([xx if xx > 0 else xx+360 for xx in
ang_converted])
# Converting back to radians, if applicable
elif unit == 'deg':
ang_final = ang_converted
# Converting to float
ang_final = ang_final.astype(float)
except:
msg = '{0} ang could not be converted!'.format(file_msg)
raise LSSUtils_Error(msg)

return ang_final

## Calculates the Angular distance between 2 points
[docs]def Ang_Distance(ra1, ra2, dec1, dec2, unit='deg', method='haversine'):
"""
Calculates angular separation between two sets of points with given
right ascensions and declinations.

Taken from: https://en.wikipedia.org/wiki/Haversine_formula

Parameters
-----------
ra1, ra2 : float
Right Ascension of the 1st and 2nd points. Units in degrees
by default.

dec1, dec2 : float
Declination of the 1st and 2nd points. Units in degrees by default.

Unit of ra1, ra2, dec1, and dec2.
This will also determine the final unit that outputs this function.

method : {'haversine', 'astropy'} str, optional
Method to use in order to calculate angular separation.
This variable is to by default to the haversine method.
If astropy, it will use the astropy framework to determine the
angular separation.

Returns
-----------
ang_sep : float
Angular separation between 1st and 2nd point.
In units of degrees.

Notes
-----------
A = 90. - dec2
B = 90. - dec1
D = ra1 - ra2
c = Angle between two points
"""
file_msg = fd.Program_Msg(__file__)
## Checking input arguments
# Valid options
methods_valid = ['haversine', 'astropy']
# Units
if not (unit in units_valid):
msg = '{0} unit ({1}) is not a valid argument'.format(
file_msg, unit)
raise LSSUtils_Error(msg)
# Method
if not (method in methods_valid):
msg = '{0} method ({1}) is not a valid argument'.format(
file_msg, method)
raise LSSUtils_Error(msg)
##
## Flipping angles
ra1 = flip_angles(ra1, unit=unit, )
ra2 = flip_angles(ra2, unit=unit)
##
## Haversine Method
if method == 'haversine':
D = np.radians(ra1 - ra2 )
# Distance
ang_sep  = (np.sin((A - B) * .5))**2.
ang_sep += np.sin(A) * np.sin(B) * (np.sin(D * .5))**2.
ang_sep  = np.degrees(2 * np.arcsin(ang_sep**0.5))
##
## Astropy Method
if method == 'astropy':
# Imports
from astropy.coordinates import SkyCoord
from astropy import units as u
# Converting to units
if unit == 'deg':
unit_opt = u.degree
# Calculations
p1 = SkyCoord(ra=ra1, dec=dec1, unit=(unit_opt, unit_opt))
p2 = SkyCoord(ra=ra2, dec=dec2, unit=(unit_opt, unit_opt))
# Angular Separation
ang_sep = p1.separation(p2)
# Converting to final units
if unit == 'deg':
ang_sep = ang_sep.degree

return ang_sep

## Coordinate Transformation function
[docs]def Coord_Transformation(ra, dec, dist, ra_cen, dec_cen, dist_cen,
trans_opt=4, return_dict=False, unit='deg'):
"""
Transforms spherical coordinates (ra, dec, dist) into cartesian
coordinates.

Parameters
-----------
ra, dec, dist : array_like, shape (N,)
Arrays of Right Ascension, declination, and distance.
Units are ['degrees', 'degrees', 'distance_units']

ra_cen, dec_cen, dist_cen : float, int
Right Ascension, declination, and distance for the center of
the coordinates. These correspond to where the corodinates
ra, dec, and dist will be centered.

trans_opt : {1, 2, 3, 4} int, optional
Option for cartesian translation/transformation for elements.
This variable ist set to 4 by default.

Options:
- 1 : No translation involved
- 2 : Translation to the center point.
- 3 : Translation and rotation to the center point.
- 4 : Translation and 2 rotaitons about the center point

return_dict : {True, False}, bool, optional
If True, this functions returns 2 dictionaries with spherical
and cartesian coordinates.
If False, it returns a pandas.DataFrame with the columns.
This variable is set to False by default.

Unit of ra1, ra2, dec1, and dec2.
This will also determine the final unit that outputs this function.
This variable is set to deg by default.

Returns
-----------
coord_dict (coord_pd) : python dictionary
Dictionary with spherical and cartesian dictionary of elements
based on trans_opt value. This value is returned if
return_dict is set to True. If not, a pandas.DataFrame is
return.
"""
file_msg = fd.Program_Msg(__file__)
## Check types of elements
# Units
if not (unit in unit_arr):
'{0} unit ({1}) is not a valid input!'.format(
file_msg, unit)
# Valid types
valid_types = (float, int, np.ndarray, list)
# Right Ascension
if not (isinstance(ra, valid_types)):
msg = '{0} ra ({1}) is not a valid type!'.format(
file_msg, type(ra))
raise LSSUtils_Error(msg)
# Declination
if not (isinstance(dec, valid_types)):
msg = '{0} dec ({1}) is not a valid type!'.format(
file_msg, type(dec))
raise LSSUtils_Error(msg)
# Distance
if not (isinstance(dist, valid_types)):
msg = '{0} dist ({1}) is not a valid type!'.format(
file_msg, type(dist))
raise LSSUtils_Error(msg)
# trans_opt
if not (trans_opt in list(range(1, 5))):
msg = '{0} trans_opt ({1}) is not within valid range!'.format(
file_msg, trans_opt)
raise LSSUtils_Error(msg)
##
## Centre's RA, DEC, DIST
# Right Ascension
if (isinstance(ra_cen, float)):
ra_cen = flip_angles(ra_cen)
else:
msg = '{0} ra_cen ({1}) is not a float!'.format(
file_msg, type(ra_cen))
raise LSSUtils_Error(msg)
# Declination
if (isinstance(dec_cen, float)):
dec_cen = flip_angles(dec_cen)
else:
msg = '{0} dec_cen ({1}) is not a float!'.format(
file_msg, type(dec_cen))
raise LSSUtils_Error(msg)
# Distance
if (isinstance(dist_cen, float)):
dist_cen = float(dist_cen)
else:
msg = '{0} dist_cen ({1}) is not a float!'.format(
file_msg, type(dist_cen))
raise LSSUtils_Error(msg)
##
## Check type of elements
# Right Ascension
if (isinstance(ra, float) or isinstance(ra, int)):
ra = np.array([ra])
else:
ra = np.array(ra)
# Declination
if (isinstance(dec, float) or isinstance(dec, int)):
dec = np.array([dec])
else:
dec = np.array(dec)
# Distance
if (isinstance(dist, float) or isinstance(dist, int)):
dist = np.array([dist])
else:
dist = np.array(dist)
##
## Converting to desired units
# Right Ascension
ra_cen_deg = np.degrees(ra_cen)
# Declination
dec_cen_deg = np.degrees(dec_cen)
elif unit == 'deg':
ra_cen_deg = ra_cen
# Declination
dec_cen_deg = dec_cen
##
## Initializing pandas DataFrame
dict_keys  = ['ra', 'dec', 'dist']
coord_dict = dict(zip(dict_keys, np.vstack([ra, dec, dist])))
##
## Spherical to Cartesian transformation
## 1st tranformation
# Centre
# All galaxies
##
## Rotations about z- and x-axis by tau and Omega
## Rotating the points, not the axes
x1 = x - x_cen
y1 = y - y_cen
z1 = z - z_cen
# Angles
omega = np.radians(90. - ra_cen_deg )
# Rotations about z-axis by omega
x2 = x1 * np.cos(omega) - y1 * np.sin(omega)
y2 = x1 * np.sin(omega) + y1 * np.cos(omega)
z2 = z1.copy()
# Rotations about X-axis by tau
x3 = x2.copy()
y3 = y2 * np.cos(tau) - z2 * np.sin(tau)
z3 = z2 * np.sin(tau) + z2 * np.cos(tau)
##
## Definining which variables to return
# No Translation
if trans_opt == 1:
coord_dict['x'] = x
coord_dict['y'] = y
coord_dict['z'] = z
# Translation
if trans_opt == 2:
coord_dict['x'] = x1
coord_dict['y'] = y1
coord_dict['z'] = z1
# Translation + Rotation
if trans_opt == 3:
coord_dict['x'] = x2
coord_dict['y'] = y2
coord_dict['z'] = z2
# Translation + 2 Rotation (centered about the centre)
if trans_opt == 4:
coord_dict['x'] = x3
coord_dict['y'] = y3
coord_dict['z'] = z3
##
## Checking what object to return, i.e. DataFrame or python dictionary
if return_dict:
return coord_dict
else:
coord_pd = pd.DataFrame(coord_dict)
return coord_pd

## Coordinate Transformation - Cartesian to Spherical coordinates
def cart_to_sph_coords(dist, x_arr, y_arr, z_arr, return_dict=False,
unit='deg'):
"""
Transforms cartesian coordinates (x, y, z) into spherical
coordinates.

Parameters
-----------
dist : np.ndarray
Array of the the distance to the N-number of objects.
Shape (N,).

x_arr, y_arr, z_arr : np.ndarray
Arrays of the cartesian coordinates (x, y, z) of the
N objects. Shape (N,). Units are in unit.

return_dict : {True, False}, bool, optional
If True, this functions returns 2 dictionaries with spherical
and cartesian coordinates.
If False, it returns a pandas.DataFrame with the columns.
This variable is set to False by default.

unit : {'dec','rad'} str, optional
Unit of the output ra, dec coordinates.
This variable is set to deg by default.

Returns
-----------
coord_dict (coord_pd) : dict
Dictionary with spherical coordinates (ra, dec) and
distance (dist) of N elements. This value is returned if
return_dict is set to True. If not, a pd.DataFrame is
returned.
"""
file_msg = fd.Program_Msg(__file__)
## Check types of elements
# Units
if not (unit in unit_arr):
'{0} unit ({1}) is not a valid input!'.format(
file_msg, unit)
# Valid types
valid_types = (float, int, np.ndarray, list)
# X-coordinate
if not (isinstance(x_arr, valid_types)):
msg = '{0} x_arr ({1}) is not a valid type!'.format(
file_msg, type(x_arr))
raise LSSUtils_Error(msg)
# Y-coordinate
if not (isinstance(y_arr, valid_types)):
msg = '{0} y_arr ({1}) is not a valid type!'.format(
file_msg, type(y_arr))
raise LSSUtils_Error(msg)
# Z-coordinate
if not (isinstance(z_arr, valid_types)):
msg = '{0} z_arr ({1}) is not a valid type!'.format(
file_msg, type(z_arr))
raise LSSUtils_Error(msg)
# Distance
if not (isinstance(dist, valid_types)):
msg = '{0} dist ({1}) is not a valid type!'.format(
file_msg, type(dist))
raise LSSUtils_Error(msg)
##
## Check type of elements
# X-Coordinate
if (isinstance(x_arr, float) or isinstance(x_arr, int)):
x_arr = np.array([x_arr])
else:
x_arr = np.array(x_arr)
# Y-Coordinate
if (isinstance(y_arr, float) or isinstance(y_arr, int)):
y_arr = np.array([y_arr])
else:
y_arr = np.array(y_arr)
# Z-Coordinate
if (isinstance(z_arr, float) or isinstance(z_arr, int)):
z_arr = np.array([z_arr])
else:
z_arr = np.array(z_arr)
# Distance
if (isinstance(dist, float) or isinstance(dist, int)):
dist = np.array([dist])
else:
dist = np.array(dist)
##
## Initializing pandas DataFrame
dict_keys  = ['x', 'y', 'z', 'dist']
coord_dict = dict(zip(dict_keys, np.vstack([x_arr, y_arr, z_arr, dist])))
##
## Cartesian to Spherical coordinates
##
## - Normalized cartesian coordinates
cart_arr = np.column_stack([x_arr, y_arr, z_arr])

(   x_val,
y_val,
z_val) = cart_arr/float(dist)
# Distance to object
dist = float(dist)
## Declination
dec_val = 90. - num.degrees(num.arccos(z_val))
## Right ascension
if x_val == 0:
if y_val > 0.:
ra_val = 90.
elif y_val < 0.:
ra_val = -90.
else:
ra_val = num.degrees(num.arctan(y_val/x_val))
##
## Seeing on which quadrant the point is at
if x_val < 0.:
ra_val += 180.
elif (x_val >= 0.) and (y_val < 0.):
ra_val += 360.

return ra_val, dec_val