# -*- coding: utf-8 -*-
"""classical astronomy functions dealing with times, dates, and coordinates.

This package is based on the packages
ManGord,    the radio astronomy package by Manchester and Gordon,
coco,       coordinate conversions by P.T. Wallace of the Starlink project,
formats,    utilities to convert coordinates to and from strings,
jplephem,   JPL ephemeris support by Rick Fisher of NRAO.

This package contains classical astronomy functions having to do with
times, dates, and coordinates. It is based on the packages ManGord, the 
radio astronomy package by Manchester and Gordon, coco, coordinate 
conversions by P.T. Wallace of the Starlink project, formats, utilities 
to convert coordinates to and from strings, jplephem, JPL ephemeris 
support by Rick Fisher of NRAO.

It also defines AU and pc in meters.

While it is possible to use functions from these modules directly, the 
Astronomy package is designed to wrap the key functions with a 
consistent input/output format. Angles are input in decimal natural 
units, i.e., hours and degrees. Single valued outputs are variables, 
arrays are lists, multiple variables are tuples. The Astronomy module 
contains functions designed to achieve this consistency as well as some 
functions not provided elsewhere.

Details are in the doc strings for the modules and functions.

If module jplephem is to be used:
 # update the IERS table with 'refresh_ierstab()'
 # load the DSN station coordinates with something like
   'station_coordinates = get_cartesian_coordinates()'
 # set the ephemeris directory if the default is inappropriate; the default is
   /usr/local/lib/python[python version]/site-packages/Astronomy/jplephem/.
"""
import urllib
import datetime

import coco
try:
  import jplephem
  have_jplephem = True
except ImportError:
  print "Module jplephem is not available."
  have_jplephem = False
import formats
import ManGord
import mysql

import math
import glob
import os
import os.path
import time
import sys
import Physics as P
import Scientific.Physics.PhysicalQuantities as PQ
pq = PQ.PhysicalQuantity

diag = False

c = P.c # in m/s
PQ._addUnit('pc', 'm*3.0856e16')
PQ._addUnit('AU', 'm*1.495985e11')
pc = pq(1,'pc').inBaseUnits().value # = 3.0856e16 m
AU = pq(1,'AU').inBaseUnits().value # = 1.495985e11 m

Months = ['January','February','March','April','May','June',
          'July','August','September','October','November','December'] 

default_astro_dir = "/usr/local/pkg/astro/lib/astro"
cat_path          = default_astro_dir+"/sources/"

#eph_url = "http://data.iers.org/products/177/11221/orig/eopc04_05_IAU2000.62-now"
eph_url = "http://data.iers.org/products/214/14443/orig/eopc04_08_IAU2000.62-now"
pyver=sys.version[0:3]

if have_jplephem:
  jplephem.set_ephemeris_dir(
  '/usr/local/lib/python' + pyver +
  '/site-packages/Astronomy/jplephem/',
  'DEc403')

def short_month(mon):
  """
  Short month name

  @param mon : int
    January = 1, ..., December = 12

  @return: 4-character string
    Ends with period
  """
  if len(Months[mon-1]) > 4:
    return Months[mon-1][0:3]+'.'
  else:
    return Months[mon-1]
  
def long_month(mon):
  """
  Full month name

  @param mon : int
    January = 1, ..., December = 12

  @return: string
    Full month name
  """
  return Months[mon-1]

def parse_iers_line(line):
  mjd = int(line[12:19])
  x = float(line[19:30])
  y = float(line[30:41])
  UT1_UTC = float(line[41:52])
  return mjd,x,y,UT1_UTC
  
def refresh_ierstab(url=eph_url,force=False):
  """
  Updates the UT1-UTC offset table from 'url'.

  Notes
  =====
  Should be invoked if working with current or future dates.
  Only 800 days are stored

  @param url : strings
    The recommended URL is
    http://data.iers.org/products/177/11221/orig/eopc04_05_IAU2000.62-now

  @return: True if iers.tab was updated
  """
  delta_t_file = "/usr/local/lib/python" + pyver + \
                 "/site-packages/Astronomy/jplephem/iers.tab"
  
  # Is it needed? Is the file more than a day old?
  if time.time() - os.path.getmtime(delta_t_file) > 24*60*60 or force == True:
    f = urllib.urlopen(url)
    iers_data = f.readlines()
    f.close()
  
    t = open(delta_t_file,"w")

    for line in iers_data[-801:-1]:
      mjd,x,y,UT1_UTC = parse_iers_line(line)
      t.write("%6d %7.4f %7.4f %8.5f\n" % (mjd,x,y,UT1_UTC))
    t.close()
    return True
  else:
    return False

