#!/usr/bin/python
#
# Example code to query Simbad and Vizier catalogs with Python
# Requires AstroPy and AstroQuery to be installed (put "install astropy" and
# "install astrowuery" in google for instructions).
# Note: This code is written for Python 2.7. At least the print statements will not
# work in Python 3 ...

# Import the libraries for Simbad and Vizier queries:
from astroquery.simbad import Simbad
from astroquery.vizier import Vizier
from astropy import units as u # Units is a handy tool to work with actual units ...
import numpy as np

# First, let's say we have a couple of targets with random IDs, but need HD numbers (the last one is a 'typo').
# Make a list of strings containing the target IDs and an empty list to be filled with the HD numbers:
random_IDs = ['GJ 9049', 'SAO 100038', 'HIP17439', 'HIPz']
HD_numbers = []

# We query Simbad for the IDs we have:
for ID_old in random_IDs:
  IDs = Simbad.query_objectids(ID_old)
  
  # What if some target is not found?
  if IDs == None:
    HD_numbers.append('NOT FOUND') # Write into the list of results something that says your query didn't give any results.
    print '' # just for a blank line ...
    print 'Warning: ' + ID_old + ' not found in SIMBAD!' # Also, be nice and print a warning here.
    print ''
    break # Break the loop here, since in this case there is nothing else to be done
  
  # For the ones found, we get a list of all IDs:
  print 'Identifiers found for ' + ID_old + ':'
  print IDs
  print ''
   
  # Now search for the HD entry and append it to the list of new IDs:
  for ID_new in IDs:
    if ID_new[0][0:2] == 'HD': # If the first two characters in the ID are 'HD'
      HD_numbers.append(ID_new[0].replace(" ", "")) # the .replace(" ", "") removes all spaces in the string

# Now a break in the presentation, wait for user input:
raw_input("Press Enter to continue...")
print ''

# Print the results:
for i in range(0, len(HD_numbers)):
  print random_IDs[i] + ' = ' + HD_numbers[i]
print ''

raw_input("Press Enter to continue...")
print ''

# Now, for simplicity, we just continue with one target, HD7788.
# We want to observe it with the VLTI (PIONIER) and need calibrators. We find those
# in the catalogue of Merand et al. (2006). They are supposed to be similar
# to HD7788 in H magnitude and position on sky. Additional relevant parameters
# are their spectral types, diameters, and uncertainties on the diameters.

# First, lets get the H magnitude of HD7788 from SIMBAD (not the most accurate
# value, probably, but good enough):

Simbad.add_votable_fields('flux(H)') # Set the output to return the flux, error, and unit in K band
# note: A list of all possible VOTable fileds can be found here: http://simbad.u-strasbg.fr/simbad/sim-help?Page=sim-fscript#VotableFields
HD7788_props = Simbad.query_object("HD7788") # Do the query

print 'Columns retrieved:'
print HD7788_props.colnames # Print names of all columns retrieved
print ''

print 'Brightness in H: ', HD7788_props['FLUX_H'][0], HD7788_props['FLUX_H'].unit # You can access the columns of the output by their name. '.unit' after a column gives you the unit.
print ''

raw_input("Press Enter to continue...")
print ''

# Now query the Merand et al. catalog for calibrators

calibrators = Vizier.query_object(object_name='HD7788', catalog='J/A+A/433/1155', radius=20.0 * u.deg)

print ''
print 'Note: I have no idea what this exception means. Seems not to make any problems ...'
print ''

raw_input("Press Enter to continue...")
print ''

print 'Columns retrieved:'
print calibrators[0].colnames
print ''

raw_input("Press Enter to continue...")
print ''

print 'calibrators found:'
for i in range(0, len(calibrators[0]['HD'])):
  if (calibrators[0]['Hmag'][i] > HD7788_props['FLUX_H'][0] - 2.0) and (calibrators[0]['Hmag'][i] < HD7788_props['FLUX_H'][0] + 2.0): # If the H mag is within a range of +/- 2 mag from HD7788
    print 'HD', calibrators[0]['HD'][i], '; sep = ', calibrators[0]['_r'][i], '; SpT = ', calibrators[0]['SpType'][i], '; Vmag = ', calibrators[0]['Vmag'][i], '; Hmag = ', calibrators[0]['Hmag'][i], '; diam = ', calibrators[0]['UDdiamH'][i], '+/-', calibrators[0]['e_UDdiam'][i]