File:Colours and temperatures of main sequence stars 1.svg

1. 2. star spectral type and temperature interpolate
   3. view star colours in filled circles
   4. needs predefined table of staller tspectral types, masses, temperatures ...
   5. python 3 code
   7. 13.2.2023 v 0000.0002

import numpy as np

1. 1. import scipy as sp

import matplotlib.pyplot as plt import math import pandas as pd

from sklearn import metrics from sklearn.linear_model import LinearRegression from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import r2_score from pygam import GAM, s, te from pygam import LinearGAM from pygam import PoissonGAM

1. 1. from pygam import LogisticGAM

from scipy import interpolate

from scipy.constants import h, c, k from matplotlib.patches import Circle

from colour_system import cs_hdtv cs = cs_hdtv

def planck(lam, T):

   lam_m = lam / 1.e9
   fac = h*c/lam_m/k/T
   B = 2*h*c**2/lam_m**5 / (np.exp(fac) - 1)
   return B

def linear_log_predict1(x,y): logx=np.log10(x).reshape(-1, 1) logy=np.log10(y).reshape(-1, 1) logmodel= LinearRegression().fit(logx, logy) logk=logmodel.coef_ logb=logmodel.intercept_ #print(logmodel) logpre=logmodel.predict(logx) pre=np.power(10, logpre) return(pre, logk, logb)

def linear_log_predict2(x,y, x2): logx=np.log10(x).reshape(-1, 1) logy=np.log10(y).reshape(-1, 1) logx2=np.log10(x2).reshape(-1, 1) logmodel= LinearRegression().fit(logx, logy) logk=logmodel.coef_ logb=logmodel.intercept_ #print(loglums_model) logpre=logmodel.predict(logx2) pre=np.power(10, logpre) return(pre)

def linear_log_predict_point(x,y, px): logx=np.log10(x).reshape(-1, 1) logy=np.log10(y).reshape(-1, 1) logmodel= LinearRegression().fit(logx, logy) logk=logmodel.coef_ logb=logmodel.intercept_ logpx=math.log10(px) logpy=logb+logpx*logk py=math.pow(10,logpy) return(py)

def gamma_log_predict2(x,y, x2): logx=np.log10(x).reshape(-1, 1) logy=np.log10(y).reshape(-1, 1) logx2=np.log10(x2).reshape(-1, 1) gam = LinearGAM(s(0, n_splines=10)).fit(logx,logy) #gam = PoissonGAM(s(0, n_splines=10)).fit(logx,logy) logpre = gam.predict(logx2) pre=np.power(10, logpre) return(pre)

def gamma_predict3(x,y, x2): ex=x.reshape(-1, 1) ey=y.reshape(-1, 1) ex2=x2.reshape(-1, 1) gam = LogisticGAM(s(0, n_splines=10)).fit(ex,ey) pre = gam.predict(ex2) return(pre)

def rf_log_predict2(x,y, x2): logx=np.log10(x).reshape(-1, 1) logy=np.log10(y).reshape(-1, 1) logx2=np.log10(x2).reshape(-1, 1) regressor = RandomForestRegressor(n_estimators=10000, random_state=0) regressor.fit(logx,logy) logpre = regressor.predict(logx2) pre=np.power(10, logpre) return(pre)

def interpolate1(x,y, x2): f = interpolate.interp1d(x, y) y2 = f(x2) return(y2)

def interpolate_point(x,y, px): f = interpolate.interp1d(x, y) y2 = f(px) return(y2)

def interpolate2(x,y, x2): logx=np.log10(x).reshape(-1, 1) logy=np.log10(y).reshape(-1, 1) logx2=np.log10(x2).reshape(-1, 1) f = interpolate.interp1d(logx, logy) y2 = f(logx2) pre=np.power(10,y2) return(y2)

def speknum(spt): spe=spt spe=spe.replace('O', '1') spe=spe.replace('B', '2') spe=spe.replace('A', '3') spe=spe.replace('F', '4') spe=spe.replace('G', '5') spe=spe.replace('K', '6') spe=spe.replace('M', '7') spe=spe.replace('L', '8') return(int(spe))

1. 1. load main sequence star data, edited from wiki 2023
   2. https://en.wikipedia.org/wiki/Main_sequence

arrname1="mainseq1.txt"

