#!/usr/bin/python2
'''Provides functions to compile and use star catalogs, including routines
to query Vizier and SIMBAD.
'''

from math import *
import sys, os, re
import xml.sax
obsprepDir=os.getenv('OBSPREPDIR')
if obsprepDir == None:
    raise EnvironmentError('*** environment variable OBSPREPDIR not defined!')

import novas
import pdb

DEBUG_CATALOGS = False   # True  -> print out extra debug messages
                         # False -> do not print out extra debug messages


knownDiameters = {}


def getKnownDiameters():
    ''' Read the BSC file compiled by Christian Hummel.

    We will now use a diameter.bsc file from 2013-04-18 with HR5774
    flagged as "B" and some newer entries.

    This file was created via the "compile_diameters" command in oyster,
    and was taken from
    the oyster source directory on April 19, 2004.  It gives diameters,
    binarity, and suitability as calibrators (based on diameters and
    binarity) for all stars in the Bright Star Catalog with V < 5.5
    and dec > -30.  See the comments at the beginning of the file for a
    description of how the file was compiled.

    The data is returned as a dictionary in the global variable
    "knownDiameters", where the keys are the BSC (HR) number of the star,
    and the values are tuples containing 1) the diameter, in mas, and
    2) the calibrator status of the star, with values: "C" = calibrator
    (not binary and diameter less than 1 mas), "B" = binary/multiple,
    and "." = anything else.
    '''

    # No need to read in the file if we already have
    if len(knownDiameters) > 0: return

    # Here's the file to read
    diameterFile = obsprepDir+'/local/npoi/diameter.bsc' 
    if not os.path.exists(diameterFile):
        if DEBUG_CATALOGS:
            print 'looking for a local copy of diameter.bsc'
        diameterFile = obsprepDir+'/diameter.bsc'
        if not os.path.exists(diameterFile):
            raise ValueErr, "couldn't find diameter file '%s' in either '%s/local/npoi' or '.'" % (diameterFile, obsprepDir)

    # Read the file and create the dictionary
    expr = re.compile('(irfm|dm1|dm2|tn)')
    f = open(diameterFile)
    for line in f:
        if line[0] == '!': continue  # Skip comment lines
        tokens = line.split(None,7)
        source = 'ri'
        if len(tokens) > 7:
            result = expr.search(tokens[7])
            if result != None: source = result.group(0)
        knownDiameters[int(tokens[0])] = (float(tokens[1]), tokens[6], source)
    f.close()
        


def simbadQuery(starName):
    '''Run SIMBAD query and return dictionary with constellation, Hipparcos,
    FK5, HR, and common names.

    ARGS:
        starName (string):   name of star to query for

    RETURNS:
        (dictionary)   catalog name versus star name in that catalog
    '''

    import urllib

    # Query the SIMBAD web page with this star name
    simbadURL='http://simbad.u-strasbg.fr/simbad/sim-id'
    params='Ident='+starName.replace(' ','+')
    myQueryText=simbadURL+'?'+params
    print 'SIMBAD query text = ',myQueryText#@@@
    url = urllib.urlopen(myQueryText)

    # Fetch the constellation and Hipparcos names for this star from the HTML
    d = {}
    constellation = re.compile(r'>\*</A> ([a-zA-Z]+) ([a-zA-Z]+)')
    common= re.compile(r'>\*</A> ([a-zA-Z0-9]+) ([a-zA-Z]+)') #@@@ used to be 2nd constellation
    hipparcos = re.compile(r'>HIP (\d+)<')
    fk5 = re.compile(r'>FK5 (\d+)<')
    hr = re.compile(r'>HR (\d+)<')
    hd = re.compile(r'>HD (\d+)<')
    variablestar = re.compile('>V\*</A> ([a-zA-Z]+) ([a-zA-Z]+)')
    # for next re:
    # note that >NAME</A> ALGOL ABC matches excluding the ABC or
    # whatever else follows Algol as we have one group of (\w+)
    common = re.compile(r'>NAME</A> (\w+)')
    line = url.readline()
    while line:
        result = constellation.search(line)
        if result != None:
            d['constellation'] = \
           '%s %s' % (result.group(1), result.group(2))
        result = hipparcos.search(line)
        if result != None: d['hipparcos'] = int(result.group(1))
        result = fk5.search(line)
        if result != None: d['fk5'] = int(result.group(1))
        result = hr.search(line)
        if result != None: d['hr'] = int(result.group(1))
        result = hd.search(line)
        if result != None: d['hd'] = int(result.group(1))
        result = variablestar.search(line)
        if result != None:
            d['variablestar'] = \
            '%s %s' % (result.group(1), result.group(2))
        if d.has_key('common') == False:
            result = common.search(line)
            if result != None:
                d['common'] = result.group(1)
        line = url.readline()
    # Return the dictionary of names to whatever called simbadQuery()
    return d


