ThinFilm.h

Go to the documentation of this file.

/*
  Implementation of Fresnel equations for thin film optical properties

  Copyright (C) 2010 Darko Veberic, darko.veberic@kit.edu

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program 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
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef _ThinFilm_h_
#define _ThinFilm_h_

#include <complex>


namespace TelescopeSimulatorLX {

  class ThinFilm {

  private:
    typedef std::complex<double> dcomplex;

    struct Etas {
      double fEta1;
      dcomplex fEta2;
      dcomplex fEta3;
    };

  public:
    static
    double
    GetReflectance(const double nIndex1,
                   const dcomplex& nIndex2,
                   const double nIndex3,
                   const double sinTheta1,
                   const double thickness,
                   const double wavelength)
    {
      return GetXtance<GetReflectance2>(nIndex1, nIndex2, nIndex3,
                                        sinTheta1, thickness, wavelength);
    }

    static
    double
    GetTransmittance(const double nIndex1,
                     const dcomplex& nIndex2,
                     const double nIndex3,
                     const double sinTheta1,
                     const double thickness,
                     const double wavelength)
    {
      return GetXtance<GetTransmittance2>(nIndex1, nIndex2, nIndex3,
                                          sinTheta1, thickness, wavelength);
    }

    static
    double
    GetAbsorptance(const double nIndex1,
                   const dcomplex& nIndex2,
                   const double nIndex3,
                   const double sinTheta1,
                   const double thickness,
                   const double wavelength)
    {
      return GetXtance<GetAbsorptance2>(nIndex1, nIndex2, nIndex3,
                                        sinTheta1, thickness, wavelength);
    }

    static
    std::pair<double, double>
    GetReflectanceAbsorptance(const double nIndex1,
                              const dcomplex& nIndex2,
                              const double nIndex3,
                              const double sinTheta1,
                              const double thickness,
                              const double wavelength)
    {
      const double n1sinTheta1 = nIndex1 * sinTheta1;
      const dcomplex phase2 =
        GetPhase(nIndex2, n1sinTheta1, thickness, wavelength);
      const Etas etasS = GetAllEtasS(nIndex1, nIndex2, nIndex3, n1sinTheta1);
      const std::pair<dcomplex, dcomplex> bcS =
        GetNormalizedFields(etasS, phase2);
      const Etas etasP = GetAllEtasP(nIndex1, nIndex2, nIndex3, n1sinTheta1);
      const std::pair<dcomplex, dcomplex> bcP =
        GetNormalizedFields(etasP, phase2);
      return
        std::make_pair(
          0.5 * (GetReflectance3(etasS, bcS) + GetReflectance3(etasP, bcP)),
          0.5 * (GetAbsorptance2(etasS, phase2) + GetAbsorptance2(etasP, phase2))
        );
    }

  private:
    template<double Xtance(const Etas&, const dcomplex&)>
    static
    double
    GetXtance(const double nIndex1,
              const dcomplex& nIndex2,
              const double nIndex3,
              const double sinTheta1,
              const double thickness,
              const double wavelength)
    {
      const double n1sinTheta1 = nIndex1 * sinTheta1;
      const dcomplex phase2 =
        GetPhase(nIndex2, n1sinTheta1, thickness, wavelength);
      const Etas etasS = GetAllEtasS(nIndex1, nIndex2, nIndex3, n1sinTheta1);
      const Etas etasP = GetAllEtasP(nIndex1, nIndex2, nIndex3, n1sinTheta1);
      return 0.5 * (Xtance(etasS, phase2) + Xtance(etasP, phase2));
    }

    template<typename T>
    static T GetSnellNCosTheta(const T& nIndex, const double n1sinTheta1)
    { return std::sqrt(nIndex*nIndex - n1sinTheta1*n1sinTheta1); }

    template<typename T>
    static T GetEtaS(const T& nIndex, const double n1sinTheta1)
    { return GetSnellNCosTheta(nIndex, n1sinTheta1); }

    template<typename T>
    static T GetEtaP(const T& nIndex, const double n1sinTheta1)
    { return nIndex*nIndex / GetEtaS(nIndex, n1sinTheta1); }

    static
    Etas
    GetAllEtasS(const double nIndex1,
                const dcomplex& nIndex2,
                const double nIndex3,
                const double n1sinTheta1)
    {
      const Etas etas = {
        GetEtaS<double>(nIndex1, n1sinTheta1),
        GetEtaS<dcomplex>(nIndex2, n1sinTheta1),
        GetEtaS<dcomplex>(nIndex3, n1sinTheta1)
      };
      return etas;
    }

    static
    Etas
    GetAllEtasP(const double nIndex1,
                const dcomplex& nIndex2,
                const double nIndex3,
                const double n1sinTheta1)
    {
      const Etas etas = {
        GetEtaP<double>(nIndex1, n1sinTheta1),
        GetEtaP<dcomplex>(nIndex2, n1sinTheta1),
        GetEtaP<dcomplex>(nIndex3, n1sinTheta1)
      };
      return etas;
    }

    static
    dcomplex
    GetPhase(const dcomplex& nIndex, const double n1sinTheta1,
             const double thickness, const double wavelength)
    {
      return (2*M_PI) * thickness / wavelength *
        GetSnellNCosTheta(nIndex, n1sinTheta1);
    }

    static
    std::pair<dcomplex, dcomplex>
    GetNormalizedFields(const Etas& etas, const dcomplex& phase)
    {
      const dcomplex cosPhase = std::cos(phase);
      const dcomplex isinPhase = dcomplex(0, 1) * std::sin(phase);
      return std::make_pair(cosPhase + isinPhase * etas.fEta3 / etas.fEta2,
                            isinPhase * etas.fEta2 + cosPhase * etas.fEta3);
    }

    static
    double
    GetReflectance2(const Etas& etas, const dcomplex& phase)
    {
      const std::pair<dcomplex, dcomplex> bc = GetNormalizedFields(etas, phase);
      return GetReflectance3(etas, bc);
    }

    static
    double
    GetReflectance3(const Etas& etas, const std::pair<dcomplex, dcomplex>& bc)
    {
      const dcomplex eta1b = etas.fEta1 * bc.first;
      const dcomplex rho = (eta1b - bc.second) / (eta1b + bc.second);
      return std::norm(rho);
    }

    static
    double
    GetTransmittance2(const Etas& etas, const dcomplex& phase)
    {
      const std::pair<dcomplex, dcomplex> bc = GetNormalizedFields(etas, phase);
      const dcomplex x = etas.fEta1 * bc.first + bc.second;
      return 4 * etas.fEta1 * std::real(etas.fEta3) / std::norm(x);
    }

    static
    double
    GetAbsorptance2(const Etas& etas, const dcomplex& phase)
    {
      const std::pair<dcomplex, dcomplex> bc = GetNormalizedFields(etas, phase);
      return GetAbsorptance3(etas, bc);
    }

    static
    double
    GetAbsorptance3(const Etas& etas, const std::pair<dcomplex, dcomplex>& bc)
    {
      const dcomplex x = etas.fEta1 * bc.first + bc.second;
      return 4 * etas.fEta1 *
        std::real(bc.first * conj(bc.second) - etas.fEta3) / std::norm(x);
    }

  };

}


#endif