tmatrix.cc

Go to the documentation of this file.

/* Copyright (C) 2013 Oliver Lemke
 *
 * 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/>.
 *
 */
 
#include "tmatrix.h"
#include <cmath>
#include <cstring>
#include <stdexcept>
#include "complex.h"
#include "math_funcs.h"
#include "matpackI.h"
#include "messages.h"
#include "optproperties.h"
 
void calc_phamat(Matrix& z,
                 const Index& nmax,
                 const Numeric& lam,
                 const Numeric& thet0,
                 const Numeric& thet,
                 const Numeric& phi0,
                 const Numeric& phi,
                 const Numeric& beta,
                 const Numeric& alpha);
 
void ampmat_to_phamat(Matrix& z,
                      const Complex& s11,
                      const Complex& s12,
                      const Complex& s21,
                      const Complex& s22);
 
void integrate_phamat_alpha10(Matrix& phamat,
                              const Index& nmax,
                              const Numeric& lam,
                              const Numeric& thet0,
                              const Numeric& thet,
                              const Numeric& phi0,
                              const Numeric& phi,
                              const Numeric& beta,
                              const Numeric& alpha_1,
                              const Numeric& alpha_2);
 
void integrate_phamat_alpha6(Matrix& phamat,
                             const Index& nmax,
                             const Numeric& lam,
                             const Numeric& thet0,
                             const Numeric& thet,
                             const Numeric& phi0,
                             const Numeric& phi,
                             const Numeric& beta,
                             const Numeric& alpha_1,
                             const Numeric& alpha_2);
 
void integrate_phamat_theta0_phi10(Matrix& phamat,
                                   const Index& nmax,
                                   const Numeric& lam,
                                   const Numeric& thet0_1,
                                   const Numeric& thet0_2,
                                   const Numeric& thet,
                                   const Numeric& phi0,
                                   const Numeric& phi_1,
                                   const Numeric& phi_2,
                                   const Numeric& beta,
                                   const Numeric& alpha);
 
void integrate_phamat_theta0_phi_alpha6(Matrix& phamat,
                                        const Index& nmax,
                                        const Numeric& lam,
                                        const Numeric& thet0_1,
                                        const Numeric& thet0_2,
                                        const Numeric& thet,
                                        const Numeric& phi0,
                                        const Numeric& phi_1,
                                        const Numeric& phi_2,
                                        const Numeric& beta,
                                        const Numeric& alpha_1,
                                        const Numeric& alpha_2);
 
#ifdef ENABLE_TMATRIX
extern "C" {
#endif
 
void tmd_(const Numeric& rat,
          const Index& ndistr,
          const Numeric& axmax,
          const Index& npnax,
          const Numeric& b,
          const Numeric& gam,
          const Index& nkmax,
          const Numeric& eps,
          const Index& np,
          const Numeric& lam,
          const Numeric& mrr,
          const Numeric& mri,
          const Numeric& ddelt,
          const Index& npna,
          const Index& ndgs,
          const Numeric& r1rat,
          const Numeric& r2rat,
          const Index& quiet,
          Numeric& reff,   // OUT
          Numeric& veff,   // OUT
          Numeric& cext,   // OUT
          Numeric& csca,   // OUT
          Numeric& walb,   // OUT
          Numeric& asymm,  // OUT
          Numeric* f11,    // Array npna
          Numeric* f22,    // Array npna
          Numeric* f33,    // Array npna
          Numeric* f44,    // Array npna
          Numeric* f12,    // Array npna
          Numeric* f34,    // Array npna
          char* errmsg);
 
void tmatrix_(const Numeric& rat,
              const Numeric& axi,
              const Index& np,
              const Numeric& lam,
              const Numeric& eps,
              const Numeric& mrr,
              const Numeric& mri,
              const Numeric& ddelt,
              const Index& quiet,
              Index& nmax,
              Numeric& csca,
              Numeric& cext,
              char* errmsg);
 
void ampl_(const Index& nmax,
           const Numeric& lam,
           const Numeric& thet0,
           const Numeric& thet,
           const Numeric& phi0,
           const Numeric& phi,
           const Numeric& alpha,
           const Numeric& beta,
           Complex& s11,
           Complex& s12,
           Complex& s21,
           Complex& s22);
 
void avgtmatrix_(const Index& nmax);
#ifdef ENABLE_TMATRIX
}
#endif
 
// Define dummy functions that throw runtime errors if ARTS is
// compiled without T-Matrix support.
#ifndef ENABLE_TMATRIX
 