starr1 = np.loadtxt(arrname1,skiprows=1,dtype="float",delimiter=",",usecols =(1,2,3,4)) print(starr1)

masses1=starr1[:,1] lums1=starr1[:,2] temps1=starr1[:,3] radiuses1=starr1[:,0]

print (masses1) print (lums1) print (temps1)

df00 = pd.read_csv(arrname1, sep=",")

1. print(df00)

spek1=df00['spectral']

1. print(spek1)

classes="OBAFGKML"

spek2=spek1.tolist()

1. print(spek2)

spek3 = [sub.replace('O', '1') for sub in spek2] spek3 = [sub.replace('B', '2') for sub in spek3] spek3 = [sub.replace('A', '3') for sub in spek3] spek3 = [sub.replace('F', '4') for sub in spek3] spek3 = [sub.replace('G', '5') for sub in spek3] spek3 = [sub.replace('K', '6') for sub in spek3] spek3 = [sub.replace('M', '7') for sub in spek3] spek3 = [sub.replace('L', '8') for sub in spek3]

1. print(spek3)

spek4=np.array(spek3).astype(float)

1. print(spek4)

speg1=np.arange(12,82,1)

1. print(speg1)

speg2=speg1.astype(int)

1. print(speg2)

speg3=speg2.tolist()

1. print(speg3)

speg4 = [str(x) for x in speg3]

1. print(speg4)

speg5=[]

first1=0

for sx in speg4:

sy=sx print(sy) chari1=sy[0] if(sy[0]=='1'): sy=sy.replace(sy[0],"O",1) if(sy[0]=='2'): sy=sy.replace(sy[0],"B",1) if(sy[0]=='3'): sy=sy.replace(sy[0],"A",1) if(sy[0]=='4'): sy=sy.replace(sy[0],"F",1) if(sy[0]=='5'): sy=sy.replace(sy[0],"G",1) if(sy[0]=='6'): sy=sy.replace(sy[0],"K",1) if(sy[0]=='7'): sy=sy.replace(sy[0],"M",1) if(sy[0]=='8'): sy=sy.replace(sy[0],"L",1) speg5.append(sy)

1. print("speg5")

1. print(spek4)
2. print(speg1)
3. print(speg5)

spektr1=spek4 spektr2=speg1

sally1=spek2 sally2=speg5

print(sally1) print(sally2)

print(spektr1) print(spektr2)

masses2=interpolate1(spektr1, masses1, spektr2) radiuses2=interpolate1(spektr1, radiuses1, spektr2) lums2=interpolate1(spektr1, lums1, spektr2) temps2=interpolate1(spektr1, temps1, spektr2)

spt1="G2"

spnum1=speknum(spt1)

teff1=interpolate_point(spektr1,temps1, spnum1)

print(spt1, spnum1,teff1)

1. lum1=interpolate_point(masses1, lums1, m1)

1. plt.xlim(0.1, 1.5)
2. plt.ylim(0.0, 3)

font = {'family' : 'DejaVu Sans',

       'weight' : 'normal',
       'size'   : 18}

1. plt.rc('font', **font)

1. plt.title ("Main sequence spectral class - surface temperature", size=20)
2. plt.plot(spektr1, temps1, "o", lw=6)

1. plt.plot(spektr2, temps2, "-", color="r", lw=10)

1. plt.grid(color='b', linestyle=':', linewidth=1)

1. plt.xticks(spektr1, sally1)

1. plt.plot(spnum1, teff1, "x")

1. plt.show()

fig, ax = plt.subplots()

1. The grid of visible wavelengths corresponding to the grid of colour-matching
2. functions used by the ColourSystem instance.

lam = np.arange(380., 781., 5)