def get_cartesian_coordinates(station):
  """
  Create Cartesian coordinate dictionary

  This creates a dictionary with the Cartesian coordinates in meters
  of the DSN stations using the ITRF1993 (assuming subreflector-fixed
  configuration).
    
  Notes
  =====
  For an explanation of the coordinate systems see
  http://dsnra.jpl.nasa.gov/Antennas/Antennas.html#anchor950381
  Some other stations have been added.
  An example to get a baseline length:
  
  >>> x1,y1,z1 = get_cartesian_coordinates('DSS 24')
  >>> x2,y2,z2 = get_cartesian_coordinates('DSS 13')
  >>> print math.sqrt(math.pow(x2-x1,2) + math.pow(y2-y1,2)+math.pow(z2-z1,2))
  12621.4825356

  @param station : string
    For example 'Venus' or 'DSS 13'.

  @return: tuple
    Cartesian coordinates in meters
  """

  coordinates = \
    {"GB OVLBI": ( 884084.2636, -4924578.7481, 3943734.3354),
     "DSS 12": (-2350443.812, -4651980.837, +3665630.988), \
     "Echo":   (-2350443.812, -4651980.837, +3665630.988), \
     "DSS 13": (-2351112.491, -4655530.714, +3660912.787), \
     "Venus":  (-2351112.491, -4655530.714, +3660912.787), \
     "DSS 14": (-2353621.251, -4641341.542, +3677052.370), \
     "Mars":   (-2353621.251, -4641341.542, +3677052.370), \
     "DSS 15": (-2353538.790, -4641649.507, +3676670.043), \
     "DSS 16": (-2354763.158, -4646787.462, +3669387.069), \
     "DSS 17": (-2354730.357, -4646751.776, +3669440.659), \
     "DSS 23": (-2354757.567, -4646934.675, +3669207.824), \
     "DSS 24": (-2354906.495, -4646840.128, +3669242.317), \
     "DSS 25": (-2355022.066, -4646953.636, +3669040.895), \
     "DSS 26": (-2354890.967,-4647166.925,+3668872.212), \
     "DSS 27": (-2349915.260,-4656756.484,+3660096.529), \
     "DSS 28": (-2350101.849,-4656673.447,+3660103.577), \
     "DSS 33": (-4461083.514,+2682281.745,-3674570.392), \
     "DSS 34": (-4461146.756,+2682439.293,-3674393.542), \
     "DSS 42": (-4460981.016,+2682413.525,-3674582.072), \
     "DSS 43": (-4460894.585,+2682361.554,-3674748.580), \
     "DSS 45": (-4460935.250,+2682765.710,-3674381.402), \
     "DSS 46": (-4460828.619,+2682129.556,-3674975.508), \
     "Parkes": (-4554231.843,+2816758.983,-3454036.065), \
     "DSS 53": (+4849330.129,-0360338.092,+4114758.766), \
     "DSS 54": (+4849434.555,-0360724.108,+4114618.643), \
     "DSS 55": (+4849525.318,-0360606.299,+4114494.905), \
     "DSS 61": (+4849245.211,-0360278.166,+4114884.445), \
     "DSS 63": (+4849092.647,-0360180.569,+4115109.113), \
     "DSS 65": (+4849336.730,-0360488.859,+4114748.775), \
     "DSS 66": (+4849148.543,-0360474.842,+4114995.021) }
  try:
    return coordinates[station]
  except:
    return None

def B_epoch_to_J(ra50,dec50,format):
  """
  Convert B1950 coordinates to J2000
  
  This is a convenience to avoid having to format angles before calling
  the coordinate conversion routine.  The need often comes up when an
  operator asks for coordinates in J2000 or decimal format.
    
  Notes
  =====
  If format = 'decimal', then:

  >>> from Astronomy import *
  >>> B_epoch_to_J('''00h02m29.056400s''', '''54d11'43.187000"''', 'decimal')

  ('0.0845454', '54.4735905')

  If format = 'formatted', then it returns:
  ('00h05m04.364s', '+54d28m24.926"')
    
  Otherwise it returns
  ([0, 5, '4.364'], [54, 28, '24.926'])

  Tests and validations
  =====================

  For the above example, compare to the VLA
  Calibrator List:
  J2000 A 00h05m04.363531s 54d28'24.926230"
  B1950 A 00h02m29.056400s 54d11'43.187000"
  
  Crossing the midnight boundary
  ------------------------------
  
  >>> B_epoch_to_J('''23h58m34.865400s''', '''18d57'51.753000"''', 'none')
  ([0, 1, '8.622'], [19, 14, '33.802'])
  
  Compare to the VLA Calibrator List:
  J2000 A 00h01m08.621563s 19d14'33.801860"
  B1950 A 23h58m34.865400s 18d57'51.753000"
  
  Negative declination
  --------------------
  
  >>> B_epochs_to_J('''00h00m48.4200s''', '''-17d43'54.000"''', 'formatted')
  ('00h03m21.997s', '-17d27m11.782"')

  @param ra50 : string
    For example: '00h02m29.056400s'

  @param dec50: string
    For example : '''54d11'43.187000"'''

  @param format : string
    'decimal', 'formatted', other

  @return: tuple of strings
    See notes for details.
  """
  if ( ra50.find('h') > -1 ):
    ra1950 = formats.hms_delimited_angle_to_rads(ra50)
    dec1950 = formats.dms_delimited_angle_to_rads(dec50)
  elif ( ra50.find(':') > -1 ):
    p1950 = formats.parse_colon_delimited_angles(ra50,dec50)
    ra1950 = p1950[0]*math.pi/12
    dec1950 = p1950[1]*math.pi/180
  # The following function wants angles in radians
  p2000   = coco.B1950_to_J2000(ra1950,dec1950,1980.0)
  ra2000  = coco.range_0_2pi(p2000[0])*12/math.pi
  dec2000 = coco.range_p_m_pi(p2000[1])*180/math.pi
  if ( len(format) == 0 ):
    format = 'None'
  if ( format[0] == "d" ):
    return (ra2000, dec2000 )
  else:
    ra_hms = coco.rads_to_HMS(3,ra2000*math.pi/12)
    ra_sign = ra_hms[0]
    ra = ra_hms[1]
    dec_dms = coco.rads_to_DMS(3,dec2000*math.pi/180)
    dec_sign = dec_dms[0]
    dec = dec_dms[1]
    if ( format[0] == "f"):
      ra_string = "%02dh%02dm%06.3fs" % \
                  (ra_sign*ra[0], ra[1], \
                  (float(ra[2]) +float(ra[3])/pow(10,3)))
      dec_string = """%+02dd%02d'%06.3f\"""" % \
                  (dec_sign*dec[0], dec[1], \
                  (float(dec[2]) +float(dec[3])/pow(10,3)))
      return (ra_string,dec_string)
    else:
      ra_list = [ra_sign*ra[0], ra[1], \
                  float(("%.3f" % (float(ra[2]) +float(ra[3])/pow(10,3))))]
      dec_list = [dec_sign*dec[0],dec[1],\
                  float(("%.3f" % (float(dec[2])+float(dec[3])/pow(10,3))))]
      return (ra_list,dec_list)