// T-matrix code for randomly oriented nonspherical particles.
void tmd_(const Numeric&,
          const Index&,
          const Numeric&,
          const Index&,
          const Numeric&,
          const Numeric&,
          const Index&,
          const Numeric&,
          const Index&,
          const Numeric&,
          const Numeric&,
          const Numeric&,
          const Numeric&,
          const Index&,
          const Index&,
          const Numeric&,
          const Numeric&,
          const Index&,
          Numeric&,
          Numeric&,
          Numeric&,
          Numeric&,
          Numeric&,
          Numeric&,
          Numeric*,
          Numeric*,
          Numeric*,
          Numeric*,
          Numeric*,
          Numeric*,
          char*) {
  throw std::runtime_error(
      "This version of ARTS was compiled without T-Matrix support.");
}
 
// T-matrix code for nonspherical particles in a fixed orientation
void tmatrix_(const Numeric&,
              const Numeric&,
              const Index&,
              const Numeric&,
              const Numeric&,
              const Numeric&,
              const Numeric&,
              const Numeric&,
              const Index&,
              Index&,
              Numeric&,
              Numeric&,
              char*) {
  throw std::runtime_error(
      "This version of ARTS was compiled without T-Matrix support.");
}
 
// T-matrix code for nonspherical particles in a fixed orientation
void ampl_(const Index&,
           const Numeric&,
           const Numeric&,
           const Numeric&,
           const Numeric&,
           const Numeric&,
           const Numeric&,
           const Numeric&,
           Complex&,
           Complex&,
           Complex&,
           Complex&) {
  throw std::runtime_error(
      "This version of ARTS was compiled without T-Matrix support.");
}
 
void avgtmatrix_(const Index&) {
  throw std::runtime_error(
      "This version of ARTS was compiled without T-Matrix support.");
}
 
#endif
 
void calc_phamat(Matrix& z,
                 const Index& nmax,
                 const Numeric& lam,
                 const Numeric& thet0,
                 const Numeric& thet,
                 const Numeric& phi0,
                 const Numeric& phi,
                 const Numeric& beta,
                 const Numeric& alpha) {
  Complex s11;
  Complex s12;
  Complex s21;
  Complex s22;
  ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
 
  ampmat_to_phamat(z, s11, s12, s21, s22);
}
 
void ampmat_to_phamat(Matrix& z,
                      const Complex& s11,
                      const Complex& s12,
                      const Complex& s21,
                      const Complex& s22) {
  z.resize(4, 4);
  z(0, 0) = 0.5 * (s11 * conj(s11) + s12 * conj(s12) + s21 * conj(s21) +
                   s22 * conj(s22))
                      .real();
  z(0, 1) = 0.5 * (s11 * conj(s11) - s12 * conj(s12) + s21 * conj(s21) -
                   s22 * conj(s22))
                      .real();
  z(0, 2) = (-s11 * conj(s12) - s22 * conj(s21)).real();
  z(0, 3) = (Complex(0., 1.) * (s11 * conj(s12) - s22 * conj(s21))).real();
 
  z(1, 0) = 0.5 * (s11 * conj(s11) + s12 * conj(s12) - s21 * conj(s21) -
                   s22 * conj(s22))
                      .real();
  z(1, 1) = 0.5 * (s11 * conj(s11) - s12 * conj(s12) - s21 * conj(s21) +
                   s22 * conj(s22))
                      .real();
  z(1, 2) = (-s11 * conj(s12) + s22 * conj(s21)).real();
  z(1, 3) = (Complex(0., 1.) * (s11 * conj(s12) + s22 * conj(s21))).real();
 
  z(2, 0) = (-s11 * conj(s21) - s22 * conj(s12)).real();
  z(2, 1) = (-s11 * conj(s21) + s22 * conj(s12)).real();
  z(2, 2) = (s11 * conj(s22) + s12 * conj(s21)).real();
  z(2, 3) = (Complex(0., -1.) * (s11 * conj(s22) + s21 * conj(s12))).real();
 
  z(3, 0) = (Complex(0., 1.) * (s21 * conj(s11) + s22 * conj(s12))).real();
  z(3, 1) = (Complex(0., 1.) * (s21 * conj(s11) - s22 * conj(s12))).real();
  z(3, 2) = (Complex(0., -1.) * (s22 * conj(s11) - s12 * conj(s21))).real();
  z(3, 3) = (s22 * conj(s11) - s12 * conj(s21)).real();
}
 
static const Numeric GaussLeg6[][3] = {{0.23861918, 0.66120939, 0.93246951},
                                       {0.46791393, 0.36076157, 0.17132449}};
 
