1"""Module for coordinate manipulation (conversions, calculations etc.).
2
3(c) 2007-2012 Matt Hilton
4
5(c) 2013-2016 Matt Hilton & Steven Boada
6
7"""
8
9import sys
10import math
11import numpy
12from PyWCSTools import wcscon
13#import IPython
14
15#-----------------------------------------------------------------------------
16def hms2decimal(RAString, delimiter):
17    """Converts a delimited string of Hours:Minutes:Seconds format into decimal
18    degrees.
19
20    @type RAString: string
21    @param RAString: coordinate string in H:M:S format
22    @type delimiter: string
23    @param delimiter: delimiter character in RAString
24    @rtype: float
25    @return: coordinate in decimal degrees
26
27    """
28    # is it in HH:MM:SS format?
29    if delimiter == "":
30        RABits = str(RAString).split()
31    else:
32        RABits = str(RAString).split(delimiter)
33    if len(RABits) > 1:
34        RAHDecimal = float(RABits[0])
35        if len(RABits) > 1:
36            RAHDecimal = RAHDecimal+(float(RABits[1])/60.0)
37        if len(RABits) > 2:
38            RAHDecimal = RAHDecimal+(float(RABits[2])/3600.0)
39        RADeg = (RAHDecimal/24.0)*360.0
40    else:
41        RADeg = float(RAString)
42
43    return RADeg
44
45#-----------------------------------------------------------------------------
46def dms2decimal(decString, delimiter):
47    """Converts a delimited string of Degrees:Minutes:Seconds format into
48    decimal degrees.
49
50    @type decString: string
51    @param decString: coordinate string in D:M:S format
52    @type delimiter: string
53    @param delimiter: delimiter character in decString
54    @rtype: float
55    @return: coordinate in decimal degrees
56
57    """
58    # is it in DD:MM:SS format?
59    if delimiter == "":
60        decBits = str(decString).split()
61    else:
62        decBits = str(decString).split(delimiter)
63    if len(decBits) > 1:
64        decDeg = float(decBits[0])
65        if decBits[0].find("-") != -1:
66            if len(decBits) > 1:
67                decDeg = decDeg-(float(decBits[1])/60.0)
68            if len(decBits) > 2:
69                decDeg = decDeg-(float(decBits[2])/3600.0)
70        else:
71            if len(decBits) > 1:
72                decDeg = decDeg+(float(decBits[1])/60.0)
73            if len(decBits) > 2:
74                decDeg = decDeg+(float(decBits[2])/3600.0)
75    else:
76        decDeg = float(decString)
77
78    return decDeg
79
80#-----------------------------------------------------------------------------
81def decimal2hms(RADeg, delimiter):
82    """Converts decimal degrees to string in Hours:Minutes:Seconds format with
83    user specified delimiter.
84
85    @type RADeg: float
86    @param RADeg: coordinate in decimal degrees
87    @type delimiter: string
88    @param delimiter: delimiter character in returned string
89    @rtype: string
90    @return: coordinate string in H:M:S format
91
92    """
93    hours = (RADeg/360.0)*24
94    #if hours < 10 and hours >= 1:
95    if 1 <= hours < 10:
96        sHours = "0"+str(hours)[0]
97    elif hours >= 10:
98        sHours = str(hours)[:2]
99    elif hours < 1:
100        sHours = "00"
101
102    if str(hours).find(".") == -1:
103        mins = float(hours)*60.0
104    else:
105        mins = float(str(hours)[str(hours).index("."):])*60.0
106    #if mins<10 and mins>=1:
107    if 1 <= mins<10:
108        sMins = "0"+str(mins)[:1]
109    elif mins >= 10:
110        sMins = str(mins)[:2]
111    elif mins < 1:
112        sMins = "00"
113
114    secs = (hours-(float(sHours)+float(sMins)/60.0))*3600.0
115    #if secs < 10 and secs>0.001:
116    if 0.001 < secs < 10:
117        sSecs = "0"+str(secs)[:str(secs).find(".")+4]
118    elif secs < 0.0001:
119        sSecs = "00.000"
120    else:
121        sSecs = str(secs)[:str(secs).find(".")+4]
122    if len(sSecs) < 5:
123        sSecs = sSecs+"00"	# So all to 3dp
124
125    if float(sSecs) == 60.000:
126        sSecs = "00.00"
127        sMins = str(int(sMins)+1)
128    if int(sMins) == 60:
129        sMins = "00"
130        sHours = str(int(sHours)+1)
131
132    return sHours+delimiter+sMins+delimiter+sSecs
133
134#------------------------------------------------------------------------------
135def decimal2dms(decDeg, delimiter):
136    """Converts decimal degrees to string in Degrees:Minutes:Seconds format
137    with user specified delimiter.
138
139    @type decDeg: float
140    @param decDeg: coordinate in decimal degrees
141    @type delimiter: string
142    @param delimiter: delimiter character in returned string
143    @rtype: string
144    @return: coordinate string in D:M:S format
145
146    """
147    # Positive
148    if decDeg > 0:
149        #if decDeg < 10 and decDeg>=1:
150        if 1 <= decDeg < 10:
151            sDeg = "0"+str(decDeg)[0]
152        elif decDeg >= 10:
153            sDeg = str(decDeg)[:2]
154        elif decDeg < 1:
155            sDeg = "00"
156
157        if str(decDeg).find(".") == -1:
158            mins = float(decDeg)*60.0
159        else:
160            mins = float(str(decDeg)[str(decDeg).index("."):])*60
161        #if mins<10 and mins>=1:
162        if 1 <= mins < 10:
163            sMins = "0"+str(mins)[:1]
164        elif mins >= 10:
165            sMins = str(mins)[:2]
166        elif mins < 1:
167            sMins = "00"
168
169        secs = (decDeg-(float(sDeg)+float(sMins)/60.0))*3600.0
170        #if secs<10 and secs>0:
171        if 0 < secs < 10:
172            sSecs = "0"+str(secs)[:str(secs).find(".")+3]
173        elif secs < 0.001:
174            sSecs = "00.00"
175        else:
176            sSecs = str(secs)[:str(secs).find(".")+3]
177        if len(sSecs) < 5:
178            sSecs = sSecs+"0"	# So all to 2dp
179
180        if float(sSecs) == 60.00:
181            sSecs = "00.00"
182            sMins = str(int(sMins)+1)
183        if int(sMins) == 60:
184            sMins = "00"
185            sDeg = str(int(sDeg)+1)
186
187        return "+"+sDeg+delimiter+sMins+delimiter+sSecs
188
189    else:
190        #if decDeg>-10 and decDeg<=-1:
191        if -10 < decDeg <= -1:
192            sDeg = "-0"+str(decDeg)[1]
193        elif decDeg <= -10:
194            sDeg = str(decDeg)[:3]
195        elif decDeg > -1:
196            sDeg = "-00"
197
198        if str(decDeg).find(".") == -1:
199            mins = float(decDeg)*-60.0
200        else:
201            mins = float(str(decDeg)[str(decDeg).index("."):])*60
202        #if mins<10 and mins>=1:
203        if 1 <= mins < 10:
204            sMins = "0"+str(mins)[:1]
205        elif mins >= 10:
206            sMins = str(mins)[:2]
207        elif mins < 1:
208            sMins = "00"
209
210        secs = (decDeg-(float(sDeg)-float(sMins)/60.0))*3600.0
211        #if secs>-10 and secs<0:
212        # so don't get minus sign
213        if -10 < secs < 0:
214            sSecs = "0"+str(secs)[1:str(secs).find(".")+3]
215        elif secs > -0.001:
216            sSecs = "00.00"
217        else:
218            sSecs = str(secs)[1:str(secs).find(".")+3]
219        if len(sSecs) < 5:
220            sSecs = sSecs+"0"	# So all to 2dp
221
222        if float(sSecs) == 60.00:
223            sSecs = "00.00"
224            sMins = str(int(sMins)+1)
225        if int(sMins) == 60:
226            sMins = "00"
227            sDeg = str(int(sDeg)-1)
228
229        return sDeg+delimiter+sMins+delimiter+sSecs
230
231#-----------------------------------------------------------------------------
232def calcAngSepDeg(RADeg1, decDeg1, RADeg2, decDeg2):
233    """Calculates the angular separation of two positions on the sky (specified
234    in decimal degrees) in decimal degrees. Note that RADeg2, decDeg2 can be numpy
235    arrays.
236
237    @type RADeg1: float
238    @param RADeg1: R.A. in decimal degrees for position 1
239    @type decDeg1: float
240    @param decDeg1: dec. in decimal degrees for position 1
241    @type RADeg2: float or numpy array
242    @param RADeg2: R.A. in decimal degrees for position 2
243    @type decDeg2: float or numpy array
244    @param decDeg2: dec. in decimal degrees for position 2
245    @rtype: float or numpy array, depending upon type of RADeg2, decDeg2
246    @return: angular separation in decimal degrees
247
248    """
249
250    a=numpy.sin(numpy.radians(decDeg1))*numpy.sin(numpy.radians(decDeg2))+numpy.cos(numpy.radians(decDeg1))*numpy.cos(numpy.radians(decDeg2))*numpy.cos(numpy.radians(RADeg1-RADeg2))
251    mask=numpy.greater(a, 1.0)
252    if mask.sum() > 0:
253        if type(a) == numpy.ndarray:
254            a[mask]=1.0
255        else:
256            a=1.0
257    mask=numpy.less(a, -1.0)
258    if mask.sum() > 0:
259        if type(a) == numpy.ndarray:
260            a[mask]=-1.0
261        else:
262            a=-1.0
263    r=numpy.degrees(numpy.arccos(a))
264
265    # Above gives nan when RADeg1, decDeg1 == RADeg1, decDeg2
266    indexList=numpy.where(numpy.isnan(r) == True)[0]
267    tolerance=1e-6
268    if len(indexList) > 0:
269        for index in indexList:
270            if type(r) == numpy.ndarray:
271                if type(RADeg2) == numpy.ndarray:
272                    if abs(RADeg1 - RADeg2[index]) < tolerance and abs(decDeg1 -decDeg2[index]) < tolerance:
273                        r[index]=0.0
274                    else:
275                        raise Exception("astCoords: calcAngSepDeg - encountered nan not due to equal RADeg, decDeg coords")
276                elif type(RADeg1) == numpy.ndarray:
277                    if abs(RADeg2 - RADeg1[index]) < tolerance and abs(decDeg2 -decDeg1[index]) < tolerance:
278                        r[index]=0.0
279                    else:
280                        raise Exception("astCoords: calcAngSepDeg - encountered nan not due to equal RADeg, decDeg coords")
281            else:
282                r=0.0
283
284    return r
285
286#-----------------------------------------------------------------------------
287def shiftRADec(ra1, dec1, deltaRA, deltaDec):
288    """Computes new right ascension and declination shifted from the original
289    by some delta RA and delta DEC. Input position is decimal degrees. Shifts
290    (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on
291    an IDL routine of the same name.
292
293    @param ra1: float
294    @type ra1: R.A. in decimal degrees
295    @param dec1: float
296    @type dec1: dec. in decimal degrees
297    @param deltaRA: float
298    @type deltaRA: shift in R.A. in arcseconds
299    @param deltaDec: float
300    @type deltaDec: shift in dec. in arcseconds
301    @rtype: float [newRA, newDec]
302    @return: shifted R.A. and dec.
303
304    """
305
306    d2r = math.pi/180.
307    as2r = math.pi/648000.
308
309    # Convert everything to radians
310    rara1 = ra1*d2r
311    dcrad1 = dec1*d2r
312    shiftRArad = deltaRA*as2r
313    shiftDCrad = deltaDec*as2r
314
315    # Shift!
316    deldec2 = 0.0
317    sindis = math.sin(shiftRArad / 2.0)
318    sindelRA = sindis / math.cos(dcrad1)
319    delra = 2.0*math.asin(sindelRA) / d2r
320
321    # Make changes
322    ra2 = ra1+delra
323    dec2 = dec1 +deltaDec / 3600.0
324
325    return ra2, dec2
326
327#-----------------------------------------------------------------------------
328def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch):
329    """Converts specified coordinates (given in decimal degrees) between J2000,
330    B1950, and Galactic.
331
332    @type inputSystem: string
333    @param inputSystem: system of the input coordinates (either "J2000",
334        "B1950" or "GALACTIC")
335    @type outputSystem: string
336    @param outputSystem: system of the returned coordinates (either "J2000",
337        "B1950" or "GALACTIC")
338    @type coordX: float
339    @param coordX: longitude coordinate in decimal degrees, e.g. R. A.
340    @type coordY: float
341    @param coordY: latitude coordinate in decimal degrees, e.g. dec.
342    @type epoch: float
343    @param epoch: epoch of the input coordinates
344    @rtype: list
345    @return: coordinates in decimal degrees in requested output system
346
347    """
348
349    if inputSystem=="J2000" or inputSystem=="B1950" or inputSystem=="GALACTIC":
350        if outputSystem=="J2000" or outputSystem=="B1950" or \
351            outputSystem=="GALACTIC":
352
353            outCoords=wcscon.wcscon(wcscon.wcscsys(inputSystem),
354                wcscon.wcscsys(outputSystem), 0, 0, float(coordX), float(coordY), epoch)
355
356            return outCoords
357
358    raise Exception("inputSystem and outputSystem must be 'J2000', 'B1950'"
359                    "or 'GALACTIC'")
360
361#-----------------------------------------------------------------------------
362def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg):
363    """Calculates minimum and maximum RA, dec coords needed to define a box
364    enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg
365    coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses
366    L{calcAngSepDeg}, so has the same limitations.
367
368    @type RADeg: float
369    @param RADeg: RA coordinate of centre of search region
370    @type decDeg: float
371    @param decDeg: dec coordinate of centre of search region
372    @type radiusSkyDeg: float
373    @param radiusSkyDeg: radius in degrees on the sky used to define search
374        region
375    @rtype: list
376    @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees
377        defining search box
378
379    """
380
381    tolerance = 1e-5  # in degrees on sky
382    targetHalfSizeSkyDeg = radiusSkyDeg
383    funcCalls = ["calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
384               "calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
385               "calcAngSepDeg(RADeg, decDeg, RADeg, guess)",
386               "calcAngSepDeg(RADeg, decDeg, RADeg, guess)"]
387    coords = [RADeg, RADeg, decDeg, decDeg]
388    signs = [1.0, -1.0, 1.0, -1.0]
389    results = []
390    for f, c, sign in zip(funcCalls, coords, signs):
391        # Initial guess range
392        maxGuess = sign*targetHalfSizeSkyDeg*80.0
393        minGuess = sign*targetHalfSizeSkyDeg/80.0
394        #guessStep = (maxGuess-minGuess)/10.0
395        guesses = numpy.linspace(minGuess+c, maxGuess+c, 1000)
396        converged=False
397        for i in range(50):
398            minSizeDiff = 1e6
399            bestGuess = None
400            for guess in guesses:
401                sizeDiff = abs(eval(f)-targetHalfSizeSkyDeg)
402                if sizeDiff < minSizeDiff:
403                    minSizeDiff = sizeDiff
404                    bestGuess = guess
405            if minSizeDiff < tolerance:
406                converged=True
407                break
408            else:
409                #print sizeDiff, bestGuess, bestGuess-minGuess, bestGuess-maxGuess
410                if bestGuess == None:
411                    raise Exception("bestGuess is None")
412                guessRange = abs((maxGuess-minGuess))
413                maxGuess = bestGuess+guessRange/4.0
414                minGuess = bestGuess-guessRange/4.0
415                # Stop us from searching the wrong side of the coordinate
416                if sign == 1:
417                    if minGuess < c:
418                        minGuess=c+tolerance
419                if sign == -1:
420                    if maxGuess > c:
421                        maxGuess=c-tolerance
422                #guessStep = (maxGuess-minGuess)/20.0
423                guesses = numpy.linspace(minGuess, maxGuess, 1000)
424        if converged == False:
425            raise Exception("calcRADecSearchBox failed to converge")
426        results.append(bestGuess)
427
428    RAMax = results[0]
429    RAMin = results[1]
430    decMax = results[2]
431    decMin = results[3]
432
433    # Sanity check
434    if (RAMax-RAMin)+(2*tolerance) < 2*targetHalfSizeSkyDeg:
435        raise Exception("calcRADecSearchBox failed sanity check")
436
437    return [RAMin, RAMax, decMin, decMax]
438
439