import csv
import colorsys
import numpy as np

# Compute the color of a star from its color index:
# see https://stackoverflow.com/questions/21977786/star-b-v-color-index-to-apparent-rgb-color
def bv2rgb(bv):
    r,g,b = 0,0,0
    #extrema and out of range values
    if (bv<-0.4):
        bv=-0.4
    if (bv> 2.0):
        bv= 2.0
    #red
    if bv>=-0.40 and bv<0.00:
        t=(bv+0.40)/(0.00+0.40)
        r=0.61+(0.11*t)+(0.1*t*t)
    elif bv>= 0.00 and bv<0.40:
        t=(bv-0.00)/(0.40-0.00)
        r=0.83+(0.17*t)
    elif bv>= 0.40 and bv<2.10:
        t=(bv-0.40)/(2.10-0.40)
        r=1.00
    #green
    if bv>=-0.40 and bv<0.00:
        t=(bv+0.40)/(0.00+0.40)
        g=0.70+(0.07*t)+(0.1*t*t)
    elif bv>= 0.00 and bv<0.40:
        t=(bv-0.00)/(0.40-0.00)
        g=0.87+(0.11*t)
    elif bv>= 0.40 and bv<1.60:
        t=(bv-0.40)/(1.60-0.40)
        g=0.98-(0.16*t)
    elif bv>= 1.60 and bv<2.00:
        t=(bv-1.60)/(2.00-1.60)
        g=0.82-(0.5*t*t)
    #blue
    if bv>=-0.40 and bv<0.40:
        t=(bv+0.40)/(0.40+0.40)
        b=1.00
    elif bv>= 0.40 and bv<1.50:
        t=(bv-0.40)/(1.50-0.40)
        b=1.00-(0.47*t)+(0.1*t*t)
    elif bv>= 1.50 and bv<1.94:
        t=(bv-1.50)/(1.94-1.50)
        b=0.63-(0.6*t*t)
    return [r,g,b]

# Darken the rgb accordingly to magnitude
# see https://blender.stackexchange.com/questions/80034/fix-hsv-to-rgb-conversion/80047
"""
def linerRGBToHSV(r,g,b)
def linear_to_srgb(r, g, b):
    def linear(c):
        a = .055
        if c <= .0031308:
            return 12.92 * c
        else:
            return (1+a) * c**(1/2.4) - a
    return [linear(c) for c in (r, g, b)]

def srgb_to_linear(r, g, b):
    def srgb(c):
        a = .055
        if c <= .04045:
            return c / 12.92
        else:
            return ((c+a) / (1+a)) ** 2.4
    return tuple(srgb(c) for c in (r, g, b))
"""

starFile = "/home/loic/Downloads/hygdata_v3_full.csv"

reader = csv.DictReader(open(starFile))
verts, tris, cols = [], [], []
ind = 0
offsetY = 0.25*np.sqrt(3)/4 #offset to make an equilateral triangle
for i, row in enumerate(reader):
    if i>0 and i<300000:
        #Get the data
        x = float(row["x"])
        y = float(row["y"])
        z = float(row["z"])
        #mag = float(row["absmag"])
        col = [int(255*w) for w in bv2rgb(float(row["ci"]))] if row["ci"]!="" else [200,200,200]
        if 1:#(x**2 + y**2 + z**2)**0.5 < 1000:
            verts.append([x,y,z])
            #tris.append(tri)
            cols.append(col)
            #iterate and info
            print i

with open("out.ply", "w") as f:
    f.write("ply\nformat ascii 1.0\n")
    f.write("element vertex " + str(len(verts)) + "\n")
    f.write("property float x\nproperty float y\nproperty float z\n")
    f.write("property uchar red\nproperty uchar green\nproperty uchar blue\n")
    f.write("end_header\n")
    for v,c in zip(verts, cols):
        f.write(" ".join([str(x) for x in [v[0], v[1], v[2], c[0], c[1], c[2]] ]) + "\n")
    f.write("\n")