class VizierStarTableHandler(xml.sax.ContentHandler):
    '''XML content handler to parse the VOTable/DTD output from Vizier
    and return the results as a list of "Star" objects.
    Handles queries to either Hipparcos, Tycho, or Tycho-2.
    '''

    # Mappings from Vizier field names to Star attributes.  First set of
    # keys map a Hipparcos query; keys starting from 'TYC' handle those fields
    # which change names in a Tycho query; keys starting from 'VTmag' handle
    # those field which change names in a Tycho-2 query.
    map = {'_RA':'raNow', '_DE':'decNow', 'HIP':'id',
           'Proxy':'proximity', 'Vmag':'v', 'RA(ICRS)':'ra',
           'DE(ICRS)':'dec', 'Plx':'parallax', 'pmRA':'pmRa',
           'pmDE':'pmDec', 'e_Plx':'parallaxErr', 'B-V':'bv',
           'V-I':'vi', 'HvarType':'variability', 'CCDM':'ccdm', 'rho':'sep',
           'dHp':'dmag', 'Notes':'note', 'SpType':'spectralType', 'HD':'hd',
           'TYC':'id', 'VarFlag':'variability', 'm_HIP':'ccdm',
           'Remark':'note', 'VTmag':'v','prox':'proximity'}

    def __init__(self):
        self.chars=[]
        self.tables=[]
        self.nFields = 0

    def characters(self, content):
        self.chars.append(content)

    def startElement(self, name, attrs):
        if name=="TD":
            self.chars=[]
        elif name=="TR":
            self.star = Star()
            self.nFields = 0
        elif name=="TABLE":
            self.stars=[]
            self.fields=[]
            self.types=[]
            self.nFields=0
        elif name=="FIELD":
            if re.search('_RA',str(attrs.getValueByQName('name'))) != None:
                self.fields.append('_RA')
            elif re.search('_DE',str(attrs.getValueByQName('name'))) != None:
                self.fields.append('_DE')
            else:
                self.fields.append(str(attrs.getValueByQName('name')))
            self.types.append(str(attrs.getValueByQName('datatype')))
            self.nFields += 1

    def endElement(self, name):
        if name=="TD":
            # The Vizier web page adds an extra undocumented column (recid)
            # which we must handle specially
            if self.nFields >= len(self.fields): return
            if self.fields[self.nFields] == 'recno':
                self.nFields += 1
                return

            # Cast the data value and assign to the appropriate attribute
            if ''.join(self.chars) == '':
                pass
                #self.star.__dict__[self.map[self.fields[self.nFields]]] = None
            else:
                if self.types[self.nFields] == 'char':
                    self.star.__dict__[self.map[self.fields[self.nFields]]] = str(''.join(self.chars))
                elif self.types[self.nFields] == 'int' or \
                     self.types[self.nFields] == 'short' or \
                     self.types[self.nFields] == 'unsignedByte':
                    self.star.__dict__[self.map[self.fields[self.nFields]]] = int(''.join(self.chars))
                elif self.types[self.nFields] == 'float' or \
                     self.types[self.nFields] == 'double':
                    self.star.__dict__[self.map[self.fields[self.nFields]]] = float(''.join(self.chars))
                else:
                    raise ValueError, 'unrecognized type "%s"' % \
                          self.types[self.nFields]
            self.nFields += 1
        elif name=="TR":
            self.stars.append(self.star)
        elif name=="TABLE":
            self.tables.append(self.stars)



class VizierDictTableHandler(xml.sax.ContentHandler):
    '''XML content handler to parse the VOTable/DTD output from Vizier
    and return the results as a list of dictionaries, one dictionary per star.
    '''

    def __init__(self):
        self.chars=[]
        self.tables=[]

    def characters(self, content):
        self.chars.append(content)

    def startElement(self, name, attrs):
        if name=="TD":
            self.chars=[]
        elif name=="TR":
            self.star = {}
            self.nFields = 0
        elif name=="TABLE":
            self.stars=[]
            self.fields=[]
            self.types=[]
            self.nFields=0
        elif name=="FIELD":
            self.fields.append(str(attrs.getValueByQName('name')))
            self.types.append(str(attrs.getValueByQName('datatype')))
            self.nFields += 1

    def endElement(self, name):
        if name=="TD":
            # The Vizier web page adds an extra undocumented column (recid)
            # which we must handle specially
            if self.nFields >= len(self.fields): return

            # Cast the data value and assign to the appropriate attribute
            if ''.join(self.chars) == '':
                self.star[self.fields[self.nFields]] = None
            else:
                if self.types[self.nFields] == 'char':
                    self.star[self.fields[self.nFields]] = str(''.join(self.chars))
                elif self.types[self.nFields] == 'int' or \
                     self.types[self.nFields] == 'short' or \
                     self.types[self.nFields] == 'unsignedByte':
                    self.star[self.fields[self.nFields]] = int(''.join(self.chars))
                elif self.types[self.nFields] == 'float' or \
                     self.types[self.nFields] == 'double':
                    self.star[self.fields[self.nFields]] = float(''.join(self.chars))
                else:
                    raise ValueError, 'unrecognized type "%s"' % \
                          self.types[self.nFields]
            self.nFields += 1
        elif name=="TR":
            self.stars.append(self.star)
        elif name=="TABLE":
            self.tables.append(self.stars)


