"""module for coordinate manipulation (conversions, calculations etc.)

(c) 2007-2012 Matt Hilton

(c) 2013 Matt Hilton & Steven Boada

U{http://astlib.sourceforge.net}

"""

import sys
import math
import numpy
import string
from PyWCSTools import wcscon

#-----------------------------------------------------------------------------
def hms2decimal(RAString, delimiter):
    """Converts a delimited string of Hours:Minutes:Seconds format into decimal
    degrees.

    @type RAString: string
    @param RAString: coordinate string in H:M:S format
    @type delimiter: string
    @param delimiter: delimiter character in RAString
    @rtype: float
    @return: coordinate in decimal degrees

    """
    # is it in HH:MM:SS format?
    if delimiter == "":
        RABits = str(RAString).split()
    else:
        RABits = str(RAString).split(delimiter)
    if len(RABits) > 1:
        RAHDecimal = float(RABits[0])
        if len(RABits) > 1:
            RAHDecimal = RAHDecimal+(float(RABits[1])/60.0)
        if len(RABits) > 2:
            RAHDecimal = RAHDecimal+(float(RABits[2])/3600.0)
        RADeg = (RAHDecimal/24.0)*360.0
    else:
        RADeg = float(RAString)

    return RADeg

#-----------------------------------------------------------------------------
def dms2decimal(decString, delimiter):
    """Converts a delimited string of Degrees:Minutes:Seconds format into
    decimal degrees.

    @type decString: string
    @param decString: coordinate string in D:M:S format
    @type delimiter: string
    @param delimiter: delimiter character in decString
    @rtype: float
    @return: coordinate in decimal degrees

    """
    # is it in DD:MM:SS format?
    if delimiter == "":
        decBits = str(decString).split()
    else:
        decBits = str(decString).split(delimiter)
    if len(decBits) > 1:
        decDeg = float(decBits[0])
        if decBits[0].find("-") != -1:
            if len(decBits) > 1:
                decDeg = decDeg-(float(decBits[1])/60.0)
            if len(decBits) > 2:
                decDeg = decDeg-(float(decBits[2])/3600.0)
        else:
            if len(decBits) > 1:
                decDeg = decDeg+(float(decBits[1])/60.0)
            if len(decBits) > 2:
                decDeg = decDeg+(float(decBits[2])/3600.0)
    else:
        decDeg = float(decString)

    return decDeg

#-----------------------------------------------------------------------------
def decimal2hms(RADeg, delimiter):
    """Converts decimal degrees to string in Hours:Minutes:Seconds format with
    user specified delimiter.

    @type RADeg: float
    @param RADeg: coordinate in decimal degrees
    @type delimiter: string
    @param delimiter: delimiter character in returned string
    @rtype: string
    @return: coordinate string in H:M:S format

    """
    hours = (RADeg/360.0)*24
    #if hours < 10 and hours >= 1:
    if 1 <= hours < 10:
        sHours = "0"+str(hours)[0]
    elif hours >= 10:
        sHours = str(hours)[:2]
    elif hours < 1:
        sHours = "00"

    if str(hours).find(".") == -1:
        mins = float(hours)*60.0
    else:
        mins = float(str(hours)[str(hours).index("."):])*60.0
    #if mins<10 and mins>=1:
    if 1 <= mins<10:
        sMins = "0"+str(mins)[:1]
    elif mins >= 10:
        sMins = str(mins)[:2]
    elif mins < 1:
        sMins = "00"

    secs = (hours-(float(sHours)+float(sMins)/60.0))*3600.0
    #if secs < 10 and secs>0.001:
    if 0.001 < secs < 10:
        sSecs = "0"+str(secs)[:str(secs).find(".")+4]
    elif secs < 0.0001:
        sSecs = "00.001"
    else:
        sSecs = str(secs)[:str(secs).find(".")+4]
    if len(sSecs) < 5:
        sSecs = sSecs+"00"	# So all to 3dp

    if float(sSecs) == 60.000:
        sSecs = "00.00"
        sMins = str(int(sMins)+1)
    if int(sMins) == 60:
        sMins = "00"
        sDeg = str(int(sDeg)+1)

    return sHours+delimiter+sMins+delimiter+sSecs

