StarPro star file generation scripts

flat_to_bin.py

#!/usr/bin/python

# Convert flat text file of star data into a packed binary format for StarPro.
#
# Input: one star per line, in the text format
# RA DEC MAG RED GREEN BLUE RANGE
#
# where RA/DEC are in degrees, RED GREEN BLUE are floating-point with the greatest
# of the three normalized to 1.0. RANGE is light-years from sun, or -1 for unknown (infinite) range.
#
# Output: (all fields LITTLE-ENDIAN!)
# struct star_rec {
#     double    ra;           8 bytes   (degrees)
#     double   dec;           8 bytes   (degrees)
#     float    lum;           4 bytes   (linear luminance, not magnitude!)
#     u8   red, green, blue;  3 bytes   (color, normalized to 255, sRGB)
#     u8   pad[1];            1 byte
#     double  range;          8 bytes   (range from sun, light years)
# };
#             total:         32 bytes

import sys, struct, math, random, getopt, os

random.seed(42) # seed the random number generator consistently

use_range = 1 # include range data in output.


# linear luminance <-> sRGB conversions
def sRGB_encode(f):
    if f <= 0.000304:
        return f * 12.92
    else:
        return 1.055*math.pow(f, 1.0/2.4) - 0.055

def float_to_sRGB(f):
    return int(sRGB_encode(f)*255.0 + 0.5)

# get the latitude and longitude of a random point
# uniformly distributed on the unit sphere
# returns (theta, phi)
# where 0 < theta < 2PI
#       -PI/2 < phi < PI/2
def sphere_point():
    # from http://mathworld.wolfram.com/SpherePointPicking.html
    theta = 2 * math.pi * random.random()
    phi = math.acos(2*random.random()-1) - math.pi/2.0
    return (theta, phi)

# return a random (right ascension, declination) tuple
def rand_ra_dec():
    (theta, phi) = sphere_point()
    ra = theta * 180.0/math.pi
    dec = phi * 180.0/math.pi
    return (ra, dec)

# utility to replace a file on disk atomically
class AtomicFileWrite:
    def __init__(self, filename, mode):
        self.filename = filename
        self.mode = mode
        self.temp_filename = filename + '.inprogress'
        self.fd = None
    def __enter__(self):
        self.fd = open(self.temp_filename, self.mode)
        return self
    def complete(self, fsync = True):
        self.fd.flush()
        if fsync:
            os.fsync(self.fd.fileno())
        os.rename(self.temp_filename, self.filename)
        self.fd.close()
        self.fd = None
    def __exit__(self, type, value, traceback):
        if self.fd is not None:
            try:
                os.unlink(self.temp_filename)
            except:
                pass
            self.fd.close()
            self.fd = None

if __name__ == '__main__':
    quiet = False

    # number of copies of the star database to include.
    copies = 1

    opts, args = getopt.gnu_getopt(sys.argv[1:], 'q', ['copies=','quiet','seed='])

    for key, val in opts:
        if key == '--copies': copies = int(val)
        elif key in ('-q', '--quiet'): quiet = True
        elif key == '--seed': random.seed(int(val))

    if len(args) < 1:
        sys.stderr.write('usage: %s input.flat output.stars\n' % sys.argv[0])
        sys.exit(1)

    # write output file atomically, so it gets deleted in case of error instead of leaving a partial file
    with AtomicFileWrite(args[1], 'wb') as out_atom:
        out_fd = out_atom.fd

        # write out the header (32 bytes)
        if use_range:
            hdr = "STARPRO2"
        else:
            hdr = "STARPRO1"

        hdr += " " * (32 - len(hdr)) # pad to 32 bytes
        out_fd.write(hdr)

        total = 0

        for copynum in xrange(copies):
            if args[0] == '-':
                in_fd = sys.stdin
            else:
                in_fd = open(args[0], 'r') # read flat data from this file

            while 1:
                line = in_fd.readline()

                if not line: break
                if line[0] == '#': continue

                (ra, dec, mag, red, green, blue, range) = map(float, line.strip().split())

                # if using multiple copies, randomize the RA/DEC of the star
                if copies > 1:
                    ra, dec = rand_ra_dec()

                # convert magnitude to linear luminance
                lum = math.pow(10.0, -(mag/2.5))

                # convert red, green, blue to 0-255
                red = float_to_sRGB(red)
                green = float_to_sRGB(green)
                blue = float_to_sRGB(blue)

                # write out the binary struct (little-endian)
                out_fd.write(struct.pack("<ddfBBBBd", ra, dec, lum, red, green, blue, 0, range))
                total += 1

        if not quiet:
            sys.stderr.write('wrote %d total stars\n' % total)

        out_atom.complete()

genstars.py

#!/usr/bin/python
#
# Generate random stars
# prints "flat" text format to stdout:
# RA DEC MAG RED GREEN BLUE

from random import random
import math, sys

def sphere_point():
    # get the latitude and longitude of a random point
    # uniformly distributed on the unit sphere
    # returns (theta, phi)
    # where 0 < theta < 2PI
    #       -PI/2 < phi < PI/2

    # from http://mathworld.wolfram.com/SpherePointPicking.html
    theta = 2 * math.pi * random()
    phi = math.acos(2*random()-1) - math.pi/2.0
    return (theta, phi)

def rand_mag_dan():
    # choose magnitude based on a heuristic
    # function of the form 1-exp(-x)
    # dimmest star magnitude possible
    max_mag = 8.0
    
    dim_pop = 0.2
    bright_pop = 0.39
    return max_mag*(1-math.exp(-(1.0/dim_pop)*math.pow(random(),bright_pop)))

def rand_mag_pavel():
    l = -1.0
    u = 8.5
    # c2 higher -> fewer bright stars
    c2 = 1.7
    return math.log((math.exp(c2*u)-math.exp(c2*l))*random() + math.exp(c2*l))/c2

def rand_star():
    # returns (ra, dec, mag, r, g, b)
    (theta, phi) = sphere_point()

    # conv to degrees
    ra = theta * 180.0/math.pi
    dec = phi * 180.0/math.pi

    #mag = rand_mag_dan()
    mag = rand_mag_pavel()

    n = random()
    if random() > 0.3:
        # reddish
        r = 1.0
        g = 1.0 - 0.4*n
        b = 1.0 - 0.99*n
    else:
        # bluish
        b = 1.0
        r = 1.0 - 0.9*n
        g = 1.0 - 0.6*n

    return (ra, dec, mag, r, g, b)

if __name__ == "__main__":
    if len(sys.argv) != 2:
        print "usage: genstars <n> (number of stars to generate)"
        sys.exit(-1)

    for i in range(int(sys.argv[1])):
        print "%f %f %f %f %f %f" % rand_star()