static const Numeric GaussLeg10[][5] = {
    {0.14887434, 0.43339539, 0.67940957, 0.86506337, 0.97390653},
    {0.29552422, 0.26926672, 0.21908636, 0.14945135, 0.06667134}};
 
void integrate_phamat_alpha10(Matrix& phamat,
                              const Index& nmax,
                              const Numeric& lam,
                              const Numeric& thet0,
                              const Numeric& thet,
                              const Numeric& phi0,
                              const Numeric& phi,
                              const Numeric& beta,
                              const Numeric& alpha_1,
                              const Numeric& alpha_2) {
  const Numeric c = 0.5 * (alpha_2 + alpha_1);
  const Numeric m = 0.5 * (alpha_2 - alpha_1);
 
  phamat.resize(4, 4);
  phamat = 0.;
  Matrix z;
 
  for (Index i = 0; i < 5; ++i) {
    const Numeric abscissa = GaussLeg10[0][i];
    const Numeric weight = GaussLeg10[1][i];
 
    calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m * abscissa);
    z *= weight;
    phamat += z;
 
    calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m * abscissa);
    z *= weight;
    phamat += z;
  }
  phamat *= m;
}
 
void integrate_phamat_alpha6(Matrix& phamat,
                             const Index& nmax,
                             const Numeric& lam,
                             const Numeric& thet0,
                             const Numeric& thet,
                             const Numeric& phi0,
                             const Numeric& phi,
                             const Numeric& beta,
                             const Numeric& alpha_1,
                             const Numeric& alpha_2) {
  const Numeric c = 0.5 * (alpha_2 + alpha_1);
  const Numeric m = 0.5 * (alpha_2 - alpha_1);
 
  phamat.resize(4, 4);
  phamat = 0.;
  Matrix z;
 
  for (Index i = 0; i < 3; ++i) {
    const Numeric abscissa = GaussLeg6[0][i];
    const Numeric weight = GaussLeg6[1][i];
 
    calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m * abscissa);
    z *= weight;
    phamat += z;
 
    calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m * abscissa);
    z *= weight;
    phamat += z;
  }
  phamat *= m;
}
 
void integrate_phamat_theta0_phi10(Matrix& phamat,
                                   const Index& nmax,
                                   const Numeric& lam,
                                   const Numeric& thet0_1,
                                   const Numeric& thet0_2,
                                   const Numeric& thet,
                                   const Numeric& phi0,
                                   const Numeric& phi_1,
                                   const Numeric& phi_2,
                                   const Numeric& beta,
                                   const Numeric& alpha) {
  extern const Numeric PI;
 
  phamat.resize(4, 4);
  phamat = 0.;
  Matrix z;
 
  const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
  const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
  const Numeric c_phi = 0.5 * (phi_2 + phi_1);
  const Numeric m_phi = 0.5 * (phi_2 - phi_1);
 
  for (Index t = 0; t < 5; ++t) {
    const Numeric abscissa_thet0 = GaussLeg10[0][t];
    const Numeric weight_thet0 = GaussLeg10[1][t];
 
    Matrix phamat_phi(4, 4, 0.);
 
    for (Index p = 0; p < 5; ++p) {
      const Numeric abscissa_phi = GaussLeg10[0][p];
      const Numeric weight_phi = GaussLeg10[1][p];
 
      Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
      calc_phamat(z,
                  nmax,
                  lam,
                  this_thet0,
                  thet,
                  phi0,
                  c_phi + m_phi * abscissa_phi,
                  beta,
                  alpha);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
 
      calc_phamat(z,
                  nmax,
                  lam,
                  this_thet0,
                  thet,
                  phi0,
                  c_phi - m_phi * abscissa_phi,
                  beta,
                  alpha);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
 
      this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
      calc_phamat(z,
                  nmax,
                  lam,
                  this_thet0,
                  thet,
                  phi0,
                  c_phi + m_phi * abscissa_phi,
                  beta,
                  alpha);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
 
      calc_phamat(z,
                  nmax,
                  lam,
                  this_thet0,
                  thet,
                  phi0,
                  c_phi - m_phi * abscissa_phi,
                  beta,
                  alpha);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
    }
    phamat_phi *= m_phi * weight_thet0;
    phamat += phamat_phi;
  }
 
  phamat *= m_thet0;
}
 
