c    Colours of a thin film in reflected light.
c    Illumination: blackbody radiation, colour temperature Tabs
c
       implicit real*8 (a-h,o-z)
       dimension iskalf(1000),icv(3,1000),Yuarray(3)
       real*4 colcomp(3,1000)
       common XuR,YuR,ZuR,XuG,YuG,ZuG,XuB,YuB,ZuB,
     1        RuX,GuX,BuX,RuY,GuY,BuY,RuZ,GuZ,BuZ,
     2        pi
c
       pi=3.141592653589d0
 	 Tabs=6504
       nhor=1000

c      sRGB primaries
       xr=0.64
       yr=0.33
       xg=0.30
       yg=0.60
       xb=0.15
       yb=0.06

c      white D65
       xw=0.3127
       yw=0.3290

c
       Gamma=2.2
       zr=1.d0-(xr+yr)
       zg=1.d0-(xg+yg)
       zb=1.d0-(xb+yb)

       zw=1.d0-(xw+yw)

!   Calibrate primaries to white
!   (As FORTRAN does not discriminate upper and lower case,
!   instead of X rather Xu is written, u like uppercase)
       XuW=xw/yw
       YuW=1.0d0
       ZuW=zw/yw
       call SSLE(xr,xg,xb,yr,yg,yb,zr,zg,zb,XuW,YuW,ZuW,
     *     u1,u2,u3)
       XuR=xr*u1
       YuR=yr*u1
       ZuR=zr*u1
       XuG=xg*u2 
       YuG=yg*u2
       ZuG=zg*u2
       XuB=xb*u3
       YuB=yb*u3
       ZuB=zb*u3

!  the brightnesses of the primaries stored in array Yuarray
      Yuarray(1)=YuR
      Yuarray(2)=YuG
      Yuarray(3)=YuB

!  Determination of RGB-coordinates of the primaries XYZ.
!  Solving the system of linear equations with unit matrix
!  on the right is matrix inversion.

c      superpose X from R, G, B 
       call SSLE(XuR,XuG,XuB,YuR,YuG,YuB,ZuR,ZuG,ZuB,
     *     1.d0,0.d0,0.d0,RuX,GuX,BuX)
c      superpose Y from R, G, B 
       call SSLE(XuR,XuG,XuB,YuR,YuG,YuB,ZuR,ZuG,ZuB,
     *     0.d0,1.d0,0.d0,RuY,GuY,BuY)
c      superpose Z from R, G, B 
       call SSLE(XuR,XuG,XuB,YuR,YuG,YuB,ZuR,ZuG,ZuB,
     *     0.d0,0.d0,1.d0,RuZ,GuZ,BuZ)

c	Transformation matrix is now ready.

c	To go from XYZ to RGB one computes
c     RuColor = RuX*XuColor + RuY*YuColor + RuZ*ZuColor
c	and analogously for G and B

c----------------------------------------------------------
                             
      do ic=1,1000
         ds=0.005d0*ic		
!  Lengths are measured in mu-m (micrometers)
!  ds is the optical path difference of the 
!  two interfering rays

!  compute XYZ
         call lambInt(XF,YF,ZF,ds,Tabs)
                              
!  convert XYZ to RGB
         RuF=RuX*XF+RuY*YF+RuZ*ZF
         GuF=GuX*XF+GuY*YF+GuZ*ZF
         BuF=BuX*XF+BuY*YF+BuZ*ZF
           colcomp(1,ic)=RuF
           colcomp(2,ic)=GuF
           colcomp(3,ic)=BuF
      enddo   ! ic    

!    Scale to make the maximum value of R, G or B equal to 1,
!    Gamma-correction and conversion to 8-bit-format

        themaximum=0.d0
        do j=1,nhor
          do k=1,3
           if (colcomp(k,j) .gt. themaximum) 
     *        themaximum=colcomp(k,j) 
          enddo
        enddo

        ammaG=1.0d0/2.4d0
        do j=1,nhor
        Yu=(colcomp(1,j)*YuR+colcomp(2,j)*YuG+colcomp(3,j)*YuB)
     *     /themaximum      !  brightness
        do i=1,3
         colcomp(i,j)=colcomp(i,j)/themaximum
         if (colcomp(i,j).lt.0.0) then
            Corrf=Yu/(Yu-colcomp(i,j)*Yuarray(i))
            do k=1,3
              colcomp(k,j)=colcomp(k,j)*Corrf   ! readjust brightness
            enddo
            colcomp(i,j)=0.0  
         endif
         if (colcomp(i,j).le.0.00304) then
           colcomp(i,j)=colcomp(i,j)*12.92
           else
           colcomp(i,j)=1.055*colcomp(i,j)**ammaG-0.055
         endif
         icv(i,j)=int(255.0*colcomp(i,j))
        enddo
        enddo

!   Produce a ppm-image
!   http://netpbm.sourceforge.net/doc/ppm.html
       ndepth=255
       nvert=80
       open(unit=7,file='thinfilm.ppm',status='unknown')
       write(7,*) 'P3             '
       write(7,*) '# thinfilm.ppm  '
       write(7,*) nhor,nvert
       write(7,*) ndepth
       do j=1,nvert-10
        do i=1,nhor
         write(7,*) icv(1,i),icv(2,i),icv(3,i)
        enddo
       enddo