def J_epoch_to_B(ra2000,dec2000,format):
  """
  Convert J2000 coordinates to B1950.
  
  See B_epoch_to_J for documentation

  Notes
  =====
  A test:
  
  >>> from Astronomy import *
  >>> B_epoch_to_J('''00h02m29.056400s''', '''54d11'43.187000"''', 'formatted')
  ('00h05m04.364s', '+54d28\'24.926"')
  >>> J_epoch_to_B('''00h05m04.364s''', '''+54d28\'24.926"''', 'formatted')
  ('00h02m29.068s', '+54d11m42.840"')
  
  Seems to be good to about an arcsecond.

  @param ra2000 : string

  @param dec2000 : string

  @param format : string

  @return: tuple of strings
  """
  if ( ra2000.find('h') > -1 ):
    ra2k = formats.hms_delimited_angle_to_rads(ra2000)
    dec2k = formats.dms_delimited_angle_to_rads(dec2000)
  elif ( ra2000.find(':') > -1 ):
    p2000 = formats.parse_colon_delimited_angles(ra2000,dec2000)
    ra2k = p2000[0]*math.pi/12
    dec2k = p2000[1]*math.pi/180
  # The following function wants angles in radians
  p1950 = coco.FK5_to_FK4(ra2k, dec2k)
  #print p1950
  ra1950 = (p1950[0][0]-p1950[2][0]*50)*12/math.pi
  dec1950 = (p1950[1][0]-p1950[3][0]*50)*180/math.pi
  if ( len(format) == 0 ):
    format = 'None'
  if ( format[0] == "d" ):
    return ("%.7f" % ra1950, "%.7f" % dec1950 )
  else:
    ra_hms = coco.rads_to_HMS(3,ra1950*math.pi/12)
    ra_sign = ra_hms[0]
    ra = ra_hms[1]
    dec_dms = coco.rads_to_DMS(3,dec1950*math.pi/180)
    dec_sign = dec_dms[0]
    dec = dec_dms[1]
    if ( format[0] == "f"):
      ra_string = "%02dh%02dm%06.3fs" % \
                  (ra_sign*ra[0], ra[1], \
                  (float(ra[2]) +float(ra[3])/pow(10,3)))
      dec_string = "%+02dd%02dm%06.3f\"" % \
                  (dec_sign*dec[0], dec[1], \
                  (float(dec[2]) +float(dec[3])/pow(10,3)))
      return (ra_string,dec_string)
    else:
      ra_list = [ra_sign*ra[0], ra[1], \
                ("%.3f" % (float(ra[2]) +float(ra[3])/pow(10,3)))]
      dec_list = [dec_sign*dec[0],dec[1],\
                ("%.3f" % (float(dec[2])+float(dec[3])/pow(10,3)))]
      return (ra_list,dec_list)

def MJD(*args):
  """
  Modified Julian date

  This can take input in various forms:
  2 or 4 digit year, followed by day-of-year or month and day number.

  Examples:
  In [2]: MJD(11,02,02)
  Out[2]: 55594
  In [3]: MJD(2011,02,02)
  Out[3]: 55594
  In [4]: MJD(2011,33)
  Out[4]: 55594
  In [5]: MJD(11,33)
  Out[5]: 55594

  @param args : [two or four digit year, month or day-of-year, optional day-of-month]
  @type  args : list of ints

  @return: int
  """
  if len(args) == 2:
    year, month, day = calendar_date(args[0],args[1])
  else:
    year, month, day = args
  result = coco.cal_to_MJD(year,month,day)
  if ( result[1] == 0 ):
    return int(result[0])
  else:
    return -result[1]

def YMD_to_MJD(year,month,day):
  """
  Modified Julian date from calendar date

  Notes
  =====
  This calls coco.full_ymd_to_mjd but output format is simpler.

  @param year : four digit year
  @type  year : int

  @param month : int

  @param day : int

  @return: long
  """
  return long(coco.full_YMD_to_MJD(year,month,day)[0])

def julian_date (year, doy):
  """
  Julian date

  Notes
  =====
  On the Julian calendar; Julian Days are used by astromoners to
  avoid the complexities of the various calendars.

  @param year : int

  @param doy : int
    Day of year
  
  @return: float
    Julian Day (J.D) = number of days since noon on Jan. 1, 4713 BC
  """
  prev_year = year - 1
  century = prev_year / 100
  num_leaps = int(prev_year / 4) - century + int(century / 4)
  jd = 1721425. + 365. * prev_year + num_leaps - 0.5 + doy
  return (jd)

def day_of_year (year, month, day):
  """
  Day of year

  @param year : int

  @param month : int

  @param day : int

  @return: int
    the day of the year where Jan. 1 is DOY 1
  """
  day = day + (month -1) * 30. + (int)((month + 1) * 0.61) - 2
  if (month <= 2):
    day = day + month
  else:
    if (leap_year(year) != 1):
      day = day - 1
  return (int(day))

def leap_year (year):
  """
  Leap year

  @param year : int

  @return: int
    1 if year is leap year, otherwise 0
  """
  if (year % 100 == 0 ): # Gregorian fix
    if (year % 400 == 0 ):
      return (1)
    else:
      return (0)
  else:
    if (year % 4 == 0 ):
      return (1)
    else:
      return (0)