#------------------------------------------------------------------------------
def decimal2dms(decDeg, delimiter):
    """Converts decimal degrees to string in Degrees:Minutes:Seconds format
    with user specified delimiter.

    @type decDeg: float
    @param decDeg: coordinate in decimal degrees
    @type delimiter: string
    @param delimiter: delimiter character in returned string
    @rtype: string
    @return: coordinate string in D:M:S format

    """
    # Positive
    if decDeg > 0:
        #if decDeg < 10 and decDeg>=1:
        if 1 <= decDeg < 10:
            sDeg = "0"+str(decDeg)[0]
        elif decDeg >= 10:
            sDeg = str(decDeg)[:2]
        elif decDeg < 1:
            sDeg = "00"

        if str(decDeg).find(".") == -1:
            mins = float(decDeg)*60.0
        else:
            mins = float(str(decDeg)[str(decDeg).index("."):])*60
        #if mins<10 and mins>=1:
        if 1 <= mins < 10:
            sMins = "0"+str(mins)[:1]
        elif mins >= 10:
            sMins = str(mins)[:2]
        elif mins < 1:
            sMins = "00"

        secs = (decDeg-(float(sDeg)+float(sMins)/60.0))*3600.0
        #if secs<10 and secs>0:
        if 0 < secs < 10:
            sSecs = "0"+str(secs)[:str(secs).find(".")+3]
        elif secs < 0.001:
            sSecs = "00.00"
        else:
            sSecs = str(secs)[:str(secs).find(".")+3]
        if len(sSecs) < 5:
            sSecs = sSecs+"0"	# So all to 2dp

        if float(sSecs) == 60.00:
            sSecs = "00.00"
            sMins = str(int(sMins)+1)
        if int(sMins) == 60:
            sMins = "00"
            sDeg = str(int(sDeg)+1)

        return "+"+sDeg+delimiter+sMins+delimiter+sSecs

    else:
        #if decDeg>-10 and decDeg<=-1:
        if -10 < decDeg <= -1:
            sDeg = "-0"+str(decDeg)[1]
        elif decDeg <= -10:
            sDeg = str(decDeg)[:3]
        elif decDeg > -1:
            sDeg = "-00"

        if str(decDeg).find(".") == -1:
            mins = float(decDeg)*-60.0
        else:
            mins = float(str(decDeg)[str(decDeg).index("."):])*60
        #if mins<10 and mins>=1:
        if 1 <= mins < 10:
            sMins = "0"+str(mins)[:1]
        elif mins >= 10:
            sMins = str(mins)[:2]
        elif mins < 1:
            sMins = "00"

        secs = (decDeg-(float(sDeg)-float(sMins)/60.0))*3600.0
        #if secs>-10 and secs<0:
        # so don't get minus sign
        if -10 < secs < 0:
            sSecs = "0"+str(secs)[1:str(secs).find(".")+3]
        elif secs > -0.001:
            sSecs = "00.00"
        else:
            sSecs = str(secs)[1:str(secs).find(".")+3]
        if len(sSecs) < 5:
            sSecs = sSecs+"0"	# So all to 2dp

        if float(sSecs) == 60.00:
            sSecs = "00.00"
            sMins = str(int(sMins)+1)
        if int(sMins) == 60:
            sMins = "00"
            sDeg = str(int(sDeg)-1)

        return sDeg+delimiter+sMins+delimiter+sSecs