void integrate_phamat_theta0_phi_alpha6(Matrix& phamat,
                                        const Index& nmax,
                                        const Numeric& lam,
                                        const Numeric& thet0_1,
                                        const Numeric& thet0_2,
                                        const Numeric& thet,
                                        const Numeric& phi0,
                                        const Numeric& phi_1,
                                        const Numeric& phi_2,
                                        const Numeric& beta,
                                        const Numeric& alpha_1,
                                        const Numeric& alpha_2) {
  extern const Numeric PI;
 
  Matrix z;
  Matrix phamat_phi(4, 4);
 
  const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
  const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
  const Numeric c_phi = 0.5 * (phi_2 + phi_1);
  const Numeric m_phi = 0.5 * (phi_2 - phi_1);
 
  phamat.resize(4, 4);
  phamat = 0.;
 
  for (Index t = 0; t < 3; ++t) {
    const Numeric abscissa_thet0 = GaussLeg6[0][t];
    const Numeric weight_thet0 = GaussLeg6[1][t];
 
    phamat_phi = 0.;
 
    for (Index p = 0; p < 3; ++p) {
      const Numeric abscissa_phi = GaussLeg6[0][p];
      const Numeric weight_phi = GaussLeg6[1][p];
 
      Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
      integrate_phamat_alpha6(z,
                              nmax,
                              lam,
                              this_thet0,
                              thet,
                              phi0,
                              c_phi + m_phi * abscissa_phi,
                              beta,
                              alpha_1,
                              alpha_2);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
 
      integrate_phamat_alpha6(z,
                              nmax,
                              lam,
                              this_thet0,
                              thet,
                              phi0,
                              c_phi - m_phi * abscissa_phi,
                              beta,
                              alpha_1,
                              alpha_2);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
 
      this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
      integrate_phamat_alpha6(z,
                              nmax,
                              lam,
                              this_thet0,
                              thet,
                              phi0,
                              c_phi + m_phi * abscissa_phi,
                              beta,
                              alpha_1,
                              alpha_2);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
 
      integrate_phamat_alpha6(z,
                              nmax,
                              lam,
                              this_thet0,
                              thet,
                              phi0,
                              c_phi - m_phi * abscissa_phi,
                              beta,
                              alpha_1,
                              alpha_2);
      z *= weight_phi * sin(this_thet0 * PI / 180.);
      phamat_phi += z;
    }
    phamat_phi *= m_phi * weight_thet0;
    phamat += phamat_phi;
  }
 
  phamat *= m_thet0;
}
 
void tmatrix_random_orientation(Numeric& cext,
                                Numeric& csca,
                                Vector& f11,
                                Vector& f22,
                                Vector& f33,
                                Vector& f44,
                                Vector& f12,
                                Vector& f34,
                                const Numeric equiv_radius,
                                const Numeric aspect_ratio,
                                const Index np,
                                const Numeric lam,
                                const Numeric ref_index_real,
                                const Numeric ref_index_imag,
                                const Numeric precision,
                                const Index nza,
                                const Index ndgs,
                                const Index quiet = 1) {
  Numeric reff;
  Numeric veff;
  Numeric walb;
  Numeric asymm;
 
  char errmsg[1024] = "";
 
  f11.resize(nza);
  f11 = NAN;
  f22.resize(nza);
  f22 = NAN;
  f33.resize(nza);
  f33 = NAN;
  f44.resize(nza);
  f44 = NAN;
  f12.resize(nza);
  f12 = NAN;
  f34.resize(nza);
  f34 = NAN;
 
  // It is necessary to make sure that the Fortran code is not
  // called from different threads at the same time. The common
  // blocks are not threadsafe.
#pragma omp critical(tmatrix_code)
  tmd_(1.0,
       4,
       equiv_radius,
       1,
       0.1,
       1.0,
       -1,
       aspect_ratio,
       np,
       lam,
       ref_index_real,
       ref_index_imag,
       precision,
       nza,
       ndgs,
       0.9999999,
       1.0000001,
       quiet,
       reff,
       veff,
       cext,
       csca,
       walb,
       asymm,
       f11.get_c_array(),
       f22.get_c_array(),
       f33.get_c_array(),
       f44.get_c_array(),
       f12.get_c_array(),
       f34.get_c_array(),
       errmsg);
 
  if (strlen(errmsg)) {
    std::ostringstream os;
    os << "T-Matrix code failed: " << errmsg;
    throw std::runtime_error(os.str());
  }
}
 
void tmatrix_fixed_orientation(Numeric& cext,
                               Numeric& csca,
                               Index& nmax,
                               const Numeric equiv_radius,
                               const Numeric aspect_ratio,
                               const Index np,
                               const Numeric lam,
                               const Numeric ref_index_real,
                               const Numeric ref_index_imag,
                               const Numeric precision,
                               const Index quiet = 1) {
  char errmsg[1024] = "";
 
  // It is necessary to make sure that the Fortran code is not
  // called from different threads at the same time. The common
  // blocks are not threadsafe.
#pragma omp critical(tmatrix_code)
  tmatrix_(1.,
           equiv_radius,
           np,
           lam,
           aspect_ratio,
           ref_index_real,
           ref_index_imag,
           precision,
           quiet,
           nmax,
           csca,
           cext,
           errmsg);
 
  if (strlen(errmsg)) {
    std::ostringstream os;
    os << "T-Matrix code failed: " << errmsg;
    throw std::runtime_error(os.str());
  }
}
 
