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

"""
Created by schaffer at 8/22/19

"""

import math

import numpy as np
from astropy import units as u
from astropy.wcs import WCS
from sunpy.coordinates import Helioprojective, frames


def calc_theta(x, y, r_ref=945.57):
    """Calculate the value of the heliocentric angle theta

    Args:
      x: helio-projective cartesian x posn in arcsec
      y: helio-projective cartesian y posn in arcsec
      r_ref: apparent angular disk radiusin arcsun

    Returns:
      theta: float

    """
    r_obs = math.sqrt(math.pow(x, 2) + math.pow(y, 2))
    radius_ratio = min([r_obs / r_ref, 1])
    theta_value = math.asin(radius_ratio) / math.pi * 180
    return theta_value


def calc_mu(x, y, r_ref=945.57):
    """Calculate the value of mu (cos of heliocentric angle)

    Args:
      x: helio-projective cartesian x posn in arcsec
      y: helio-projective cartesian y posn in arcsec
      r_ref: apparent angular disk radiusin arcsun

    Returns:
      mu: float, range 1 to 0; 1 means disk center

    """
    mu = math.cos(calc_theta(x, y, r_ref))
    return mu


def to_heliographic(coords, date):
    """
    Transform a coordinate tuple from Helioprojective Cartesian to Heliographic coordinates

    """
    x, y = coords
    input = Helioprojective(x * u.arcsec, y * u.arcsec, obstime=date, observer='Earth')
    output = input.transform_to(frames.HeliographicStonyhurst)
    lon, lat = output.lon.degree, output.lat.degree
    return lon, lat


def calc_wcs_footprint(header):
    """Calculate the WCS footprint of the image

    Args:
      header: FITS header

    Returns:
      NumPy array of dimensionless min and max values of x and y
      with shape [[xmin, ymin], [xmax, ymax]]

    """
    wcs = WCS(header)
    pixel_shape = wcs.pixel_shape
    n_x = pixel_shape[0]
    n_y = pixel_shape[1]
    naxis = len(pixel_shape)

    slit_rect = np.array([[0, n_x - 1] + [0] * (naxis - 2), [0, n_y - 1] + [0] * (naxis - 2)])
    slit_world = np.array(wcs.pixel_to_world_values(slit_rect))
    slit_deg = u.degree * slit_world
    slit_arcsec = slit_deg.to(u.arcsec)

    return np.array([[min(slit_arcsec[0, 0], slit_arcsec[1, 0]).value,
                      min(slit_arcsec[0, 1], slit_arcsec[1, 1]).value],
                     [max(slit_arcsec[0, 0], slit_arcsec[1, 0]).value,
                      max(slit_arcsec[0, 1], slit_arcsec[1, 1]).value]])