!   make scale-ticks in 500-nm-intervals
      do i=1,nhor
       iskalf(i)=0
      enddo
      do i=1,10
       iskalf(100*i-1)=255
       iskalf(100*i)=255
       if (i.lt.10) iskalf(100*i+1)=255
      enddo
      do j=nvert-9,nvert
       do i=1,nhor
        write(7,*) iskalf(i),iskalf(i),iskalf(i)
       enddo 
      enddo
      stop
      end 
c----------------------------------------------------------
      function Refl(s,alambd)
!   Reflection at thin film, not normalized
      implicit real*8 (a-h,o-z)
      common XuR,YuR,ZuR,XuG,YuG,ZuG,XuB,YuB,ZuB,
     1        RuX,GuX,BuX,RuY,GuY,BuY,RuZ,GuZ,BuZ,
     2        pi
      beta=0.95d0
      Refl=1.d0+beta*beta-2.d0*beta*cos(s*2.d0*pi/alambd)
      end
c----------------------------------------------------------
        subroutine lambInt(Xw,Yw,Zw,ds,Tabs)
        implicit real*8 (a-h,o-z)
        real*4 xbar(95),ybar(95),zbar(95)
        common XuR,YuR,ZuR,XuG,YuG,ZuG,XuB,YuB,ZuB,
     1        RuX,GuX,BuX,RuY,GuY,BuY,RuZ,GuZ,BuZ,
     2        pi

!  Colour matching functions, 5 nm intervals, 
!  starting at 360 nm
!  Source:  
!  http://cvrl.ioo.ucl.ac.uk/database/data/cmfs/ciexyz31.txt
!  (new ordering) 
        data (xbar(i),i=1,95)/
     *   0.129900E-03,  0.232100E-03,  0.414900E-03, 
     *   0.741600E-03,  0.136800E-02,  0.223600E-02, 
     *   0.424300E-02,  0.765000E-02,  0.143100E-01, 
     *   0.231900E-01,  0.435100E-01,  0.776300E-01, 
     *   0.134380E+00,  0.214770E+00,  0.283900E+00, 
     *   0.328500E+00,  0.348280E+00,  0.348060E+00, 
     *   0.336200E+00,  0.318700E+00,  0.290800E+00, 
     *   0.251100E+00,  0.195360E+00,  0.142100E+00, 
     *   0.956400E-01,  0.579500E-01,  0.320100E-01, 
     *   0.147000E-01,  0.490000E-02,  0.240000E-02, 
     *   0.930000E-02,  0.291000E-01,  0.632700E-01, 
     *   0.109600E+00,  0.165500E+00,  0.225750E+00, 
     *   0.290400E+00,  0.359700E+00,  0.433450E+00, 
     *   0.512050E+00,  0.594500E+00,  0.678400E+00, 
     *   0.762100E+00,  0.842500E+00,  0.916300E+00, 
     *   0.978600E+00,  0.102630E+01,  0.105670E+01, 
     *   0.106220E+01,  0.104560E+01,  0.100260E+01, 
     *   0.938400E+00,  0.854450E+00,  0.751400E+00, 
     *   0.642400E+00,  0.541900E+00,  0.447900E+00, 
     *   0.360800E+00,  0.283500E+00,  0.218700E+00, 
     *   0.164900E+00,  0.121200E+00,  0.874000E-01, 
     *   0.636000E-01,  0.467700E-01,  0.329000E-01, 
     *   0.227000E-01,  0.158400E-01,  0.113592E-01, 
     *   0.811092E-02,  0.579035E-02,  0.410946E-02, 
     *   0.289933E-02,  0.204919E-02,  0.143997E-02, 
     *   0.999949E-03,  0.690079E-03,  0.476021E-03, 
     *   0.332301E-03,  0.234826E-03,  0.166150E-03, 
     *   0.117413E-03,  0.830753E-04,  0.587065E-04, 
     *   0.415099E-04,  0.293533E-04,  0.206738E-04, 
     *   0.145598E-04,  0.102540E-04,  0.722146E-05, 
     *   0.508587E-05,  0.358165E-05,  0.252252E-05, 
     *   0.177651E-05,  0.125114E-05/ 
        data (ybar(i),i=1,95)/
     *   0.391700E-05,  0.696500E-05,  0.123900E-04, 
     *   0.220200E-04,  0.390000E-04,  0.640000E-04, 
     *   0.120000E-03,  0.217000E-03,  0.396000E-03, 
     *   0.640000E-03,  0.121000E-02,  0.218000E-02, 
     *   0.400000E-02,  0.730000E-02,  0.116000E-01, 
     *   0.168400E-01,  0.230000E-01,  0.298000E-01, 
     *   0.380000E-01,  0.480000E-01,  0.600000E-01, 
     *   0.739000E-01,  0.909800E-01,  0.112600E+00, 
     *   0.139020E+00,  0.169300E+00,  0.208020E+00, 
     *   0.258600E+00,  0.323000E+00,  0.407300E+00, 
     *   0.503000E+00,  0.608200E+00,  0.710000E+00, 
     *   0.793200E+00,  0.862000E+00,  0.914850E+00, 
     *   0.954000E+00,  0.980300E+00,  0.994950E+00, 
     *   0.100000E+01,  0.995000E+00,  0.978600E+00, 
     *   0.952000E+00,  0.915400E+00,  0.870000E+00, 
     *   0.816300E+00,  0.757000E+00,  0.694900E+00, 
     *   0.631000E+00,  0.566800E+00,  0.503000E+00, 
     *   0.441200E+00,  0.381000E+00,  0.321000E+00, 
     *   0.265000E+00,  0.217000E+00,  0.175000E+00, 
     *   0.138200E+00,  0.107000E+00,  0.816000E-01, 
     *   0.610000E-01,  0.445800E-01,  0.320000E-01, 
     *   0.232000E-01,  0.170000E-01,  0.119200E-01, 
     *   0.821000E-02,  0.572300E-02,  0.410200E-02, 
     *   0.292900E-02,  0.209100E-02,  0.148400E-02, 
     *   0.104700E-02,  0.740000E-03,  0.520000E-03, 
     *   0.361100E-03,  0.249200E-03,  0.171900E-03, 
     *   0.120000E-03,  0.848000E-04,  0.600000E-04, 
     *   0.424000E-04,  0.300000E-04,  0.212000E-04, 
     *   0.149900E-04,  0.106000E-04,  0.746570E-05, 
     *   0.525780E-05,  0.370290E-05,  0.260780E-05, 
     *   0.183660E-05,  0.129340E-05,  0.910930E-06, 
     *   0.641530E-06,  0.451810E-06/ 
        data (zbar(i),i=1,95)/
     *   0.606100E-03,  0.108600E-02,  0.194600E-02, 
     *   0.348600E-02,  0.645000E-02,  0.105500E-01, 
     *   0.200500E-01,  0.362100E-01,  0.678500E-01, 
     *   0.110200E+00,  0.207400E+00,  0.371300E+00, 
     *   0.645600E+00,  0.103905E+01,  0.138560E+01, 
     *   0.162296E+01,  0.174706E+01,  0.178260E+01, 
     *   0.177211E+01,  0.174410E+01,  0.166920E+01, 
     *   0.152810E+01,  0.128764E+01,  0.104190E+01, 
     *   0.812950E+00,  0.616200E+00,  0.465180E+00, 
     *   0.353300E+00,  0.272000E+00,  0.212300E+00, 
     *   0.158200E+00,  0.111700E+00,  0.782500E-01, 
     *   0.572500E-01,  0.421600E-01,  0.298400E-01, 
     *   0.203000E-01,  0.134000E-01,  0.875000E-02, 
     *   0.575000E-02,  0.390000E-02,  0.275000E-02, 
     *   0.210000E-02,  0.180000E-02,  0.165000E-02, 
     *   0.140000E-02,  0.110000E-02,  0.100000E-02, 
     *   0.800000E-03,  0.600000E-03,  0.340000E-03, 
     *   0.240000E-03,  0.190000E-03,  0.100000E-03, 
     *   0.500000E-04,  0.300000E-04,  0.200000E-04, 
     *   0.100000E-04,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00/ 

      alf=14.388e3      