def day_of_week (doy, year):
  """
  Numeric value for the day of week

  @param doy : int
    Day of year

  @param year : int

  @return: int
  1 - sunday,
  2 - monday,
  3 - tuesday,
  4 - wednesday,
  5 - thursday,
  6 - friday,
  7 - saturday,
  """
  day = julian_date( year, (doy + 0.5) ) + 2
  return int((day - 7 * (int(day - 1) / 7)))

def v_sun(mjd,ra_source,dec_source):
  """
  LSR velocity of the sun

  Notes
  =====
  A routine for the vsun in ManGord.
  
  The velocity of the Sun with respect to the Local Standard of Rest is
  the conventional value of 20.0 km/s toward
  RA, dec = 18h03m50.29s, +30d00'16.8 (J2000)

  @param mjd : int
    mean Julian day

  @param ra_source : float
    Right ascension in decimal hours

  @param dec_source : float
    Declination in decimal degrees

  @return: float
    the projection onto the line of sight of the velocity of the
    Sun with respect to the Local Standard of Rest on the date given by
    mjd at UT = 0h.
  """
  # precession routine takes decimal angles
  ra_SM = formats.hms_delimited_angle_to_rads("""18h03m50.29s""")
  dec_SM = formats.dms_delimited_angle_to_rads('''+30d00'16.8"''')
  ra_SM = ra_SM*12/math.pi
  dec_SM = dec_SM*180/math.pi
  # Precess solar motion from J2000 to epoch
  SM_epoch = jplephem.j2000_to_epoch(mjd, 0.0, ra_SM, dec_SM)
  ra_epoch = SM_epoch['ra']*math.pi/12
  dec_epoch = SM_epoch['dec']*math.pi/180
  # Cartesian coordinates
  v_x = 20.0 * math.cos(ra_epoch) * math.cos(dec_epoch)
  v_y = 20.0 * math.sin(ra_epoch) * math.cos(dec_epoch)
  v_z = 20.0 * math.sin(dec_epoch)
  # Source direction in Cartesian coordinates
  ra = ra_source*math.pi/12
  dec = dec_source*math.pi/180
  x = math.cos(dec) * math.cos(ra)
  y = math.cos(dec) * math.sin(ra)
  z = math.sin(dec)
  # Dot product
  v_sun = -v_x * x - v_y * y - v_z * z
  return v_sun

def HaDec_to_AzEl(HourAngle,Declination,Latitude):
  """
  Celestial to horizon coordinates

  @param HourAngle : float
    decimal hours

  @param Declination : float
    decimal degrees

  @param Latitude : float
    decimal degrees

  @return: tuple
    Azimuth and elevation in degrees.
  """
  ha = HourAngle*math.pi/12.
  dec = Declination*math.pi/180.
  lat = Latitude*math.pi/180.
  azel = ManGord.coord(math.pi,math.pi/2-lat,0.,lat,ha,dec)
  az = coco.range_0_2pi(float(azel[0]))*180/math.pi
  el = coco.range_p_m_pi(float(azel[1]))*180/math.pi
  return (az,el)

def AzEl_to_HaDec(Azimuth,Elevation,Latitude):
  """
  Horizon coordinates to celestial

  @param Azimuth : float
    degrees, clockwise from north

  @param Elevation : float
    degrees

  @param Latitude : float
    degrees

  @return: tuple
    Hour angle in float hours and declination in float degrees
  """
  az = Azimuth*math.pi/180
  el = Elevation*math.pi/180
  lat = Latitude*math.pi/180
  HaDec = ManGord.coord(math.pi,math.pi/2-lat,0.,lat,az,el)
  HA = coco.range_0_2pi(float(HaDec[0]))*12/math.pi
  dec = coco.range_p_m_pi(float(HaDec[1]))*180/math.pi
  return (HA,dec)

def ha_rise_set(el_limit,lat,dec):
  """
  Hour angle from transit for rising and setting.

  Returns pi for a source that never sets and 0 for a source always below
  the horizon.

  @param el_limit : the elevation limit in radians
  @type  el_limit : float

  @param lat : the observatory latitude in radians
  @type lat : float

  @param dec : the source declination in radians
  @type  dec : float

  @return: hour angle from transit in radians
  """
  cos_ha = (math.sin(el_limit) - math.sin(lat)*math.sin(dec)) \
           /(math.cos(lat)*math.cos(dec))
  if cos_ha <= -1:
    # never sets
    return math.pi
  elif cos_ha >= 1:
    # never visible
    return 0
  else:
    return math.acos(cos_ha)

def show_catalogs():
  """
  List source catalogs

  Prints all the catalogs in the default catalog directory
  and returns them as a list.
  """
  files = glob.glob(default_astro_dir+'/sources/*.cat')
  print "Index of cats[]"
  cats = []
  for index in range(len(files)):
    cats.append(os.path.basename(files[index]))
    print index,cats[index]
  return cats

def get_sources(catalog):
  """
  Get all the sources from a catalog

  Notes
  =====
  Example:
  
  >>> source_data = get_sources('Methanol-II') source_data.keys()
  ['W51 E1/E2', 'NGC6334F', 'Orion KL', 'W31(1)', '34.26+0.16 ?', \\
  '345.01+1.79', 'NGC 7538', 'Cepheus A', 'W3(OH)', '0.54-0.86']
  >>> ra,dec,Vlsr = source_data['Orion KL'] print ra,dec,Vlsr
  05:32:47.0 -05:24:23 7.7

  @param catalog : string

  @return: dictionary of lists
    RA, decl. and Vlsr keyed to the source name
  """
  sources = open(default_astro_dir+'/sources/'+catalog+".cat",'r')
  source_data = sources.readlines()
  sources.close()
  names = []
  result={}
  for line in source_data[5:]:
    name = line.strip()[0:17].strip()
    names.append(name)
    data = line.strip()[18:].split()
    ra = data[0]
    dec = data[1]
    vel = data[2]
    result[name] = ra,dec,float(vel)
  return result

