#!/usr/bin/python
# makebinphot.py 1.0
# Jure Skvarc, May 2000
#
# Make binary version of photometric standard files and create
# accelerator files
# Take advantage of the fact that files are already sorted by RA

import fileinput
from string import *
import struct

#  Class which holds all data about the star
class Star:
    pass

#  Compare two stars on declination
def cmpdec(a1, a2):
    if a1.dec > a2.dec:
        return 1
    elif a1.dec < a2.dec:
        return -1
    else:
        return 0

#  Compare two stars on RA
def cmpra(a1, a2):
    if a1.ra > a2.ra:
        return 1
    elif a1.ra < a2.ra:
        return -1
    else:
        return 0



#  In a sorted array finds an index i where dec[i] >= d and dec[i - 1] < d 
def findex(stars, d):

  first = 0
  nstar = len(stars)
  last = nstar

  #  0 is a special case
  if d < stars[0].dec: return 0
  #  Check if d is out of range
  if d > stars[nstar - 1].dec: return -nstar
  
  while last - first > 1:
    k = (first + last) / 2
    if d < stars[k].dec:
      if d < stars[k - 1].dec:
	last = k + 1
      else:
	return k
    elif d > stars[k].dec:
      first = k
    else: return k
  return k - 1



usnofile = 'ubvri_all.cat'  #  Name of the USNO file
i = 0
stars = []
for line in fileinput.input(usnofile):
    # Skip lines starting with #
    if line[0] == '#': continue
    cols = split(line)
    star = Star()
    i = i + 1
    star.name = '%s-%d' % (cols[0], i)
    star.ra = float(cols[1])
    star.sra = float(cols[2])
    star.dec = float(cols[3])
    star.sdec = float(cols[4])
    star.u = float(cols[6])
    star.su = float(cols[7])
    star.b = float(cols[8])
    star.sb = float(cols[9])
    star.v = float(cols[10])
    star.sv = float(cols[11])
    star.r = float(cols[12])
    star.sr = float(cols[13])
    star.i = float(cols[14])
    star.si = float(cols[15])
    stars.append(star)  #  Add star to the list

loneosfile = 'loneos.cat'
for line in fileinput.input(loneosfile):
    # Skip lines starting with #
    if line[0] == '#': continue
    cols = split(line)
    star = Star()
    i = i + 1
    star.name = ('%s-%d' % (cols[0], i))[0:19]
    star.ra = float(cols[1])
    star.dec = float(cols[2])
    star.sra = 1
    star.sdec = 1
    star.u = 99.9
    star.su = 99.9
    star.b = float(cols[6])
    star.v = float(cols[7])
    star.r = float(cols[8])
    star.i = float(cols[9])
    if cols[5] == 'y':
        dm = 0.01
    else:
        dm = 0.1
    if star.b != 99.0:
        star.sb = dm
    else: 
        star.sb = 99.0
    if star.v != 99.0:
        star.sv = dm
    else:
        star.sv = 99.0
    if star.r != 99.0:
        star.sr = dm
    else:
        star.sr = 99.0
    if star.i != 99.0:
        star.si = dm
    else:
        star.si = 99.0
    stars.append(star)  #  Add star to the list


#  Sort the complete list on RA
stars.sort(cmpra)

# Find any stars that may be duplicates, print them to a special file
# and remove the one with lower accuracy
i = 0
fp = open('duplicates.dat', 'w')
while i < len(stars):
    print i
    j = i + 1
    ra = stars[i].ra
    dec = stars[i].dec
    sra = stars[i].sra
    sdec = stars[i].sdec
    while j < len(stars):
        #  Difference in coordinates
        dra = abs(ra - stars[j].ra)
        s5ra = 2.0 / 3600
        if dra > s5ra: break  #  No match is possible any more, so exit loop
        ddec = abs(dec - stars[j].dec)
        #  Measure for error (2 arcseconds)
        s5dec = 2.0 / 3600
        if ddec <= s5dec:
            #  Two stars match, print them out
            tpl = '%s %f %f %f %f %f %f %f %f %f %f '
            dstar1 = tpl % (stars[i].name, ra, dec,\
                            stars[i].b, stars[i].sb, stars[i].v, stars[i].sv,\
                            stars[i].r, stars[i].sr, stars[i].i, stars[i].si)
            dstar2 = tpl % (stars[j].name, stars[j].ra, stars[j].dec,\
                            stars[j].b, stars[j].sb, stars[j].v, stars[j].sv,\
                            stars[j].r, stars[j].sr, stars[j].i, stars[j].si)
            fp.write(dstar1 + dstar2 + '\n')
            #  Combine magnitudes and 
            if stars[i].sb > stars[j].sb:
                stars[i].b = stars[j].b
                stars[i].sb = stars[j].sb
            if stars[i].sv > stars[j].sv:
                stars[i].v = stars[j].v
                stars[i].sv = stars[j].sv
            if stars[i].sr > stars[j].sr:
                stars[i].r = stars[j].r
                stars[i].sr = stars[j].sr
            if stars[i].si > stars[j].si:
                stars[i].i = stars[j].i
                stars[i].si = stars[j].si
            #  Remove the star j
            del stars[j]
            break
        j = j + 1
    i = i + 1
fp.close()    
limit = 7.5
stars1 = []
for star in stars:
    if star.ra < limit:
        stars1.append(star)
    else:
        #  Sort the stars by declination
        stars1.sort(cmpdec)
        #  Output them to a binary file
        h, m = divmod(limit - 7.5, 15)
        m = 4 * m
        name = 'phot%02d%02d.bin' % (h, m)  # create file name
        fp = open(name, 'wb')  # open the file for binary writing
        for a in stars1:
            fp.write(struct.pack('20sdfdfffffffffff', a.name, a.ra,\
                                 a.sra, a.dec, a.sdec, \
                                 a.u, a.su, a.b, a.sb,\
                                 a.v, a.sv, a.r, a.sr,\
                                 a.i, a.si))
        fp.close()
        # Now create the accelerator file
        accname = 'phot%02d%02d.acc' % (h, m) 
        fp = open(accname, 'w')  # open the file for binary writing
        kprev = 0
        for j in range(0, 19):
            m = 10.0 * (j - 9)
            k = findex(stars1, m)
            if k > 0:
                fp.write('%10f %6d %10f\n' % (m, k, stars1[k].dec))
                kprev = k
            else:
                fp.write('%10f %6d %10f\n' % (m, kprev, m))
        fp.close()
        #  Go to the next subsection
        limit = limit + 7.5
        #  Destroy star list
        del stars1
        #  Create s new one
	stars1 = []

        
#  Sort the stars by declination
stars1.sort(cmpdec)
#  Output them to a binary file
h, m = divmod(limit - 7.5, 15)
m = 4 * m
name = 'phot%02d%02d.bin' % (h, m)  # create file name
fp = open(name, 'wb')  # open the file for binary writing
for a in stars1:
    fp.write(struct.pack('20sdfdfffffffffff', a.name, a.ra,\
                 a.sra, a.dec, a.sdec, \
                 a.u, a.su, a.b, a.sb,\
                 a.v, a.sv, a.r, a.sr,\
                 a.i, a.si))
fp.close()
# Now create the accelerator file
accname = 'phot%02d%02d.acc' % (h, m) 
fp = open(accname, 'w')  # open the file for binary writing
for j in range(0, 19):
    m = 10.0 * (j - 9)
    k = findex(stars1, m)
    if k > 0:
        fp.write('%10f %6d %10f\n' % (m, k, stars1[k].dec))
    else:
        fp.write('%10f %6d %10f\n' % (m, k, m))
fp.close()