void calcSingleScatteringDataProperties(SingleScatteringData& ssd,
                                        ConstMatrixView ref_index_real,
                                        ConstMatrixView ref_index_imag,
                                        const Numeric equiv_radius,
                                        const Index np,
                                        const Numeric aspect_ratio,
                                        const Numeric precision,
                                        const Index ndgs,
                                        const Index robust,
                                        const Index quiet) {
  const Index nf = ssd.f_grid.nelem();
  const Index nT = ssd.T_grid.nelem();
 
  extern const Numeric PI;
  extern const Numeric SPEED_OF_LIGHT;
 
  if (ref_index_real.nrows() != nf)
    throw std::runtime_error(
        "Number of rows of refractive index real part must match ssd.f_grid.");
  if (ref_index_real.ncols() != nT)
    throw std::runtime_error(
        "Number of cols of refractive index real part must match ssd.T_grid.");
  if (ref_index_imag.nrows() != nf)
    throw std::runtime_error(
        "Number of rows of refractive index imaginary part must match ssd.f_grid.");
  if (ref_index_imag.ncols() != nT)
    throw std::runtime_error(
        "Number of cols of refractive index imaginary part must match ssd.T_grid.");
 
  // The T-Matrix code needs wavelength
  Vector lam(nf, SPEED_OF_LIGHT);
  lam /= ssd.f_grid;
 
  switch (ssd.ptype) {
    case PTYPE_TOTAL_RND: {
      const Index nza = ssd.za_grid.nelem();
 
      ssd.pha_mat_data.resize(nf, nT, nza, 1, 1, 1, 6);
      ssd.ext_mat_data.resize(nf, nT, 1, 1, 1);
      ssd.abs_vec_data.resize(nf, nT, 1, 1, 1);
 
      ssd.pha_mat_data = NAN;
      ssd.ext_mat_data = NAN;
      ssd.abs_vec_data = NAN;
 
      // Output variables
      Numeric cext = NAN;
      Numeric csca = NAN;
      Vector f11;
      Vector f22;
      Vector f33;
      Vector f44;
      Vector f12;
      Vector f34;
      Matrix mono_pha_mat_data(nza, 6, NAN);
 
      ostringstream os;
      os << "Calculation of SingleScatteringData properties failed for\n\n";
      bool anyfailed = false;
#pragma omp critical(tmatrix_ssp)
      for (Index f_index = 0; f_index < nf; ++f_index)
        for (Index T_index = 0; T_index < nT; ++T_index) {
          bool thisfailed = false;
          try {
            tmatrix_random_orientation(cext,
                                       csca,
                                       f11,
                                       f22,
                                       f33,
                                       f44,
                                       f12,
                                       f34,
                                       equiv_radius,
                                       aspect_ratio,
                                       np,
                                       lam[f_index],
                                       ref_index_real(f_index, T_index),
                                       ref_index_imag(f_index, T_index),
                                       precision,
                                       nza,
                                       ndgs,
                                       quiet);
          } catch (const std::runtime_error& e) {
            //ostringstream os;
            //os << "Calculation of SingleScatteringData properties failed for\n"
            os << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
               << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index]
               << "\n\n";
            //<< e.what();
            //throw std::runtime_error(os.str());
            thisfailed = true;
            anyfailed = true;
            cout << "\n\n";
          }
 
          if (!thisfailed) {
            mono_pha_mat_data(joker, 0) = f11;
            mono_pha_mat_data(joker, 1) = f12;
            mono_pha_mat_data(joker, 2) = f22;
            mono_pha_mat_data(joker, 3) = f33;
            mono_pha_mat_data(joker, 4) = f34;
            mono_pha_mat_data(joker, 5) = f44;
 
            mono_pha_mat_data *= csca / 4. / PI;
            ssd.pha_mat_data(f_index, T_index, joker, 0, 0, 0, joker) =
                mono_pha_mat_data;
 
            ssd.ext_mat_data(f_index, T_index, 0, 0, 0) = cext;
            ssd.abs_vec_data(f_index, T_index, 0, 0, 0) = cext - csca;
          }
        }
      if (anyfailed)
        if (robust)
          cout << os.str();
        else
          throw std::runtime_error(os.str());
      else {
        os << "None\n";
      }
 
      break;
    }
    case PTYPE_AZIMUTH_RND: {
      const Index nza = ssd.za_grid.nelem();
      const Index naa = ssd.aa_grid.nelem();
 
      ssd.pha_mat_data.resize(nf, nT, nza, naa, nza, 1, 16);
      ssd.ext_mat_data.resize(nf, nT, nza, 1, 3);
      ssd.abs_vec_data.resize(nf, nT, nza, 1, 2);
 
      ssd.ext_mat_data = NAN;
      ssd.pha_mat_data = NAN;
      ssd.abs_vec_data = NAN;
 
      // Output variables
      Numeric cext = NAN;
      Numeric csca = NAN;
      Index nmax = -1;
 
      Tensor5 csca_data(nf, nT, nza, 1, 2);
 
#pragma omp critical(tmatrix_ssp)
      for (Index f_index = 0; f_index < nf; ++f_index) {
        const Numeric lam_f = lam[f_index];
 
        for (Index T_index = 0; T_index < nT; ++T_index) {
          try {
            tmatrix_fixed_orientation(cext,
                                      csca,
                                      nmax,
                                      equiv_radius,
                                      aspect_ratio,
                                      np,
                                      lam_f,
                                      ref_index_real(f_index, T_index),
                                      ref_index_imag(f_index, T_index),
                                      precision);
          } catch (const std::runtime_error& e) {
            ostringstream os;
            os << "Calculation of SingleScatteringData properties failed for\n"
               << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
               << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index] << "\n"
               << e.what();
            throw std::runtime_error(os.str());
          }
 
          Matrix phamat;
          for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index)
            for (Index aa_index = 0; aa_index < naa; ++aa_index)
              for (Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
                if (aspect_ratio < 1.0) {
                  // Phase matrix for prolate particles
                  integrate_phamat_alpha10(phamat,
                                           nmax,
                                           lam_f,
                                           ssd.za_grid[za_inc_index],
                                           ssd.za_grid[za_scat_index],
                                           0.0,
                                           ssd.aa_grid[aa_index],
                                           90.0,
                                           0.0,
                                           180.0);
                  phamat /= 180.;
                } else {
                  // Phase matrix for oblate particles
                  calc_phamat(phamat,
                              nmax,
                              lam_f,
                              ssd.za_grid[za_inc_index],
                              ssd.za_grid[za_scat_index],
                              0.0,
                              ssd.aa_grid[aa_index],
                              0.0,
                              0.0);
                }
 
                ssd.pha_mat_data(f_index,
                                 T_index,
                                 za_scat_index,
                                 aa_index,
                                 za_inc_index,
                                 0,
                                 Range(0, 4)) = phamat(0, joker);
                ssd.pha_mat_data(f_index,
                                 T_index,
                                 za_scat_index,
                                 aa_index,
                                 za_inc_index,
                                 0,
                                 Range(4, 4)) = phamat(1, joker);
                ssd.pha_mat_data(f_index,
                                 T_index,
                                 za_scat_index,
                                 aa_index,
                                 za_inc_index,
                                 0,
                                 Range(8, 4)) = phamat(2, joker);
                ssd.pha_mat_data(f_index,
                                 T_index,
                                 za_scat_index,
                                 aa_index,
                                 za_inc_index,
                                 0,
                                 Range(12, 4)) = phamat(3, joker);
              }
 
          // Csca integral
          for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index) {
            Matrix csca_integral;
            if (aspect_ratio < 1.0) {
              // Csca for prolate particles
              integrate_phamat_theta0_phi_alpha6(csca_integral,
                                                 nmax,
                                                 lam_f,
                                                 0,
                                                 180,
                                                 ssd.za_grid[za_scat_index],
                                                 0.,
                                                 0.,
                                                 180.,
                                                 90.,
                                                 0.,
                                                 180.);
              csca_integral /= 180.;
            } else {
              // Csca for oblate particles
              integrate_phamat_theta0_phi10(csca_integral,
                                            nmax,
                                            lam_f,
                                            0,
                                            180,
                                            ssd.za_grid[za_scat_index],
                                            0.,
                                            0.,
                                            180,
                                            0.,
                                            0.);
            }
            csca_data(f_index, T_index, za_scat_index, 0, joker) =
                csca_integral(Range(0, 2), 0);
          }
 
          // Extinction matrix
          if (aspect_ratio < 1.0) {
            // Average T-Matrix for prolate particles
            avgtmatrix_(nmax);
          }
 
          for (Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
            Complex s11;
            Complex s12;
            Complex s21;
            Complex s22;
            VectorView K =
                ssd.ext_mat_data(f_index, T_index, za_inc_index, 0, joker);
 
            const Numeric beta = 0.;
            const Numeric alpha = 0.;
            ampl_(nmax,
                  lam_f,
                  ssd.za_grid[za_inc_index],
                  ssd.za_grid[za_inc_index],
                  0.,
                  0.,
                  alpha,
                  beta,
                  s11,
                  s12,
                  s21,
                  s22);
 
            K[0] = (Complex(0., -1.) * (s11 + s22)).real();
            K[1] = (Complex(0., 1.) * (s22 - s11)).real();
            K[2] = (s22 - s11).real();
 
            K *= lam_f;
          }
        }
      }
 
      csca_data *= 2. * PI * PI / 32400.;
      ssd.abs_vec_data =
          ssd.ext_mat_data(joker, joker, joker, joker, Range(0, 2));
      ssd.abs_vec_data -= csca_data;
 
      break;
    }
    default: {
      std::ostringstream os;
      os << "Only particle types totally_random and azimuthally_random\n"
         << "are currently supported: " << ssd.ptype;
      throw std::runtime_error(os.str());
      break;
    }
  }
}
 
