"""Defines a 'BlackBody' function that returns an r,g,b color tuple (range 0.0
   to 1.0 for each component) for a given temp. An optional gamma argument
   specifies the output gamma function.
"""

import math

h=6.626e-34   #Joule-Sec
c=2.9979e8    #m/sec
k=1.381e-23   #Joule/K

a=8*math.pi*h*c
b=h*c/k

def rho(l=1e-10, T=6000):   #return energy density, l=lambda in metres
  try:
    return (a/(l**5)) * (1/(math.exp(b/(l*T))-1))
  except:
    return 0


def BlackBody(temp=1000, gamma=0.8):
  "Takes temp in degrees, returns r,g,b in the range 0.0 to 1.0"
  #CIE Color Matching Functions (x_bar,y_bar,z_bar)
  #for wavelenghts in 5 nm increments from 380 nm to 780 nm.
  fColorMatch=[
         (0.0014, 0.0000, 0.0065),
         (0.0022, 0.0001, 0.0105),
         (0.0042, 0.0001, 0.0201),
         (0.0076, 0.0002, 0.0362),
         (0.0143, 0.0004, 0.0679),
         (0.0232, 0.0006, 0.1102),
         (0.0435, 0.0012, 0.2074),
         (0.0776, 0.0022, 0.3713),
         (0.1344, 0.0040, 0.6456),
         (0.2148, 0.0073, 1.0391),
         (0.2839, 0.0116, 1.3856),
         (0.3285, 0.0168, 1.6230),
         (0.3483, 0.0230, 1.7471),
         (0.3481, 0.0298, 1.7826),
         (0.3362, 0.0380, 1.7721),
         (0.3187, 0.0480, 1.7441),
         (0.2908, 0.0600, 1.6692),
         (0.2511, 0.0739, 1.5281),
         (0.1954, 0.0910, 1.2876),
         (0.1421, 0.1126, 1.0419),
         (0.0956, 0.1390, 0.8130),
         (0.0580, 0.1693, 0.6162),
         (0.0320, 0.2080, 0.4652),
         (0.0147, 0.2586, 0.3533),
         (0.0049, 0.3230, 0.2720),
         (0.0024, 0.4073, 0.2123),
         (0.0093, 0.5030, 0.1582),
         (0.0291, 0.6082, 0.1117),
         (0.0633, 0.7100, 0.0782),
         (0.1096, 0.7932, 0.0573),
         (0.1655, 0.8620, 0.0422),
         (0.2257, 0.9149, 0.0298),
         (0.2904, 0.9540, 0.0203),
         (0.3597, 0.9803, 0.0134),
         (0.4334, 0.9950, 0.0087),
         (0.5121, 1.0000, 0.0057),
         (0.5945, 0.9950, 0.0039),
         (0.6784, 0.9786, 0.0027),
         (0.7621, 0.9520, 0.0021),
         (0.8425, 0.9154, 0.0018),
         (0.9163, 0.8700, 0.0017),
         (0.9786, 0.8163, 0.0014),
         (1.0263, 0.7570, 0.0011),
         (1.0567, 0.6949, 0.0010),
         (1.0622, 0.6310, 0.0008),
         (1.0456, 0.5668, 0.0006),
         (1.0026, 0.5030, 0.0003),
         (0.9384, 0.4412, 0.0002),
         (0.8544, 0.3810, 0.0002),
         (0.7514, 0.3210, 0.0001),
         (0.6424, 0.2650, 0.0000),
         (0.5419, 0.2170, 0.0000),
         (0.4479, 0.1750, 0.0000),
         (0.3608, 0.1382, 0.0000),
         (0.2835, 0.1070, 0.0000),
         (0.2187, 0.0816, 0.0000),
         (0.1649, 0.0610, 0.0000),
         (0.1212, 0.0446, 0.0000),
         (0.0874, 0.0320, 0.0000),
         (0.0636, 0.0232, 0.0000),
         (0.0468, 0.0170, 0.0000),
         (0.0329, 0.0119, 0.0000),
         (0.0227, 0.0082, 0.0000),
         (0.0158, 0.0057, 0.0000),
         (0.0114, 0.0041, 0.0000),
         (0.0081, 0.0029, 0.0000),
         (0.0058, 0.0021, 0.0000),
         (0.0041, 0.0015, 0.0000),
         (0.0029, 0.0010, 0.0000),
         (0.0020, 0.0007, 0.0000),
         (0.0014, 0.0005, 0.0000),
         (0.0010, 0.0004, 0.0000),
         (0.0007, 0.0002, 0.0000),
         (0.0005, 0.0002, 0.0000),
         (0.0003, 0.0001, 0.0000),
         (0.0002, 0.0001, 0.0000),
         (0.0002, 0.0001, 0.0000),
         (0.0001, 0.0000, 0.0000),
         (0.0001, 0.0000, 0.0000),
         (0.0001, 0.0000, 0.0000),
         (0.0000, 0.0000, 0.0000)  ]

  XX=0.0
  YY=0.0
  ZZ=0.0    # initialize accumulators

  temp=temp/1.3   #Hack to make the color look more realistic

  # loop over wavelength bands
  # integration by trapezoid method
  nbands=len(fColorMatch)
  for band in range(nbands):
    weight=1.0;
    if ((band==0) or (band==nbands-1)):
      weight=0.5    #properly weight end points
    wavelength=380.0+band*5.0    # wavelength in nm
        #  generate a black body spectrum
    con=1240.0/8.617e-5;

    dis=( 3.74183e-16*(1.0/pow(wavelength,5))/
                     (math.exp(con/(wavelength*temp))-1.0) )
    #simple integration over bands
    XX=XX+weight*dis*fColorMatch[band][0];
    YY=YY+weight*dis*fColorMatch[band][1];
    ZZ=ZZ+weight*dis*fColorMatch[band][2];

  # re-normalize the color scale
  x=XX/max(XX,YY,ZZ);
  y=YY/max(XX,YY,ZZ);
  z=ZZ/max(XX,YY,ZZ);

  rr,gg,bb=xyz2rgb(x,y,z)
  
  maxv=max(1.0e-10, rr, gg, bb)
  