def vizierQuery(starName, catalog=None, radius=-1, Vmax=1000, columns=None):
    '''Query Vizier for a set of stars.  Returns a list of dictionaries.

    Run Hipparcos or Tycho query and each catalog field as a dictionary.
    ARGS:
        starName (string):  Name of star - any name resolvable by Simbad is
                            legal
        catalog (string):   'h' = Hipparcos, 't1' = Tycho, 't2' = Tycho-2
                            'b' = BSC, 'f' = FK5
        radius (double):    find all stars within this radius of the star
                            (arcsec);  if <= 0, fetch just the single star
        Vmax (double)       faintest star, in V, to return
        columns (list)      columns to return (defaults to all)
    '''

    import urllib

    # Parse the catalog name and star number/name
    if catalog == None:
        # No catalog specified, so we're fetching it from the star name.
        # We can handle only the Hipparcos, Tycho, FK5, and HR catalogs.
        expr = re.compile("(HIP|TYC|FK5|HR|BSC)(\d+)")
        result = expr.match(starName.upper().replace(' ',''))
        if result == None: raise ValueError, "illegal star name '%s'"%starName
        catalog = result.group(1)
        if catalog == 'BSC': catalog = 'HR'
    else:
        # Catalog was specified.  Any star name OK, and hopefully will be
        # resolved by Vizier.
        catalog = catalog.upper()

    # Set the catalog in the query.
    if catalog == 'H' or catalog == 'HIP':
        cat = 'I/239/hip_main'
        vvar = 'VMag'
    elif catalog == 'F' or catalog == 'FK5':
        cat = 'I/149A/catalog'
        vvar = 'VMag'
    elif catalog == 'B' or catalog == 'HR':
        cat = 'V/50/catalog'
        vvar = 'VMag'
    elif catalog == 'T1' or catalog == 'TYC':
        cat = 'I/239/tyc_main'
        vvar = 'Vmag'
    elif catalog == 'T2':
        cat = 'I/259/tyc2'
        vvar = 'VTmag'
    else:
        raise ValueError, "catalog '%s' not supported" % catalog
    
    # Run the VizieR query
    vizierURL="http://vizier.u-strasbg.fr/viz-bin/votable/-dtd"
    params = '-c='+starName.replace(' ','+')+'&-source='+cat
    if columns == None:
        params += '&-out.add=.&-out.all'
    else:
        params += '&-out='+columns
    if radius > 0: params += '&-c.rs='+str(radius)+'&-sort=_r'
    if Vmax < 999: params += '&'+vvar+'=<'+str(Vmax)
    url = urllib.urlopen(vizierURL,params)

    # Parse the XML returned from Vizier, returning the results as a dictionary
    p = xml.sax.make_parser()
    handler = VizierDictTableHandler()
    p.setContentHandler(handler)
    line = url.readline()
    while line:
        p.feed(line)
        line = url.readline()

    # We expect to have read only one table
    if len(handler.tables) > 1:
        raise ValueError, "Vizier query didn't return exactly one table"

    # Return star if a single star search, or list of stars if a proximity
    # search
    if radius > 0:
        if len(handler.tables) == 0:
            return []
        else:
            return handler.tables[0]
    else:
        if len(handler.tables) != 1:
            raise ValueError, "Vizier query didn't return exactly one table"
        if len(handler.tables[0]) != 1:
            raise ValueError, "Vizier query didn't return exactly one star"
        return handler.tables[0][0]


