SITools2

"""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]

#!/usr/bin/python

# ******************************************************************************

#    Copyright 2015 IAS - IDOC

#

#    This program is free software: you can redistribute it and/or modify

#    it under the terms of the GNU General Public License as published by

#    the Free Software Foundation, either version 3 of the License, or

#    (at your option) any later version.

#

#    This program is distributed in the hope that it will be useful,

#    but WITHOUT ANY WARRANTY; without even the implied warranty of

#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

#    GNU General Public License for more details.

#

#    You should have received a copy of the GNU General Public License

#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

#

# ******************************************************************************

'''

This script updates the content of a fits table, adding new columns

and/or rows (i.e. objects) to it, by considering as input a user-defined ASCII table.                       

The new columns (rows) defined in the ascii file are appended at the end

(bottom) of the fits table.                                             

IMPORTANT: The 1st line of the ASCII table must contain the names   

of the columns, and must be UNCOMMENTED!                            

NOTE: Ra and DEC must be in **decimal degrees**, both in FITS and                                          

ASCII tables.

The syntax is:

$ python edit_FITS.py <table>.fits <ascii_file>

@author: Alessandro NASTASI for IAS - IDOC 

@date: 21/05/2015

'''

__author__ = "Alessandro Nastasi"

__credits__ = ["Alessandro Nastasi"]

__license__ = "GPL"

__version__ = "1.0"

__date__ = "21/05/2015"

import numpy as np

import os, sys, re, time

import string

import asciidata

import pyfits

from datetime import date

#Use the provided astCoords.py file rather than the default module of astLib, 

#since the calcAngSepDeg() of the latter works only for separation <90 deg 

#(tangent plane projection approximation)

import astCoords

class bcolors:

    HEADER = '\033[95m'

    OKBLUE = '\033[94m'

    OKGREEN = '\033[92m'

    WARNING = '\033[93m'

    FAIL = '\033[91m'

    ENDC = '\033[0m'

#Dictionary containing the FORMAT and UNITS of all (or most of) the fields

_FIELDS_DICTIONARY = { 

  'INDEX': { 'format': 'I', 'unit': 'None' },

  'COORD_SOURCE': { 'format': '5A', 'unit': 'None' },

  'x':{ 'format': 'E', 'unit': 'None' },

  'y':{ 'format': 'E', 'unit': 'None' },

  'z':{ 'format': 'E', 'unit': 'None' },

  # ****** ACT ******

  'ACT_INDEX': { 'format': 'I', 'unit': 'None' },

  'INDEX_ACT': { 'format': 'I', 'unit': 'None' },

  'CATALOG': { 'format': '7A', 'unit': 'None' },

  #'NAME': { 'format': '18A', 'unit': 'None' },

  #'GLON': { 'format': 'E', 'unit': 'degrees' },

  #'GLAT': { 'format': 'E', 'unit': 'degrees' },

  #'RA': { 'format': 'E', 'unit': 'degrees' },

  #'DEC': { 'format': 'E', 'unit': 'degrees' },

  'SNR': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'ERR_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'M500': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'ERR_M500': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'ERR_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'THETA': { 'format': 'E', 'unit': 'arcmin' },

  #'PAPER': { 'format': '56A', 'unit': 'None' },      # Use PAPER in SPT as bigger format '59A'

  'ACT_CATALOG': { 'format': '7A', 'unit': 'None' },

  'ACT_NAME': { 'format': '18A', 'unit': 'None' },

  'ACT_GLON': { 'format': 'E', 'unit': 'degrees' },

  'ACT_GLAT': { 'format': 'E', 'unit': 'degrees' },

  'ACT_RA': { 'format': 'E', 'unit': 'degrees' },

  'ACT_DEC': { 'format': 'E', 'unit': 'degrees' },

  'ACT_SNR': { 'format': 'E', 'unit': 'None' },

  'ACT_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'ACT_ERR_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'ACT_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'ACT_REDSHIFT_REF': { 'format': '19A', 'unit': 'None' },

  'ACT_M500': { 'format': 'E', 'unit': '10^14 h^-1 solar mass' },

  'ACT_ERR_M500': { 'format': 'E', 'unit': '10^14 h^-1 solar mass' },

  'ACT_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'ACT_ERR_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'ACT_THETA': { 'format': 'E', 'unit': 'arcmin' },

  'ACT_PAPER': { 'format': '56A', 'unit': 'None' },

  # ****** AMI ******

  'INDEX_AMI': { 'format': 'I', 'unit': 'None' },

  'AMI_INDEX': { 'format': 'I', 'unit': 'None' },

  #'NAME': { 'format': '18A', 'unit': 'None' },

  #'RA': { 'format': 'E', 'unit': 'Degrees' },

  #'DEC': { 'format': 'E', 'unit': 'Degrees' },

  #'GLON': { 'format': 'E', 'unit': 'Degrees' },

  #'GLAT': { 'format': 'E', 'unit': 'Degrees' },

  #'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  #'REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  #'ALT_NAME': { 'format': '60A', 'unit': 'None' },

  #'COORD_SOURCE': { 'format': '5A', 'unit': 'None' },

  'AMI_NAME': { 'format': '18A', 'unit': 'None' },

  'AMI_RA': { 'format': 'E', 'unit': 'Degrees' },

  'AMI_DEC': { 'format': 'E', 'unit': 'Degrees' },

  'AMI_GLON': { 'format': 'E', 'unit': 'Degrees' },

  'AMI_GLAT': { 'format': 'E', 'unit': 'Degrees' },

  'AMI_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'AMI_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'AMI_REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  'AMI_ALT_NAME': { 'format': '60A', 'unit': 'None' },

  # ****** CARMA ******

  'INDEX_CARMA': { 'format': 'I', 'unit': 'None' },

  'CARMA_INDEX': { 'format': 'I', 'unit': 'None' },

  #'NAME': { 'format': '18A', 'unit': 'None' },

  #'RA': { 'format': 'E', 'unit': 'Degrees' },

  #'DEC': { 'format': 'E', 'unit': 'Degrees' },

  #'GLON': { 'format': 'E', 'unit': 'Degrees' },

  #'GLAT': { 'format': 'E', 'unit': 'Degrees' },

  #'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  #'REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  #'COORD_SOURCE': { 'format': '5A', 'unit': 'None' },  

  'CARMA_NAME': { 'format': '18A', 'unit': 'None' },

  'CARMA_RA': { 'format': 'E', 'unit': 'Degrees' },

  'CARMA_DEC': { 'format': 'E', 'unit': 'Degrees' },

  'CARMA_GLON': { 'format': 'E', 'unit': 'Degrees' },

  'CARMA_GLAT': { 'format': 'E', 'unit': 'Degrees' },

  'CARMA_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'CARMA_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'CARMA_REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  'CARMA_M500': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'CARMA_ERR_M500': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  #****** PSZ1 ******

  'PSZ1_INDEX': { 'format': 'I', 'unit': 'None' },

  'INDEX_PSZ1': { 'format': 'I', 'unit': 'None' },

  'NAME': { 'format': '18A', 'unit': 'None' },

  'GLON': { 'format': 'D', 'unit': 'degrees' },

  'GLAT': { 'format': 'D', 'unit': 'degrees' },

  'RA': { 'format': 'D', 'unit': 'degrees' },

  'DEC': { 'format': 'D', 'unit': 'degrees' },

  'RA_MCXC': { 'format': 'E', 'unit': 'degrees' },

  'DEC_MCXC': { 'format': 'E', 'unit': 'degrees' },

  'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'REDSHIFT_SOURCE': { 'format': 'I', 'unit': 'None' },

  'REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  'ALT_NAME': { 'format': '66A', 'unit': 'None' },

  'YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'ERRP_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'ERRM_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'ERRP_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'ERRM_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  'ERR_S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  'Y_PSX_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'SN_PSX': { 'format': 'E', 'unit': 'None' },

  'PIPELINE': { 'format': 'I', 'unit': 'None' },

  'PIPE_DET': { 'format': 'I', 'unit': 'None' },

  'PCCS': { 'format': 'L', 'unit': 'None' },

  'VALIDATION': { 'format': 'I', 'unit': 'None' },

  'ID_EXT': { 'format': '25A', 'unit': 'None' },

  'POS_ERR': { 'format': 'E', 'unit': 'arcmin' },

  #'SNR': { 'format': 'E', 'unit': 'None' },

  'COSMO': { 'format': 'L', 'unit': 'None' },

  'COMMENT': { 'format': 'L', 'unit': 'None' },

  'QN': { 'format': 'E', 'unit': 'None' }, 

  'PSZ1_NAME': { 'format': '18A', 'unit': 'None' },

  'PSZ1_GLON': { 'format': 'D', 'unit': 'degrees' },

  'PSZ1_GLAT': { 'format': 'D', 'unit': 'degrees' },

  'PSZ1_RA': { 'format': 'D', 'unit': 'degrees' },

  'PSZ1_DEC': { 'format': 'D', 'unit': 'degrees' },

  'PSZ1_RA_MCXC': { 'format': 'E', 'unit': 'degrees' },

  'PSZ1_DEC_MCXC': { 'format': 'E', 'unit': 'degrees' },

  'PSZ1_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'PSZ1_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'PSZ1_REDSHIFT_SOURCE': { 'format': 'I', 'unit': 'None' },

  'PSZ1_REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  'PSZ1_ALT_NAME': { 'format': '66A', 'unit': 'None' },

  'PSZ1_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PSZ1_ERRP_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PSZ1_ERRM_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PSZ1_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'PSZ1_ERRP_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'PSZ1_ERRM_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'PSZ1_S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  'PSZ1_ERR_S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  'PSZ1_Y_PSX_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PSZ1_SN_PSX': { 'format': 'E', 'unit': 'None' },

  'PSZ1_PIPELINE': { 'format': 'I', 'unit': 'None' },

  'PSZ1_PIPE_DET': { 'format': 'I', 'unit': 'None' },

  'PSZ1_PCCS': { 'format': 'L', 'unit': 'None' },

  'PSZ1_VALIDATION': { 'format': 'I', 'unit': 'None' },

  'PSZ1_ID_EXT': { 'format': '25A', 'unit': 'None' },

  'PSZ1_POS_ERR': { 'format': 'E', 'unit': 'arcmin' },

  'PSZ1_SNR': { 'format': 'E', 'unit': 'None' },

  'PSZ1_COSMO': { 'format': 'L', 'unit': 'None' },

  'PSZ1_COMMENT': { 'format': 'L', 'unit': 'None' },

  'PSZ1_QN': { 'format': 'E', 'unit': 'None' }, 

  # ****** PSZ2 ******

  'PSZ2_INDEX': { 'format': 'I', 'unit': 'None' },

  'INDEX_PSZ2': { 'format': 'I', 'unit': 'None' },

  #'NAME': { 'format': '18A', 'unit': 'None' },

  #'GLON': { 'format': 'D', 'unit': 'degrees' },

  #'GLAT': { 'format': 'D', 'unit': 'degrees' },

  #'RA': { 'format': 'D', 'unit': 'degrees' },

  #'DEC': { 'format': 'D', 'unit': 'degrees' },

  #'POS_ERR': { 'format': 'E', 'unit': 'arcmin' },

  #'SNR': { 'format': 'E', 'unit': 'None' },

  #'PIPELINE': { 'format': 'I', 'unit': 'None' },

  #'PIPE_DET': { 'format': 'I', 'unit': 'None' },

  'PCCS2': { 'format': 'L', 'unit': 'None' },

  'PSZ': { 'format': 'I', 'unit': 'None' },

  'IR_FLAG': { 'format': 'I', 'unit': 'None' },

  'Q_NEURAL': { 'format': 'E', 'unit': 'None' },

  'Y5R500': { 'format': 'E', 'unit': '10^-3 arcmin^2' },

  'Y5R500_ERR': { 'format': 'E', 'unit': '10^-3 arcmin^2' },

  'PSZ2_VALIDATION': { 'format': 'I', 'unit': 'None' },

  'REDSHIFT_ID': { 'format': '25A', 'unit': 'None' },

  #'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'MSZ': { 'format': 'E', 'unit': '10^14 Msol' },

  'MSZ_ERR_UP': { 'format': 'E', 'unit': '10^14 Msol' },

  'MSZ_ERR_LOW': { 'format': 'E', 'unit': '10^14 Msol' },

  'MCXC': { 'format': '25A', 'unit': 'None' },

  'REDMAPPER': { 'format': '25A', 'unit': 'None' },

  'ACT': { 'format': '25A', 'unit': 'None' },

  'SPT': { 'format': '25A', 'unit': 'None' },

  'WISE_SIGNF': { 'format': 'E', 'unit': 'None' },

  'WISE_FLAG': { 'format': 'I', 'unit': 'None' },

  'AMI_EVIDENCE': { 'format': 'E', 'unit': 'None' },

  #'COSMO': { 'format': 'L', 'unit': 'None' },

  'PSZ2_COMMENT': { 'format': '128A', 'unit': 'None' },

  'PSZ2_NAME': { 'format': '18A', 'unit': 'None' },

  'PSZ2_GLON': { 'format': 'D', 'unit': 'degrees' },

  'PSZ2_GLAT': { 'format': 'D', 'unit': 'degrees' },

  'PSZ2_RA': { 'format': 'D', 'unit': 'degrees' },

  'PSZ2_DEC': { 'format': 'D', 'unit': 'degrees' },

  'PSZ2_POS_ERR': { 'format': 'E', 'unit': 'arcmin' },

  'PSZ2_SNR': { 'format': 'E', 'unit': 'None' },

  'PSZ2_PIPELINE': { 'format': 'I', 'unit': 'None' },

  'PSZ2_PIPE_DET': { 'format': 'I', 'unit': 'None' },

  'PSZ2_PCCS2': { 'format': 'L', 'unit': 'None' },

  'PSZ2_PSZ': { 'format': 'I', 'unit': 'None' },

  'PSZ2_IR_FLAG': { 'format': 'I', 'unit': 'None' },

  'PSZ2_Q_NEURAL': { 'format': 'E', 'unit': 'None' },

  'PSZ2_Y5R500': { 'format': 'E', 'unit': '10^-3 arcmin^2' },

  'PSZ2_Y5R500_ERR': { 'format': 'E', 'unit': '10^-3 arcmin^2' },

  #'PSZ2_VALIDATION': { 'format': 'I', 'unit': 'None' },

  'PSZ2_REDSHIFT_ID': { 'format': '25A', 'unit': 'None' },

  'PSZ2_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'PSZ2_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'PSZ2_MSZ': { 'format': 'E', 'unit': '10^14 Msol' },

  'PSZ2_MSZ_ERR_UP': { 'format': 'E', 'unit': '10^14 Msol' },

  'PSZ2_MSZ_ERR_LOW': { 'format': 'E', 'unit': '10^14 Msol' },

  'PSZ2_MCXC': { 'format': '25A', 'unit': 'None' },

  'PSZ2_REDMAPPER': { 'format': '25A', 'unit': 'None' },

  'PSZ2_ACT': { 'format': '25A', 'unit': 'None' },

  'PSZ2_SPT': { 'format': '25A', 'unit': 'None' },

  'PSZ2_WISE_SIGNF': { 'format': 'E', 'unit': 'None' },

  'PSZ2_WISE_FLAG': { 'format': 'I', 'unit': 'None' },

  'PSZ2_AMI_EVIDENCE': { 'format': 'E', 'unit': 'None' },

  'PSZ2_COSMO': { 'format': 'L', 'unit': 'None' },

  #'PSZ2_COMMENT': { 'format': '128A', 'unit': 'None' },

  # ****** PLCK ******

  'PLCK_INDEX': { 'format': 'I', 'unit': 'None' },

  'INDEX_PLCK': { 'format': 'I', 'unit': 'None' },

  #'NAME': { 'format': '18A', 'unit': 'None' },

  #'GLON': { 'format': 'D', 'unit': 'degrees' },

  #'GLAT': { 'format': 'D', 'unit': 'degrees' },

  #'RA': { 'format': 'D', 'unit': 'degrees' },

  #'DEC': { 'format': 'D', 'unit': 'degrees' },

  #'RA_MCXC': { 'format': 'E', 'unit': 'degrees' },

  #'DEC_MCXC': { 'format': 'E', 'unit': 'degrees' },

  #'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  #'REDSHIFT_SOURCE': { 'format': 'I', 'unit': 'None' },

  #'REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  #'ALT_NAME': { 'format': '66A', 'unit': 'None' },

  #'YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  #'ERRP_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  #'ERRM_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  #'M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  #'ERRP_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  #'ERRM_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  #'S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  #'ERR_S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  #'Y_PSX_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  #'SN_PSX': { 'format': 'E', 'unit': 'None' },

  #'PIPELINE': { 'format': 'I', 'unit': 'None' },

  #'PIPE_DET': { 'format': 'I', 'unit': 'None' },

  #'PCCS': { 'format': 'L', 'unit': 'None' },

  #'VALIDATION': { 'format': 'I', 'unit': 'None' },

  #'ID_EXT': { 'format': '25A', 'unit': 'None' },

  #'POS_ERR': { 'format': 'E', 'unit': 'arcmin' },

  #'SNR': { 'format': 'E', 'unit': 'None' },

  #'COSMO': { 'format': 'L', 'unit': 'None' },

  #'COMMENT': { 'format': 'L', 'unit': 'None' },

  #'QN': { 'format': 'E', 'unit': 'None' }, 

  'PLCK_NAME': { 'format': '18A', 'unit': 'None' },

  'PLCK_GLON': { 'format': 'D', 'unit': 'degrees' },

  'PLCK_GLAT': { 'format': 'D', 'unit': 'degrees' },

  'PLCK_RA': { 'format': 'D', 'unit': 'degrees' },

  'PLCK_DEC': { 'format': 'D', 'unit': 'degrees' },

  'PLCK_RA_MCXC': { 'format': 'E', 'unit': 'degrees' },

  'PLCK_DEC_MCXC': { 'format': 'E', 'unit': 'degrees' },

  'PLCK_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'PLCK_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'PLCK_REDSHIFT_SOURCE': { 'format': 'I', 'unit': 'None' },

  'PLCK_REDSHIFT_REF': { 'format': '36A', 'unit': 'None' },

  'PLCK_ALT_NAME': { 'format': '66A', 'unit': 'None' },

  'PLCK_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PLCK_ERRP_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PLCK_ERRM_YZ_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PLCK_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'PLCK_ERRP_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'PLCK_ERRM_M_YZ_500': { 'format': 'E', 'unit': '10^14 solar mass' },

  'PLCK_S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  'PLCK_ERR_S_X': { 'format': 'E', 'unit': 'erg/s/cm2' },

  'PLCK_Y_PSX_500': { 'format': 'E', 'unit': '10^-4 arcmin squared' },

  'PLCK_SN_PSX': { 'format': 'E', 'unit': 'None' },

  'PLCK_PIPELINE': { 'format': 'I', 'unit': 'None' },

  'PLCK_PIPE_DET': { 'format': 'I', 'unit': 'None' },

  'PLCK_PCCS': { 'format': 'L', 'unit': 'None' },

  'PLCK_VALIDATION': { 'format': 'I', 'unit': 'None' },

  'PLCK_ID_EXT': { 'format': '25A', 'unit': 'None' },

  'PLCK_POS_ERR': { 'format': 'E', 'unit': 'arcmin' },

  'PLCK_SNR': { 'format': 'E', 'unit': 'None' },

  'PLCK_COSMO': { 'format': 'L', 'unit': 'None' },

  'PLCK_COMMENT': { 'format': 'L', 'unit': 'None' },

  'PLCK_QN': { 'format': 'E', 'unit': 'None' }, 

  # ****** SPT ******

  'SPT_INDEX': { 'format': 'I', 'unit': 'None' },

  'INDEX_SPT': { 'format': 'I', 'unit': 'None' },

  #'CATALOG': { 'format': '7A', 'unit': 'None' },

  #'NAME': { 'format': '16A', 'unit': 'None' },

  #'GLON': { 'format': 'E', 'unit': 'degrees' },

  #'GLAT': { 'format': 'E', 'unit': 'degrees' },

  #'RA': { 'format': 'E', 'unit': 'degrees' },

  #'DEC': { 'format': 'E', 'unit': 'degrees' },

  #'SNR': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT': { 'format': 'E', 'unit': 'None' },

  #'ERR_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  #'REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'REDSHIFT_LIMIT': { 'format': 'E', 'unit': 'None' },

  'M500_fidCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'ERR_M500_fidCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'M500_PlanckCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'ERR_M500_PlanckCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'ERR_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'LX': { 'format': 'E', 'unit': '10^44 erg/s' },

  'YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'ERR_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  #'THETA': { 'format': 'E', 'unit': 'arcmin' },

  'PAPER': { 'format': '59A', 'unit': 'None' },

  'XRAY': { 'format': 'L', 'unit': 'None' },

  'STRONG_LENS': { 'format': 'L', 'unit': 'None' },

  'SPT_CATALOG': { 'format': '7A', 'unit': 'None' },

  'SPT_NAME': { 'format': '16A', 'unit': 'None' },

  'SPT_GLON': { 'format': 'E', 'unit': 'degrees' },

  'SPT_GLAT': { 'format': 'E', 'unit': 'degrees' },

  'SPT_RA': { 'format': 'E', 'unit': 'degrees' },

  'SPT_DEC': { 'format': 'E', 'unit': 'degrees' },

  'SPT_SNR': { 'format': 'E', 'unit': 'None' },

  'SPT_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'SPT_ERR_REDSHIFT': { 'format': 'E', 'unit': 'None' },

  'SPT_REDSHIFT_TYPE': { 'format': '5A', 'unit': 'None' },

  'SPT_REDSHIFT_REF': { 'format': '19A', 'unit': 'None' },

  'SPT_REDSHIFT_LIMIT': { 'format': 'E', 'unit': 'None' },

  'SPT_XRAY': { 'format': 'L', 'unit': 'None' },

  'SPT_STRONG_LENS': { 'format': 'L', 'unit': 'None' },

  'SPT_M500_fidCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'SPT_ERR_M500_fidCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'SPT_M500_PlanckCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'SPT_ERR_M500_PlanckCosmo': { 'format': 'E', 'unit': '10^14 h70^-1 solar mass' },

  'SPT_LX': { 'format': 'E', 'unit': '10^44 erg/s' },

  'SPT_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'SPT_ERR_YSZ': { 'format': 'E', 'unit': '10^-6 arcmin squared' },

  'SPT_THETA': { 'format': 'E', 'unit': 'arcmin' },

  'SPT_PAPER': { 'format': '59A', 'unit': 'None' }

  }

#Name of fields (in FITS/ASCII) sometimes called individually in the script

#For Mass and Err_Mass, more than one label could be defined just as elements of the arrays

name_mass_key = ['M500']

name_errMass_key = ['ERR_M500']

name_ra_key = 'RA'

name_dec_key = 'DEC'

name_coordinates_keys =  ['RA_MCXC', 'DEC_MCXC', name_ra_key, name_dec_key]

name_Name_key = 'NAME'  

name_index_key = 'INDEX'

name_catalog_key = 'CATALOG'

name_redshift_key = 'REDSHIFT'

name_zLimit_key = 'REDSHIFT_LIMIT'

name_zErr_key = 'ERR_REDSHIFT'

name_zType_key = 'REDSHIFT_TYPE'

name_zRef_key = 'REDSHIFT_REF'

name_altName_key = 'ALT_NAME'

name_paper_key = 'PAPER'

#Undef values for some kind of fields

_UNDEF_VALUES_ = {

  'FLOAT' : {np.nan},

  'INT' : {-1},

  'STRING' : {'NULL'},

  name_zType_key : {'undef'},

  'PIPELINE' : {0},

  'PIPE_DET' : {0}

  }

def remove_duplicated_names(string):

  '''

  This function removes duplicated names of a string, assuming they are separated by ';'

  In addition, it takes out 'NULL', 'NaN', 'False' from the final, composite string.

  It is used for the creation of ALT_NAME field.

  '''

  string = string.replace('; ',';')

  tmp = [item for item in string.split(';') if item.upper() not in ["-", "NULL", "NAN", "NONE", "FALSE"] and len(item)>0 ]

  # *** This lines of code help preserving the order ot the names ***

  tmp_uniq = []

  set_tmp = set()

  for item in tmp:

    if item not in set_tmp:

      tmp_uniq.append(item)

      set_tmp.add(item)

  # ******************************************************************

  if len(tmp)==0: new_string = '-'

  else: new_string =  "; ".join(tmp_uniq)

  return new_string

def set_undef_values(fits_data):

  '''

  Set the proper 'undef' values according to the format/name of the field

  '''

  print "\n\t>> Checking/setting undefined values for the different fields ..."

  for i, name in enumerate(fits_data.names):

    sys.stdout.write('\t%i/%i > %s (format %s) : Done                                        \r' % (i+1, len(fits_data.names), name, fits_data.formats[i]))

    sys.stdout.flush()

    for j in range(fits_data.size):

      if name == name_index_key:

	if fits_data[name_Name_key][j] <= 0: fits_data[j][i] = -1

      elif name == name_redshift_key and fits_data[name][j] == -1.0:

	fits_data[j][i] = np.nan

      elif name.find(name_zType_key) >= 0 and str(fits_data[name][j]) == 'Null':

	fits_data[j][i] = "undef"

      elif fits_data.formats[i] in 'EDJ':

	if str(fits_data[j][i]) in ['-1.6375e+30','-1.63750e+30', '-1.6375E+30', '-1.63750E+30', 'None', 'NULL']:

	  fits_data[j][i] = np.nan

      elif fits_data.formats[i].find('A') >= 0:

	fits_data[j][i] = remove_duplicated_names(fits_data[j][i])

	if str(fits_data[j][i]).upper() in ["", "0.0", "NULL", "NAN", "NONE", "FALSE"]  or str(fits_data[j][i]) == 'False':

	  fits_data[j][i] = "-"

      elif name in ['PIPELINE','PIPE_DET']:

	if fits_data[j][i] <= 0: fits_data[j][i] = 0

  print '\n'

  return fits_data

def recreate_reformatted_column(hdulist, field_name, new_format, new_vector):

  '''

  Update the length (format) of a 'STRING' (format = 'xA') FIELD. 

  The only way, though, is to re-create the column with the new format.

  It is used during the creation of NAME, ALT_NAME or REDSHIFT_REF.

  '''

  name_vec = []

  format_vec = []

  unit_vec = []

  fits_keywds = hdulist.data.names

  coldefs = pyfits.ColDefs(hdulist.columns)

  #Store attributes of the kewyords after FIELD

  for j in range(fits_keywds.index(field_name)+1, len(fits_keywds)):

    name_vec.append(coldefs.names[j])

    format_vec.append(coldefs.formats[j])

    unit_vec.append(coldefs.units[j])

  #Delete the kewyords after FIELD

  tmp = 0

  for j in range(fits_keywds.index(field_name)+1, len(fits_keywds)):

    coldefs.del_col(name_vec[tmp])

    tmp+=1

  #Delete FIELD

  coldefs.del_col(field_name)

  #Re-create FIELD with the new format

  col_tmp = pyfits.Column(name = field_name, format = new_format, unit = 'None', array = new_vector)

  coldefs.add_col(col_tmp)

  hdulist.columns = coldefs

  #Re-create all the kewyords after FIELD, with their attributes

  tmp = 0

  data_vec_tmp = []

  for j in range(fits_keywds.index(field_name)+1, len(fits_keywds)):

    data_vec_tmp = hdulist.data[name_vec[tmp]]

    col_tmp = pyfits.Column(name = name_vec[tmp], format = format_vec[tmp], unit = unit_vec[tmp], array = data_vec_tmp)

    coldefs.add_col(col_tmp)

    tmp +=1

    data_vec_tmp = []

  hdulist = pyfits.new_table(coldefs)

  return hdulist

'''

  *** >> START << ***

'''

if (len(sys.argv) > 1):

    fits_file = sys.argv[1] 

    ascii_file = sys.argv[2] 

else:

    print bcolors.WARNING +  "\n\tSintax:\t$ python edit_FITS.py <fits_file> <ascii_file>\n" + bcolors.ENDC

    os._exit(0)

#Open the output file

file_report_name = 'summary_updates.tab'

file_report = open(file_report_name, 'w')

question = bcolors.OKBLUE+ "[Q]" + bcolors.ENDC

info = bcolors.WARNING+ "[I]" + bcolors.ENDC

error = bcolors.FAIL+ "[ERR]" + bcolors.ENDC

#User can define the columns delimiter in the ASCII table.

delim=raw_input("\n%s Please enter the column delimiter of the ASCII table (default is ','):\t" % question)

if not delim:

# Read the ascii table, with the structure: ascii_table[COLUMNS][ROWS]

  ascii_table=asciidata.open(ascii_file, 'r', delimiter=',')

else:

  ascii_table=asciidata.open(ascii_file, 'r', delimiter=delim)

Ncol_ascii = ascii_table.ncols

Nrows_ascii = (ascii_table.nrows) - 1 #1st excluded because of the header

print "\n\t\t **** ASCII table details ****"

print "\t\t Number of columns: %s" % (Ncol_ascii)

print "\t\t Number of rows: %s" % (Nrows_ascii)

print "\t\t **** **** **** **** **** ****"

ascii_keywds=[]

keys_form_unit = {}

for i in range(ascii_table.ncols):

  tmpKey = str(ascii_table[i][0]).strip()

  ascii_keywds.append(tmpKey)

  if tmpKey in _FIELDS_DICTIONARY:

    keys_form_unit[tmpKey] = {}

    keys_form_unit[tmpKey]['TFORM'] = _FIELDS_DICTIONARY[tmpKey]['format']

    keys_form_unit[tmpKey]['TUNIT'] = _FIELDS_DICTIONARY[tmpKey]['unit']

#Read the fits table

hdulist = pyfits.open(fits_file)

fits_header = hdulist[1].header	# HEADER 

fits_data = hdulist[1].data 		# DATA 

#Number of columns in FITS table

Ncol_fits = int(fits_header['TFIELDS'])

#Number of rows in FITS table

Nrows_fits = fits_header['NAXIS2']

print "\n\t\t **** FITS table details ****"

print "\t\t Number of columns: %s" % (Ncol_fits)

print "\t\t Number of rows: %s" % (Nrows_fits)

print "\t\t **** *** *** *** *** *** ***"

#Fits keywords read from the header

fits_keywds=[]

original_fits_keywds = []

for i in range(Ncol_fits):

  original_fits_keywds.append(fits_data.names[i])

  fits_keywds.append(fits_data.names[i])

#Find the keywords written in the ASCII file and the corresponding columns in the FITS table...

common_keywds=[]

commonKeywds_index=[]

keywds_to_update=[]

for j in range(Ncol_fits):

  if fits_keywds[j] in ascii_keywds:

    common_keywds.append(fits_keywds[j])

    commonKeywds_index.append(j+1)

    #Only the fields with new values will be updated. Also NAME, RA and DEC are allowed to change

    keywds_to_update.append(fits_keywds[j])

print "\n\t%s The following keyword(s) will be updated in the FITS table: " % info , keywds_to_update

#...also selecting the NEW keywords defined in the ASCII file...

keywds_to_add=[item for item in ascii_keywds if item not in fits_keywds]

print "\n\t%s The following new keyword(s) will be added to the FITS table: " % info , keywds_to_add

#To associate TFORM and TUNIT to each field, first look into _FIELDS_DICTIONARY

#If nothing is found ther, ask the user to enter them manually

for i in range(len(keywds_to_add)):

  if keywds_to_add[i] not in _FIELDS_DICTIONARY:

    keys_form_unit[keywds_to_add[i]] = {}

    message = "\n%s Please enter the format (\'TFORM\') of the new field \"%s\" (e.g.: 5A, E, L, ...): " % (question, keywds_to_add[i])

    keys_form_unit[keywds_to_add[i]]['TFORM'] = raw_input(message)

    message = "\n%s Please enter the unit (\'TUNIT\') of the new field \"%s\" (e.g.: None, arcmin, ...): " % (question, keywds_to_add[i])

    keys_form_unit[keywds_to_add[i]]['TUNIT'] = raw_input(message)

  else:

    keys_form_unit[keywds_to_add[i]] = {}

    keys_form_unit[keywds_to_add[i]]['TFORM'] = _FIELDS_DICTIONARY[keywds_to_add[i]]['format']

    keys_form_unit[keywds_to_add[i]]['TUNIT'] = _FIELDS_DICTIONARY[keywds_to_add[i]]['unit']

  # ...to be appended into the 'fits_keywds' array

  fits_keywds.append(keywds_to_add[i])

'''

  *** Add the NEW COLUMNS to FITS table ***

'''

#Initialize new columns

a_tmp = []

coldefs = pyfits.ColDefs(hdulist[1].columns)

columns = []

for keys in keywds_to_add:

  if keys_form_unit[keys]['TFORM'] == 'E' or keys_form_unit[keys]['TFORM'] == 'D':

    a_tmp = [-1.6375E+30] * Nrows_fits # Initialize Float with empty array

  elif keys_form_unit[keys]['TFORM'] == 'I':

    a_tmp = [-1] * Nrows_fits # Initialize Integer with -1 array

  elif keys_form_unit[keys]['TFORM'] == 'L':

    a_tmp = [False] * Nrows_fits # Initialize logical with True array

  elif keys_form_unit[keys]['TFORM'].find('A') >= 0:

    a_tmp = ['Null'] * Nrows_fits

  while True:

    #Check between field format and values.

    try:

      col_tmp = pyfits.Column(name=keys, format=keys_form_unit[keys]['TFORM'], unit=keys_form_unit[keys]['TUNIT'], array=a_tmp)

      columns.append(col_tmp)

      break

    except ValueError:

      print bcolors.FAIL+ "\n\t\t*** FORMAT INCONSISTENT WITH DATA ***" + bcolors.ENDC

      keys_form_unit[keys]['TFORM'] = raw_input("\n%s Please, enter again the format (\'TFORM\') of the new field \"%s\": " % (question, keys))

'''

  *** 1st data UPDATE: new fields added as new columns ***

'''

#New for Pyfits > 2.3

for i in columns: coldefs.add_col(i)

hdulist = pyfits.new_table(coldefs)

#Old Python version

#hdulist = pyfits.BinTableHDU.from_columns(coldefs)

fits_data = hdulist.data

'''

  *** Object identification via POSITION matching, NAME or INDEX  ***

'''

match_option = False

match_radius = 300.0 # default = 5 arcmin

name_index_fits = ''

name_index_ascii = ''

print '\n%s Which method do you want to use for the object matching: by POSITION (1) by NAME (2) or by INDEX (3)?' % question

while match_option == False:

  message = "\n\t-> Please enter 1, 2 or 3:   "

  method = raw_input(message)

  if method == '1': 

    #Check if RA & DEC are actually in FITS and ASCII tables

    if name_ra_key not in fits_data.names or name_dec_key not in fits_data.names or name_ra_key not in ascii_keywds or name_dec_key not in ascii_keywds:

      print bcolors.FAIL+ "\n\t>> NO %s and %s found in FITS and ASCII tables: POSITION matching not possible <<" % (name_ra_key, name_dec_key) + bcolors.ENDC

    else:

      match_option = method

      match_radius = float(raw_input('\n\t%s Please enter the match radius (in arcsec): ' % question))

  elif method == '2' : match_option = method

  elif method == '3' :

    check_name_index_fits = False

    while check_name_index_fits == False:

      name_index_fits = raw_input('\n\t-> Please enter the column name of the INDEX in the FITS file: ')

      if name_index_fits not in fits_keywds:

	print bcolors.FAIL+ "\n\t*** '%s' NOT in FITS Keywords ***" % name_index_fits+ bcolors.ENDC

      else:

	check_name_index_fits = True

	index_fits = np.array( fits_data[name_index_fits] )

    check_name_index_ascii = False

    while check_name_index_ascii == False:

      name_index_ascii = raw_input('\n\t-> Please enter the column name of the INDEX in the ASCII file: ')

      if name_index_ascii not in ascii_keywds:

	print bcolors.FAIL+ "\n\t*** '%s' NOT in ASCII Keywords ***" % name_index_ascii+ bcolors.ENDC

      else:

	check_name_index_ascii = True

	index_ascii = [ (ascii_table[k][j]) for k in range(ascii_table.ncols) if ascii_table[k][0] == name_index_ascii for j in range(1,ascii_table.nrows) ]

    match_option = method

  else: print bcolors.FAIL+ "\n\t*** Wrong option ***"+ bcolors.ENDC

name_fits = np.array(fits_data[name_Name_key])

ra_fits = np.array(fits_data[ name_ra_key ])

dec_fits = np.array(fits_data[ name_dec_key ])

name_ascii = []

ra_ascii = []

dec_ascii = []

for k in range(ascii_table.ncols):

  if ascii_keywds[k]==name_Name_key:

    for j in range(ascii_table.nrows -1): name_ascii.append((ascii_table[k][j+1]).strip())

  if ascii_keywds[k]==name_ra_key:

    for j in range(ascii_table.nrows -1): ra_ascii.append(float(ascii_table[k][j+1]))

  elif ascii_keywds[k]==name_dec_key:

    for j in range(ascii_table.nrows -1): dec_ascii.append(float(ascii_table[k][j+1]))

dist_asec = []

#Two arrays with the indexes of the matching objects

rowAscii_match = []	# ASCII rows

rowFits_match = []	# FITS rows

#Array with the indexes of the NEW objects found in the ASCII file (if any)

rowAscii_new = []

method_dict = {

  '1' : 'POSITION (dist < %.1f")' % match_radius,

  '2' : 'NAME',

  '3' : 'INDEX'

  }

print "\n\t>> Matching ASCII/FITS tables by %s ...\n" % method_dict[method]

num_tot_matches = 0

for j in range(Nrows_fits):

  num_multiple_matches = 0

  id_matches = []

  ra_dec_matches = []

  if match_option == '1':

    tmp_idxs_matches = []

    tmp_dist_matches = []

    for i in range(Nrows_ascii):

      dist_tmp = 3600. * astCoords.calcAngSepDeg(float(ra_fits[j]), float(dec_fits[j]), ra_ascii[i], dec_ascii[i])

      if dist_tmp <= match_radius:

	tmp_idxs_matches.append(i)

	tmp_dist_matches.append(round(dist_tmp,1))

	num_tot_matches += 1

	num_multiple_matches += 1

    idx_match = 0

    if len( tmp_idxs_matches ) > 1:

      print bcolors.WARNING+ "\n\t! WARNING ! %i objects found within %.1f arcsec from %s \n" % ( len(tmp_idxs_matches), match_radius, name_fits[j]) + bcolors.ENDC

      for idx in range( len(tmp_idxs_matches) ): print '\t%i: %s (dist = %s")' % ( (idx+1, name_ascii[ tmp_idxs_matches[idx]], tmp_dist_matches[idx] ) )

      tmp_check = False

      while tmp_check == False:

	tmp_entry = int(raw_input('\t-> Please enter the number of the matching object: '))

	if tmp_entry in range(1, len(tmp_idxs_matches)+1 ): 

	  tmp_check = True

	  idx_match = tmp_idxs_matches[ tmp_entry - 1 ]

	else:

	  print bcolors.FAIL+ "\n\t*** Wrong option ***\n"+ bcolors.ENDC

      id_matches.append((name_ascii[idx_match]).strip())

      ra_dec_matches.append(ra_ascii[idx_match])

      ra_dec_matches.append(dec_ascii[idx_match])

      # NOTE: When the ascii_table is called, the corresponding index is (rowAscii_match + 1) because of the additional line for the HEADER

      rowAscii_match.append(idx_match)

      rowFits_match.append(j)

    elif len( tmp_idxs_matches ) == 1:

      idx_match = tmp_idxs_matches[0]

      id_matches.append((name_ascii[idx_match]).strip())

      ra_dec_matches.append(ra_ascii[idx_match])

      ra_dec_matches.append(dec_ascii[idx_match])

      # NOTE: When the ascii_table is called, the corresponding index is (rowAscii_match + 1) because of the additional line for the HEADER

      rowAscii_match.append(idx_match)

      rowFits_match.append(j)

  elif match_option == '2':

    for i in range(Nrows_ascii):

      if (name_fits[j]).strip() == (name_ascii[i]).strip():

	num_multiple_matches += 1

	num_tot_matches += 1

	if num_multiple_matches > 1: 

	  print '%s Found %i objects with the same name : %s\nAborted.\n' % (error, num_multiple_matches, name_fits[j]); os._exit(0)

	# NOTE: When the ascii_table is called, the corresponding index is (rowAscii_match + 1) because of the additional line for the HEADER  

	rowAscii_match.append(i)

	rowFits_match.append(j)

  elif match_option == '3':

    for i in range(Nrows_ascii):

      if int(index_fits[j]) == int(index_ascii[i]) and (int(index_fits[j]) >= 0 and int(index_ascii[i]) >= 0):

	num_tot_matches += 1

	# NOTE: When the ascii_table is called, the corresponding index is (rowAscii_match + 1) because of the additional line for the HEADER

	rowAscii_match.append(i)

	rowFits_match.append(j)

	break

for i in range(Nrows_ascii):

  if i not in rowAscii_match: rowAscii_new.append(i) # Rows numbers of the NEW clusters, in the ASCII file

print "\n\t%s Found %s matching clusters between FITS/ASCII table to be UPDATED in the FITS table" % (info, len(rowAscii_match))

print "\n\t%s Found %s NEW clusters in the ASCII table to be ADDED to the FITS table" % (info, len(rowAscii_new))

#Store the names of the common/new clusters

idx_name = fits_keywds.index(name_Name_key)

clName_fits=[]

for k in range(Nrows_fits):

  clName_fits.append(fits_data[k][idx_name])

common_clNames=[]

new_clNames=[]

for i, idx in enumerate(rowAscii_match):

  common_clNames.append(clName_fits[rowFits_match[i]])

for idx in rowAscii_new:

  idx_name = ascii_keywds.index(name_Name_key)

  new_clNames.append( ascii_table[idx_name][idx+1] )

#Define the MASS conversion factor, only if it is found in ASCII table:

h_factor = 1.0

tmp_check = False

mass_in_ascii = set(name_mass_key) & set(ascii_keywds)

if mass_in_ascii:

  print "\n%s Concerning %s, do you want to:\n\t1) Convert from h70^-1 -> h100^-1\n\t2) Convert from h100^-1 -> h70^-1\n\t3) Keep the original values of the ASCII table" % (question, mass_in_ascii.pop())

  while tmp_check == False:

    message = "\n\t-> Please enter 1, 2 or 3:   "

    h_opt = raw_input(message)

    if h_opt == '1': h_factor = 0.7; tmp_check = True

    elif h_opt == '2': h_factor = 1./0.7; tmp_check = True

    elif h_opt == '3': h_factor = 1.; tmp_check = True

    else: print bcolors.FAIL+ "\n\t*** Wrong option ***"+ bcolors.ENDC

newRow_num = Nrows_fits + len(rowAscii_new)

'''

  *** 2nd data UPDATE: add the new clusters as new (initially empty) rows ***

'''

hdulist = pyfits.new_table(hdulist, nrows=newRow_num)

#Add 'CATALOG' to the new clusters (if any)

if name_catalog_key in fits_keywds and name_catalog_key not in ascii_keywds and len(rowAscii_new) > 0:

  new_catalog = raw_input("\n%s Please enter the value of %s for the new cluster(s): " % (question, name_catalog_key))

'''

  *** Update the PAPER column ***

'''

paper_flag = False

updated_paper_vec = []

max_length_paper = 0

cnt = 0

#If 'PAPER' is defined in the ASCII table, but it is not in the FITS and there are NO NEW objects, the latter is updated with the former

if name_paper_key in ascii_keywds and name_paper_key not in fits_keywds and len(rowAscii_new) == 0:

  for j in range(Nrows_fits):

    #Update only those clusters specified in the ASCII table

    if j in rowFits_match:

      paper_tmp = ascii_table[ascii_keywds.index(name_paper_key)][rowAscii_match[cnt]+1]

      cnt += 1

    else:

      paper_tmp = "Null"

    paper_tmp = remove_duplicated_names(paper_tmp)

    updated_paper_vec.append(paper_tmp)

    if len(paper_tmp) > max_length_paper: max_length_paper = len(paper_tmp)

#If 'PAPER' is defined in the FITS table, it is updated for the common clusters (and created for the new clusters) with the one defined in the ASCII table.

#If no 'PAPER' is found in ASCII, user is asked to enter it manually.

elif name_paper_key in fits_keywds:

  new_paper_vec = []

  col_paper_fits = fits_keywds.index(name_paper_key)

  if name_paper_key in ascii_keywds:

    paper_flag = True

    for i in range(Nrows_ascii):

      new_paper_vec.append( ascii_table[ascii_keywds.index(name_paper_key)][i+1].strip() )

  else:

    #The new reference is asked to be added manually only if new clusters are found

    if len(new_clNames)>0:

      tmp_new_paper = raw_input("\n%s Please insert the new reference to add: " % question)

      new_paper_vec=[tmp_new_paper for x in range( Nrows_ascii ) ]

      paper_flag = True

    else:

      new_paper_vec=['' for x in range( Nrows_ascii ) ]

  #Update those clusters in common with ASCII and FITS table

  for j in range(Nrows_fits):

    paper_old = (fits_data[j][col_paper_fits]).strip()

    if j in rowFits_match:

      if paper_old == "Null":

	paper_tmp = new_paper_vec[ rowAscii_match[cnt] ]	#Here the '+1' correction is not necessary because also new_paper_vec[] contains the header line

	cnt+=1

      else:

	paper_tmp = paper_old+"; "+new_paper_vec[ rowAscii_match[cnt] ]	#Here the '+1' correction is not necessary because also new_paper_vec[] contains the header line

	cnt += 1

    else:

      paper_tmp = paper_old

    paper_tmp = remove_duplicated_names(paper_tmp)

    updated_paper_vec.append(paper_tmp)

    if len(paper_tmp) > max_length_paper: max_length_paper = len(paper_tmp)

#Delete the old 'PAPER' column and update it with a new one defined according to the above case.

if name_paper_key in fits_keywds and paper_flag: 

  hdulist.columns.del_col(name_paper_key)

  #Add the new PAPER field, as last column

  col_tmp = pyfits.Column(name=name_paper_key, format=str(max_length_paper)+'A', unit = 'None', array=updated_paper_vec)

  paper_flag = True

if paper_flag:

  coldefs = pyfits.ColDefs(hdulist.columns)

  coldefs.add_col(col_tmp)

  hdulist = pyfits.new_table(coldefs)

#Update PAPER for common cluster

len_ALT_NAME = []

new_altName_vec = []

old_altName_vec = []

new_altName = ""

cnt = 0

altName_flag = False

name_in_altName = False

replace_altName = False

#Handle the NAME/ALT_NAME update in case of position/index matching:

if len(common_clNames) > 0:

  if name_altName_key not in ascii_keywds and name_altName_key in fits_keywds:

    if name_Name_key in fits_keywds and name_Name_key in ascii_keywds:

      answer_check = False

      tmp = raw_input("\n\t%s Do you want to add the old clusters' %s listed in FITS table to %s? [y/n]: " % (question, name_Name_key, name_altName_key) )

      while answer_check == False:

	if tmp in 'yesYES1' and tmp != '': 

	  name_in_altName = True

	  answer_check = True

	elif tmp in 'nN' and tmp != '': answer_check = True

	else: tmp = raw_input(bcolors.FAIL+ "\n\t\t*** Please enter a valid answer ***" + bcolors.ENDC + ' [y/n] : ')

    col_altName_fits = fits_keywds.index( name_altName_key )

    for j in range(Nrows_fits):