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) 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 ! 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 do i=1,3 colcomp(i,j)=colcomp(i,j)/themaximum if (colcomp(i,j).lt.0.0) colcomp(i,j)=0.0 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----------------------------------------------------------