optics.f90

Go to the documentation of this file.

!> \file optics.f90  Procedures to do computations in optics
 
 
!  Copyright (c) 2002-2025  Marc van der Sluys - Nikhef/Utrecht University - marc.vandersluys.nl
!   
!  This file is part of the libSUFR package, 
!  see: http://libsufr.sourceforge.net/
!   
!  This is free software: you can redistribute it and/or modify it under the terms of the European Union
!  Public Licence 1.2 (EUPL 1.2).  This software is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
!  PURPOSE.  See the EU Public Licence for more details.  You should have received a copy of the European
!  Union Public Licence along with this code.  If not, see <https://www.eupl.eu/1.2/en/>.
!  
!  
 
 
 
 
!***********************************************************************************************************************************
!> \brief  Procedures to do computations in optics
 
module sufr_optics
  implicit none
  save
  
contains
  
  
  !***********************************************************************************************************************************
  !> \brief  Compute the reflectance for the transition from a medium with refractive index Nref1 to one with
  !!         Nref2, under an incident angle ang.  Optionally, compute the transmittance, and the polarised
  !!         components.  The media are assumed to be non-magnetic.
  !!
  !! \param angI   Angle of incidence (rad)
  !! \param Nref1  Refractive index of first medium, incoming ray
  !! \param Nref2  Refractive index of second medium, transmitted ray
  !! 
  !! \param Runp   Unpolarised reflectance (output)
  !! 
  !! \param Tunp   Unpolarised transmittance (output, optional)
  !! 
  !! \param Rprp   Perpendicular polarised reflectance (output, optional)
  !! \param Rpar   Parallel polarised reflectance (output, optional)
  !! \param Tprp   Perpendicular polarised transmittance (output, optional)
  !! \param Tpar   Parallel polarised transmittance (output, optional)
  !! 
  !! \param angT   Angle of transmittance (rad; optional) (output, optional)
  !!
  !! \see
  !! - Hecht, Optics, 3rd Ed. (1998), p.113ff
  !! - https://en.wikipedia.org/wiki/Fresnel_equations#Power_or_intensity_equations
  
  elemental subroutine reflectance_transmittance(angI, Nref1,Nref2,  Runp, Tunp, Rprp,Rpar, Tprp,Tpar, angT)
    use sufr_kinds, only: double
    use sufr_constants, only: pio2
    
    implicit none
    real(double), intent(in) :: angi, nref1,nref2
    real(double), intent(out) :: runp
    real(double), intent(out), optional :: rprp,rpar, tunp,tprp,tpar, angt
    real(double) :: var, langt, cosangi,cosangt, lrprp,lrpar
    
    var = nref1/nref2 * sin(angi)  ! Argument for Snell's law
    if(var.gt.1.d0 .or. abs(angi).gt.pio2) then  ! Total internal reflection or an impossible input value
       lrprp = 1.d0
       lrpar = 1.d0
       langt = 0.d0
    else
       langt = asin(var)  ! Angle of transmittance - Snell's law - local variable
       
       ! Reused variables:
       cosangi = cos(angi)
       cosangt = cos(langt)
       
       lrprp = ( (nref1 * cosangi - nref2 * cosangt) / (nref1 * cosangi + nref2 * cosangt) )**2
       lrpar = ( (nref1 * cosangt - nref2 * cosangi) / (nref1 * cosangt + nref2 * cosangi) )**2
    end if
    
    ! Unpolarised reflectance:
    runp = 0.5d0 * (lrprp + lrpar)
    
    ! Assign optional return values:
    if(present(rprp))  rprp = lrprp
    if(present(rpar))  rpar = lrpar
    
    if(present(tunp))  tunp = 1.d0 -  runp
    if(present(tprp))  tprp = 1.d0 - lrprp
    if(present(tpar))  tpar = 1.d0 - lrpar
    
    if(present(angt))  angt = langt
    
  elemental subroutine reflectance_transmittance(angI, Nref1,Nref2,  Runp, Tunp, Rprp,Rpar, Tprp,Tpar, angT) …
  end subroutine reflectance_transmittance
  !***********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief Computes real index of refraction of air given wavelength, temperature, pressure and relative humidity
  !!
  !! \param wavelength   Wavelength (nanometers)
  !! \param temperature  Temperature (degrees Celsius)
  !! \param pressure     Pressure (millibar)
  !! \param humidity     Relative humidity (%)
  !!
  !! \retval refractive_index_air  Refractive index of air
  !!
  !! \see http://emtoolbox.nist.gov/Wavelength/Documentation.asp#AppendixA
  
  pure function refractive_index_air(wavelength, temperature, pressure, humidity)
    use sufr_kinds, only: double
    implicit none
    
    real(double), intent(in) :: wavelength, temperature, pressure, humidity
    real(double) :: refractive_index_air
    real(double) :: s,t, omega, a,b,c,x, a1,a2,theta,y, psv, pv
    real(double) :: d,e,f,g, ns,ntp
    
    t = temperature + 273.15d0  ! deg Celsius -> Kelvin
    
    if(temperature.ge.0.d0) then  ! Above 0deg Celsius;  use IAPWS model:
       omega = t - 2.38555575678d-1/(t-6.50175348448d2)
       a =                  omega**2 + 1.16705214528d3*omega - 7.24213167032d5
       b = -1.70738469401d1*omega**2 + 1.20208247025d4*omega - 3.23255503223d6
       c =  1.49151086135d1*omega**2 - 4.82326573616d3*omega + 4.05113405421d5
       x = -b + sqrt(b**2 - 4*a*c)
       
       psv = 1.d6*(2*c/x)**4      ! saturation vapor pressure
       
    else                          ! Freezing
       
       a1 = -13.928169d0
       a2 =  34.7078238d0
       theta = t/273.16d0
       y = a1*(1.d0 - theta**(-1.5d0)) + a2*(1.d0 - theta**(-1.25d0))
       psv = 611.657d0 * exp(y)   ! saturation vapor pressure
    end if
    
    pv = (humidity/100.d0)*psv  ! Water-vapor partial pressure
    
    !  Modified Edlen Equation
    a = 8342.54d0
    b = 2406147.d0
    c = 15998.d0
    d = 96095.43d0
    e = 0.601d0
    f = 0.00972d0
    g = 0.003661d0
    
    s = 1.d6/wavelength**2
    ns = 1.d0 + 1.d-8 * (a + b/(130.d0-s) + c/(38.9d0-s))
    x = (1.d0 + 1.d-8 * (e - f*temperature) * (100*pressure)) / (1.d0 + g*temperature)
    ntp = 1.d0 + (100*pressure) * (ns-1.d0) * x/d
    
    refractive_index_air = ntp - 1.d-10*(292.75d0/t) * (3.7345d0 - 0.0401d0*s) * pv
    
  pure function refractive_index_air(wavelength, temperature, pressure, humidity) …
  end function refractive_index_air
  !*********************************************************************************************************************************
 
 
  !*********************************************************************************************************************************
  !> \brief  Refractive index of PMMA as a function of wavelength
  !!
  !! \param wavelength  Wavelength in nanometres
  !!
  !! \retval refractive_index_pmma  Refractive index of PMMA
  !!
  !! \note
  !! - valid for 23 degrees C
  !! - valid between 404.7 and 1083 nm
  !! 
  !! \see Marcin Szczurowski, http://refractiveindex.info/?shelf=3d&book=plastics&page=pmma
  
  pure function refractive_index_pmma(wavelength)
    use sufr_kinds, only: double
    implicit none
    real(double), intent(in) :: wavelength
    real(double) :: refractive_index_pmma, wl
    
    wl = wavelength/1000.d0  ! nm -> microns
    refractive_index_pmma = sqrt(1.d0 + 0.99654d0/(1.d0-0.00787d0/wl**2) + 0.18964d0/(1.d0-0.02191d0/wl**2) + &
         0.00411d0/(1.d0-3.85727d0/wl**2))
    
  pure function refractive_index_pmma(wavelength) …
  end function refractive_index_pmma
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Return an RGB-value representing the colour corresponding to a light ray with a given wavelength
  !!
  !! \param wl  Wavelength in nm (390-770 nm)
  !! \param df  Dimming factor 0-1 to scale RGB values with (0-1; 0: black, 1: full colour)
  !!
  !! \retval  wavelength2rgb  RGB values: (0-1, 0-1, 0-1)
  !!
  !!
  !! \note 
  !! - Assumed pure colour at centre of their wavelength bands:
  !!   Colours: Violet,   blue,   green,    yellow,    orange,     red:
  !!            (1,0,1)  (0,0,1)  (0,1,0)   (1,1,0)   (1,0.5,0)  (1,0,0)
  !! 
  !! - Colours between 390 and 420nm and 696 and 770nm are darker due to 'fading in' and 'fading out' effects
  
  function wavelength2rgb(wl, df)
    use sufr_kinds, only: double
    use sufr_system, only: warn
    
    implicit none
    real(double), intent(in) :: wl,df
    integer, parameter :: nc = 6  ! Number of colours in the spectrum (violet, blue, green, yellow, orange, red)
    integer :: ic
    real(double) :: wavelength2rgb(3), xipol,dxipol(nc),xipol0(nc), cbbnd(nc+1),cbctr(0:nc),cbdst(nc+1), rgb(3)
    
    
    ! Set the boundaries of the colour bands and compute their centres and mutual distances:
    cbbnd = dble([390,450,492,577,597,622,770])  ! Colour-band boundaries: 1:UV-V, 2:VB, 3:BG, 4:GY, 5:YO, 6:OR, 7:R-IR (nc+1)
    do ic=1,nc
       cbctr(ic)  = (cbbnd(ic) + cbbnd(ic+1))/2.d0      ! Centre of a colour band
       if(ic.gt.1) cbdst(ic) = cbctr(ic) - cbctr(ic-1)  ! Distance between colour bands
    end do
    cbdst(1) =  cbctr(1)-cbbnd(1)
    cbdst(nc+1) = cbbnd(nc+1)-cbctr(nc)
    
    
    ! Convert wavelength to 0 <= x <= 4.5:
    dxipol = dble([0.5,1.0,1.0,0.5,0.5,0.5])
    xipol0 = dble([0.5,1.0,2.0,3.0,3.5,4.0])
    xipol = (wl-390.d0)/30.d0*0.5d0                                         ! Black to violet    xIpol: 0.0 - 0.5
    do ic=1,nc
       if(wl.gt.cbctr(ic))  xipol = (wl-cbctr(ic)) / cbdst(ic+1) * dxipol(ic) + xipol0(ic)
    end do
    xipol = min(4.5d0, max(0.d0, xipol) )  ! xIpol should be 0 <= xIpol <= 4.5
    
    
    ! Convert xIpol to RGB.  Use black as 'invisible', i.e. UV or IR:
    rgb = [0.d0, 0.d0, 0.d0]  ! Ensure always defined
    if(xipol.ge.0.d0 .and. xipol.le.0.5d0) then
       rgb =                         [xipol,           0.d0,       xipol]       ! Black to violet    xIpol: 0.0 - 0.5
    else if(xipol.le.1.0d0) then
       rgb =                         [1.d0-xipol,      0.d0,       xipol]       ! Violet to blue     xIpol: 0.5 - 1.0
    else if(xipol.le.2.0d0) then
       rgb =                         [0.d0,            xipol-1.d0, 2.d0-xipol]  ! Blue to green      xIpol: 1.0 - 2.0
    else if(xipol.le.3.0d0) then
       rgb =                         [xipol-2.d0,      1.d0,       0.d0]        ! Green to yellow    xIpol: 2.0 - 3.0
    else if(xipol.le.3.5d0) then
       rgb =                         [1.d0,            4.d0-xipol, 0.d0]        ! Yellow to orange   xIpol: 3.0 - 3.5
    else if(xipol.le.4.0d0) then
       rgb =                         [1.d0,            4.d0-xipol, 0.d0]        ! Orange to red      xIpol: 3.5 - 4.0
    else if(xipol.le.4.5d0) then
       rgb =                         [(4.5d0-xipol)*2, 0.d0,       0.d0]        ! Red to black       xIpol: 4.0 - 4.5
    else
       call warn('SUFR_optics/wavelength2RGB(): xIpol out of bounds')
    end if
    
    
    ! Apply dimming and return:
    wavelength2rgb = rgb*df
    
  function wavelength2rgb(wl, df) …
  end function wavelength2rgb
  !*********************************************************************************************************************************
  
  
end module sufr_optics
!***********************************************************************************************************************************