def get_coords(source,timestamp):
  """
  Returns the J2000.0 position of the source.
   
  For a planet the position is returned for the time given by
  a UNIX timestamp.

  Note
  ====
  Sidereal sources have not yet been coded.

  @param source : string
    source name

  @param timestamp : float
    UNIX time stamp for planetary position
  """
  planets = ['Moon', 'Mercury', 'Venus', 'Mars', 'Jupiter', 'Saturn',
             'Uranus', 'Neptune', 'Pluto']
  timetuple = time.gmtime(timestamp)
  year,month,day = timetuple[:3]
  mjd = MJD(year,month,day)
  hr,mn,sc = timetuple[3:6]
  UT = (hr + (mn + sc/60.)/60.)/24.
  try:
    ptr = planets.index(source.capitalize())
  except KeyError:
    # Could be a reference position:
    ptr = planets.index(source.capitalize()[:-2])
  except ValueError:
    # May be a sidereal object. Use database or catalogues
    return 0.,0.,0.
  track = jplephem.object_track(source,mjd,UT,1,0)
  ra_obs = track['ra'][0]
  dec_obs = track['dec'][0]
  coords2000 = jplephem.epoch_to_j2000(mjd,UT,ra_obs,dec_obs)
  epoch = coords2000['epoch']
  return coords2000['ra'],coords2000['dec'],coords2000['epoch']

def get_geodetic_coords(dss=0,observatory=None):
  """
  Observatory's geodetic coordinates
  
  An observatory may be specified with a DSS number or a name.

  If no observatory is specified, this returns a dictionary in which
  the station IDs are the keys and the values are tuples.

  If the observatory is specified it is used as a key and only that
  observatory's coordinates are returned as a tuple.

  If the observatory is not found, an empty dictionary is returned.
    
  Notes
  =====
  For an explanation of Earth coordinate systems so
  http://dsnra.jpl.nasa.gov/Antennas/Antennas.html#anchor950381
  
  Example:

  >>> station_data = get_geodedic_coords()
  >>> long,lat,alt,tz,name, diam = station_data[13]
  >>> print long,lat
  116.794009 35.2477189

  @param dss : int
    optional DSS number

  @param observatory : string
    optional observatory name

  @return: dictionary of tuples
    (longitude, latitude, altitude(m), time zone, name, diameter)
    keyed to the station name(s), if any
  """
  stations = open(default_astro_dir+"/stations.txt")
  stn_data = stations.readlines()
  stations.close()
  ids = []
  result={}
  for line in stn_data:
    stn_data = line.strip().split()
    # print stn_data
    id = int(stn_data[0])
    ids.append(id)
    long=float(stn_data[1])
    lat =float(stn_data[2])
    elev=float(stn_data[3])
    tz=int(stn_data[4])
    diam = int(stn_data[5])
    name=line.strip()[64:]
    if ( observatory != None \
      and observatory.strip.tolower() == name.strip().tolower()) \
    or dss == id:
      return     long,lat,elev,tz,name,diam
    result[id] = long,lat,elev,tz,name,diam
  return result

def decimal_day_to_HMS(day):
  """
  Formatted time to 0.01 seconds

  @param day : float

  @return: string
    'HH:MM:SS.SS"
  """
  return coco.dec_days_to_HMS(2,day)[1]

def time_aliases(year,UTdoy,obs_long):
  """
  Time as days since 1900, centuries since 1900 and LST

  Notes
  =====
  This simplifies the output from ManGord.atime().

  @param year : int

  @param UTdoy : float
    Day of year and niversal time

  @param obs_long : float
    Observatory west longitude in degrees

  @return: tuple of floats
    (number of days since 1900.5, the same
    in Julian centuries, and the local sidereal time in decimal days
  """
  days_since_1900, lst = ManGord.atime(year,UTdoy,obs_long)
  julian_centuries_since_1900 = days_since_1900/36525.
  return days_since_1900, julian_centuries_since_1900, lst

def current_to_ecliptic(julian_centuries_since_1900,\
                        current_ra,\
                        current_dec):
  """
  Apparent ecliptic coordinates from apparent RA and dec

  Notes
  =====
  Reformats output from ManGord.eclip()

  @param julian_centuries_since_1900 : float

  @param current_ra : float
    radians

  @param current_dec : float
    radians

  @return: tuple of floats
    (eclip_long, eclip_lat) in radians
  """
  eclip_long, eclip_lat = ManGord.eclip(julian_centuries_since_1900,\
                                        current_ra,\
                                        current_dec)
  return eclip_long[0], eclip_lat[0]

def J2000_to_apparent(MJD,UT,ra2000,dec2000):
  """
  Apparent right ascension and declination

  Notes
  =====
  This reformats input and output for jplephem.j2000_to_epoch()

  @param MJD : int
    mean Julian day

  @param UT : float
    decimal hours

  @param ra2000 : float
    radians

  @param dec2000 : float
    radians

  @return: tuple
  """
  p_apparent = jplephem.j2000_to_epoch(MJD,UT/24.,ra2000,dec2000)
  ra = float(p_apparent['ra'])
  dec = float(p_apparent['dec'])
  return float(p_apparent['ra']),float(p_apparent['dec'])

def decimal_hrs_to_hms(hours,num_dec_places):
  """
  Format decimal hours
  
  This function can process negative angles and does not put
  them in the range of 0 - 24.

  Notes
  =====
  rads_to_hms has a bug.  num_dec_places doesn't work)

  @param hours : float

  @param num_dec_places : int

  @return: string
  """
  list = formats.rads_to_hms(num_dec_places,hours*math.pi/12)
  if num_dec_places > 0:
    format = "%d:%02d:%0"+str(3+num_dec_places)+"."+str(num_dec_places)+"f"
    return format % (list[0],list[1],float(list[2]))
  else:
    format = "%d:%02d:%02d"
    return format % (list[0],list[1],int(round(float(list[2]))))