class Star(object):
    '''A Hipparcos star.'''
    
    def __init__(self,starName=None):
        '''Create a new star.  If a star name is specified, then Vizier
        is queried to initialize the star.

        ARGS:
            starName (string):   name of the star
        '''

        self.constellation = None  # Constellation name of the star
        self.fk5 = None            # FK5 number of the star
        self.hr = None             # Bright star catalog number of the star
        self.hd = None             # HD number of the star
        self.common = None         # Common name of the star
        self.variablestar = None   # Variable star name of the star

        # The following attributes are read from the Hipparcos catalog.
        # The comments indicate which field in the Hipparcos catalog each
        # attribute corresponds to.
        self.id = None             # Catalog ID number of the star (H1)
        self.proximity = None      # Is there another object within 10 arcsec
                                   # (H2)?
                                   # ''=none, 'H'=Hipparcos, 'T'=Tycho
        self.v = None              # Johnson V magnitude (H5)
        self.ra = None             # ICRS J1991.25 right ascension (deg, H8)
        self.dec = None            # ICRS J1991.25 declination (deg, H9)
        self.parallax = None       # Trigonometric parallax (mas, H11)
        self.pmRa = None           # Proper motion in ra
                                   # (mas/yr, epoch J1991.25, H12)
        self.pmDec = None          # Proper motion in dec
                                   # (mas/yr, epoch J1991.25, H13)
        self.parallaxErr = None    # Error in the parallax (mas, H16)
        self.bv = None             # Johnson B-V color (H37)
        self.vi = None             # Johnson V-I color (H40)
        self.variability = None    # Variability flag (H52)
                                   #   'C' = not variable, anything
                                   #         else = variable
        self.ccdm = None           # Entry in the Catalogue of Components of
                                   # Double and Multiple Stars (H55):
                                   #    0=not a multiple system
        self.sep = None            # Angular separation for multiple systems
                                   # (arcsec, H64)
        self.dmag = None           # Magnitude difference of components in
                                   # multiple systems (H66)
        self.note = None           # Is there a note on this star? (H70)
                                   #   '' =no note
                                   #   'D'=double & multiple systems note only
                                   #   'G'=general note only
                                   #   'P'=photometric notes only
                                   #   'W'='D'+'P' only
                                   #   'X'='D'+'G' only
                                   #   'Y'='G'+'P' only
                                   #   'Z'='D'+'P'+'G'
        self.spectralType = None   # Spectral type (H76)
        self.raNow = None          # Right ascension at today's equinox and
                                   # epoch (deg)
        self.decNow = None         # Declination at today's equinox and
                                   # epoch (deg)

        # Here's a piece of nearby blank sky
        self.skyRa = None          # Right ascension at today's equinox and
                                   # epoch of a piece of blank sky
        self.skyDec = None         # Declination at today's equinox and
                                   # epoch of a piece of blank sky

        # Information from Oyster's 'diameter.bsc' file
        self.diameter = None       # Apparent diameter (mas)
        self.diameterSource = None # Source of apparent diameter
                                   #   'ri'   = R-I plus m_V calibration
                                   #   'irfm' = infrared flux method
                                   #   'dm1'  = Mark III, first pub
                                   #   'dm2'  = Mark III, second pub
                                   #   'tn'   = NPOI
        self.status = None         # 'C' = calibrator (d < 1 mas, not binary)
                                   # 'B' = binary/multiple
                                   # '.' = anything else

        # If no star name was specified, we're done.
        if starName == None: return

        # If star name was specified, query Vizier for this star in Hipparcos
        # and set fields accordingly.
        target = self.vizierQuery(starName,catalog='h')
        for var in target.__dict__:
            self.__dict__[var] = target.__dict__[var]

        # Query SIMBAD for its various names
        s = simbadQuery(starName)
        self.constellation = s.get('constellation')
        self.fk5 = s.get('fk5')
        self.hr = s.get('hr')
        self.common = s.get('common')

        # Fetch its diameter from Oyster's list of diameters
        getKnownDiameters()
        if knownDiameters.has_key(self.hr):
            tuple = knownDiameters[self.hr]
            self.diameter = tuple[0]
            self.status = tuple[1]
            self.diameterSource = tuple[2]

        # Find blank sky for it
        self.blankSky()

    def novas(self):
        '''Return the star as a NOVAS "cat_entry" structure.'''

        # Create an entry for the Hipparcos star
        hip = novas.Cat()
        hip.catalog = 'HIP'
        hip.starname = self.npoiName()
        hip.starnumber = self.id
        hip.ra = self.ra
        hip.dec = self.dec
        hip.promora = self.pmRa
        hip.promodec = self.pmDec
        hip.parallax = self.parallax
        hip.radialvelocity = 0

        # Must transform from Hipparcos units and epoch to NOVAS units and
        # epoch (FK5).
        fk5 = novas.Cat()
        novas.TransformHip(hip,fk5)
        return fk5

    def vizierQuery(starName, catalog=None, radius=-1, Vmax=1000,predicate=''):
        '''Query Vizier for the specified star, returning a single Star
        if querying only that star, or a list of Stars if doing a proximity
        query around the star.

        ARGS:
           starName (string):  Name of star - any name resolvable by Simbad
                               is legal
           catalog (string):   Catalog to query.
                                 "h"  = Hipparcos
                                 "t1" = Tycho,
                                 "t2" = Tycho-2
           radius (double):    Find all stars within this radius of the star
                               (arcsec); if <= 0 (default), fetch just the
                               single star
           Vmax (double):      Faintest star, in V, to return
           predicate (string): Search predicate if not searching around a star
        '''

        import urllib

        # Parse the catalog name
        if catalog == None:
            # No catalog specified, so we're fetching it from the star name.
            # We can handle only the Hipparcos and Tycho catalogs.
            if starName != '':
                expr = re.compile("(HIP|TYC)(\d+)")
                result = expr.match(starName.upper().replace(' ',''))
                if result == None:
                    raise ValueError, "illegal star name '%s'" % starName
                catalog = result.group(1)
        else:
            # Catalog was specified.  Any star name OK, and hopefully will be
            # resolved by Vizier.
            catalog = catalog.upper()

        # Here's the Vizier query web site
        vizierURL="http://vizier.u-strasbg.fr/viz-bin/votable/-dtd"

        # Set parameters based on which catalog
        params = '-out=*ID_MAIN,RA(ICRS),DE(ICRS),pmRA,pmDE,_RAtoday,_DEtoday'
        if catalog == 'H' or catalog == 'HIP':
            # Searching Hipparcos
            cat = 'I/239/hip_main'
            params += ',Proxy,Vmag,B-V,Plx,e_Plx,V-I,HvarType,CCDM,rho,dHp,Notes,SpType,HD'
            vvar = 'VMag'
        elif catalog == 'T1' or catalog == 'TYC':
            # Searching Tycho
            cat = 'I/239/tyc_main'
            params += ',Proxy,Vmag,B-V,Plx,e_Plx,VarFlag,m_HIP,Remark'
            vvar = 'Vmag'
        elif catalog == 'T2':
            # Searching Tycho-2
            cat = 'I/259/tyc2'
            params = '-out=HIP,RA(ICRS),DE(ICRS),CCDM,VTmag,pmRA,pmDE'
            vvar = 'VTmag'
        elif catalog == 'HR' or catalog == 'BSC':
            # Searching 5th bright star catalog
            cat = 'V/50/catalog'
            params = '-out=Vmag,B-V'
        else:
            raise ValueError, "illegal catalog '%s'" % catalog
        params += '&-source='+cat
        if starName != '': params += '&-c='+starName.replace(' ', '+')
        if radius > 0: params += '&-c.rs='+str(radius)+'&-sort=_r'
        if Vmax < 999: params += '&'+vvar+'=<'+str(Vmax)
        if predicate != '':
            params += '&%s&-out.max=9999' % predicate
        #@@@ print 'vizier query will be: ',vizierURL,params#@@@
        url = urllib.urlopen(vizierURL,params)
               
        # Parse the XML returned from Vizier.
        p = xml.sax.make_parser()
        handler = VizierStarTableHandler()
        p.setContentHandler(handler)
        line = url.readline()
        while line:
            p.feed(line)
            line = url.readline()

        # We expect to have read only one table
        if len(handler.tables) > 1:
            raise ValueError, "Vizier query didn't return exactly one table"

        # Return star if a single star search, or list of stars if
        # a proximity search
        if radius > 0 or starName == '':
            if len(handler.tables) == 0:
                return []
            else:
                return handler.tables[0]
        else:
            if len(handler.tables) != 1:
                raise ValueError,"Vizier query didn't return exactly one table"
            if len(handler.tables[0]) == 0:
                raise ValueError, "Vizier query didn't find the star"
            if len(handler.tables[0]) > 1:
                # For some reason returned more than one star.  See if one of
                # them is the correct one.
                nicknames = simbadQuery(starName)
                if not nicknames.has_key('hipparcos'):
                    raise ValueError, "Vizier query found more than one star"
                hipId = nicknames['hipparcos']
                for idx in xrange(0, len(handler.tables[0])):
                    if handler.tables[0][idx].id == hipId:
                        return handler.tables[0][idx]
                raise ValueError, "Vizier query found more than one star"
            return handler.tables[0][0]
    vizierQuery = staticmethod(vizierQuery)

    def show(self):
        '''Show its attributes in columnar format on stdout.'''

        for var in ('constellation', 'common', 'fk5', 'hr', 'hd', 'id', 'ra', 'dec', 'pmRa', 'pmDec', 'parallax', 'parallaxErr', 'raNow', 'decNow', 'v', 'bv', 'vi', 'proximity', 'variability', 'ccdm', 'sep', 'dmag', 'spectralType', 'note', 'skyRa', 'skyDec', 'diameter', 'diameterSource', 'status'):
            print '%15s:' % var, self.__dict__[var]


    def commentFileWrite(self, fd=None, comment=None):
        '''Append an entry for this star to the specified observing comment
        file.

        ARGS:
            fd (file descriptor):  file descriptor for the observing file to
                                   append to; if None, write to stdout
            comment (string):      additional comment for the file
        '''

        # Write to stdout if no file descriptor specified
        if fd == None: fd = sys.stdout

        # Print comments
        delimiter = '!-------------------------------------------------------------------------------\n'
        fd.write(delimiter)
        if comment != None: fd.write('! %s\n' % comment)
        fd.write('! HR %d, HIP %d' % (self.hr, self.id))
        if self.fk5 != None: fd.write(', FK5 %s' % self.fk5)
        if self.common != None: fd.write(', %s' % self.common)
        if self.constellation != None:
            fd.write(', %s' % self.constellation)
        fd.write('\n')
        fd.write('!    Hipparcos data:\n')
        fd.write('!       RAJ1991.5 = %9.5f deg, DECJ1991.5 = %8.5f deg\n' % \
                 (self.ra, self.dec))
        fd.write('!       muRA = %9.2f mas/yr, muDEC = %9.2f mas/yr\n' % \
                 (self.pmRa, self.pmDec))
        fd.write('!       parallax = %7.2f +- %5.2f mas\n' % \
                 (self.parallax, self.parallaxErr))
        fd.write('!       V = %5.2f, B-V = %5.2f, V-I = %5.2f, SpType = %s\n' \
                 % (self.v, self.bv, self.vi, self.spectralType))
        fd.write('!       ')
        if self.proximity.upper() == 'H':
            fd.write('Hipparcos')
        elif self.proximity.upper() == 'T':
            fd.write('Tycho')
        else:
            assert self.proximity == ' '
            fd.write('No Hipparcos or Tycho')
        fd.write(' stars within 10 arcsec\n')
        fd.write('!       Variable: ')
        if self.variability == 'C' or self.variability == ' ':
            fd.write(' No\n')
        else:
            fd.write(' Yes (%s)\n' % self.variability)
        fd.write('!       Multiple system: ')
        if self.ccdm == 0 or self.ccdm == None:
            fd.write(' No\n')
        else:
            fd.write(' Yes (CCDM = %s' % self.ccdm)
            if self.dmag != None: fd.write(', deltaMag = %.2f' % self.dmag)
            if self.sep != None:
                fd.write(', separation = %.2f arcsec' % self.sep)
            fd.write(')\n')
        fd.write('!       Notes: %s\n' % self.note)
        fd.write('!   Diameter source: ')
        if self.diameterSource == 'ri': fd.write('R-I\n')
        elif self.diameterSource == 'irfm': fd.write('IRFM\n')
        elif self.diameterSource == 'dm1': fd.write('MarkIIIa\n')
        elif self.diameterSource == 'dm2': fd.write('MarkIIIb\n')
        elif self.diameterSource == 'tn': fd.write('NPOI\n')
        else: raise AssertionError, 'illegal diameter source'
        fd.write(delimiter)

    def blankSky(self,inner=10,outer=30,avoidance=10):
        '''Find a blank area of sky around a star.

        The blank area is a point in the sky, located within an annulus around
        the star, that does not have a Tycho-2 star within a specified distance
        of it.

        ARGS:
            inner (double):     inner radius of the search annulus (arcsec)
            outer (double):     outer radius of the search annulus (arcsec)
            avoidance (double): minimum separation from any Tycho-2 star
                                (arcsec)
        '''

        # Verify that the ra and dec for this star are set
        assert self.raNow != None
        assert self.decNow != None

        # Simply search along the middle of the annulus in 30 degree steps and
        # query the Tycho-2 catalog till we find a blank spot
        radius = (inner + outer) / 2.
        nSteps = 12
        cosDec = cos(pi * self.dec / 180.)
        for i in range(nSteps):
            angle = 2 * pi * i / nSteps
            dra = radius * cos(angle) / cosDec
            ddec = radius * sin(angle)
            ra = self.ra + dra / 3600.
            dec = self.dec + ddec / 3600.
            stars = self.vizierQuery('%.6f%+.6f' % (ra,dec), 't2', \
                                     radius=avoidance)
            if len(stars) == 0: break
        else:
            raise ValueError, 'no blank sky found'

        # Found some blank sky.  Set it.
        self.skyRa = ra
        self.skyDec = dec

        #@@@ for a quick copy
        #@@@ print 'blank ra =: ', ra
        #@@@ print 'blank dec =: ', dec
    
    def separation(self,ra,dec):
        '''Return the separation in degrees of this star from the specified
        coordinates.  Use the current coordinates for this star.

        Taken from the SLALIB routine slaDsep().

        ARGS:
            ra (double):  right ascension (deg)
            dec (double): declination (deg)
        '''

        # Convert coordinates from spherical to Cartesian
        deg2rad = pi / 180.
        ra1 = self.raNow * deg2rad
        dec1 = self.decNow * deg2rad
        cosb = cos(dec1)
        v1 = [cos(ra1) * cosb, sin(ra1) * cosb, sin(dec1)]
        ra2 = ra * deg2rad
        dec2 = dec * deg2rad
        cosb = cos(dec2)
        v2 = [cos(ra2) * cosb, sin(ra2) * cosb, sin(dec2)]

        # Modulus squared of half the difference vector
        s2 = 0.0
        for i in range(3):
            d = v1[i] - v2[i]
            s2 += d * d
        s2 /= 4.0

        # Angle between the vectors
        c2 = 1.0 - s2
        if c2 < 0: c2 = 0.
        return (2.0 * atan2(sqrt(s2), sqrt(c2))) / deg2rad
        

    def findCalibrators(self,catalog):
        '''Find all possible calibrator stars for this star from the specified
        catalog.

        Returns a list of tuples, where the first entry is the calibrating
        star, and the second is its distance from the target star.

        ARGS:
            catalog (StarCatalog):    Catalog of potential calibrating stars
        '''

        # Find all calibrators within 20 degrees
        sep = [(star,star.separation(self.ra,self.dec)) \
               for star in catalog.targets if star.status == 'C']
        calibrators = [star for star in sep if star[1] < 20.]

        # Sort by increasing separation
        seps = [(calibrators[i][1], i) for i in range(len(calibrators))]
        seps.sort()
        return [calibrators[i] for junk,i in seps]

    def npoiName(self):
        '''Return the preferred NPOI name for the star.

        If the star has an FK5 entry, use that name [with FKV rather than FK5],
        else use the HR name [with BSC rather than HR]), else use the
        Hipparcos name.

        This will be called for each target star and each calibrator of that
        target.

        npoiName will also read in (once only) the cross index (cri) list
        Bob Zavala obtained from Jim Benson. This list contains FKV numbers
        which have HIC numbers. Some FK5 numbers do not have HIC numbers
        and the embedded system needs the BSC number in this case. 
        '''
        
        # setup dictionary for cross-index list

        criList = {}

        if len(criList) == 0: #only read file and fill dict once
            criFile = '%s/local/npoi/fkv.cri.bz2.ascii' % obsprepDir
            f = open(criFile)
            for line in f:
                tokens = line.split(None,1)
                criList[int(tokens[0])] = tokens[1]
            f.close()
        
        if self.fk5 != None:
            if self.fk5 in criList:
                return 'FKV%04d' % self.fk5
            elif self.hr != None:
                return 'BSC%04d' % self.hr
            else:
                return 'HIP%d' % self.id
        elif self.hr != None:
            return 'BSC%04d' % self.hr
        else:
            return 'HIP%d' % self.id
    


