339#ifndef ENABLE_TMATRIX
373 throw std::runtime_error(
374 "This version of ARTS was compiled without T-Matrix support.");
391 throw std::runtime_error(
392 "This version of ARTS was compiled without T-Matrix support.");
408 throw std::runtime_error(
409 "This version of ARTS was compiled without T-Matrix support.");
413 throw std::runtime_error(
414 "This version of ARTS was compiled without T-Matrix support.");
432 ampl_(nmax, lam, thet0, thet, phi0, phi, alpha,
beta, s11, s12, s21, s22);
456 z(0, 0) = 0.5 * (s11 *
conj(s11) + s12 *
conj(s12) + s21 *
conj(s21) +
459 z(0, 1) = 0.5 * (s11 *
conj(s11) - s12 *
conj(s12) + s21 *
conj(s21) -
463 z(0, 3) = (
Complex(0., 1.) * (s11 *
conj(s12) - s22 *
conj(s21))).real();
465 z(1, 0) = 0.5 * (s11 *
conj(s11) + s12 *
conj(s12) - s21 *
conj(s21) -
468 z(1, 1) = 0.5 * (s11 *
conj(s11) - s12 *
conj(s12) - s21 *
conj(s21) +
472 z(1, 3) = (
Complex(0., 1.) * (s11 *
conj(s12) + s22 *
conj(s21))).real();
477 z(2, 3) = (
Complex(0., -1.) * (s11 *
conj(s22) + s21 *
conj(s12))).real();
479 z(3, 0) = (
Complex(0., 1.) * (s21 *
conj(s11) + s22 *
conj(s12))).real();
480 z(3, 1) = (
Complex(0., 1.) * (s21 *
conj(s11) - s22 *
conj(s12))).real();
481 z(3, 2) = (
Complex(0., -1.) * (s22 *
conj(s11) - s12 *
conj(s21))).real();
485static const Numeric GaussLeg6[][3] = {{0.23861918, 0.66120939, 0.93246951},
486 {0.46791393, 0.36076157, 0.17132449}};
488static const Numeric GaussLeg10[][5] = {
489 {0.14887434, 0.43339539, 0.67940957, 0.86506337, 0.97390653},
490 {0.29552422, 0.26926672, 0.21908636, 0.14945135, 0.06667134}};
520 const Numeric c = 0.5 * (alpha_2 + alpha_1);
521 const Numeric m = 0.5 * (alpha_2 - alpha_1);
527 for (
Index i = 0; i < 5; ++i) {
528 const Numeric abscissa = GaussLeg10[0][i];
529 const Numeric weight = GaussLeg10[1][i];
570 const Numeric c = 0.5 * (alpha_2 + alpha_1);
571 const Numeric m = 0.5 * (alpha_2 - alpha_1);
577 for (
Index i = 0; i < 3; ++i) {
578 const Numeric abscissa = GaussLeg6[0][i];
579 const Numeric weight = GaussLeg6[1][i];
629 const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
630 const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
631 const Numeric c_phi = 0.5 * (phi_2 + phi_1);
632 const Numeric m_phi = 0.5 * (phi_2 - phi_1);
634 for (
Index t = 0; t < 5; ++t) {
635 const Numeric abscissa_thet0 = GaussLeg10[0][t];
636 const Numeric weight_thet0 = GaussLeg10[1][t];
638 Matrix phamat_phi(4, 4, 0.);
640 for (
Index p = 0; p < 5; ++p) {
641 const Numeric abscissa_phi = GaussLeg10[0][p];
642 const Numeric weight_phi = GaussLeg10[1][p];
644 Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
651 c_phi + m_phi * abscissa_phi,
654 z *= weight_phi * sin(this_thet0 *
PI / 180.);
663 c_phi - m_phi * abscissa_phi,
666 z *= weight_phi * sin(this_thet0 *
PI / 180.);
669 this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
676 c_phi + m_phi * abscissa_phi,
679 z *= weight_phi * sin(this_thet0 *
PI / 180.);
688 c_phi - m_phi * abscissa_phi,
691 z *= weight_phi * sin(this_thet0 *
PI / 180.);
694 phamat_phi *= m_phi * weight_thet0;
695 phamat += phamat_phi;
740 const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
741 const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
742 const Numeric c_phi = 0.5 * (phi_2 + phi_1);
743 const Numeric m_phi = 0.5 * (phi_2 - phi_1);
748 for (
Index t = 0; t < 3; ++t) {
749 const Numeric abscissa_thet0 = GaussLeg6[0][t];
750 const Numeric weight_thet0 = GaussLeg6[1][t];
754 for (
Index p = 0; p < 3; ++p) {
755 const Numeric abscissa_phi = GaussLeg6[0][p];
756 const Numeric weight_phi = GaussLeg6[1][p];
758 Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
765 c_phi + m_phi * abscissa_phi,
769 z *= weight_phi * sin(this_thet0 *
PI / 180.);
778 c_phi - m_phi * abscissa_phi,
782 z *= weight_phi * sin(this_thet0 *
PI / 180.);
785 this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
792 c_phi + m_phi * abscissa_phi,
796 z *= weight_phi * sin(this_thet0 *
PI / 180.);
805 c_phi - m_phi * abscissa_phi,
809 z *= weight_phi * sin(this_thet0 *
PI / 180.);
812 phamat_phi *= m_phi * weight_thet0;
813 phamat += phamat_phi;
860 const Index quiet = 1) {
866 char errmsg[1024] =
"";
884#pragma omp critical(tmatrix_code)
917 if (strlen(errmsg)) {
918 std::ostringstream os;
919 os <<
"T-Matrix code failed: " << errmsg;
920 throw std::runtime_error(os.str());
951 const Index quiet = 1) {
952 char errmsg[1024] =
"";
957#pragma omp critical(tmatrix_code)
972 if (strlen(errmsg)) {
973 std::ostringstream os;
974 os <<
"T-Matrix code failed: " << errmsg;
975 throw std::runtime_error(os.str());
995 if (ref_index_real.
nrows() != nf)
996 throw std::runtime_error(
997 "Number of rows of refractive index real part must match ssd.f_grid.");
998 if (ref_index_real.
ncols() != nT)
999 throw std::runtime_error(
1000 "Number of cols of refractive index real part must match ssd.T_grid.");
1001 if (ref_index_imag.
nrows() != nf)
1002 throw std::runtime_error(
1003 "Number of rows of refractive index imaginary part must match ssd.f_grid.");
1004 if (ref_index_imag.
ncols() != nT)
1005 throw std::runtime_error(
1006 "Number of cols of refractive index imaginary part must match ssd.T_grid.");
1012 switch (ssd.
ptype) {
1033 Matrix mono_pha_mat_data(nza, 6, NAN);
1036 os <<
"Calculation of SingleScatteringData properties failed for\n\n";
1037 bool anyfailed =
false;
1038#pragma omp critical(tmatrix_ssp)
1039 for (
Index f_index = 0; f_index < nf; ++f_index)
1040 for (
Index T_index = 0; T_index < nT; ++T_index) {
1041 bool thisfailed =
false;
1055 ref_index_real(f_index, T_index),
1056 ref_index_imag(f_index, T_index),
1061 }
catch (
const std::runtime_error& e) {
1064 os <<
"f_grid[" << f_index <<
"] = " << ssd.
f_grid[f_index] <<
"\n"
1065 <<
"T_grid[" << T_index <<
"] = " << ssd.
T_grid[T_index]
1075 mono_pha_mat_data(
joker, 0) = f11;
1076 mono_pha_mat_data(
joker, 1) = f12;
1077 mono_pha_mat_data(
joker, 2) = f22;
1078 mono_pha_mat_data(
joker, 3) = f33;
1079 mono_pha_mat_data(
joker, 4) = f34;
1080 mono_pha_mat_data(
joker, 5) = f44;
1082 mono_pha_mat_data *= csca / 4. /
PI;
1087 ssd.
abs_vec_data(f_index, T_index, 0, 0, 0) = cext - csca;
1094 throw std::runtime_error(os.str());
1118 Tensor5 csca_data(nf, nT, nza, 1, 2);
1120#pragma omp critical(tmatrix_ssp)
1121 for (
Index f_index = 0; f_index < nf; ++f_index) {
1122 const Numeric lam_f = lam[f_index];
1124 for (
Index T_index = 0; T_index < nT; ++T_index) {
1133 ref_index_real(f_index, T_index),
1134 ref_index_imag(f_index, T_index),
1136 }
catch (
const std::runtime_error& e) {
1138 os <<
"Calculation of SingleScatteringData properties failed for\n"
1139 <<
"f_grid[" << f_index <<
"] = " << ssd.
f_grid[f_index] <<
"\n"
1140 <<
"T_grid[" << T_index <<
"] = " << ssd.
T_grid[T_index] <<
"\n"
1142 throw std::runtime_error(os.str());
1146 for (
Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index)
1147 for (
Index aa_index = 0; aa_index < naa; ++aa_index)
1148 for (
Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
1149 if (aspect_ratio < 1.0) {
1206 for (
Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index) {
1208 if (aspect_ratio < 1.0) {
1222 csca_integral /= 180.;
1237 csca_data(f_index, T_index, za_scat_index, 0,
joker) =
1238 csca_integral(
Range(0, 2), 0);
1242 if (aspect_ratio < 1.0) {
1247 for (
Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
1270 K[0] = (
Complex(0., -1.) * (s11 + s22)).real();
1271 K[1] = (
Complex(0., 1.) * (s22 - s11)).real();
1272 K[2] = (s22 - s11).
real();
1279 csca_data *= 2. *
PI *
PI / 32400.;
1287 std::ostringstream os;
1288 os <<
"Only particle types totally_random and azimuthally_random\n"
1289 <<
"are currently supported: " << ssd.
ptype;
1290 throw std::runtime_error(os.str());
1300 out0 <<
"======================================================\n";
1301 out0 <<
"Test for nonspherical particles in a fixed orientation\n";
1302 out0 <<
"Output should match 3rdparty/tmatrix/tmatrix_ampld.ref\n";
1303 out0 <<
"======================================================\n";
1321 char errmsg[1024] =
"";
1324 rat, axi, np, lam, eps, mrr, mri, ddelt, quiet, nmax, csca, cext, errmsg);
1326 out0 <<
"nmax: " << nmax <<
"\n";
1327 out0 <<
"csca: " << csca <<
" um2\n";
1328 out0 <<
"cext: " << cext <<
" um2\n";
1330 out0 <<
"Error message: " << (strlen(errmsg) ? errmsg :
"None") <<
"\n";
1345 ampl_(nmax, lam, thet0, thet, phi0, phi, alpha,
beta, s11, s12, s21, s22);
1347 out0 <<
"AMPLITUDE MATRIX (all in [um]): \n";
1348 out0 <<
"s11: " << s11 <<
"\n";
1349 out0 <<
"s12: " << s12 <<
"\n";
1350 out0 <<
"s21: " << s21 <<
"\n";
1351 out0 <<
"s22: " << s22 <<
"\n";
1357 out0 <<
"PHASE MATRIX (all un [um2]): \n" << z <<
"\n";
1364 out0 <<
"======================================================\n";
1365 out0 <<
"Test for randomly oriented nonspherical particles\n";
1366 out0 <<
"Output should match 3rdparty/tmatrix/tmatrix_tmd.ref\n";
1367 out0 <<
"======================================================\n";
1403 char errmsg[1024] =
"";
1437 out0 <<
"reff: " << reff <<
" um\n";
1438 out0 <<
"veff: " << veff <<
"\n";
1439 out0 <<
"cext: " << cext <<
" um2\n";
1440 out0 <<
"csca: " << csca <<
" um2\n";
1441 out0 <<
"walb: " << walb <<
"\n";
1442 out0 <<
"asymm: " << asymm <<
"\n";
1443 out0 <<
"f11: " << f11 <<
"\n";
1444 out0 <<
"f22: " << f22 <<
"\n";
1445 out0 <<
"f33: " << f33 <<
"\n";
1446 out0 <<
"f44: " << f44 <<
"\n";
1447 out0 <<
"f12: " << f12 <<
"\n";
1448 out0 <<
"f34: " << f34 <<
"\n";
1450 out0 <<
"Error message: " << (strlen(errmsg) ? errmsg :
"None") <<
"\n";
1456 out0 <<
"======================================================\n";
1457 out0 <<
"Test calculation of single scattering data\n";
1458 out0 <<
"for randomly oriented, oblate particles\n";
1459 out0 <<
"======================================================\n";
1464 ssd.
f_grid = {230e9, 240e9};
1474 mrr(0, 0) = 1.78031135;
1475 mrr(0, 1) = 1.78150475;
1476 mrr(1, 0) = 1.78037238;
1477 mrr(1, 1) = 1.78147686;
1479 mri(0, 0) = 0.00278706;
1480 mri(0, 1) = 0.00507565;
1481 mri(1, 0) = 0.00287245;
1482 mri(1, 1) = 0.00523012;
1486 out0 <<
"ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1489 out0 <<
"ssd.ext_mat_data:\n" << ssd.
ext_mat_data <<
"\n\n";
1490 out0 <<
"ssd.abs_vec_data:\n" << ssd.
abs_vec_data <<
"\n\n";
1492 out0 <<
"======================================================\n";
1493 out0 <<
"Test calculation of single scattering data\n";
1494 out0 <<
"for randomly oriented, prolate particles\n";
1495 out0 <<
"======================================================\n";
1499 out0 <<
"ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1502 out0 <<
"ssd.ext_mat_data:\n" << ssd.
ext_mat_data <<
"\n\n";
1503 out0 <<
"ssd.abs_vec_data:\n" << ssd.
abs_vec_data <<
"\n\n";
1509 out0 <<
"======================================================\n";
1510 out0 <<
"Test calculation of single scattering data\n";
1511 out0 <<
"for oblate particles with fixed orientation\n";
1512 out0 <<
"======================================================\n";
1517 ssd.
f_grid = {230e9, 240e9};
1527 mrr(0, 0) = 1.78031135;
1528 mrr(0, 1) = 1.78150475;
1529 mrr(1, 0) = 1.78037238;
1530 mrr(1, 1) = 1.78147686;
1532 mri(0, 0) = 0.00278706;
1533 mri(0, 1) = 0.00507565;
1534 mri(1, 0) = 0.00287245;
1535 mri(1, 1) = 0.00523012;
1539 out0 <<
"ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1542 out0 <<
"ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1545 out0 <<
"abs_vec_data:\n" << ssd.
abs_vec_data <<
"\n\n";
1547 out0 <<
"======================================================\n";
1548 out0 <<
"Test calculation of single scattering data\n";
1549 out0 <<
"for prolate particles with fixed orientation\n";
1550 out0 <<
"======================================================\n";
1553 out0 <<
"ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1556 out0 <<
"ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1559 out0 <<
"abs_vec_data:\n" << ssd.
abs_vec_data <<
"\n\n";
const Numeric SPEED_OF_LIGHT
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
Index nelem() const noexcept
Returns the number of elements.
void resize(Index r, Index c)
Resize function.
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
void resize(Index n)
Resize function.
const Numeric PI
Global constant, pi.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
constexpr Complex conj(Complex c) noexcept
the conjugate of c
std::complex< Numeric > Complex
constexpr Numeric real(Complex c) noexcept
real
Declarations having to do with the four output streams.
Scattering database structure and functions.
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
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)
Calculate SingleScatteringData properties.
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)
Integrate phase matrix over particle orientation angle.
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)
Calculate properties for particles in a fixed orientation.
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 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)
Calculate properties for randomly oriented particles.
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
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)
T-matrix code for nonspherical particles in a fixed orientation.
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
void avgtmatrix_(const Index &nmax)
Perform orientation averaging.
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
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)
T-matrix code for nonspherical particles in a fixed orientation.
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, Numeric &veff, Numeric &cext, Numeric &csca, Numeric &walb, Numeric &asymm, Numeric *f11, Numeric *f22, Numeric *f33, Numeric *f44, Numeric *f12, Numeric *f34, char *errmsg)
T-matrix code for randomly oriented nonspherical particles.
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)
Integrate phase matrix over particle orientation angle.
void ampmat_to_phamat(Matrix &z, const Complex &s11, const Complex &s12, const Complex &s21, const Complex &s22)
Calculate phase matrix from amplitude matrix.
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)
Integrate phase matrix over angles thet0 and phi.
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)
Integrate phase matrix over angles thet0, phi and alpha.
Declarations for the T-Matrix interface.