Friday, September 24, 2021

Star Structure GAIA Sat ESA ER3 Dataset 5000 stars of least parallax error







import OpenGL
import csv
import math
import numpy as np
import locale
import random
from locale import atof
from time import sleep
from OpenGL.GL import *

from OpenGL.GLUT import *

from OpenGL.GLU import *

def wavelength_to_rgb(wavelength, gamma=0.8):

    '''This converts a given wavelength of light to an 
    approximate RGB color value. The wavelength must be given
    in nanometers in the range from 380 nm through 750 nm
    (789 THz through 400 THz).

    Based on code by Dan Bruton
    http://www.physics.sfasu.edu/astro/color/spectra.html
    '''

    wavelength = 380.0+370.0*float(wavelength)
    if wavelength >= 380 and wavelength <= 440:
        attenuation = 0.3 + 0.7 * (wavelength - 380) / (440 - 380)
        R = ((-(wavelength - 440) / (440 - 380)) * attenuation) ** gamma
        G = 0.0
        B = (1.0 * attenuation) ** gamma
    elif wavelength >= 440 and wavelength <= 490:
        R = 0.0
        G = ((wavelength - 440) / (490 - 440)) ** gamma
        B = 1.0
    elif wavelength >= 490 and wavelength <= 510:
        R = 0.0
        G = 1.0
        B = (-(wavelength - 510) / (510 - 490)) ** gamma
    elif wavelength >= 510 and wavelength <= 580:
        R = ((wavelength - 510) / (580 - 510)) ** gamma
        G = 1.0
        B = 0.0
    elif wavelength >= 580 and wavelength <= 645:
        R = 1.0
        G = (-(wavelength - 645) / (645 - 580)) ** gamma
        B = 0.0
    elif wavelength >= 645 and wavelength <= 750:
        attenuation = 0.3 + 0.7 * (750 - wavelength) / (750 - 645)
        R = (1.0 * attenuation) ** gamma
        G = 0.0
        B = 0.0
    else:
        R = 0.0
        G = 0.0
        B = 0.0
    return [R,G, B]

def init():
global Z
global triangles
global angle
global myscale
global colorbprp
global cbprp_rgb
global between_brightest
global between_smallest
random.seed()
cbprp_rgb=[0.0,0.0,0.0]
Z=0.0
angle=0.0
triangles= []
myscale=[0.0,0.0,0.0] 
average=[0.0,0.0,0.0]
maximum=[0.0,0.0,0.0]
minimum=[1000000.0,1000000.0,100000.0]
line_count =0.0
step=random.randint(80,120)
with open('/home/pamela/Documents/gaia-sample-1.csv') as csv_file:
csv_reader = csv.reader(csv_file, delimiter=',')

for row in csv_reader:
if line_count==0.0 :
     line_count=0.0001
     print row[69],row[9],row[5],row[7],row[86]
else :
if int(line_count)<10000 :
step=random.randint(80,120)
c1=9
c2=5
c3=7
c5=86
row[c1]=''.join(c for c in row[c1] if c.isdigit() or '.')
row[c2]=''.join(c for c in row[c2] if c.isdigit() or '.')
row[c3]=''.join(c for c in row[c3] if c.isdigit() or '.' )
row[44]=''.join(c for c in row[44] if c.isdigit()or '.')
row[c5]=''.join(c for c in row[c5] if c.isdigit() or '.')
if row[c1]!='' and row[c2]!='' and row[c3]!='' and row[44]!=''and row[c5]!=''   :
parallax=1.0/atof(row[c1])*10000000000.0
RA=(atof(row[c2])-180.0)
DEC=atof(row[c3])
x=parallax*np.cos(DEC)*np.cos(RA)
y=parallax*np.cos(DEC)*np.sin(RA)
z=parallax*np.sin(DEC)
if x>maximum[0] :
maximum[0]=x
if y>maximum[1] :
maximum[1]=y
if z>maximum[2] :
maximum[2]=z
if x<minimum[0] :
minimum[0]=x
if y<minimum[1] :
minimum[1]=y
if z<minimum[2] :
minimum[2]=z
average[0] += x
average[1] += y
average[2] += z
myrgb= wavelength_to_rgb((atof(row[c5])-0.6)/1.17)
triangles.append([[x,y,z],0.50+math.log(atof(row[69])/atof(row[c1]),10.0),myrgb])
print myrgb,(atof(row[c5])-0.6)/1.17
line_count += 2.0
print(line_count)
average[0] = average[0]/line_count
average[1] = average[1]/line_count
average[2] = average[2]/line_count
print(maximum[0],",",maximum[1],",",maximum[2])
print(minimum[0],",",minimum[1],",",minimum[2])
myscale[0]=math.fabs(maximum[0]-minimum[0])
myscale[1]=math.fabs(maximum[1]-minimum[1])
myscale[2]=math.fabs(maximum[2]-minimum[2])
medianx=[]
mediany=[]
medianz=[]
for v,cbprp,myrgb in triangles:
      medianx.append(v[0])