!   alf = h.c/kB (Planck constant times velocity of light                 
!   divided by Boltzmann constant)
      alamo=0.560       
!   the spectral distribution sl of the illumination is 
!   normalized to 1 at this wavelength
        Xw=0.   
        Yw=0.   
        Zw=0.   
        do i=1,95
        wlambda=0.355+0.005*i
        sl=(alamo/wlambda)**5*(exp(alf/(Tabs*alamo))-1.0)
     +      /(exp(alf/(Tabs*wlambda))-1.0)
        hil=Refl(ds,wlambda)*sl
        Xw=Xw+xbar(i)*hil
        Yw=Yw+ybar(i)*hil
        Zw=Zw+zbar(i)*hil
        enddo
        return
        end
c----------------------------------------------------------
      subroutine SSLE(a11,a12,a13,a21,a22,a23,a31,a32,a33,
     1     b1,b2,b3,x,y,z)
!   solve system of 3 linear equations
!   A_ik . x_k = b_i
      implicit real*8 (a-h,o-z)
      det=a11*a22*a33+a12*a23*a31+a13*a21*a32
     1  -a13*a22*a31-a12*a21*a33-a11*a23*a32
      x =(b1*a22*a33+a12*a23*b3+a13*b2*a32
     1  -a13*a22*b3-a12*b2*a33-b1*a23*a32)/det
      y =(a11*b2*a33+b1*a23*a31+a13*a21*b3
     1  -a13*b2*a31-b1*a21*a33-a11*a23*b3)/det
      z =(a11*a22*b3+a12*b2*a31+b1*a21*a32
     1  -b1*a22*a31-a12*a21*b3-a11*b2*a32)/det
      return
      end
c----------------------------------------------------------