def decimal_deg_to_dms(degrees,num_dec_places):
  """
  Format decimal degrees

  Given an angle in decimal degrees, returns a string with
  SDD:MM:SS.SS... with the specified number of decimal places
  
  @param degrees : float

  @param num_dec_places : int

  @return: string
  """
  list = formats.rads_to_dms(num_dec_places,degrees*math.pi/180)
  if num_dec_places > 0:
    format = "%d:%02d:%0"+str(3+num_dec_places)+"."+str(num_dec_places)+"f"
    return format % (list[0],list[1],float(list[2]))
  else:
    format = "%d:%02d:%02d"
    return format % (list[0],list[1],int(round(float(list[2]))))

def blueshift_fraction(ra,dec,MJD,UT):
  """
  v/c of the observer relative to the solar system barycenter

  Notes
  =====
  This is the JPL ephemeris function doppler_fraction() with
  different input and output formats.

  @param ra : float
    RA of epoch in radians

  @param dec : float
    decl of epoch in radians

  @param MJD : int
    modified Julian date

  @param UT : float

  @return: float
    v/c of the observer relative to the solar system barycenter
    where v is positive for approaching objects.
  """
  dopfac = jplephem.doppler_fraction(ra*12/math.pi,dec*180/math.pi, \
                                     MJD,UT/24.,1,1)
  return dopfac['frac'][0]

def UTdatetime_from_ISO(ISO_string):
  """
  ISO time to time tuple

  @param ISO_string : string
    yyyy-mm-ddThh:mm:ss

  @return: tuple
    (year, day of year, time in decimal hours)
  """
  T = time.strptime(ISO_string,'%Y-%m-%dT%H:%M:%S')
  year = T[0]
  doy = T[7]
  decimal_time = T[3]+(T[4]+T[5]/60.)/60.
  return year,doy,decimal_time

def ecliptic_to_J2000(long,lat,mjd):
  """
  ecliptic longitude and latitude to celestial coordinates

  Notes
  =====
  Depends on module coco.
   
  @param long : float
    observatory longitude in degrees

  @param lat : float
    observatory latitude in degrees

  @param mjd : int
    modified Julian data

  @return: tuple of floats
    (J2000 right ascension and declination) in decimal hours and degrees.
  """
  longr = long*math.pi/180.
  latr = lat*math.pi/180.
  result = coco.eclip_to_J2000(longr,latr,mjd)
  rah = result[0]*12/math.pi
  decd = result[1]*180./math.pi
  return rah,decd

def great_circle_angle(lat1,long1,lat2,long2):
  """
  Great circle angle between two points on a sphere

  Given two positions on a sphere, return the angle between them
  in degrees.
  
  @param lat1 : float
    initial latitude in degrees

  @param long1 : float
    initial longitude in degrees

  @param lat2 : float
    final latitude in degrees

  @param long2 : float
    final longitude in degrees

  @return: float
  """
  dlongr = (long1 - long2)*math.pi/180.
  lat1r = lat1*math.pi/180.
  long1r = long1*math.pi/180.
  lat2r = lat2*math.pi/180.
  long2r = long2*math.pi/180.
  term1 = math.cos(lat2r)*math.sin(dlongr)
  term2 = math.cos(lat1r)*math.sin(lat2r)
  term3 = math.sin(lat1r)*math.cos(lat2r)*math.cos(dlongr)
  numer = math.sqrt(math.pow(term1,2)+math.pow(term2-term3,2))
  term4 = math.sin(lat1r)*math.sin(lat2r)
  term5 = math.cos(lat1r)*math.cos(lat2r)*math.cos(dlongr)
  denom = term4+term5
  cangler = math.atan2(numer,denom)
  return cangler*180/math.pi

def great_circle_distance(central_angle):
  """
  Great circle distance on Earth from great circle angle

  Notes
  =====
  Assumes a spherical earth of radius 6372.795 km.
  This is accurate to about 0.5%

  @param central_angle : float
    angle subtended at the center of the earth

  @return: float
    the distance in km
  """
  return 6372.795*math.pi*central_angle/180.

def isotime_to_time_struct(iso):
  """
  ISO time string to time tuple
  
  Converts an ISO timestamp
  (e.g. 2005-08-12T23:22:10 or 20050812T2322)
  to a time tuple as defined in module time
  (e.g. (2005, 8, 12, 23, 22, 10))
  
  @param iso : string

  @return: time tuple
  """
  try:
    result = time.strptime(iso, "%Y-%m-%dT%H:%M:%S")[0:6]
  except ValueError:
    try:
      result = time.strptime(iso,"%Y%m%dT%H%M")[0:6]
    except:
      result = (0, 0, 0, 0, 0, 0)
  return result

def isotime_to_datetime(iso):
  """
  ISO time to datetime.datetime object

  Converts an ISO timestamp
  (e.g. 2005-08-12T23:22:10 or 20050812T2322)
  to a datetime object, as defined in module datetime
  (e.g. datetime.datetime(2008, 5, 19, 13, 53, 7)

  @param iso : string

  @return: datetime:datetime object
  """
  time_str = isotime_to_time_struct(iso)
  result = datetime.datetime(*time_str)
  return result