// Documentation in header file.
void tmatrix_ampld_test(const Verbosity& verbosity) {
  CREATE_OUT0;
 
  out0 << "======================================================\n";
  out0 << "Test for nonspherical particles in a fixed orientation\n";
  out0 << "Output should match 3rdparty/tmatrix/tmatrix_ampld.ref\n";
  out0 << "======================================================\n";
 
  // Same inputs as in example included in original ampld.lp.f
  Numeric rat = 1.;
  Numeric axi = 10.;  //[um]
  Index np = -1;
  Numeric lam = acos(-1.) * 2.;  //[um]
  Numeric eps = 0.5;
  Numeric mrr = 1.5;
  Numeric mri = 0.02;
  Numeric ddelt = 0.001;
 
  Index quiet = 1;
 
  // Output variables
  Index nmax;
  Numeric csca;
  Numeric cext;
  char errmsg[1024] = "";
 
  tmatrix_(
      rat, axi, np, lam, eps, mrr, mri, ddelt, quiet, nmax, csca, cext, errmsg);
 
  out0 << "nmax: " << nmax << "\n";
  out0 << "csca: " << csca << " um2\n";
  out0 << "cext: " << cext << " um2\n";
 
  out0 << "Error message: " << (strlen(errmsg) ? errmsg : "None") << "\n";
 
  // Same inputs as in example included in original ampld.lp.f
  Numeric alpha = 145.;
  Numeric beta = 52.;
  Numeric thet0 = 56.;
  Numeric thet = 65.;
  Numeric phi0 = 114.;
  Numeric phi = 128.;
 
  // Output variables
  Complex s11;
  Complex s12;
  Complex s21;
  Complex s22;
  ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
 
  out0 << "AMPLITUDE MATRIX (all in [um]): \n";
  out0 << "s11: " << s11 << "\n";
  out0 << "s12: " << s12 << "\n";
  out0 << "s21: " << s21 << "\n";
  out0 << "s22: " << s22 << "\n";
 
  Matrix z;
  ampmat_to_phamat(z, s11, s12, s21, s22);
  //z *= 1e12; // meter^2 to micron^2 for comparison with original results
 
  out0 << "PHASE MATRIX (all un [um2]): \n" << z << "\n";
}
 