#  print "x=",x," y=",y," z=",z
#  print "XX=",XX," YY=",YY," ZZ=",ZZ
#  print "rr=",rr," gg=",gg," bb=",bb
  
  return (rr/maxv)**gamma, (gg/maxv)**gamma, (bb/maxv)**gamma
  
  
def xyz2rgb(xColor,yColor,zColor):
  xr=0.64
  yr=0.33
  zr=1.0-(xr+yr)
  xg=0.29
  yg=0.60
  zg=1.0-(xg+yg)
  xb=0.15
  yb=0.06
  zb=1.0-(xb+yb)
  xw=0.3127
  yw=0.3291
  zw=1.0-(xw+yw)

  denominator=(xr*yg-xg*yr)*zb+(xb*yr-xr*yb)*zg+(xg*yb-xb*yg)*zr
  
  red=((xColor*yg-xg*yColor)*zb+(xg*yb-xb*yg)*zColor
        +(xb*yColor-xColor*yb)*zg)/denominator;
  green=((xr*yColor-xColor*yr)*zb+(xb*yr-xr*yb)*zColor
          +(xColor*yb-xb*yColor)*zr)/denominator
  blue=((xr*yg-xg*yr)*zColor+(xColor*yr-xr*yColor)*zg
         +(xg*yColor-xColor*yg)*zr)/denominator

  if(red<0.0):
    red=0.0
  elif(red>1.0):
    red=1.0
  if(green<0.0):
    green=0.0
  elif(green>1.0):
    green=1.0
  if(blue<0.0):
    blue=0.0
  elif(blue>1.0):
    blue=1.0

  return red,green,blue