BASIC Code

This approximation is based on data provided by Mitchell Charity.
Function StarColor(T, RR, GG, BB)
   'T = Temperature, RR=Red, GG=Green, BB=Blue
   RR = 56100000 * T ^ (-3 / 2) + 148
   GG = 100.04 * Log(T) - 623.6
   If T > 6500 Then GG = 35200000 * T ^ (-3 / 2) + 184
   BB = 194.18 * Log(T) - 1448.6
   If RR > 255 Then RR = 255
   If GG > 255 Then GG = 255
   If BB > 255 Then BB = 255
   If RR < 0 Then RR = 0
   If GG < 0 Then GG = 0
   If BB < 0 Then BB = 0
End Function

FORTRAN Code

c
c      RGB VALUES FOR HOT OBJECTS   by Dan Bruton (astro@sfasu.edu)
c
c      This program can be found at 
c      http://www.physics.sfasu.edu/astro/color.html
c      and was last updated on March 19, 1996.
c
c      Reference Book
c      Color Science : Concepts and Methods, Quantitative Data and Formula 
c                      by Gunter Wyszecki, W. S. Stiles 
c                      John Wiley & Sons; ISBN: 0471021067 
c
c      This program will calculate the RGB values for a given
c      energy distribution as a function of wavelength of light.
c      A black body is used and an example.
c
c      NetPBM's ppmtogif can be used to convert the ppm image generated
c      to a gif.  The red, green and blue values (RGB) are
c      assumed to vary linearly with wavelength (for GAM=1).  
c      NetPBM Software: ftp://ftp.cs.ubc.ca/ftp/archive/netpbm/
c
       IMPLICIT REAL*8 (a-h,o-z)
       DIMENSION ICV(600,100,3)
c
c      Chromaticity Coordinates for Red, Green, Blue and White
c
       XR=0.64
       YR=0.33
       XG=0.29
       YG=0.60
       XB=0.15
       YB=0.06
       XW=0.3127
       YW=0.3291
       ZR=1.D0-(XR+YR)
       ZG=1.D0-(XG+YG)
       ZB=1.D0-(XB+YB)
       ZW=1.D0-(XW+YW)
c
c      IMAGE INFO - WIDTH, HEIGHT, DEPTH, GAMMA
c
       M=600
       N=50
       MX=255
       GAM=0.8
c
c      FIND COLOR VALUE, ICV, OF EACH PIXEL 
c
       RMX=0.
       DO I=1,M
        T=1000.0 + REAL((I-1) * 10000./M)
        ITYPE=1
        CALL BLKBOD(T,X,Y,Z)
        CALL XYZTORGB(XR,YR,ZR,XG,YG,ZG,XB,YB,ZB,X,Y,Z,R,G,B) 
        DO J=1,N
          IF (J.GT.40) CALL MARK(T,ITYPE,R,G,B)
          RMAX=1.D-10
          IF (R.GT.RMAX) RMAX=R
          IF (G.GT.RMAX) RMAX=G
          IF (B.GT.RMAX) RMAX=B
         ICV(I,J,1)=INT(REAL(MX)*(R/RMAX)**GAM)
         ICV(I,J,2)=INT(REAL(MX)*(G/RMAX)**GAM)
         ICV(I,J,3)=INT(REAL(MX)*(B/RMAX)**GAM)
        ENDDO
        IF (DMAX1(R,G,B).GT.RMX) THEN 
          RMX=DMAX1(R,G,B)
        ENDIF
       ENDDO
c
c*****************************************************************************
c
c      WRITE OUTPUT TO PPM FILE
c
c      Output color values are normalized to color depth, MX.
c
       OPEN(UNIT=20,FILE='temp.ppm',STATUS='UNKNOWN')
1      FORMAT (A10)
       WRITE(20,1) 'P3        '
       WRITE(20,1) '# temp.ppm'
       WRITE(20,*) M,N
       WRITE(20,*) MX
       DO J=1,N
        DO I=1,M
         WRITE(20,*) ICV(I,J,1),ICV(I,J,2),ICV(I,J,3)
        ENDDO
       ENDDO
       STOP
       END 
c
c ****************************************************************************
c
c      XYZ VALUES FROM TEMPERATURE OF OBJECT
c
c      A black body approximation is used where the temperature,
c      T, is given in Kelvin.  The XYZ values are determined by
c      "integrating" the product of the wavelength distribution of
c      energy and the XYZ functions for a uniform source.  
c
       SUBROUTINE BLKBOD(T,X,Y,Z)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION  WXYZ(3,81)
c
c      CIE Color Matching Functions (x_bar,y_bar,z_bar)
c      for wavelenghts in 5 nm increments from 380 nm to 780 nm.
c
       DATA ((WXYZ(I,J),I=1,3),J=1,81)/
     *  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/
c
       XX=0.D0
       YY=0.D0
       ZZ=0.D0
       DO K=1,81
         WL=380.D0+REAL(K-1)*5.D0
         CON=1240./8.617D-5
         DIS=3.74183D-16*(1./WL**5)/(DEXP(CON/(WL*T))-1.)
         XX=XX+DIS*WXYZ(1,K)
         YY=YY+DIS*WXYZ(2,K)
         ZZ=ZZ+DIS*WXYZ(3,K)
       ENDDO
       X=XX/DMAX1(XX,YY,ZZ) 
       Y=YY/DMAX1(XX,YY,ZZ)
       Z=ZZ/DMAX1(XX,YY,ZZ)
       RETURN
       END
c
c ****************************************************************************
c
c      PLACE MARKERS ON IMAGE
c
       SUBROUTINE MARK(T,ITYPE,R,G,B)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c
c      ITYPE=1 - PLAIN IMAGE
c      ITYPE=2 - MARK IMAGE AT 1000 K INTEVALS
c
         IF (ITYPE.EQ.2) THEN
            DO K=1000,10000,1000
              IF (ABS(INT(T)-K).LE.1) THEN
               R=0.
               G=0.
               B=0.
              ENDIF
            ENDDO
         ENDIF
      RETURN
      END
c
c*********************************************************************
c
      SUBROUTINE XYZTORGB(xr,yr,zr,xg,yg,zg,xb,yb,zb,xc,yc,zc,R,G,B)
      IMPLICIT REAL*8 (a-h,o-z)
       r=(-xg*yc*zb+xc*yg*zb+xg*yb*zc-xb*yg*zc-xc*yb*zg+xb*yc*zg)/
     *   (+xr*yg*zb-xg*yr*zb-xr*yb*zg+xb*yr*zg+xg*yb*zr-xb*yg*zr)
       g=(+xr*yc*zb-xc*yr*zb-xr*yb*zc+xb*yr*zc+xc*yb*zr-xb*yc*zr)/
     *   (+xr*yg*zb-xg*yr*zb-xr*yb*zg+xb*yr*zg+xg*yb*zr-xb*yg*zr)
       b=(+xr*yg*zc-xg*yr*zc-xr*yc*zg+xc*yr*zg+xg*yc*zr-xc*yg*zr)/
     *   (+xr*yg*zb-xg*yr*zb-xr*yb*zg+xb*yr*zg+xg*yb*zr-xb*yg*zr)
       IF (R.LT.0.) R=0.
       IF (G.LT.0.) G=0.
       IF (B.LT.0.) B=0.
       IF (R.GT.1.) R=1.
       IF (G.GT.1.) G=1.
       IF (B.GT.1.) B=1.
      RETURN
      END
c *************************************************************************