ARTS 2.5.4 (git: 31ce4f0e)
tmatrix.cc
Go to the documentation of this file.
1/* Copyright (C) 2013 Oliver Lemke
2 *
3 * This program is free software: you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation, either version 3 of the License, or
6 * (at your option) any later version.
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with this program. If not, see <http://www.gnu.org/licenses/>.
15 *
16 */
17
26#include "tmatrix.h"
27#include <cmath>
28#include <cstring>
29#include <stdexcept>
30#include "arts_constants.h"
31#include "matpack_complex.h"
32#include "math_funcs.h"
33#include "matpackI.h"
34#include "messages.h"
35#include "optproperties.h"
36
37void calc_phamat(Matrix& z,
38 const Index& nmax,
39 const Numeric& lam,
40 const Numeric& thet0,
41 const Numeric& thet,
42 const Numeric& phi0,
43 const Numeric& phi,
44 const Numeric& beta,
45 const Numeric& alpha);
46
48 const Complex& s11,
49 const Complex& s12,
50 const Complex& s21,
51 const Complex& s22);
52
54 const Index& nmax,
55 const Numeric& lam,
56 const Numeric& thet0,
57 const Numeric& thet,
58 const Numeric& phi0,
59 const Numeric& phi,
60 const Numeric& beta,
61 const Numeric& alpha_1,
62 const Numeric& alpha_2);
63
65 const Index& nmax,
66 const Numeric& lam,
67 const Numeric& thet0,
68 const Numeric& thet,
69 const Numeric& phi0,
70 const Numeric& phi,
71 const Numeric& beta,
72 const Numeric& alpha_1,
73 const Numeric& alpha_2);
74
76 const Index& nmax,
77 const Numeric& lam,
78 const Numeric& thet0_1,
79 const Numeric& thet0_2,
80 const Numeric& thet,
81 const Numeric& phi0,
82 const Numeric& phi_1,
83 const Numeric& phi_2,
84 const Numeric& beta,
85 const Numeric& alpha);
86
88 const Index& nmax,
89 const Numeric& lam,
90 const Numeric& thet0_1,
91 const Numeric& thet0_2,
92 const Numeric& thet,
93 const Numeric& phi0,
94 const Numeric& phi_1,
95 const Numeric& phi_2,
96 const Numeric& beta,
97 const Numeric& alpha_1,
98 const Numeric& alpha_2);
99
100#ifdef ENABLE_TMATRIX
101extern "C" {
102#endif
210void tmd_(const Numeric& rat,
211 const Index& ndistr,
212 const Numeric& axmax,
213 const Index& npnax,
214 const Numeric& b,
215 const Numeric& gam,
216 const Index& nkmax,
217 const Numeric& eps,
218 const Index& np,
219 const Numeric& lam,
220 const Numeric& mrr,
221 const Numeric& mri,
222 const Numeric& ddelt,
223 const Index& npna,
224 const Index& ndgs,
225 const Numeric& r1rat,
226 const Numeric& r2rat,
227 const Index& quiet,
228 Numeric& reff, // OUT
229 Numeric& veff, // OUT
230 Numeric& cext, // OUT
231 Numeric& csca, // OUT
232 Numeric& walb, // OUT
233 Numeric& asymm, // OUT
234 Numeric* f11, // Array npna
235 Numeric* f22, // Array npna
236 Numeric* f33, // Array npna
237 Numeric* f44, // Array npna
238 Numeric* f12, // Array npna
239 Numeric* f34, // Array npna
240 char* errmsg);
241
275void tmatrix_(const Numeric& rat,
276 const Numeric& axi,
277 const Index& np,
278 const Numeric& lam,
279 const Numeric& eps,
280 const Numeric& mrr,
281 const Numeric& mri,
282 const Numeric& ddelt,
283 const Index& quiet,
284 Index& nmax,
285 Numeric& csca,
286 Numeric& cext,
287 char* errmsg);
288
312void ampl_(const Index& nmax,
313 const Numeric& lam,
314 const Numeric& thet0,
315 const Numeric& thet,
316 const Numeric& phi0,
317 const Numeric& phi,
318 const Numeric& alpha,
319 const Numeric& beta,
320 Complex& s11,
321 Complex& s12,
322 Complex& s21,
323 Complex& s22);
324
333void avgtmatrix_(const Index& nmax);
334#ifdef ENABLE_TMATRIX
335}
336#endif
337
338// Define dummy functions that throw runtime errors if ARTS is
339// compiled without T-Matrix support.
340#ifndef ENABLE_TMATRIX
341
342// T-matrix code for randomly oriented nonspherical particles.
343void tmd_(const Numeric&,
344 const Index&,
345 const Numeric&,
346 const Index&,
347 const Numeric&,
348 const Numeric&,
349 const Index&,
350 const Numeric&,
351 const Index&,
352 const Numeric&,
353 const Numeric&,
354 const Numeric&,
355 const Numeric&,
356 const Index&,
357 const Index&,
358 const Numeric&,
359 const Numeric&,
360 const Index&,
361 Numeric&,
362 Numeric&,
363 Numeric&,
364 Numeric&,
365 Numeric&,
366 Numeric&,
367 Numeric*,
368 Numeric*,
369 Numeric*,
370 Numeric*,
371 Numeric*,
372 Numeric*,
373 char*) {
374 throw std::runtime_error(
375 "This version of ARTS was compiled without T-Matrix support.");
376}
377
378// T-matrix code for nonspherical particles in a fixed orientation
379void tmatrix_(const Numeric&,
380 const Numeric&,
381 const Index&,
382 const Numeric&,
383 const Numeric&,
384 const Numeric&,
385 const Numeric&,
386 const Numeric&,
387 const Index&,
388 Index&,
389 Numeric&,
390 Numeric&,
391 char*) {
392 throw std::runtime_error(
393 "This version of ARTS was compiled without T-Matrix support.");
394}
395
396// T-matrix code for nonspherical particles in a fixed orientation
397void ampl_(const Index&,
398 const Numeric&,
399 const Numeric&,
400 const Numeric&,
401 const Numeric&,
402 const Numeric&,
403 const Numeric&,
404 const Numeric&,
405 Complex&,
406 Complex&,
407 Complex&,
408 Complex&) {
409 throw std::runtime_error(
410 "This version of ARTS was compiled without T-Matrix support.");
411}
412
413void avgtmatrix_(const Index&) {
414 throw std::runtime_error(
415 "This version of ARTS was compiled without T-Matrix support.");
416}
417
418#endif
419
421 const Index& nmax,
422 const Numeric& lam,
423 const Numeric& thet0,
424 const Numeric& thet,
425 const Numeric& phi0,
426 const Numeric& phi,
427 const Numeric& beta,
428 const Numeric& alpha) {
429 Complex s11;
430 Complex s12;
431 Complex s21;
432 Complex s22;
433 ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
434
435 ampmat_to_phamat(z, s11, s12, s21, s22);
436}
437
452 const Complex& s11,
453 const Complex& s12,
454 const Complex& s21,
455 const Complex& s22) {
456 z.resize(4, 4);
457 z(0, 0) = 0.5 * (s11 * conj(s11) + s12 * conj(s12) + s21 * conj(s21) +
458 s22 * conj(s22))
459 .real();
460 z(0, 1) = 0.5 * (s11 * conj(s11) - s12 * conj(s12) + s21 * conj(s21) -
461 s22 * conj(s22))
462 .real();
463 z(0, 2) = (-s11 * conj(s12) - s22 * conj(s21)).real();
464 z(0, 3) = (Complex(0., 1.) * (s11 * conj(s12) - s22 * conj(s21))).real();
465
466 z(1, 0) = 0.5 * (s11 * conj(s11) + s12 * conj(s12) - s21 * conj(s21) -
467 s22 * conj(s22))
468 .real();
469 z(1, 1) = 0.5 * (s11 * conj(s11) - s12 * conj(s12) - s21 * conj(s21) +
470 s22 * conj(s22))
471 .real();
472 z(1, 2) = (-s11 * conj(s12) + s22 * conj(s21)).real();
473 z(1, 3) = (Complex(0., 1.) * (s11 * conj(s12) + s22 * conj(s21))).real();
474
475 z(2, 0) = (-s11 * conj(s21) - s22 * conj(s12)).real();
476 z(2, 1) = (-s11 * conj(s21) + s22 * conj(s12)).real();
477 z(2, 2) = (s11 * conj(s22) + s12 * conj(s21)).real();
478 z(2, 3) = (Complex(0., -1.) * (s11 * conj(s22) + s21 * conj(s12))).real();
479
480 z(3, 0) = (Complex(0., 1.) * (s21 * conj(s11) + s22 * conj(s12))).real();
481 z(3, 1) = (Complex(0., 1.) * (s21 * conj(s11) - s22 * conj(s12))).real();
482 z(3, 2) = (Complex(0., -1.) * (s22 * conj(s11) - s12 * conj(s21))).real();
483 z(3, 3) = (s22 * conj(s11) - s12 * conj(s21)).real();
484}
485
486static const Numeric GaussLeg6[][3] = {{0.23861918, 0.66120939, 0.93246951},
487 {0.46791393, 0.36076157, 0.17132449}};
488
489static const Numeric GaussLeg10[][5] = {
490 {0.14887434, 0.43339539, 0.67940957, 0.86506337, 0.97390653},
491 {0.29552422, 0.26926672, 0.21908636, 0.14945135, 0.06667134}};
492
512 const Index& nmax,
513 const Numeric& lam,
514 const Numeric& thet0,
515 const Numeric& thet,
516 const Numeric& phi0,
517 const Numeric& phi,
518 const Numeric& beta,
519 const Numeric& alpha_1,
520 const Numeric& alpha_2) {
521 const Numeric c = 0.5 * (alpha_2 + alpha_1);
522 const Numeric m = 0.5 * (alpha_2 - alpha_1);
523
524 phamat.resize(4, 4);
525 phamat = 0.;
526 Matrix z;
527
528 for (Index i = 0; i < 5; ++i) {
529 const Numeric abscissa = GaussLeg10[0][i];
530 const Numeric weight = GaussLeg10[1][i];
531
532 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m * abscissa);
533 z *= weight;
534 phamat += z;
535
536 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m * abscissa);
537 z *= weight;
538 phamat += z;
539 }
540 phamat *= m;
541}
542
562 const Index& nmax,
563 const Numeric& lam,
564 const Numeric& thet0,
565 const Numeric& thet,
566 const Numeric& phi0,
567 const Numeric& phi,
568 const Numeric& beta,
569 const Numeric& alpha_1,
570 const Numeric& alpha_2) {
571 const Numeric c = 0.5 * (alpha_2 + alpha_1);
572 const Numeric m = 0.5 * (alpha_2 - alpha_1);
573
574 phamat.resize(4, 4);
575 phamat = 0.;
576 Matrix z;
577
578 for (Index i = 0; i < 3; ++i) {
579 const Numeric abscissa = GaussLeg6[0][i];
580 const Numeric weight = GaussLeg6[1][i];
581
582 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m * abscissa);
583 z *= weight;
584 phamat += z;
585
586 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m * abscissa);
587 z *= weight;
588 phamat += z;
589 }
590 phamat *= m;
591}
592
614 const Index& nmax,
615 const Numeric& lam,
616 const Numeric& thet0_1,
617 const Numeric& thet0_2,
618 const Numeric& thet,
619 const Numeric& phi0,
620 const Numeric& phi_1,
621 const Numeric& phi_2,
622 const Numeric& beta,
623 const Numeric& alpha) {
624 static constexpr Numeric PI=Constant::pi;
625
626 phamat.resize(4, 4);
627 phamat = 0.;
628 Matrix z;
629
630 const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
631 const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
632 const Numeric c_phi = 0.5 * (phi_2 + phi_1);
633 const Numeric m_phi = 0.5 * (phi_2 - phi_1);
634
635 for (Index t = 0; t < 5; ++t) {
636 const Numeric abscissa_thet0 = GaussLeg10[0][t];
637 const Numeric weight_thet0 = GaussLeg10[1][t];
638
639 Matrix phamat_phi(4, 4, 0.);
640
641 for (Index p = 0; p < 5; ++p) {
642 const Numeric abscissa_phi = GaussLeg10[0][p];
643 const Numeric weight_phi = GaussLeg10[1][p];
644
645 Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
646 calc_phamat(z,
647 nmax,
648 lam,
649 this_thet0,
650 thet,
651 phi0,
652 c_phi + m_phi * abscissa_phi,
653 beta,
654 alpha);
655 z *= weight_phi * sin(this_thet0 * PI / 180.);
656 phamat_phi += z;
657
658 calc_phamat(z,
659 nmax,
660 lam,
661 this_thet0,
662 thet,
663 phi0,
664 c_phi - m_phi * abscissa_phi,
665 beta,
666 alpha);
667 z *= weight_phi * sin(this_thet0 * PI / 180.);
668 phamat_phi += z;
669
670 this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
671 calc_phamat(z,
672 nmax,
673 lam,
674 this_thet0,
675 thet,
676 phi0,
677 c_phi + m_phi * abscissa_phi,
678 beta,
679 alpha);
680 z *= weight_phi * sin(this_thet0 * PI / 180.);
681 phamat_phi += z;
682
683 calc_phamat(z,
684 nmax,
685 lam,
686 this_thet0,
687 thet,
688 phi0,
689 c_phi - m_phi * abscissa_phi,
690 beta,
691 alpha);
692 z *= weight_phi * sin(this_thet0 * PI / 180.);
693 phamat_phi += z;
694 }
695 phamat_phi *= m_phi * weight_thet0;
696 phamat += phamat_phi;
697 }
698
699 phamat *= m_thet0;
700}
701
725 const Index& nmax,
726 const Numeric& lam,
727 const Numeric& thet0_1,
728 const Numeric& thet0_2,
729 const Numeric& thet,
730 const Numeric& phi0,
731 const Numeric& phi_1,
732 const Numeric& phi_2,
733 const Numeric& beta,
734 const Numeric& alpha_1,
735 const Numeric& alpha_2) {
736 static constexpr Numeric PI=Constant::pi;
737
738 Matrix z;
739 Matrix phamat_phi(4, 4);
740
741 const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
742 const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
743 const Numeric c_phi = 0.5 * (phi_2 + phi_1);
744 const Numeric m_phi = 0.5 * (phi_2 - phi_1);
745
746 phamat.resize(4, 4);
747 phamat = 0.;
748
749 for (Index t = 0; t < 3; ++t) {
750 const Numeric abscissa_thet0 = GaussLeg6[0][t];
751 const Numeric weight_thet0 = GaussLeg6[1][t];
752
753 phamat_phi = 0.;
754
755 for (Index p = 0; p < 3; ++p) {
756 const Numeric abscissa_phi = GaussLeg6[0][p];
757 const Numeric weight_phi = GaussLeg6[1][p];
758
759 Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
761 nmax,
762 lam,
763 this_thet0,
764 thet,
765 phi0,
766 c_phi + m_phi * abscissa_phi,
767 beta,
768 alpha_1,
769 alpha_2);
770 z *= weight_phi * sin(this_thet0 * PI / 180.);
771 phamat_phi += z;
772
774 nmax,
775 lam,
776 this_thet0,
777 thet,
778 phi0,
779 c_phi - m_phi * abscissa_phi,
780 beta,
781 alpha_1,
782 alpha_2);
783 z *= weight_phi * sin(this_thet0 * PI / 180.);
784 phamat_phi += z;
785
786 this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
788 nmax,
789 lam,
790 this_thet0,
791 thet,
792 phi0,
793 c_phi + m_phi * abscissa_phi,
794 beta,
795 alpha_1,
796 alpha_2);
797 z *= weight_phi * sin(this_thet0 * PI / 180.);
798 phamat_phi += z;
799
801 nmax,
802 lam,
803 this_thet0,
804 thet,
805 phi0,
806 c_phi - m_phi * abscissa_phi,
807 beta,
808 alpha_1,
809 alpha_2);
810 z *= weight_phi * sin(this_thet0 * PI / 180.);
811 phamat_phi += z;
812 }
813 phamat_phi *= m_phi * weight_thet0;
814 phamat += phamat_phi;
815 }
816
817 phamat *= m_thet0;
818}
819
845 Numeric& csca,
846 Vector& f11,
847 Vector& f22,
848 Vector& f33,
849 Vector& f44,
850 Vector& f12,
851 Vector& f34,
852 const Numeric equiv_radius,
853 const Numeric aspect_ratio,
854 const Index np,
855 const Numeric lam,
856 const Numeric ref_index_real,
857 const Numeric ref_index_imag,
858 const Numeric precision,
859 const Index nza,
860 const Index ndgs,
861 const Index quiet = 1) {
862 Numeric reff;
863 Numeric veff;
864 Numeric walb;
865 Numeric asymm;
866
867 char errmsg[1024] = "";
868
869 f11.resize(nza);
870 f11 = NAN;
871 f22.resize(nza);
872 f22 = NAN;
873 f33.resize(nza);
874 f33 = NAN;
875 f44.resize(nza);
876 f44 = NAN;
877 f12.resize(nza);
878 f12 = NAN;
879 f34.resize(nza);
880 f34 = NAN;
881
882 // It is necessary to make sure that the Fortran code is not
883 // called from different threads at the same time. The common
884 // blocks are not threadsafe.
885#pragma omp critical(tmatrix_code)
886 tmd_(1.0,
887 4,
888 equiv_radius,
889 1,
890 0.1,
891 1.0,
892 -1,
893 aspect_ratio,
894 np,
895 lam,
896 ref_index_real,
897 ref_index_imag,
898 precision,
899 nza,
900 ndgs,
901 0.9999999,
902 1.0000001,
903 quiet,
904 reff,
905 veff,
906 cext,
907 csca,
908 walb,
909 asymm,
910 f11.get_c_array(),
911 f22.get_c_array(),
912 f33.get_c_array(),
913 f44.get_c_array(),
914 f12.get_c_array(),
915 f34.get_c_array(),
916 errmsg);
917
918 if (strlen(errmsg)) {
919 std::ostringstream os;
920 os << "T-Matrix code failed: " << errmsg;
921 throw std::runtime_error(os.str());
922 }
923}
924
943 Numeric& csca,
944 Index& nmax,
945 const Numeric equiv_radius,
946 const Numeric aspect_ratio,
947 const Index np,
948 const Numeric lam,
949 const Numeric ref_index_real,
950 const Numeric ref_index_imag,
951 const Numeric precision,
952 const Index quiet = 1) {
953 char errmsg[1024] = "";
954
955 // It is necessary to make sure that the Fortran code is not
956 // called from different threads at the same time. The common
957 // blocks are not threadsafe.
958#pragma omp critical(tmatrix_code)
959 tmatrix_(1.,
960 equiv_radius,
961 np,
962 lam,
963 aspect_ratio,
964 ref_index_real,
965 ref_index_imag,
966 precision,
967 quiet,
968 nmax,
969 csca,
970 cext,
971 errmsg);
972
973 if (strlen(errmsg)) {
974 std::ostringstream os;
975 os << "T-Matrix code failed: " << errmsg;
976 throw std::runtime_error(os.str());
977 }
978}
979
981 ConstMatrixView ref_index_real,
982 ConstMatrixView ref_index_imag,
983 const Numeric equiv_radius,
984 const Index np,
985 const Numeric aspect_ratio,
986 const Numeric precision,
987 const Index ndgs,
988 const Index robust,
989 const Index quiet) {
990 const Index nf = ssd.f_grid.nelem();
991 const Index nT = ssd.T_grid.nelem();
992
993 static constexpr Numeric PI=Constant::pi;
995
996 if (ref_index_real.nrows() != nf)
997 throw std::runtime_error(
998 "Number of rows of refractive index real part must match ssd.f_grid.");
999 if (ref_index_real.ncols() != nT)
1000 throw std::runtime_error(
1001 "Number of cols of refractive index real part must match ssd.T_grid.");
1002 if (ref_index_imag.nrows() != nf)
1003 throw std::runtime_error(
1004 "Number of rows of refractive index imaginary part must match ssd.f_grid.");
1005 if (ref_index_imag.ncols() != nT)
1006 throw std::runtime_error(
1007 "Number of cols of refractive index imaginary part must match ssd.T_grid.");
1008
1009 // The T-Matrix code needs wavelength
1010 Vector lam(nf, SPEED_OF_LIGHT);
1011 lam /= ssd.f_grid;
1012
1013 switch (ssd.ptype) {
1014 case PTYPE_TOTAL_RND: {
1015 const Index nza = ssd.za_grid.nelem();
1016
1017 ssd.pha_mat_data.resize(nf, nT, nza, 1, 1, 1, 6);
1018 ssd.ext_mat_data.resize(nf, nT, 1, 1, 1);
1019 ssd.abs_vec_data.resize(nf, nT, 1, 1, 1);
1020
1021 ssd.pha_mat_data = NAN;
1022 ssd.ext_mat_data = NAN;
1023 ssd.abs_vec_data = NAN;
1024
1025 // Output variables
1026 Numeric cext = NAN;
1027 Numeric csca = NAN;
1028 Vector f11;
1029 Vector f22;
1030 Vector f33;
1031 Vector f44;
1032 Vector f12;
1033 Vector f34;
1034 Matrix mono_pha_mat_data(nza, 6, NAN);
1035
1036 ostringstream os;
1037 os << "Calculation of SingleScatteringData properties failed for\n\n";
1038 bool anyfailed = false;
1039#pragma omp critical(tmatrix_ssp)
1040 for (Index f_index = 0; f_index < nf; ++f_index)
1041 for (Index T_index = 0; T_index < nT; ++T_index) {
1042 bool thisfailed = false;
1043 try {
1045 csca,
1046 f11,
1047 f22,
1048 f33,
1049 f44,
1050 f12,
1051 f34,
1052 equiv_radius,
1053 aspect_ratio,
1054 np,
1055 lam[f_index],
1056 ref_index_real(f_index, T_index),
1057 ref_index_imag(f_index, T_index),
1058 precision,
1059 nza,
1060 ndgs,
1061 quiet);
1062 } catch (const std::runtime_error& e) {
1063 //ostringstream os;
1064 //os << "Calculation of SingleScatteringData properties failed for\n"
1065 os << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
1066 << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index]
1067 << "\n\n";
1068 //<< e.what();
1069 //throw std::runtime_error(os.str());
1070 thisfailed = true;
1071 anyfailed = true;
1072 cout << "\n\n";
1073 }
1074
1075 if (!thisfailed) {
1076 mono_pha_mat_data(joker, 0) = f11;
1077 mono_pha_mat_data(joker, 1) = f12;
1078 mono_pha_mat_data(joker, 2) = f22;
1079 mono_pha_mat_data(joker, 3) = f33;
1080 mono_pha_mat_data(joker, 4) = f34;
1081 mono_pha_mat_data(joker, 5) = f44;
1082
1083 mono_pha_mat_data *= csca / 4. / PI;
1084 ssd.pha_mat_data(f_index, T_index, joker, 0, 0, 0, joker) =
1085 mono_pha_mat_data;
1086
1087 ssd.ext_mat_data(f_index, T_index, 0, 0, 0) = cext;
1088 ssd.abs_vec_data(f_index, T_index, 0, 0, 0) = cext - csca;
1089 }
1090 }
1091 if (anyfailed)
1092 if (robust)
1093 cout << os.str();
1094 else
1095 throw std::runtime_error(os.str());
1096 else {
1097 os << "None\n";
1098 }
1099
1100 break;
1101 }
1102 case PTYPE_AZIMUTH_RND: {
1103 const Index nza = ssd.za_grid.nelem();
1104 const Index naa = ssd.aa_grid.nelem();
1105
1106 ssd.pha_mat_data.resize(nf, nT, nza, naa, nza, 1, 16);
1107 ssd.ext_mat_data.resize(nf, nT, nza, 1, 3);
1108 ssd.abs_vec_data.resize(nf, nT, nza, 1, 2);
1109
1110 ssd.ext_mat_data = NAN;
1111 ssd.pha_mat_data = NAN;
1112 ssd.abs_vec_data = NAN;
1113
1114 // Output variables
1115 Numeric cext = NAN;
1116 Numeric csca = NAN;
1117 Index nmax = -1;
1118
1119 Tensor5 csca_data(nf, nT, nza, 1, 2);
1120
1121#pragma omp critical(tmatrix_ssp)
1122 for (Index f_index = 0; f_index < nf; ++f_index) {
1123 const Numeric lam_f = lam[f_index];
1124
1125 for (Index T_index = 0; T_index < nT; ++T_index) {
1126 try {
1128 csca,
1129 nmax,
1130 equiv_radius,
1131 aspect_ratio,
1132 np,
1133 lam_f,
1134 ref_index_real(f_index, T_index),
1135 ref_index_imag(f_index, T_index),
1136 precision);
1137 } catch (const std::runtime_error& e) {
1138 ostringstream os;
1139 os << "Calculation of SingleScatteringData properties failed for\n"
1140 << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
1141 << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index] << "\n"
1142 << e.what();
1143 throw std::runtime_error(os.str());
1144 }
1145
1146 Matrix phamat;
1147 for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index)
1148 for (Index aa_index = 0; aa_index < naa; ++aa_index)
1149 for (Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
1150 if (aspect_ratio < 1.0) {
1151 // Phase matrix for prolate particles
1153 nmax,
1154 lam_f,
1155 ssd.za_grid[za_inc_index],
1156 ssd.za_grid[za_scat_index],
1157 0.0,
1158 ssd.aa_grid[aa_index],
1159 90.0,
1160 0.0,
1161 180.0);
1162 phamat /= 180.;
1163 } else {
1164 // Phase matrix for oblate particles
1165 calc_phamat(phamat,
1166 nmax,
1167 lam_f,
1168 ssd.za_grid[za_inc_index],
1169 ssd.za_grid[za_scat_index],
1170 0.0,
1171 ssd.aa_grid[aa_index],
1172 0.0,
1173 0.0);
1174 }
1175
1176 ssd.pha_mat_data(f_index,
1177 T_index,
1178 za_scat_index,
1179 aa_index,
1180 za_inc_index,
1181 0,
1182 Range(0, 4)) = phamat(0, joker);
1183 ssd.pha_mat_data(f_index,
1184 T_index,
1185 za_scat_index,
1186 aa_index,
1187 za_inc_index,
1188 0,
1189 Range(4, 4)) = phamat(1, joker);
1190 ssd.pha_mat_data(f_index,
1191 T_index,
1192 za_scat_index,
1193 aa_index,
1194 za_inc_index,
1195 0,
1196 Range(8, 4)) = phamat(2, joker);
1197 ssd.pha_mat_data(f_index,
1198 T_index,
1199 za_scat_index,
1200 aa_index,
1201 za_inc_index,
1202 0,
1203 Range(12, 4)) = phamat(3, joker);
1204 }
1205
1206 // Csca integral
1207 for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index) {
1208 Matrix csca_integral;
1209 if (aspect_ratio < 1.0) {
1210 // Csca for prolate particles
1212 nmax,
1213 lam_f,
1214 0,
1215 180,
1216 ssd.za_grid[za_scat_index],
1217 0.,
1218 0.,
1219 180.,
1220 90.,
1221 0.,
1222 180.);
1223 csca_integral /= 180.;
1224 } else {
1225 // Csca for oblate particles
1226 integrate_phamat_theta0_phi10(csca_integral,
1227 nmax,
1228 lam_f,
1229 0,
1230 180,
1231 ssd.za_grid[za_scat_index],
1232 0.,
1233 0.,
1234 180,
1235 0.,
1236 0.);
1237 }
1238 csca_data(f_index, T_index, za_scat_index, 0, joker) =
1239 csca_integral(Range(0, 2), 0);
1240 }
1241
1242 // Extinction matrix
1243 if (aspect_ratio < 1.0) {
1244 // Average T-Matrix for prolate particles
1245 avgtmatrix_(nmax);
1246 }
1247
1248 for (Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
1249 Complex s11;
1250 Complex s12;
1251 Complex s21;
1252 Complex s22;
1253 VectorView K =
1254 ssd.ext_mat_data(f_index, T_index, za_inc_index, 0, joker);
1255
1256 const Numeric beta = 0.;
1257 const Numeric alpha = 0.;
1258 ampl_(nmax,
1259 lam_f,
1260 ssd.za_grid[za_inc_index],
1261 ssd.za_grid[za_inc_index],
1262 0.,
1263 0.,
1264 alpha,
1265 beta,
1266 s11,
1267 s12,
1268 s21,
1269 s22);
1270
1271 K[0] = (Complex(0., -1.) * (s11 + s22)).real();
1272 K[1] = (Complex(0., 1.) * (s22 - s11)).real();
1273 K[2] = (s22 - s11).real();
1274
1275 K *= lam_f;
1276 }
1277 }
1278 }
1279
1280 csca_data *= 2. * PI * PI / 32400.;
1281 ssd.abs_vec_data =
1282 ssd.ext_mat_data(joker, joker, joker, joker, Range(0, 2));
1283 ssd.abs_vec_data -= csca_data;
1284
1285 break;
1286 }
1287 default: {
1288 std::ostringstream os;
1289 os << "Only particle types totally_random and azimuthally_random\n"
1290 << "are currently supported: " << ssd.ptype;
1291 throw std::runtime_error(os.str());
1292 break;
1293 }
1294 }
1295}
1296
1297// Documentation in header file.
1298void tmatrix_ampld_test(const Verbosity& verbosity) {
1300
1301 out0 << "======================================================\n";
1302 out0 << "Test for nonspherical particles in a fixed orientation\n";
1303 out0 << "Output should match 3rdparty/tmatrix/tmatrix_ampld.ref\n";
1304 out0 << "======================================================\n";
1305
1306 // Same inputs as in example included in original ampld.lp.f
1307 Numeric rat = 1.;
1308 Numeric axi = 10.; //[um]
1309 Index np = -1;
1310 Numeric lam = acos(-1.) * 2.; //[um]
1311 Numeric eps = 0.5;
1312 Numeric mrr = 1.5;
1313 Numeric mri = 0.02;
1314 Numeric ddelt = 0.001;
1315
1316 Index quiet = 1;
1317
1318 // Output variables
1319 Index nmax;
1320 Numeric csca;
1321 Numeric cext;
1322 char errmsg[1024] = "";
1323
1324 tmatrix_(
1325 rat, axi, np, lam, eps, mrr, mri, ddelt, quiet, nmax, csca, cext, errmsg);
1326
1327 out0 << "nmax: " << nmax << "\n";
1328 out0 << "csca: " << csca << " um2\n";
1329 out0 << "cext: " << cext << " um2\n";
1330
1331 out0 << "Error message: " << (strlen(errmsg) ? errmsg : "None") << "\n";
1332
1333 // Same inputs as in example included in original ampld.lp.f
1334 Numeric alpha = 145.;
1335 Numeric beta = 52.;
1336 Numeric thet0 = 56.;
1337 Numeric thet = 65.;
1338 Numeric phi0 = 114.;
1339 Numeric phi = 128.;
1340
1341 // Output variables
1342 Complex s11;
1343 Complex s12;
1344 Complex s21;
1345 Complex s22;
1346 ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
1347
1348 out0 << "AMPLITUDE MATRIX (all in [um]): \n";
1349 out0 << "s11: " << s11 << "\n";
1350 out0 << "s12: " << s12 << "\n";
1351 out0 << "s21: " << s21 << "\n";
1352 out0 << "s22: " << s22 << "\n";
1353
1354 Matrix z;
1355 ampmat_to_phamat(z, s11, s12, s21, s22);
1356 //z *= 1e12; // meter^2 to micron^2 for comparison with original results
1357
1358 out0 << "PHASE MATRIX (all un [um2]): \n" << z << "\n";
1359}
1360
1361// Documentation in header file.
1362void tmatrix_tmd_test(const Verbosity& verbosity) {
1364
1365 out0 << "======================================================\n";
1366 out0 << "Test for randomly oriented nonspherical particles\n";
1367 out0 << "Output should match 3rdparty/tmatrix/tmatrix_tmd.ref\n";
1368 out0 << "======================================================\n";
1369
1370 // Same inputs as in example included in original tmd.lp.f
1371 Numeric rat = 0.5;
1372 Index ndistr = 3;
1373 Numeric axmax = 1.; //[um]
1374 Index npnax = 2;
1375 Numeric b = 0.1;
1376 Numeric gam = 0.5;
1377 Index nkmax = 5;
1378 Numeric eps = 2;
1379 Index np = -1;
1380 Numeric lam = 0.5; //[um]
1381 Numeric mrr = 1.53;
1382 Numeric mri = 0.008;
1383 Numeric ddelt = 0.001;
1384 Index npna = 19;
1385 Index ndgs = 2;
1386 Numeric r1rat = 0.89031;
1387 Numeric r2rat = 1.56538;
1388
1389 Index quiet = 1;
1390
1391 // Output variables
1392 Numeric reff;
1393 Numeric veff;
1394 Numeric cext;
1395 Numeric csca;
1396 Numeric walb;
1397 Numeric asymm;
1398 Vector f11(npna, 0.);
1399 Vector f22(npna, 0.);
1400 Vector f33(npna, 0.);
1401 Vector f44(npna, 0.);
1402 Vector f12(npna, 0.);
1403 Vector f34(npna, 0.);
1404 char errmsg[1024] = "";
1405
1406 tmd_(rat,
1407 ndistr,
1408 axmax,
1409 npnax,
1410 b,
1411 gam,
1412 nkmax,
1413 eps,
1414 np,
1415 lam,
1416 mrr,
1417 mri,
1418 ddelt,
1419 npna,
1420 ndgs,
1421 r1rat,
1422 r2rat,
1423 quiet,
1424 reff,
1425 veff,
1426 cext,
1427 csca,
1428 walb,
1429 asymm,
1430 f11.get_c_array(),
1431 f22.get_c_array(),
1432 f33.get_c_array(),
1433 f44.get_c_array(),
1434 f12.get_c_array(),
1435 f34.get_c_array(),
1436 errmsg);
1437
1438 out0 << "reff: " << reff << " um\n";
1439 out0 << "veff: " << veff << "\n";
1440 out0 << "cext: " << cext << " um2\n";
1441 out0 << "csca: " << csca << " um2\n";
1442 out0 << "walb: " << walb << "\n";
1443 out0 << "asymm: " << asymm << "\n";
1444 out0 << "f11: " << f11 << "\n";
1445 out0 << "f22: " << f22 << "\n";
1446 out0 << "f33: " << f33 << "\n";
1447 out0 << "f44: " << f44 << "\n";
1448 out0 << "f12: " << f12 << "\n";
1449 out0 << "f34: " << f34 << "\n";
1450
1451 out0 << "Error message: " << (strlen(errmsg) ? errmsg : "None") << "\n";
1452}
1453
1454// Documentation in header file.
1455void calc_ssp_random_test(const Verbosity& verbosity) {
1457 out0 << "======================================================\n";
1458 out0 << "Test calculation of single scattering data\n";
1459 out0 << "for randomly oriented, oblate particles\n";
1460 out0 << "======================================================\n";
1461
1463
1464 ssd.ptype = PTYPE_TOTAL_RND;
1465 ssd.f_grid = {230e9, 240e9};
1466 ssd.T_grid = {220, 250};
1467 nlinspace(ssd.za_grid, 0, 180, 19);
1468 nlinspace(ssd.aa_grid, 0, 180, 19);
1469
1470 // Refractive index real and imagenary parts
1471 // Dimensions: [nf, nT];
1472 Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
1473 Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
1474
1475 mrr(0, 0) = 1.78031135;
1476 mrr(0, 1) = 1.78150475;
1477 mrr(1, 0) = 1.78037238;
1478 mrr(1, 1) = 1.78147686;
1479
1480 mri(0, 0) = 0.00278706;
1481 mri(0, 1) = 0.00507565;
1482 mri(1, 0) = 0.00287245;
1483 mri(1, 1) = 0.00523012;
1484
1485 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 1.5);
1486
1487 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1488 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1489
1490 out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
1491 out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1492
1493 out0 << "======================================================\n";
1494 out0 << "Test calculation of single scattering data\n";
1495 out0 << "for randomly oriented, prolate particles\n";
1496 out0 << "======================================================\n";
1497
1498 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 0.7);
1499
1500 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1501 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1502
1503 out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
1504 out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1505}
1506
1507// Documentation in header file.
1508void calc_ssp_fixed_test(const Verbosity& verbosity) {
1510 out0 << "======================================================\n";
1511 out0 << "Test calculation of single scattering data\n";
1512 out0 << "for oblate particles with fixed orientation\n";
1513 out0 << "======================================================\n";
1514
1516
1518 ssd.f_grid = {230e9, 240e9};
1519 ssd.T_grid = {220, 250};
1520 nlinspace(ssd.za_grid, 0, 180, 19);
1521 nlinspace(ssd.aa_grid, 0, 180, 19);
1522
1523 // Refractive index real and imagenary parts
1524 // Dimensions: [nf, nT];
1525 Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
1526 Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
1527
1528 mrr(0, 0) = 1.78031135;
1529 mrr(0, 1) = 1.78150475;
1530 mrr(1, 0) = 1.78037238;
1531 mrr(1, 1) = 1.78147686;
1532
1533 mri(0, 0) = 0.00278706;
1534 mri(0, 1) = 0.00507565;
1535 mri(1, 0) = 0.00287245;
1536 mri(1, 1) = 0.00523012;
1537
1538 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 1.5);
1539
1540 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1541 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1542
1543 out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1544 << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
1545
1546 out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1547
1548 out0 << "======================================================\n";
1549 out0 << "Test calculation of single scattering data\n";
1550 out0 << "for prolate particles with fixed orientation\n";
1551 out0 << "======================================================\n";
1552 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 0.7);
1553
1554 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1555 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1556
1557 out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1558 << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
1559
1560 out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1561}
Constants of physical expressions as constexpr.
constexpr Numeric SPEED_OF_LIGHT
Definition: cia.cc:40
A constant view of a Matrix.
Definition: matpackI.h:1043
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Numeric * get_c_array() const noexcept
Conversion to plain C-array, const-version.
Definition: matpackI.h:622
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
The range class.
Definition: matpackI.h:159
The Tensor5 class.
Definition: matpackV.h:516
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1741
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5488
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
constexpr Numeric PI
Definition: disort.cc:56
#define beta
#define precision
Definition: logic.cc:46
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
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.
#define CREATE_OUT0
Definition: messages.h:203
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Numeric e
Elementary charge convenience name [C].
constexpr Numeric alpha
Fine structure constant convenience name [-].
Scattering database structure and functions.
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:54
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
Definition: tmatrix.cc:1508
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.
Definition: tmatrix.cc:980
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.
Definition: tmatrix.cc:561
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.
Definition: tmatrix.cc:942
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)
Definition: tmatrix.cc:420
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.
Definition: tmatrix.cc:844
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1362
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.
Definition: tmatrix.cc:379
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition: tmatrix.cc:1455
void avgtmatrix_(const Index &nmax)
Perform orientation averaging.
Definition: tmatrix.cc:413
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1298
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.
Definition: tmatrix.cc:397
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.
Definition: tmatrix.cc:343
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.
Definition: tmatrix.cc:511
void ampmat_to_phamat(Matrix &z, const Complex &s11, const Complex &s12, const Complex &s21, const Complex &s22)
Calculate phase matrix from amplitude matrix.
Definition: tmatrix.cc:451
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.
Definition: tmatrix.cc:613
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.
Definition: tmatrix.cc:724
Declarations for the T-Matrix interface.
#define c
#define b