#-----------------------------------------------------------------------------
def calcAngSepDeg(RA1, dec1, RA2, dec2):
    """Calculates the angular separation of two positions on the sky (specified
    in decimal degrees) in decimal degrees, assuming a tangent plane projection
    (so separation has to be <90 deg). Note that RADeg2, decDeg2 can be numpy
    arrays.

    @type RADeg1: float
    @param RADeg1: R.A. in decimal degrees for position 1
    @type decDeg1: float
    @param decDeg1: dec. in decimal degrees for position 1
    @type RADeg2: float or numpy array
    @param RADeg2: R.A. in decimal degrees for position 2
    @type decDeg2: float or numpy array
    @param decDeg2: dec. in decimal degrees for position 2
    @rtype: float or numpy array, depending upon type of RADeg2, decDeg2
    @return: angular separation in decimal degrees

    ***** Updates added by Alessandro Nastasi - 26/06/2014 *****

    Now RA and DEC can be entered indifferently as decimal degrees or as 
    hh:mm:ss.s or hh mm ss.s (provided that they are passed as STRING
    in the latter two cases). An internal routine recognizes the correct format 
    and converts them into decimal degrees.
 
    In addition, the tangent plane approximation restriction (i.e., dist < 90 deg) 
    has been removed and the complete formula is now implemented. 
    Pythagoras approximation is applied only for dist < 1 arcsec.
    """
    # ** ORIGINAL CODE **
    #
    #cRA = numpy.radians(RADeg1)
    #cDec = numpy.radians(decDeg1)

    #gRA = numpy.radians(RADeg2)
    #gDec = numpy.radians(decDeg2)

    #dRA = cRA-gRA
    #dDec = gDec-cDec
    #cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) +
        #(numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA)))
    #x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC
    #y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) *
        #numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC)
    #r = numpy.degrees(numpy.sqrt(x*x+y*y))

    RA1 = str(RA1).strip()
    dec1 = str(dec1).strip()
    RA2 = str(RA2).strip()
    dec2 = str(dec2).strip()
    
    inp_param = [str(RA1), dec1, RA2, dec2]
    newinp = []
    for x in inp_param: 
      newinp.append(string.replace(x, ":", " "))
      
    if str(newinp[0]).find(' ') >= 0:
      RADeg1 = hms2decimal(newinp[0], ' ')
      decDeg1 = dms2decimal(newinp[1], ' ')
    else:
      RADeg1 = float(newinp[0])
      decDeg1 = float(newinp[1])

    if str(newinp[2]).find(' ') >= 0:
      RADeg2 = hms2decimal(newinp[2], ' ')
      decDeg2 = dms2decimal(newinp[3], ' ')
    else:
      RADeg2 = float(newinp[2])
      decDeg2 = float(newinp[3])
	
    cRA = numpy.radians(RADeg1)
    cDec = numpy.radians(decDeg1)

    gRA = numpy.radians(RADeg2)
    gDec = numpy.radians(decDeg2)
      
    x=numpy.cos(cRA)*numpy.cos(cDec)*numpy.cos(gRA)*numpy.cos(gDec)
    y=numpy.sin(cRA)*numpy.cos(cDec)*numpy.sin(gRA)*numpy.cos(gDec)
    z=numpy.sin(cDec)*numpy.sin(gDec)

    rad=numpy.degrees(numpy.arccos(round(x+y+z,10)))

    dRA = cRA-gRA
    dDec = gDec-cDec
    #cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) +
        #(numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA)))
    #x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC
    #y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) *
        #numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC)
    #r = numpy.degrees(numpy.sqrt(x*x+y*y))
    
    # use Pythagoras approximation if rad < 1 arcsec
    if rad<0.000004848:
      rad=numpy.degrees(numpy.sqrt((numpy.cos(cDec)*dRA)**2+dDec**2))

    return rad
    
#-----------------------------------------------------------------------------
def calcAngSepDegPythagoras(RADeg1, decDeg1, RADeg2, decDeg2):
    # ** ORIGINAL CODE **
    #
    cRA = numpy.radians(RADeg1)
    cDec = numpy.radians(decDeg1)

    gRA = numpy.radians(RADeg2)
    gDec = numpy.radians(decDeg2)

    dRA = cRA-gRA
    dDec = gDec-cDec
    cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) +
        (numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA)))
    x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC
    y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) *
        numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC)
    r = numpy.degrees(numpy.sqrt(x*x+y*y))
  
    return r

