FORTRAN Code

c
c      CHROMATCITY DIAGRAM by Dan Bruton (astro@tamu.edu)
c   
c      This program can be found at 
c      http://www.physics.sfasu.edu/astro/color.html
c      and was last updated on February 20, 1996.
c
c      This program will create a ppm (portable pixmap) 
c      image of an approximate chromaticity diagram using 
c      equations from the Color Equations FAQ at
c      ftp://ftp.wmin.ac.uk/pub/itrg/coloureq.txt
c      or the Color Space FAQ.
c      The unix program ppmtogif is used to convert 
c      the ppm image to a gif.  NetPBM's ppmtogif can be found at
c      ftp://ftp.cs.ubc.ca/ftp/archive/netpbm/
c
       IMPLICIT REAL*8 (a-h,o-z)
       REAL*8 CV(500,500,3)
c      
c      Width, height and color depth of the ppm image
c
       M=500
       N=M
       L=255
c
c      Gamma Correction
c
       GAM=1.0
c
c      Luminance=YY
c
       YY=1.0
c
c      Initialize Scaling Factor
c
       RMAX=0.0
c
c      Calculate RGB Values for x and y coordinates
c
       do J=1,N
        do I=1,M
            X=REAL(1.*I/M)
            Y=REAL(1.*(N-J)/N)+0.00001
            XX=X*YY/Y
            ZZ=(1.-X-Y)*YY/Y
            R=(1.910*XX)-(0.532*YY)-(0.288*ZZ)
            G=-(0.985*XX)+(1.999*YY)-(0.028*ZZ)
            B=(0.058*XX)-(0.118*YY)+(0.898*ZZ)
            SSS=(R+G+B)/3.0
            R=R/SSS
            G=G/SSS
            B=B/SSS
            IF ((R.LT.0.).OR.(G.LT.0.).OR.(B.LT.0.)) then
               R=0.
               G=0.
               B=0.
            ENDIF
            R=R**GAM
            G=G**GAM
            B=B**GAM
            IF (R.GT.RMAX) then
              RMAX=R
            ENDIF
            IF (G.GT.RMAX) then
              RMAX=G
            ENDIF
            IF (B.GT.RMAX) then
              RMAX=B
            ENDIF
            CV(I,J,1)=R
            CV(I,J,2)=G
            CV(I,J,3)=B
         enddo
       enddo
c
c      Write to PPM File
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,*) L
       do J=1,N
        do I=1,M
100      FORMAT(I4,2X,I4,2X,I4,2X)
         IR=INT(255.*CV(I,J,1)/RMAX)
         IG=INT(255.*CV(I,J,2)/RMAX)
         IB=INT(255.*CV(I,J,3)/RMAX)
         write(20,*) IR,IG,IB
        enddo
       enddo
       stop
       end