From 8c0af0e467b0ea4b689d45cc4d5142463c663352 Mon Sep 17 00:00:00 2001
From: hut <hut@hut.pm>
Date: Fri, 19 Apr 2024 03:39:12 +0200
Subject: [PATCH] add star radius calculation based on stefan-boltzman law

---
 doc/scripts/generate_starchart.py | 32 ++++++++++++++++++++++++-------
 1 file changed, 25 insertions(+), 7 deletions(-)

diff --git a/doc/scripts/generate_starchart.py b/doc/scripts/generate_starchart.py
index 55f0fc5..e16ec51 100755
--- a/doc/scripts/generate_starchart.py
+++ b/doc/scripts/generate_starchart.py
@@ -34,7 +34,7 @@ MAX_APPARENT_MAGNITUDE = 7.0
 SOL_ABSMAG = 4.83
 SOL_RADIUS = 696.3e6
 SOL_LUMINOSITY = 3.828e26
-STEFAN_BOLZMANN_CONSTANT = 5.670374419e-8
+STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8
 
 print("// This file was autogenerated by \"genrate_starchart.py\" using data from the")
 print("// HYG database: https://github.com/astronexus/HYG-Database/tree/main/hyg")
@@ -43,7 +43,7 @@ print("// Each star's values: (x, y, z, magnitude, absolute magnitude, color ind
 print("[")
 
 
-def render(x, y, z, ra, dec, mag, absmag, ci, dist, name):
+def render(i, x, y, z, ra, dec, mag, absmag, ci, dist, name):
     # Takes ra/deg in degrees
     ra = float(ra)
     dec = float(dec)
@@ -53,6 +53,8 @@ def render(x, y, z, ra, dec, mag, absmag, ci, dist, name):
     if name == 'Sol':
         return
 
+    #radius = star_radius(float(ci), float(absmag))
+
 #    distance = 1.0
 #    ra_radians = math.radians(ra * 15)  # ra is in [0, 24], multiplying by 15 gives degrees in [0, 360]
 #    dec_radians = math.radians(dec)
@@ -73,17 +75,33 @@ def clean_name(string):
     return "".join([c for c in string if c.isalnum() or c in ' -_'])
 
 
+def star_radius(color_index, absolute_magnitude):
+    # Convert color index to temperature (using Ballesteros' formula approximation for B-V)
+    temp = 4600 * ((1 / (0.92 * color_index + 1.7)) + (1 / (0.92 * color_index + 0.62)))
+
+    # Convert absolute magnitude to luminosity (in solar luminosities)
+    lum = 10 ** (0.4 * (SOL_ABSMAG - absolute_magnitude))
+
+    # Calculate the radius using Stefan-Boltzmann law (output should be in meters)
+    radius = math.sqrt(lum * SOL_LUMINOSITY / (4 * math.pi * STEFAN_BOLTZMANN_CONSTANT * temp ** 4))  # Luminosity in watts for the Sun ~ 3.828e26
+
+    # Convert radius from meters to solar radii (1 solar radius ≈ 6.96e8 meters)
+    radius_solar = radius / SOL_RADIUS
+
+    return radius_solar
+
+
 total = 0
 count = 0
 count_extra = 0
 entries = []
 with open(PATH, "r", encoding="utf-8") as f:
-    for entry in csv.DictReader(f):
-        entries.append(entry)
+    for index, entry in enumerate(csv.DictReader(f)):
+        entries.append((index, entry))
 
-entries.sort(key=lambda entry: float(entry['mag']))
+entries.sort(key=lambda entry: float(entry[1]['mag']))
 
-for entry in entries:
+for index, entry in entries:
     total += 1
     ra = entry['ra']
     dec = entry['dec']
@@ -104,7 +122,7 @@ for entry in entries:
         continue
     if float(mag) > MAX_APPARENT_MAGNITUDE:
         continue
-    render(x, y, z, ra, dec, mag, absmag, ci, dist, name)
+    render(index, x, y, z, ra, dec, mag, absmag, ci, dist, name)
     count += 1
 #for entry in CUSTOM_ENTRIES:
 #    render(entry[0], entry[1], entry[2], entry[3], entry[4])