// Documentation in header file.
void tmatrix_tmd_test(const Verbosity& verbosity) {
  CREATE_OUT0;
 
  out0 << "======================================================\n";
  out0 << "Test for randomly oriented nonspherical particles\n";
  out0 << "Output should match 3rdparty/tmatrix/tmatrix_tmd.ref\n";
  out0 << "======================================================\n";
 
  // Same inputs as in example included in original tmd.lp.f
  Numeric rat = 0.5;
  Index ndistr = 3;
  Numeric axmax = 1.;  //[um]
  Index npnax = 2;
  Numeric b = 0.1;
  Numeric gam = 0.5;
  Index nkmax = 5;
  Numeric eps = 2;
  Index np = -1;
  Numeric lam = 0.5;  //[um]
  Numeric mrr = 1.53;
  Numeric mri = 0.008;
  Numeric ddelt = 0.001;
  Index npna = 19;
  Index ndgs = 2;
  Numeric r1rat = 0.89031;
  Numeric r2rat = 1.56538;
 
  Index quiet = 1;
 
  // Output variables
  Numeric reff;
  Numeric veff;
  Numeric cext;
  Numeric csca;
  Numeric walb;
  Numeric asymm;
  Vector f11(npna, 0.);
  Vector f22(npna, 0.);
  Vector f33(npna, 0.);
  Vector f44(npna, 0.);
  Vector f12(npna, 0.);
  Vector f34(npna, 0.);
  char errmsg[1024] = "";
 
  tmd_(rat,
       ndistr,
       axmax,
       npnax,
       b,
       gam,
       nkmax,
       eps,
       np,
       lam,
       mrr,
       mri,
       ddelt,
       npna,
       ndgs,
       r1rat,
       r2rat,
       quiet,
       reff,
       veff,
       cext,
       csca,
       walb,
       asymm,
       f11.get_c_array(),
       f22.get_c_array(),
       f33.get_c_array(),
       f44.get_c_array(),
       f12.get_c_array(),
       f34.get_c_array(),
       errmsg);
 
  out0 << "reff: " << reff << " um\n";
  out0 << "veff: " << veff << "\n";
  out0 << "cext: " << cext << " um2\n";
  out0 << "csca: " << csca << " um2\n";
  out0 << "walb: " << walb << "\n";
  out0 << "asymm: " << asymm << "\n";
  out0 << "f11: " << f11 << "\n";
  out0 << "f22: " << f22 << "\n";
  out0 << "f33: " << f33 << "\n";
  out0 << "f44: " << f44 << "\n";
  out0 << "f12: " << f12 << "\n";
  out0 << "f34: " << f34 << "\n";
 
  out0 << "Error message: " << (strlen(errmsg) ? errmsg : "None") << "\n";
}
 