#-----------------------------------------------------------------------------
def shiftRADec(ra1, dec1, deltaRA, deltaDec):
    """Computes new right ascension and declination shifted from the original
    by some delta RA and delta DEC. Input position is decimal degrees. Shifts
    (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on
    an IDL routine of the same name.

    @param ra1: float
    @type ra1: R.A. in decimal degrees
    @param dec1: float
    @type dec1: dec. in decimal degrees
    @param deltaRA: float
    @type deltaRA: shift in R.A. in arcseconds
    @param deltaDec: float
    @type deltaDec: shift in dec. in arcseconds
    @rtype: float [newRA, newDec]
    @return: shifted R.A. and dec.

    """

    d2r = math.pi/180.
    as2r = math.pi/648000.

    # Convert everything to radians
    rara1 = ra1*d2r
    dcrad1 = dec1*d2r
    shiftRArad = deltaRA*as2r
    shiftDCrad = deltaDec*as2r

    # Shift!
    deldec2 = 0.0
    sindis = math.sin(shiftRArad / 2.0)
    sindelRA = sindis / math.cos(dcrad1)
    delra = 2.0*math.asin(sindelRA) / d2r

    # Make changes
    ra2 = ra1+delra
    dec2 = dec1 +deltaDec / 3600.0

    return ra2, dec2

#-----------------------------------------------------------------------------
def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch):
    """Converts specified coordinates (given in decimal degrees) between J2000,
    B1950, and Galactic.

    @type inputSystem: string
    @param inputSystem: system of the input coordinates (either "J2000",
        "B1950" or "GALACTIC")
    @type outputSystem: string
    @param outputSystem: system of the returned coordinates (either "J2000",
        "B1950" or "GALACTIC")
    @type coordX: float
    @param coordX: longitude coordinate in decimal degrees, e.g. R. A.
    @type coordY: float
    @param coordY: latitude coordinate in decimal degrees, e.g. dec.
    @type epoch: float
    @param epoch: epoch of the input coordinates
    @rtype: list
    @return: coordinates in decimal degrees in requested output system

    """

    if inputSystem=="J2000" or inputSystem=="B1950" or inputSystem=="GALACTIC":
        if outputSystem=="J2000" or outputSystem=="B1950" or \
            outputSystem=="GALACTIC":

            outCoords=wcscon.wcscon(wcscon.wcscsys(inputSystem),
                wcscon.wcscsys(outputSystem), 0, 0, coordX, coordY, epoch)

            return outCoords

    raise Exception("inputSystem and outputSystem must be 'J2000', 'B1950'"
                    "or 'GALACTIC'")

#-----------------------------------------------------------------------------
def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg):
    """Calculates minimum and maximum RA, dec coords needed to define a box
    enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg
    coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses
    L{calcAngSepDeg}, so has the same limitations.

    @type RADeg: float
    @param RADeg: RA coordinate of centre of search region
    @type decDeg: float
    @param decDeg: dec coordinate of centre of search region
    @type radiusSkyDeg: float
    @param radiusSkyDeg: radius in degrees on the sky used to define search
        region
    @rtype: list
    @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees
        defining search box

    """

    tolerance = 1e-5  # in degrees on sky
    targetHalfSizeSkyDeg = radiusSkyDeg
    funcCalls = ["calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
               "calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
               "calcAngSepDeg(RADeg, decDeg, RADeg, guess)",
               "calcAngSepDeg(RADeg, decDeg, RADeg, guess)"]
    coords = [RADeg, RADeg, decDeg, decDeg]
    signs = [1.0, -1.0, 1.0, -1.0]
    results = []
    for f, c, sign in zip(funcCalls, coords, signs):
        # Initial guess range
        maxGuess = sign*targetHalfSizeSkyDeg*10.0
        minGuess = sign*targetHalfSizeSkyDeg/10.0
        guessStep = (maxGuess-minGuess)/10.0
        guesses = numpy.arange(minGuess+c, maxGuess+c, guessStep)
        for i in range(50):
            minSizeDiff = 1e6
            bestGuess = None
            for guess in guesses:
                sizeDiff = abs(eval(f)-targetHalfSizeSkyDeg)
                if sizeDiff < minSizeDiff:
                    minSizeDiff = sizeDiff
                    bestGuess = guess
            if minSizeDiff < tolerance:
                break
            else:
                guessRange = abs((maxGuess-minGuess))
                maxGuess = bestGuess+guessRange/4.0
                minGuess = bestGuess-guessRange/4.0
                guessStep = (maxGuess-minGuess)/10.0
                guesses = numpy.arange(minGuess, maxGuess, guessStep)
        results.append(bestGuess)

    RAMax = results[0]
    RAMin = results[1]
    decMax = results[2]
    decMin = results[3]

    return [RAMin, RAMax, decMin, decMax]