class StarCatalog(object):
    '''A catalog of Star objects.'''

    def __init__(self):
        '''Create a new, empty catalog.'''

        # Empty list of targets
        self.targets = []

        # Set up dictionaries of names versus star
        self.hip = {}
        self.hr = {}
        self.fk5 = {}
        self.hd = {}
        self.constellation = {}
        self.common = {}
        
    def build(Vmax=7.001,decMin=-30):
        '''Create a new catalog based around a Vizier query.

        The catalog will consist of all stars
        from the Hipparcos catalog brighter than the specified V and
        above the specified declination limit.

        ARGS:
            Vmax (double):   Faintest object (in V) to include in the catalog.
            decMin (double): Lower declination limit (deg) for the catalog.
        '''

        # Create the catalog
        catalog = StarCatalog()

        # Query Vizier for all Hipparcos stars meeting the specified V
        # and declination limits.
        print 'Querying Vizier with constraints:'
        Vmin=6.501 #@@@
        print 'Vmag=<%s & DE(ICRS)=>%s' % (Vmax,decMin)
        if Vmin > -0.1:
            print 'Vmin = %s' % (Vmin)
        catalog.targets = Star.vizierQuery('','h',
                                           predicate='Vmag=<%s&Vmag=>=%s&DE(ICRS)=>%s' \
                                           % (Vmax,Vmin,decMin))

        # Set dictionary of Hipparcos number versus star
        for star in catalog.targets: catalog.hip[star.id] = star

        # Make sure raNow and decNow are in decimal degree format
        catalog.checkNowCoords()

        # Find various other names for each star from SIMBAD
        print 'Fetching other names from SIMBAD'
        print 'This takes some time, recall that patience is a virtue.'
        catalog.otherNames()
        
        # Find diameters for each star from Oyster's list.
        print 'Updating diameters'
        catalog.diameters()

        # Find blank sky for each star
        print 'Updating blank sky' #@@@
        catalog.skys()             #@@@

        # Return the catalog
        return catalog
    build = staticmethod(build)

    def dumpDict(self,fileName):
        '''Dump the catalog as an ascii dictionary file.

        ARGS:
            fileName (string):  name of the file to dump the catalog to
        '''

        fd = open(fileName, 'w')
        for star in self.targets:
            print >>fd, star.__dict__
        fd.close()

    def loadDict(fileName):
        '''Load a catalog from an ascii dictionary file.

        ARGS:
            fileName (string):  name of the file to load the catalog from
        '''

        # Create the catalog
        catalog = StarCatalog()

        # Read in each star from the file
        fd = open(fileName, 'r')
        starCount = 0
        for line in fd:
            star = Star()
            star.__dict__ = eval(line)
            if star.id != None:
                catalog.hip[star.id] = star
                starCount+=1            
            for attr in ('hr', 'fk5', 'hd', 'constellation', 'common'):
                if star.__dict__[attr] != None:
                    catalog.__dict__[attr][star.__dict__[attr]] = star
            print 'Added star HIP',star.id,' to the catalog'
            catalog.targets.append(star)
        fd.close()
        print 'The catalog has '+str(starCount)+' stars.'
        return catalog
    loadDict = staticmethod(loadDict)

    def dump(self,fileName):
        '''Dump the catalog to a Pickle file.

        ARGS:
            fileName (string):  name of the file to dump the catalog to
        '''

        import cPickle

        fd = open(fileName, 'w.b')
        cPickle.dump(self,fd,True)
        fd.close()
        os.system('date')

    def load(fileName):
        '''Load a catalog from a Pickle file.

        INPUT:
            fileName (string):  name of the file to load the catalog from

        RETURNS:
            (StarCatalog) The catalog.
        '''

        import cPickle

        if DEBUG_CATALOGS:
           print 'Ready to open the pickled file: %s' % fileName

        fd = open(fileName, 'r.b')
        catalog = cPickle.load(fd)
        
        if DEBUG_CATALOGS: 
           print 'Sucessfully opened the pickled file'

        fd.close()
        return catalog
    load = staticmethod(load)


    def otherNames(self):
        '''Query SIMBAD using simbadQuery for the HR, FK5, constellation, and
        common names for each star in the catalog, and update those entries
        in the catalog.'''

        noOtherName = 0
        foundAName = 0
        print 'Looking for aliases for ',len(self.targets),' stars.'
        for star in self.targets:
            otherNamesDict = simbadQuery('HIP %d' % star.id)
            # Update the alias for each star. Update count
            # assuming we'll find other name(s). If we don't
            # we'll decrement below in the else: block
            #pdb.set_trace()
            if otherNamesDict.has_key('hr'):
                star.hr = otherNamesDict['hr']
                foundAName+=1
            if otherNamesDict.has_key('fk5'):
                star.fk5 = otherNamesDict['fk5']
                foundAName+=1
            if otherNamesDict.has_key('hd'):
                star.hd = otherNamesDict['hd']
                foundAName+=1
            if otherNamesDict.has_key('constellation'):
                star.constellation = otherNamesDict['constellation']
                foundAName+=1
            if otherNamesDict.has_key('variablestar'):
                star.variablestar = otherNamesDict['variablestar']
                foundAName+=1
            if otherNamesDict.has_key('common'):
                star.common = otherNamesDict['common']
                foundAName+=1
            if foundAName == 0:
                print 'No otherName found for HIP ',star.id
                noOtherName+=1
        print 'Found names for', len(self.targets)-noOtherName, 'of', len(self.targets),' stars.'            
 
        # Create dictionary versus HD names
        for star in self.targets:
            if star.hd != None: self.hd[star.hd] = star

         
    def skys(self):
        '''Set the blank sky for all stars in the catalog.'''
        idx = 0
        for star in self.targets:
            try:
                star.blankSky()
            except:
                print 'Trying blankSky() for', idx, ' HIP ',star.id,' again'
                star.blankSky()
            idx += 1
            if idx % 25 == 0: print idx, 'of', len(self.targets), 'finished'

    def diameters(self):
        '''Fetch diameters from Oysters list of diameters for all stars.'''

        getKnownDiameters()
        ndone = 0
        for star in self.targets:
            if knownDiameters.has_key(star.hr):
                ndone += 1
                tuple = knownDiameters[star.hr]
                star.diameter = tuple[0]
                star.status = tuple[1]
                star.diameterSource = tuple[2]
        print 'Found diameters for', ndone, 'of', len(self.targets)

    def checkNowCoords(self):
        '''Check if raNow or decNow are in sexagesimal format.
        If they are convert them to decimal degrees.
        I do watch for a negative zero degrees declination.
        '''
        for star in self.targets:
            if type(star.raNow) == type('a'):
                raParts=star.raNow.split()
                raDecimal=float(raParts[0])+(float(raParts[1])+float(raParts[2])/60.0)/60.0
                star.raNow=15.0*raDecimal
            if type(star.decNow) == type('a'):
                decParts=star.decNow.split()
                decimalPart=(float(decParts[1])+float(decParts[2])/60.0)/60.0
                if decParts[0][0] != '-':
                    decDecimal=float(decParts[0])+decimalPart
                else:
                    decDecimal=-1.0*(float(decParts[0][1:len(decParts[0])])+decimalPart)
                star.decNow=decDecimal                    
                
    def get(self,starName):
        '''Return the catalog star with the specified name.

        The name may
        be any name resolvable by SIMBAD, however Hipparcos, HR, HD, FK5,
        constellation, and some common names do not require a SIMBAD query.

        ARGS:
            starName (string):    name of the star to fetch

        RETURNS:
            (StarCatalog) The star, or None if the star is not in the catalog
        '''

        # First search the built in dictionaries
        starName = starName.upper()
        reg = re.compile(r'^(HIP|FK5|FKV|HR|BSC|HD) *(\d+)$')
        result = reg.match(starName)
 
        if DEBUG_CATALOGS:
            print 'starName is :',starName
            print 'starName in catalogs.py is :',starName #@@@
            
        if result != None:
            id = int(result.group(2).lstrip('0'))
            
            if DEBUG_CATALOGS:
                print 'id = %d' % id
            
            if result.group(1) == 'HIP':
                if DEBUG_CATALOGS:
                   print 'Searching catalog with hip number' 
                star = self.hip.get(id)
            elif result.group(1) == 'FK5' or result.group(1) == 'FKV':
                if DEBUG_CATALOGS:
                    print 'Searching catalog with fk5 number'
                star = self.fk5.get(id)
            elif result.group(1) == 'HR' or result.group(1) == 'BSC':
                star = self.hr.get(id)
            elif result.group(1) == 'HD':
                star = self.hd.get(id)
        else:
            star = self.constellation.get(starName)
            if star == None: star = self.common.get(starName)
        
        if star != None: return star
       
        # Not found.  Query SIMBAD for the corresponding Hipparcos name
        # and look up in the catalog using the Hipparcos name.
        print 'Will use simbad in StarCatalog.get() as star '+starName+' not recognized.'
        s = simbadQuery(starName)
        if not s.has_key('hipparcos'): 
            return None                
        if self.hip.has_key(s['hipparcos']):
            print 'Simbad query worked, found star in the catalog'
            return self.hip.get(s['hipparcos'])

        print 'Star not in catalog, perform Vizier query.'
        # Still not found, therefore not in our standard catalog.  Query
        # Vizier for the star and add it to the catalog (for this session
        # only). Note that Vizier query is contained within class Star
        star = Star(starName)
        if star == None: return None
        self.targets.append(star)
        if not star.id == None: self.hip[star.id] = star
        if not star.hr == None: self.hr[star.hr] = star
        if not star.fk5 == None: self.fk5[star.fk5] = star
        if not star.common == None: self.common[star.common] = star
        if not star.constellation == None:
            self.constellation[star.constellation] = star
        # Check for and fix if needed the raNow and decNow formats
        # pdb.set_trace()
        if type(star.raNow) == type('a'):
            raNowParts=star.raNow.split()
            raHours=float(raNowParts[0])+(float(raNowParts[1])+float(raNowParts[2])/60.0)/60.0
            star.raNow=raHours*15.0
        if type(star.decNow) == type('a'):
            decNowParts=star.decNow.split()
            decDegrees=float(decNowParts[0])+(float(decNowParts[1])+float(decNowParts[2])/60.0)/60.0
            star.decNow=decDegrees
        return star