for i in range(len(spektr1)):

   # T = 500 to 12000 K
   #T = 500*i + 500
   T=int(temps1[i])
   sp1=sally1[i]

   # Calculate the black body spectrum and the HTML hex RGB colour string
   # it looks like
   spec = planck(lam, T)
   html_rgb = cs.spec_to_rgb(spec, out_fmt='html')

   # Place and label a circle with the colour of a black body at temperature T
   x, y = i % 5, -(i // 5)
   circle = Circle(xy=(x, y*1.2), radius=0.4, fc=html_rgb)
   ax.add_patch(circle)
   ax.annotate('{:4d} K'.format(T), xy=(x, y*1.2-0.6), va='center',
               ha='center', color=html_rgb, fontsize=15)
   ax.annotate('{:4s}'.format(sp1), xy=(x, y*1.2), va='center',
               ha='center', color='k', fontsize=18)   
               

1. Set the limits and background colour; remove the ticks

ax.set_xlim(-0.5,5.0) ax.set_ylim(-4.35, 0.5) ax.set_xticks([]) ax.set_yticks([]) ax.set_facecolor('k')

1. Make sure our circles are circular!

ax.set_aspect("equal") plt.title("Colours and temperatures of main sequence stars", size=22) plt.show()

1. colour_system.py

import numpy as np

def xyz_from_xy(x, y):

   """Return the vector (x, y, 1-x-y)."""
   return np.array((x, y, 1-x-y))

class ColourSystem:

   """A class representing a colour system.

   A colour system defined by the CIE x, y and z=1-x-y coordinates of
   its three primary illuminants and its "white point".

   TODO: Implement gamma correction

   """

   # The CIE colour matching function for 380 - 780 nm in 5 nm intervals
   cmf = np.loadtxt('cie-cmf.txt', usecols=(1,2,3))

   def __init__(self, red, green, blue, white):
       """Initialise the ColourSystem object.

       Pass vectors (ie NumPy arrays of shape (3,)) for each of the
       red, green, blue  chromaticities and the white illuminant
       defining the colour system.

       """

       # Chromaticities
       self.red, self.green, self.blue = red, green, blue
       self.white = white
       # The chromaticity matrix (rgb -> xyz) and its inverse
       self.M = np.vstack((self.red, self.green, self.blue)).T 
       self.MI = np.linalg.inv(self.M)
       # White scaling array
       self.wscale = self.MI.dot(self.white)
       # xyz -> rgb transformation matrix
       self.T = self.MI / self.wscale[:, np.newaxis]

   def xyz_to_rgb(self, xyz, out_fmt=None):
       """Transform from xyz to rgb representation of colour.

       The output rgb components are normalized on their maximum
       value. If xyz is out the rgb gamut, it is desaturated until it
       comes into gamut.

       By default, fractional rgb components are returned; if
       out_fmt='html', the HTML hex string '#rrggbb' is returned.

       """

       rgb = self.T.dot(xyz)
       if np.any(rgb < 0):
           # We're not in the RGB gamut: approximate by desaturating
           w = - np.min(rgb)
           rgb += w
       if not np.all(rgb==0):
           # Normalize the rgb vector
           rgb /= np.max(rgb)

       if out_fmt == 'html':
           return self.rgb_to_hex(rgb)
       return rgb

   def rgb_to_hex(self, rgb):
       """Convert from fractional rgb values to HTML-style hex string."""

       hex_rgb = (255 * rgb).astype(int)
       return '#{:02x}{:02x}{:02x}'.format(*hex_rgb)

   def spec_to_xyz(self, spec):
       """Convert a spectrum to an xyz point.

       The spectrum must be on the same grid of points as the colour-matching
       function, self.cmf: 380-780 nm in 5 nm steps.

       """

       XYZ = np.sum(spec[:, np.newaxis] * self.cmf, axis=0)
       den = np.sum(XYZ)
       if den == 0.:
           return XYZ
       return XYZ / den

   def spec_to_rgb(self, spec, out_fmt=None):
       """Convert a spectrum to an rgb value."""

       xyz = self.spec_to_xyz(spec)
       return self.xyz_to_rgb(xyz, out_fmt)

illuminant_D65 = xyz_from_xy(0.3127, 0.3291) cs_hdtv = ColourSystem(red=xyz_from_xy(0.67, 0.33),

                      green=xyz_from_xy(0.21, 0.71),
                      blue=xyz_from_xy(0.15, 0.06),
                      white=illuminant_D65)

cs_smpte = ColourSystem(red=xyz_from_xy(0.63, 0.34),

                       green=xyz_from_xy(0.31, 0.595),
                       blue=xyz_from_xy(0.155, 0.070),
                       white=illuminant_D65)

cs_srgb = ColourSystem(red=xyz_from_xy(0.64, 0.33),

                      green=xyz_from_xy(0.30, 0.60),
                      blue=xyz_from_xy(0.15, 0.06),
                      white=illuminant_D65)