average[0]=np.median(medianx)
for v,cbprp,myrgb in triangles:
      mediany.append(v[1])
average[1]=np.median(mediany)
for v,cbprp,myrgb in triangles:
      medianz.append(v[2])
average[2]=np.median(medianz)
print(average[0],average[1],average[2])
print(myscale[0],myscale[1],myscale[2])

triangles2=[]
for v,cbprp,myrgb in triangles:
v[0]=(v[0]-average[0])/math.fabs(myscale[0])
v[1]=(v[1]-average[1])/math.fabs(myscale[1])
v[2]=(v[2]-average[2])/math.fabs(myscale[2])
triangles2.append([v,cbprp,myrgb])
triangles=triangles2
bright_stars=[]
small_stars=[]
for v,cbprp,myrgb in triangles:
if cbprp>2.0 :
bright_stars.append([v,cbprp])
else :
small_stars.append([v,cbprp])
between_brightest=[]
between_smallest=[]
for v1,cbprp1 in bright_stars:
for v2,cbprp2 in bright_stars:
dist=pow((v1[0]-v2[0])*(v1[0]-v2[0])+(v1[1]-v2[1])*(v1[1]-v2[1])+(v1[2]-v2[2])*(v1[2]-v2[2]),0.5)
if v1[0]!=v2[0] and v1[1]!=v2[1]  and v1[2]!=v2[2] and cbprp1>cbprp2 and dist<0.06:
between_brightest.append([v1,v2])
for v1,cbprp1 in small_stars:
for v2,cbprp2 in small_stars:
dist=pow((v1[0]-v2[0])*(v1[0]-v2[0])+(v1[1]-v2[1])*(v1[1]-v2[1])+(v1[2]-v2[2])*(v1[2]-v2[2]),0.5)
if v1[0]!=v2[0] and v1[1]!=v2[1]  and v1[2]!=v2[2] and cbprp1>cbprp2 and dist<0.0045:
between_smallest.append([v1,v2])
glClearColor(0.0, 0.0, 0.0, 1.0)
glViewport(0,0 ,640, 480)
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
gluPerspective(120.0,640.0/480.0,0.0000001,50.0)
glPointSize(5.0)
glColor3f(0.5, 0.5, 0.0)

def normalize(v):
    norm = np.linalg.norm(v)
    if norm == 0: 
       return v
    return v / norm

def mouse( x, y):
global Z
Z=(y-480.0/2.0)/600.0
Z=Z*Z
def plotpoints():
global triangles
global angle
global colorbprp
if 0==0 :
ViewD=Z
la=[0.0,0.0,0.0]
la[0]=ViewD*np.cos(angle/180.0*math.pi)
la[1]=0.0
la[2]=ViewD*np.sin(angle/180.0*math.pi);
normalize(la)
angle +=0.1
glClear(GL_COLOR_BUFFER_BIT)
glMatrixMode(GL_MODELVIEW) 
glLoadIdentity()
gluLookAt(la[0],la[1],la[2],0.0,0.0,0.0,0.0,1.0,0.0)
glColor3f(0.5, 0.5, 0.0)
glPointSize(5.0)
for v1,cbprp1,myrgb1 in triangles:
glColor3f(myrgb1[0]*255.0, myrgb1[1]*255.0,myrgb1[2]*255.0)
glPointSize(pow(cbprp1,2.0))
glBegin(GL_POINTS)
glVertex3f(v1[0], v1[1],v1[2])
glEnd()
glColor3f(0.5, 0.5, 0.5)
for v1,v2 in between_brightest:
glBegin(GL_LINES)
glVertex3f(v1[0],v1[1],v1[2])
glVertex3f(v2[0],v2[1],v2[2])
glEnd()
glColor3f(0.5, 0.0, 0.5)
for v1,v2 in between_smallest:
glBegin(GL_LINES)
glVertex3f(v1[0],v1[1],v1[2])
glVertex3f(v2[0],v2[1],v2[2])
glEnd()
glFlush()
glutPostRedisplay()



if __name__ == "__main__":


glutInit()
glutInitDisplayMode(GLUT_RGB |GLUT_MULTISAMPLE  |GLUT_DEPTH)
glutInitWindowSize(640, 480)
glutCreateWindow("PLot CSV")
glutDisplayFunc(plotpoints)
glutPassiveMotionFunc( mouse )
init()
glutMainLoop()