import fitsio
import numpy as np
import time
import math
import sys

def euler(ra,dec,select=1,fk4=False,radian=False):
  """
  Copying euler.pro from IDL, found here: http://idlastro.gsfc.nasa.gov/ftp/pro/astro/euler.pro
  Only coded select = 1 (transform to galactic coordinates)
  From IL:
  /FK4 - If this keyword is set and non-zero, then input and output 
  celestial and ecliptic coordinates should be given in equinox 
  B1950.
       /RAIAN - if set, then all input and output angles are in radians rather
  than degrees.
       SELECT  - The coordinate conversion integer (1-6) may alternatively be 
   specified as a keyword (CURRENTLY ONLY 1 SUPPORTE).
  example: euler(299.590315,35.201604)
  Has not really been tested!!!  Works for splitting N/S, no other tests performed.
  Should just try it on a couple of values.
  Ok, tested on a single ra/dec combo and it works after arctan --> arctan2
  """
  assert select == 1
  assert fk4 == False

  rad_to_deg = 180.0/np.pi
  fourpi = 4.0*np.pi
  twopi = 2.0*np.pi

  if fk4 == True:
    equinox = '(B1950)'
    ## these are all radians, I think. 
    psi   = np.array([ 0.57595865315, 4.9261918136,  0.00000000000, 0.0000000000,  0.11129056012, 4.7005372834])
    stheta =np.array([ 0.88781538514,-0.88781538514,  0.39788119938,-0.39788119938,  0.86766174755,-0.86766174755])
    ctheta =np.array([ 0.46019978478, 0.46019978478,  0.91743694670, 0.91743694670,  0.49715499774, 0.49715499774])
    phi  = np.array([ 4.9261918136,  0.57595865315, 0.0000000000, 0.00000000000, 4.7005372834, 0.11129056012])
  else:
    equinox = '(J2000)'
    psi   = np.array([ 0.57477043300, 4.9368292465,  0.00000000000, 0.0000000000, 0.11142137093, 4.71279419371])
    stheta = np.array([ 0.88998808748,-0.88998808748, 0.39777715593,-0.39777715593,  0.86766622025,-0.86766622025])
    ctheta = np.array([ 0.45598377618, 0.45598377618,  0.91748206207, 0.91748206207,  0.49714719172, 0.49714719172])
    phi  = np.array([4.9368292465,  0.57477043300,  0.0000000000, 0.00000000000, 4.71279419371, 0.11142137093])

  indx = select - 1

  ## does this make sense?  i think ra/dec will be in degrees.
  
  if radian == False:
    ao = ra/rad_to_deg - phi[indx]
    bo = dec/rad_to_deg
  else:
    ao = ra - phi[indx]
    bo = dec - 0.0 ## make a copy of dec.
  
  sb = np.sin(bo)
  cb = np.cos(bo)
  cbsa = cb * np.sin(ao)
  bo  = -stheta[indx] * cbsa + ctheta[indx] * sb
  
  xx = np.where(bo > 1.0)[0]
  bo[xx] = 1.0
  xx = np.where(bo < -1.0)[0]
  bo[xx] = -1.0
  bo = np.arcsin(bo)
  if radian == False: 
    bo = bo * rad_to_deg
  ## silently gives in the wrong answer with np.arctan(x,y) !!! doh.
  ao = np.arctan2(ctheta[indx] * cbsa + stheta[indx] * sb, cb * np.cos(ao))
  ao = ((ao + psi[indx] + fourpi) % twopi)
  if radian == False:
    ao = ao * rad_to_deg
  return ao, bo

## for testing purposes.
def angsep(ra1, dec1, ra2, dec2):
  """
  Compute the angular separation of two ra/dec pairs, assuming they are input in radians.
  Returns radians.
  """
  myarg = np.cos(dec1)*np.cos(dec2)*np.cos(ra1-ra2)+np.sin(dec1)*np.sin(dec2)
  if(myarg > 1):
    #print 'ahh!', ra1, dec1, ra2, dec2
    return 0.0
  return math.acos(myarg)

def angsepdeg(ra1, dec1, ra2, dec2):
  """
  Compute the angular separation of two ra/dec pairs, assuming they are input in degrees.
  Returns degrees.
  """
  degtorad = np.pi/180.
  myarg = np.cos(dec1*degtorad)*np.cos(dec2*degtorad)*np.cos((ra1-ra2)*degtorad)+np.sin(dec1*degtorad)*np.sin(dec2*degtorad)
  if(myarg > 1):
    #print 'ahh!', ra1, dec1, ra2, dec2
    return 0.0
  return math.acos(myarg)/degtorad