// Documentation in header file.
void calc_ssp_random_test(const Verbosity& verbosity) {
  CREATE_OUT0;
  out0 << "======================================================\n";
  out0 << "Test calculation of single scattering data\n";
  out0 << "for randomly oriented, oblate particles\n";
  out0 << "======================================================\n";
 
  SingleScatteringData ssd;
 
  ssd.ptype = PTYPE_TOTAL_RND;
  ssd.f_grid = {230e9, 240e9};
  ssd.T_grid = {220, 250};
  nlinspace(ssd.za_grid, 0, 180, 19);
  nlinspace(ssd.aa_grid, 0, 180, 19);
 
  // Refractive index real and imagenary parts
  // Dimensions: [nf, nT];
  Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
  Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
 
  mrr(0, 0) = 1.78031135;
  mrr(0, 1) = 1.78150475;
  mrr(1, 0) = 1.78037238;
  mrr(1, 1) = 1.78147686;
 
  mri(0, 0) = 0.00278706;
  mri(0, 1) = 0.00507565;
  mri(1, 0) = 0.00287245;
  mri(1, 1) = 0.00523012;
 
  calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 1.5);
 
  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
       << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
 
  out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
  out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
 
  out0 << "======================================================\n";
  out0 << "Test calculation of single scattering data\n";
  out0 << "for randomly oriented, prolate particles\n";
  out0 << "======================================================\n";
 
  calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 0.7);
 
  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
       << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
 
  out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
  out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
}
 
// Documentation in header file.
void calc_ssp_fixed_test(const Verbosity& verbosity) {
  CREATE_OUT0;
  out0 << "======================================================\n";
  out0 << "Test calculation of single scattering data\n";
  out0 << "for oblate particles with fixed orientation\n";
  out0 << "======================================================\n";
 
  SingleScatteringData ssd;
 
  ssd.ptype = PTYPE_AZIMUTH_RND;
  ssd.f_grid = {230e9, 240e9};
  ssd.T_grid = {220, 250};
  nlinspace(ssd.za_grid, 0, 180, 19);
  nlinspace(ssd.aa_grid, 0, 180, 19);
 
  // Refractive index real and imagenary parts
  // Dimensions: [nf, nT];
  Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
  Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
 
  mrr(0, 0) = 1.78031135;
  mrr(0, 1) = 1.78150475;
  mrr(1, 0) = 1.78037238;
  mrr(1, 1) = 1.78147686;
 
  mri(0, 0) = 0.00278706;
  mri(0, 1) = 0.00507565;
  mri(1, 0) = 0.00287245;
  mri(1, 1) = 0.00523012;
 
  calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 1.5);
 
  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
       << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
 
  out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
       << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
 
  out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
 
  out0 << "======================================================\n";
  out0 << "Test calculation of single scattering data\n";
  out0 << "for prolate particles with fixed orientation\n";
  out0 << "======================================================\n";
  calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 0.7);
 
  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
       << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
 
  out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
       << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
 
  out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
}