#!/usr/bin/env python3
# -*- coding: utf-8 -*-

#**************************************************************************
#   Copyright (C) 2019, Paul Lutus                                        *
#                                                                         *
#   This program is free software; you can redistribute it and/or modify  *
#   it under the terms of the GNU General Public License as published by  *
#   the Free Software Foundation; either version 2 of the License, or     *
#   (at your option) any later version.                                   *
#                                                                         *
#   This program is distributed in the hope that it will be useful,       *
#   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
#   GNU General Public License for more details.                          *
#                                                                         *
#   You should have received a copy of the GNU General Public License     *
#   along with this program; if not, write to the                         *
#   Free Software Foundation, Inc.,                                       *
#   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
#**************************************************************************


import re
import sys
import math
from PIL import Image, ImageDraw

# user customizations:
  
image_width = 16000
radius_multiplier = 3
brightness_multiplier = 8
# CSV star database from http://www.astronexus.com/hyg
path_to_star_database = 'star_positions/hygdata_v3.csv'
path_to_output_graphic = 'star_field_equirectangular.png'

# end of user settings


ww = image_width
hh = int(ww/2)

img = Image.new('RGB', (ww,hh))
 
draw = ImageDraw.Draw(img)

# interpolation function

def ntrp(x,xa,xb,ya,yb):
     return (x-xa) * (yb-ya) / (xb-xa) + ya

# convert magnitude to brightness
# https://www.phys.ksu.edu/personal/wysin/astro/magnitudes.html
     
def bri(m): return math.exp(-0.4 * m)

#n = 0
found = 0
fhash = {}

with open(path_to_star_database) as f:
  for n,line in enumerate(f):
    fields = re.split(',',line)
    # field name record
    if(n == 0):
      #print(fields)
      for idx,name in enumerate(fields):
        fhash[name] = idx
      #print(fhash)
    else:
      ra = float(fields[fhash['ra']])
      dec = float(fields[fhash['dec']])
      mag = float(fields[fhash['mag']])
      # brightness
      b = bri(mag)
      # don't include the sun among the stars
      if(b < 100):
        found += 1
        r = b * radius_multiplier
        # convert ra, decl into graphics coordinates
        # ra range 0 to 24h
        # decl range -90 to 90 degrees
        x = int(ntrp(ra,24,0,0,ww))
        y = int(ntrp(dec,90,-90,0,hh))
        # color (pixel brightness)
        c = b * 255 * brightness_multiplier
        # control color range
        c = max(0,c)
        c = min(255,c)
        c = int(c)
        draw.ellipse(((x-r,y-r),(x+r,y+r)),fill=(c,c,c))
        sys.stdout.write('creating star %06d ...\r' % n)
       

print("\nCreating graphic...")
        
del(draw)

img.save(path_to_output_graphic)