def calendar_date(year, doy):
  """Calendar date from day of year

  @param year : int

  @param doy : int
    day of year

  @return: tuple of ints
    (year, month, day)
  """
  if doy < 32:
    month = 1
    day = doy
  elif doy < 60 + leap_year(year):
    month = 2
    day = doy - 31
  else:
    if leap_year(year) == 0:
      doy += 1
    month = int((doy+31.39)/30.61)
    day = doy + 2 - (month-1)*30-int((month+1)*0.61)
  return year,month,day

def v_dop(rah,decd,MJD,UT):
  """
  Radial velocity of the Local Standard of Rest
  
  Returns the radial velocity of the Local Standard of Rest
  with respect to the observer.

  Notes
  =====
  This should be more accurate than ManGord.dop() because it uses
  the JPL ephemeris.

  @param rah : float
    right ascension of epoch (observation) in hours

  @param decd : float
    declination of epoch (observation) in hours

  @param MJD : int
    modified Julian date

  @param UT : float
    hours

  @return: float
    km/s
  """
  ra = rah*math.pi/12
  dec = decd*math.pi/180
  v_sun_obs = -(c/1000.)*blueshift_fraction(ra,dec,MJD,UT)
  v_LSR_sun = v_sun(MJD,rah,decd)
  return v_sun_obs+v_LSR_sun

def east_long_to_west(deg,min,sec):
  """
  East logitude to west longitude.

  @param deg : int
    degrees

  @param min : int
    minutes

  @param sec : float
    seconds

  @return: float
    west longitude
  """
  east_long = deg+(min+sec/60.)/60.
  west_long = 360.-east_long
  return west_long


def get_dUT():
  """UT1 - UTC

  Returns a dictionary of UT1-UTC indexed by year and day-of-year

  Notes
  =====
  The data files are in
  /usr/local/lib/python2.5/site-packages/Astronomy/
  with names 'ut1-utc-*'
  """
  dUT = {}
  files = glob.glob("/usr/local/lib/python2.5/site-packages/Astronomy/ut1-utc-*")
  for f in files:
    nm = os.path.basename(f)
    yr = int(nm.split('-')[-1])
    dUT[yr] = {}
    fd = open(f,'r')
    lines = fd.readlines()
    fd.close()
    for l in lines:
      mn = int(l[4:6])
      dy = int(l[6:8])
      DOY = day_of_year(yr,mn,dy)
      the_rest = l[8:].split()
      MJD = float(the_rest[0])
      dUT[yr][DOY] = float(the_rest[1])
  return dUT

def find_source_name(ra,dec):
  """Find source name from 1950 coordinates.

  Searches the DSNRA catalogs for a source matching the coordinates.

  @param ra : float::
    Degrees, FITS compatible

  @param dec : float::
    Degrees

  @return: string::
    Source name, or source name followed by ?, or U followed by
    coordinates in IAU convention
  """
  ra_str = decimal_hrs_to_hms(ra/15,0).split('.')[0]
  dec_str = decimal_deg_to_dms(dec,0).split('.')[0]
  if diag:
    print "find_source_name: coordinates:",ra_str,dec_str
  ra_test = []
  ra_test.append(str((int(ra) - .1) % 360)+"[.]")
  ra_test.append(str((int(ra)     ) % 360)+"[.]")
  ra_test.append(str((int(ra) + .1) % 360)+"[.]")
  hits = []
  for test in ra_test:
    command = 'grep -h '+test+' '+cat_path+'*.deg'
    if diag:
      print command
    fd = os.popen(command)
    found = fd.readlines()
    fd.close()
    if diag:
      print len(found),"hits"
    hits += found
  if len(hits) > 0:
    # Find the nearest, but not more than a degree away
    delta_min = 1.
    for hit in hits:
      if diag:
        print "find_source_name: match:",hit.strip()
      parts = hit[20:].split()
      ra_found,dec_found = float(parts[0]), float(parts[1])
      if diag:
        print "find_source_name: coords:",ra_found,dec_found
      delta = great_circle_angle(ra, dec, ra_found, dec_found)
      if diag:
        print "find_source_name: delta:",delta
      if delta < delta_min:
        ra_save,dec_save = ra_found,dec_found
        delta_min = delta
        source = hit[:17].strip()
    if delta_min < 0.0333:
      # This is an exact identification
      return source
    elif delta_min < 1.:
      # This is an approximate identification
      return source+'?'
  # No hits
  ra_parts = ra_str.split(":")
  # We expect NNN
  if len(ra_parts[0]) < 2:
    # Only one digit; prefix with 0
    ra_parts[0] = "0" + ra_parts[0]
  dec_parts = dec_str.split(":")
  if len(dec_parts[0]) == 1:
    # Only one digit; must be positive; prefix with +0
    dec_parts[0] = "+0" + dec_parts[0]
  elif len(dec_parts[0]) == 2:
    # could be NN or -N
    if dec > 0.0:
      # must be NN; prefix with
      dec_parts[0] = '+' + dec_parts[0]
    else:
      dec_parts[0] = dec_parts[0][0]+'0'+dec_parts[0][1]
  return "U "+ra_parts[0]+ra_parts[1]+dec_parts[0]+dec_parts[1]


def days_since_1900_noon(MJD):
  """
  MJD to ephemeris date in fractional days since noon Jan.1, 1900

  @param MJD : modified Julian date
  @type  MJD : float

  @return: days since noon Jan.1, 1900, float
  """
  return MJD - 15019.5

def MJD_de_ephem_date(ephem_date):
  """
  Convert ephemeris date () to modified JD

  @param ephem_date : days since noon, Jan. 1, 1900
  @type  ephem_date : float

  @return: modified Julian date (float)
  """
  return ephem_date + 15019.5

def cal_date_de_ephem_date(ephem_date):
  MJD = MJD_de_ephem_date(ephem_date)
  return MJD_to_cal(1,MJD)

