Welcome to Fractal Forums

from math import exp, floor, ceil
import random
import numpy as np
from matplotlib import pyplot as plt


C = .70351417223810256  # prefactor such that peak is 1
K = 4.79924  # prefactor if units are in PHz and 10^4 K
# desired range for v is .4 to .789
# desired range for T is 7.98e-2 to 2.
# typical range would be .5 to .7, commonly .65


def Planck(v, T):  # the distribution of spectral frequencies for a given temperature
    KTv = K * v / T
    return C * KTv * KTv * KTv * v * v / (exp(KTv) - 1)


def MC(T):  # generates samples from the distribution using Markov Chain Monte Carlo method
    v0 = random.uniform(.4, .789)
    I0 = Planck(v0, T)
    while True:
        v = random.uniform(.4, .789)
        I = Planck(v, T)
        a = I/I0
        if a >= 1:
            v0 = v
            I0 = I
        else:
            p = random.random()
            if a <= p:
                v0 = v
                I0 = I
        yield v0


# a bunch of utility functions to get to proper units


def to_eV(v):
    return v * 4.135668


def to_nm(v):
    return 299.8/v


def from_K(T):
    return T * 1.e-4


def from_THz(v):
    return v * 1.e-3

# example histogram for some temperature

temp = 6500  # in Kelvin

data = [to_nm(v) for v, _ in zip(MC(from_K(temp)), range(100000))]  # in nanometers, from 100k samples.

binbot = ceil(min(data))
bintop = floor(max(data))
bincount = int(bintop-binbot)
bins = np.linspace(binbot,
                   bintop,
                   bincount)

plt.xlim([min(data), max(data)])

plt.hist(data, bins=bins, alpha=0.5)
plt.title('spectral distribution of ' + repr(temp) + 'K light source')
plt.xlabel('nm')
plt.ylabel('count')

plt.show()

from math import exp, floor, ceil
import random
import numpy as np
from matplotlib import pyplot as plt
import csv


# a bunch of utility functions to get to proper units


def to_eV(v):
    return v * 4.135668


def to_nm(v):
    return 299.8/v

def from_nm(v):
    return 299.8/v

def from_K(T):
    return T * 1.e-4


def from_THz(v):
    return v * 1.e-3

C = 77.766  # prefactor such that peak is 1
K = 4.79924  # prefactor if units are in PHz and 10^4 K
# desired range for v is .4 to .789
# desired range for T is 7.98e-2 to 2.
# typical range would be .5 to .7, commonly .65


def Planck(v, T):  # the distribution of spectral frequencies for a given temperature
    vT = v / T
    return C * vT * vT * vT / (exp(K * vT) - 1)

topPHz = from_nm(390)
botPHz = from_nm(830)


def MC(T):  # generates samples from the distribution using Markov Chain Monte Carlo method
    v0 = random.uniform(botPHz, topPHz)
    I0 = Planck(v0, T)
    while True:
        v = random.uniform(botPHz, topPHz)
        I = Planck(v, T)
        a = I/I0
        if a >= 1:
            v0 = v
            I0 = I
        else:
            p = random.random()
            if p < a:
                v0 = v
                I0 = I
        yield v0

# example histogram for some temperature

samples = 1000000
temp = 5777  # in Kelvin

data = [to_nm(v) for v, _ in zip(MC(from_K(temp)), range(samples))]  # in nanometers, from 100k samples.

specxyz = []
with open('lin2012xyz10e_fine_7sf.csv', newline='') as inputfile:
    for row in csv.reader(inputfile):
        specxyz.append(row)

col = (0, 0, 0)

r, g, b = 0., 0., 0.
for d in data:
    x = d*10 - 3900
    i = floor(x)
    y = x - i
    vals = zip(specxyz[i], specxyz[i + 1])
    ivals = []
    for v in vals:
        ivals.append(float(v[0])*(1-y) + float(v[1])*y)
    r += ivals[1]
    g += ivals[2]
    b += ivals[3]

t = r + g + b
r /= g

b /= g
g = 1

def xyztorgb(x, y, z):
    r = 3.2404542 * x - 1.5371385 * y - 0.4985314 * z
    g = -0.969266 * x + 1.8760108 * y + 1.8760108 * z
    b = 0.0556434 * x - 0.2040259 * y + 1.0572252 * z
    return r, g, b

print('XYZ')
print(r)
print(g)
print(b)
print()
r, g, b = xyztorgb(r, g, b)

r = int(min(1, max(0, r)) * 255)
g = int(min(1, max(0, g)) * 255)
b = int(min(1, max(0, b)) * 255)

print('RGB')
print(r)
print(g)
print(b)
print()

binbot = ceil(min(data))
bintop = floor(max(data))
bincount = int(bintop-binbot)
bins = np.linspace(binbot,
                   bintop,
                   bincount)

plt.xlim([min(data), max(data)])

plt.hist(data, bins=bins, alpha=0.5)
plt.title('spectral distribution of ' + repr(temp) + 'K light source')
plt.xlabel('nm')
plt.ylabel('count')

plt.show()

XYZ
0.9739410242962122
1
0.9990353465096308

RGB
255
255
231