def format_date(date_tuple):
  """
  Convert a date tuple to YYYY/MM/DD string

  @param date_tuple : calendar date, (year, month, day)
  @type  date_tuple : tuple of ints
  """
  return str(date_tuple[0])+"/"+str(date_tuple[1])+"/"+str(date_tuple[2])

def LST_de_UT(UT,LST0):
  """
  Convert UT to LST

  The inputs and outputs have the same units, whatever they are.

  @param UT : universal time
  @type  UT : float

  @param LST0 : local sidereal time at UT = 0
  @type  LST0 : float

  @return: local sidereal time in the same units as UT and LST0
  """
  return LST0 + UT

def UT_de_LST(LST,LST0):
  """
  Convert LST to UT

  The inputs and outputs have the same units, whatever they are.

  @param LST : local sidereal time
  @type  LST : float

  @param LST0 : local sidereal time at UT = 0
  @type  LST0 : float

  @return: universal time in the same units as UT and LST0
  """
  return LST - LST0

def greenwich_sidereal_time(year,doy):
  """
  Approximate sidereal time at Greenwich

  @param year : year of date
  @type  year : int

  @param doy : day of year
  @type  doy : float

  @return: Greenwich sidereal time in hours
  """
  year_from_1966 = year-1966
  dt = (year_from_1966*365 + int((year_from_1966 + 1)/4.) + int(doy)-1)/36525.
  dst = 0.278329562 + (8640184.67*dt+0.0929*dt**2)/86400
  gst0 = dst % 1 # GST on Jan. 0 of current year
  return 24*(gst0 + (doy % 1)/0.997269566) % 24

def get_source_positions(catalog,targets):
  """
  Get the positions for selected sources in the catalog

  @param catalog : catalog name without suffix
  @type  catalog : str

  @param targets : list of names from the catalog
  @type  targets : list

  @return: dictionary of decimal coordinate pairs keyed on targets
  """
  sources = open(cat_path+catalog+'.cat','r')
  source_data = sources.readlines()
  sources.close()
  P = {}
  for line in source_data[5:]:
    name = line.strip()[0:17].strip()
    if name in targets:
      data = line.strip()[18:].split()
      position_str = [data[0],data[1]]
      ra_string = position_str[0]
      dec_string = position_str[1]
      P[name] = formats.parse_colon_delimited_angles(ra_string,dec_string)
  if diag:
    print "Positions:",P
  return P

def choose_cal_source(dss,cal_sources,timestamp):
  """
  Select calibrator by number from a list that are up at a certain time

  Notes
  =====
  This needs a flux calculation. Use MySQL database. See
  /usr/local/pkg/astro/lib/astro/sources/calibrators/plot_spectra.py
  for interpolation.

  @param dss : DSN station number
  @type  dss : int

  @param cal_sources : list of sources from catalog 'Strong_cals'
  @type  cal_sources : list of ints

  @param timestamp : time for computing azimuth and elevation
  @type  timestamp : time in seconds since the epoch, as in module time

  @return: name of the selected calibrator
  """
  lng,lat,elev,tz,name = get_geodetic_coords(dss=dss)
  source_data = get_sources('Strong_cals')
  time_tuple = time.gmtime(timestamp)
  year = time_tuple[0]
  DOY = time_tuple[7]
  UT = time_tuple[3]+time_tuple[4]/60.
  days_since_1900, centuries_since_1900,LST_0h = time_aliases(year, DOY, lng)
  lst = math.fmod(UT/24.+LST_0h,1.)*24
  print "Calibrators in order of flux:"
  for index in range(len(cal_sources)):
    ra_str,dec_str,Vlsr = source_data[cal_sources[index]]
    ra,dec = formats.parse_colon_delimited_angles(ra_str, dec_str)
    HA = coco.range_p_m_pi((lst - ra)*hrs2rad)*rad2hrs
    az,el = HaDec_to_AzEl(HA, dec, lat)
    if el > 10:
      print "%d: %10s, Az=%3.0f, El=%3f" % (index,cal_sources[index],az,el)
  choice = raw_input("Select calibrator by numbers: ")
  return cal_sources[int(choice)]

def get_station_parameters(stations):
  """
  Get the geodetic and Cartesian coordinates of the DSN stations

  @param stations : a list of stations
  @type  stations : list of ints

  @return: dictionary of dictionaries.
  """
  station_pars = {}
  for station in stations:
    station_data = get_geodetic_coords(dss=station)
    station_pars[station] = {}
    station_pars[station]['long'] = float(station_data[0])
    station_pars[station]['lat']  = float(station_data[1])
    station_pars[station]['elev'] = float(station_data[2])
    station_pars[station]['TZ']   = int(station_data[3])
    station_pars[station]['name'] = station_data[4]
    x,y,z = get_cartesian_coordinates('DSS '+str(station))
    station_pars[station]['x'] = x
    station_pars[station]['y'] = y
    station_pars[station]['z'] = z
  return station_pars

def AzEl_to_RaDec(azimuth,elevation,latitude,longitude,date_time):
  """
  Convert azimuth and elevation to right ascension and declination

  @param azimuth : east from north (clockwise) in degrees
  @type  azimuth : float

  @param elevation : above the horizon in degrees
  @type  elevation :

  @param latitude : above the equator, in degrees
  @type  latitude : float

  @param longitude : west from Greenwich
  @type  longitude : float

  @param date_time : (year, DOY) tuple, with fractional day of year
  @type  date_time : (int, float)

  @return: (RA (hrs), dec (degs))
  """
  LST = greenwich_sidereal_time(*date_time)-longitude/15.
  HA,dec  = AzEl_to_HaDec(azimuth, elevation, latitude)
  RA = math.fmod(LST - HA, 24.)
  if RA < 0:
    RA += 24.
  return RA,dec