ARTS 2.5.11 (git: 725533f0)
tmatrix.cc
Go to the documentation of this file.
1
9#include "tmatrix.h"
10#include <cmath>
11#include <cstring>
12#include <stdexcept>
13#include "arts_constants.h"
14#include "matpack_complex.h"
15#include "math_funcs.h"
16#include "matpack_data.h"
17#include "messages.h"
18#include "optproperties.h"
19
20void calc_phamat(Matrix& z,
21 const Index& nmax,
22 const Numeric& lam,
23 const Numeric& thet0,
24 const Numeric& thet,
25 const Numeric& phi0,
26 const Numeric& phi,
27 const Numeric& beta,
28 const Numeric& alpha);
29
30void ampmat_to_phamat(Matrix& z,
31 const Complex& s11,
32 const Complex& s12,
33 const Complex& s21,
34 const Complex& s22);
35
36void integrate_phamat_alpha10(Matrix& phamat,
37 const Index& nmax,
38 const Numeric& lam,
39 const Numeric& thet0,
40 const Numeric& thet,
41 const Numeric& phi0,
42 const Numeric& phi,
43 const Numeric& beta,
44 const Numeric& alpha_1,
45 const Numeric& alpha_2);
46
47void integrate_phamat_alpha6(Matrix& phamat,
48 const Index& nmax,
49 const Numeric& lam,
50 const Numeric& thet0,
51 const Numeric& thet,
52 const Numeric& phi0,
53 const Numeric& phi,
54 const Numeric& beta,
55 const Numeric& alpha_1,
56 const Numeric& alpha_2);
57
58void integrate_phamat_theta0_phi10(Matrix& phamat,
59 const Index& nmax,
60 const Numeric& lam,
61 const Numeric& thet0_1,
62 const Numeric& thet0_2,
63 const Numeric& thet,
64 const Numeric& phi0,
65 const Numeric& phi_1,
66 const Numeric& phi_2,
67 const Numeric& beta,
68 const Numeric& alpha);
69
70void integrate_phamat_theta0_phi_alpha6(Matrix& phamat,
71 const Index& nmax,
72 const Numeric& lam,
73 const Numeric& thet0_1,
74 const Numeric& thet0_2,
75 const Numeric& thet,
76 const Numeric& phi0,
77 const Numeric& phi_1,
78 const Numeric& phi_2,
79 const Numeric& beta,
80 const Numeric& alpha_1,
81 const Numeric& alpha_2);
82
83#ifdef ENABLE_TMATRIX
84extern "C" {
85#endif
193void tmd_(const Numeric& rat,
194 const Index& ndistr,
195 const Numeric& axmax,
196 const Index& npnax,
197 const Numeric& b,
198 const Numeric& gam,
199 const Index& nkmax,
200 const Numeric& eps,
201 const Index& np,
202 const Numeric& lam,
203 const Numeric& mrr,
204 const Numeric& mri,
205 const Numeric& ddelt,
206 const Index& npna,
207 const Index& ndgs,
208 const Numeric& r1rat,
209 const Numeric& r2rat,
210 const Index& quiet,
211 Numeric& reff, // OUT
212 Numeric& veff, // OUT
213 Numeric& cext, // OUT
214 Numeric& csca, // OUT
215 Numeric& walb, // OUT
216 Numeric& asymm, // OUT
217 Numeric* f11, // Array npna
218 Numeric* f22, // Array npna
219 Numeric* f33, // Array npna
220 Numeric* f44, // Array npna
221 Numeric* f12, // Array npna
222 Numeric* f34, // Array npna
223 char* errmsg);
224
258void tmatrix_(const Numeric& rat,
259 const Numeric& axi,
260 const Index& np,
261 const Numeric& lam,
262 const Numeric& eps,
263 const Numeric& mrr,
264 const Numeric& mri,
265 const Numeric& ddelt,
266 const Index& quiet,
267 Index& nmax,
268 Numeric& csca,
269 Numeric& cext,
270 char* errmsg);
271
295void ampl_(const Index& nmax,
296 const Numeric& lam,
297 const Numeric& thet0,
298 const Numeric& thet,
299 const Numeric& phi0,
300 const Numeric& phi,
301 const Numeric& alpha,
302 const Numeric& beta,
303 Complex& s11,
304 Complex& s12,
305 Complex& s21,
306 Complex& s22);
307
316void avgtmatrix_(const Index& nmax);
317#ifdef ENABLE_TMATRIX
318}
319#endif
320
321// Define dummy functions that throw runtime errors if ARTS is
322// compiled without T-Matrix support.
323#ifndef ENABLE_TMATRIX
324
325// T-matrix code for randomly oriented nonspherical particles.
326void tmd_(const Numeric&,
327 const Index&,
328 const Numeric&,
329 const Index&,
330 const Numeric&,
331 const Numeric&,
332 const Index&,
333 const Numeric&,
334 const Index&,
335 const Numeric&,
336 const Numeric&,
337 const Numeric&,
338 const Numeric&,
339 const Index&,
340 const Index&,
341 const Numeric&,
342 const Numeric&,
343 const Index&,
344 Numeric&,
345 Numeric&,
346 Numeric&,
347 Numeric&,
348 Numeric&,
349 Numeric&,
350 Numeric*,
351 Numeric*,
352 Numeric*,
353 Numeric*,
354 Numeric*,
355 Numeric*,
356 char*) {
357 throw std::runtime_error(
358 "This version of ARTS was compiled without T-Matrix support.");
359}
360
361// T-matrix code for nonspherical particles in a fixed orientation
362void tmatrix_(const Numeric&,
363 const Numeric&,
364 const Index&,
365 const Numeric&,
366 const Numeric&,
367 const Numeric&,
368 const Numeric&,
369 const Numeric&,
370 const Index&,
371 Index&,
372 Numeric&,
373 Numeric&,
374 char*) {
375 throw std::runtime_error(
376 "This version of ARTS was compiled without T-Matrix support.");
377}
378
379// T-matrix code for nonspherical particles in a fixed orientation
380void ampl_(const Index&,
381 const Numeric&,
382 const Numeric&,
383 const Numeric&,
384 const Numeric&,
385 const Numeric&,
386 const Numeric&,
387 const Numeric&,
388 Complex&,
389 Complex&,
390 Complex&,
391 Complex&) {
392 throw std::runtime_error(
393 "This version of ARTS was compiled without T-Matrix support.");
394}
395
396void avgtmatrix_(const Index&) {
397 throw std::runtime_error(
398 "This version of ARTS was compiled without T-Matrix support.");
399}
400
401#endif
402
403void calc_phamat(Matrix& z,
404 const Index& nmax,
405 const Numeric& lam,
406 const Numeric& thet0,
407 const Numeric& thet,
408 const Numeric& phi0,
409 const Numeric& phi,
410 const Numeric& beta,
411 const Numeric& alpha) {
412 Complex s11;
413 Complex s12;
414 Complex s21;
415 Complex s22;
416 ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
417
418 ampmat_to_phamat(z, s11, s12, s21, s22);
419}
420
434void ampmat_to_phamat(Matrix& z,
435 const Complex& s11,
436 const Complex& s12,
437 const Complex& s21,
438 const Complex& s22) {
439 z.resize(4, 4);
440 z(0, 0) = 0.5 * (s11 * conj(s11) + s12 * conj(s12) + s21 * conj(s21) +
441 s22 * conj(s22))
442 .real();
443 z(0, 1) = 0.5 * (s11 * conj(s11) - s12 * conj(s12) + s21 * conj(s21) -
444 s22 * conj(s22))
445 .real();
446 z(0, 2) = (-s11 * conj(s12) - s22 * conj(s21)).real();
447 z(0, 3) = (Complex(0., 1.) * (s11 * conj(s12) - s22 * conj(s21))).real();
448
449 z(1, 0) = 0.5 * (s11 * conj(s11) + s12 * conj(s12) - s21 * conj(s21) -
450 s22 * conj(s22))
451 .real();
452 z(1, 1) = 0.5 * (s11 * conj(s11) - s12 * conj(s12) - s21 * conj(s21) +
453 s22 * conj(s22))
454 .real();
455 z(1, 2) = (-s11 * conj(s12) + s22 * conj(s21)).real();
456 z(1, 3) = (Complex(0., 1.) * (s11 * conj(s12) + s22 * conj(s21))).real();
457
458 z(2, 0) = (-s11 * conj(s21) - s22 * conj(s12)).real();
459 z(2, 1) = (-s11 * conj(s21) + s22 * conj(s12)).real();
460 z(2, 2) = (s11 * conj(s22) + s12 * conj(s21)).real();
461 z(2, 3) = (Complex(0., -1.) * (s11 * conj(s22) + s21 * conj(s12))).real();
462
463 z(3, 0) = (Complex(0., 1.) * (s21 * conj(s11) + s22 * conj(s12))).real();
464 z(3, 1) = (Complex(0., 1.) * (s21 * conj(s11) - s22 * conj(s12))).real();
465 z(3, 2) = (Complex(0., -1.) * (s22 * conj(s11) - s12 * conj(s21))).real();
466 z(3, 3) = (s22 * conj(s11) - s12 * conj(s21)).real();
467}
468
469static const Numeric GaussLeg6[][3] = {{0.23861918, 0.66120939, 0.93246951},
470 {0.46791393, 0.36076157, 0.17132449}};
471
472static const Numeric GaussLeg10[][5] = {
473 {0.14887434, 0.43339539, 0.67940957, 0.86506337, 0.97390653},
474 {0.29552422, 0.26926672, 0.21908636, 0.14945135, 0.06667134}};
475
494void integrate_phamat_alpha10(Matrix& phamat,
495 const Index& nmax,
496 const Numeric& lam,
497 const Numeric& thet0,
498 const Numeric& thet,
499 const Numeric& phi0,
500 const Numeric& phi,
501 const Numeric& beta,
502 const Numeric& alpha_1,
503 const Numeric& alpha_2) {
504 const Numeric c = 0.5 * (alpha_2 + alpha_1);
505 const Numeric m = 0.5 * (alpha_2 - alpha_1);
506
507 phamat.resize(4, 4);
508 phamat = 0.;
509 Matrix z;
510
511 for (Index i = 0; i < 5; ++i) {
512 const Numeric abscissa = GaussLeg10[0][i];
513 const Numeric weight = GaussLeg10[1][i];
514
515 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m * abscissa);
516 z *= weight;
517 phamat += z;
518
519 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m * abscissa);
520 z *= weight;
521 phamat += z;
522 }
523 phamat *= m;
524}
525
544void integrate_phamat_alpha6(Matrix& phamat,
545 const Index& nmax,
546 const Numeric& lam,
547 const Numeric& thet0,
548 const Numeric& thet,
549 const Numeric& phi0,
550 const Numeric& phi,
551 const Numeric& beta,
552 const Numeric& alpha_1,
553 const Numeric& alpha_2) {
554 const Numeric c = 0.5 * (alpha_2 + alpha_1);
555 const Numeric m = 0.5 * (alpha_2 - alpha_1);
556
557 phamat.resize(4, 4);
558 phamat = 0.;
559 Matrix z;
560
561 for (Index i = 0; i < 3; ++i) {
562 const Numeric abscissa = GaussLeg6[0][i];
563 const Numeric weight = GaussLeg6[1][i];
564
565 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m * abscissa);
566 z *= weight;
567 phamat += z;
568
569 calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m * abscissa);
570 z *= weight;
571 phamat += z;
572 }
573 phamat *= m;
574}
575
597 const Index& nmax,
598 const Numeric& lam,
599 const Numeric& thet0_1,
600 const Numeric& thet0_2,
601 const Numeric& thet,
602 const Numeric& phi0,
603 const Numeric& phi_1,
604 const Numeric& phi_2,
605 const Numeric& beta,
606 const Numeric& alpha) {
607 static constexpr Numeric PI=Constant::pi;
608
609 phamat.resize(4, 4);
610 phamat = 0.;
611 Matrix z;
612
613 const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
614 const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
615 const Numeric c_phi = 0.5 * (phi_2 + phi_1);
616 const Numeric m_phi = 0.5 * (phi_2 - phi_1);
617
618 for (Index t = 0; t < 5; ++t) {
619 const Numeric abscissa_thet0 = GaussLeg10[0][t];
620 const Numeric weight_thet0 = GaussLeg10[1][t];
621
622 Matrix phamat_phi(4, 4, 0.);
623
624 for (Index p = 0; p < 5; ++p) {
625 const Numeric abscissa_phi = GaussLeg10[0][p];
626 const Numeric weight_phi = GaussLeg10[1][p];
627
628 Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
629 calc_phamat(z,
630 nmax,
631 lam,
632 this_thet0,
633 thet,
634 phi0,
635 c_phi + m_phi * abscissa_phi,
636 beta,
637 alpha);
638 z *= weight_phi * sin(this_thet0 * PI / 180.);
639 phamat_phi += z;
640
641 calc_phamat(z,
642 nmax,
643 lam,
644 this_thet0,
645 thet,
646 phi0,
647 c_phi - m_phi * abscissa_phi,
648 beta,
649 alpha);
650 z *= weight_phi * sin(this_thet0 * PI / 180.);
651 phamat_phi += z;
652
653 this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
654 calc_phamat(z,
655 nmax,
656 lam,
657 this_thet0,
658 thet,
659 phi0,
660 c_phi + m_phi * abscissa_phi,
661 beta,
662 alpha);
663 z *= weight_phi * sin(this_thet0 * PI / 180.);
664 phamat_phi += z;
665
666 calc_phamat(z,
667 nmax,
668 lam,
669 this_thet0,
670 thet,
671 phi0,
672 c_phi - m_phi * abscissa_phi,
673 beta,
674 alpha);
675 z *= weight_phi * sin(this_thet0 * PI / 180.);
676 phamat_phi += z;
677 }
678 phamat_phi *= m_phi * weight_thet0;
679 phamat += phamat_phi;
680 }
681
682 phamat *= m_thet0;
683}
684
708 const Index& nmax,
709 const Numeric& lam,
710 const Numeric& thet0_1,
711 const Numeric& thet0_2,
712 const Numeric& thet,
713 const Numeric& phi0,
714 const Numeric& phi_1,
715 const Numeric& phi_2,
716 const Numeric& beta,
717 const Numeric& alpha_1,
718 const Numeric& alpha_2) {
719 static constexpr Numeric PI=Constant::pi;
720
721 Matrix z;
722 Matrix phamat_phi(4, 4);
723
724 const Numeric c_thet0 = 0.5 * (thet0_2 + thet0_1);
725 const Numeric m_thet0 = 0.5 * (thet0_2 - thet0_1);
726 const Numeric c_phi = 0.5 * (phi_2 + phi_1);
727 const Numeric m_phi = 0.5 * (phi_2 - phi_1);
728
729 phamat.resize(4, 4);
730 phamat = 0.;
731
732 for (Index t = 0; t < 3; ++t) {
733 const Numeric abscissa_thet0 = GaussLeg6[0][t];
734 const Numeric weight_thet0 = GaussLeg6[1][t];
735
736 phamat_phi = 0.;
737
738 for (Index p = 0; p < 3; ++p) {
739 const Numeric abscissa_phi = GaussLeg6[0][p];
740 const Numeric weight_phi = GaussLeg6[1][p];
741
742 Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
744 nmax,
745 lam,
746 this_thet0,
747 thet,
748 phi0,
749 c_phi + m_phi * abscissa_phi,
750 beta,
751 alpha_1,
752 alpha_2);
753 z *= weight_phi * sin(this_thet0 * PI / 180.);
754 phamat_phi += z;
755
757 nmax,
758 lam,
759 this_thet0,
760 thet,
761 phi0,
762 c_phi - m_phi * abscissa_phi,
763 beta,
764 alpha_1,
765 alpha_2);
766 z *= weight_phi * sin(this_thet0 * PI / 180.);
767 phamat_phi += z;
768
769 this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
771 nmax,
772 lam,
773 this_thet0,
774 thet,
775 phi0,
776 c_phi + m_phi * abscissa_phi,
777 beta,
778 alpha_1,
779 alpha_2);
780 z *= weight_phi * sin(this_thet0 * PI / 180.);
781 phamat_phi += z;
782
784 nmax,
785 lam,
786 this_thet0,
787 thet,
788 phi0,
789 c_phi - m_phi * abscissa_phi,
790 beta,
791 alpha_1,
792 alpha_2);
793 z *= weight_phi * sin(this_thet0 * PI / 180.);
794 phamat_phi += z;
795 }
796 phamat_phi *= m_phi * weight_thet0;
797 phamat += phamat_phi;
798 }
799
800 phamat *= m_thet0;
801}
802
827void tmatrix_random_orientation(Numeric& cext,
828 Numeric& csca,
829 Vector& f11,
830 Vector& f22,
831 Vector& f33,
832 Vector& f44,
833 Vector& f12,
834 Vector& f34,
835 const Numeric equiv_radius,
836 const Numeric aspect_ratio,
837 const Index np,
838 const Numeric lam,
839 const Numeric ref_index_real,
840 const Numeric ref_index_imag,
841 const Numeric precision,
842 const Index nza,
843 const Index ndgs,
844 const Index quiet = 1) {
845 Numeric reff;
846 Numeric veff;
847 Numeric walb;
848 Numeric asymm;
849
850 char errmsg[1024] = "";
851
852 f11.resize(nza);
853 f11 = NAN;
854 f22.resize(nza);
855 f22 = NAN;
856 f33.resize(nza);
857 f33 = NAN;
858 f44.resize(nza);
859 f44 = NAN;
860 f12.resize(nza);
861 f12 = NAN;
862 f34.resize(nza);
863 f34 = NAN;
864
865 // It is necessary to make sure that the Fortran code is not
866 // called from different threads at the same time. The common
867 // blocks are not threadsafe.
868#pragma omp critical(tmatrix_code)
869 tmd_(1.0,
870 4,
871 equiv_radius,
872 1,
873 0.1,
874 1.0,
875 -1,
876 aspect_ratio,
877 np,
878 lam,
879 ref_index_real,
880 ref_index_imag,
881 precision,
882 nza,
883 ndgs,
884 0.9999999,
885 1.0000001,
886 quiet,
887 reff,
888 veff,
889 cext,
890 csca,
891 walb,
892 asymm,
893 f11.unsafe_data_handle(),
894 f22.unsafe_data_handle(),
895 f33.unsafe_data_handle(),
896 f44.unsafe_data_handle(),
897 f12.unsafe_data_handle(),
898 f34.unsafe_data_handle(),
899 errmsg);
900
901 if (strlen(errmsg)) {
902 std::ostringstream os;
903 os << "T-Matrix code failed: " << errmsg;
904 throw std::runtime_error(os.str());
905 }
906}
907
925void tmatrix_fixed_orientation(Numeric& cext,
926 Numeric& csca,
927 Index& nmax,
928 const Numeric equiv_radius,
929 const Numeric aspect_ratio,
930 const Index np,
931 const Numeric lam,
932 const Numeric ref_index_real,
933 const Numeric ref_index_imag,
934 const Numeric precision,
935 const Index quiet = 1) {
936 char errmsg[1024] = "";
937
938 // It is necessary to make sure that the Fortran code is not
939 // called from different threads at the same time. The common
940 // blocks are not threadsafe.
941#pragma omp critical(tmatrix_code)
942 tmatrix_(1.,
943 equiv_radius,
944 np,
945 lam,
946 aspect_ratio,
947 ref_index_real,
948 ref_index_imag,
949 precision,
950 quiet,
951 nmax,
952 csca,
953 cext,
954 errmsg);
955
956 if (strlen(errmsg)) {
957 std::ostringstream os;
958 os << "T-Matrix code failed: " << errmsg;
959 throw std::runtime_error(os.str());
960 }
961}
962
964 ConstMatrixView ref_index_real,
965 ConstMatrixView ref_index_imag,
966 const Numeric equiv_radius,
967 const Index np,
968 const Numeric aspect_ratio,
969 const Numeric precision,
970 const Index ndgs,
971 const Index robust,
972 const Index quiet) {
973 const Index nf = ssd.f_grid.nelem();
974 const Index nT = ssd.T_grid.nelem();
975
976 static constexpr Numeric PI=Constant::pi;
977 static constexpr Numeric SPEED_OF_LIGHT=Constant::speed_of_light;
978
979 if (ref_index_real.nrows() != nf)
980 throw std::runtime_error(
981 "Number of rows of refractive index real part must match ssd.f_grid.");
982 if (ref_index_real.ncols() != nT)
983 throw std::runtime_error(
984 "Number of cols of refractive index real part must match ssd.T_grid.");
985 if (ref_index_imag.nrows() != nf)
986 throw std::runtime_error(
987 "Number of rows of refractive index imaginary part must match ssd.f_grid.");
988 if (ref_index_imag.ncols() != nT)
989 throw std::runtime_error(
990 "Number of cols of refractive index imaginary part must match ssd.T_grid.");
991
992 // The T-Matrix code needs wavelength
993 Vector lam(nf, SPEED_OF_LIGHT);
994 lam /= ssd.f_grid;
995
996 switch (ssd.ptype) {
997 case PTYPE_TOTAL_RND: {
998 const Index nza = ssd.za_grid.nelem();
999
1000 ssd.pha_mat_data.resize(nf, nT, nza, 1, 1, 1, 6);
1001 ssd.ext_mat_data.resize(nf, nT, 1, 1, 1);
1002 ssd.abs_vec_data.resize(nf, nT, 1, 1, 1);
1003
1004 ssd.pha_mat_data = NAN;
1005 ssd.ext_mat_data = NAN;
1006 ssd.abs_vec_data = NAN;
1007
1008 // Output variables
1009 Numeric cext = NAN;
1010 Numeric csca = NAN;
1011 Vector f11;
1012 Vector f22;
1013 Vector f33;
1014 Vector f44;
1015 Vector f12;
1016 Vector f34;
1017 Matrix mono_pha_mat_data(nza, 6, NAN);
1018
1019 ostringstream os;
1020 os << "Calculation of SingleScatteringData properties failed for\n\n";
1021 bool anyfailed = false;
1022#pragma omp critical(tmatrix_ssp)
1023 for (Index f_index = 0; f_index < nf; ++f_index)
1024 for (Index T_index = 0; T_index < nT; ++T_index) {
1025 bool thisfailed = false;
1026 try {
1028 csca,
1029 f11,
1030 f22,
1031 f33,
1032 f44,
1033 f12,
1034 f34,
1035 equiv_radius,
1036 aspect_ratio,
1037 np,
1038 lam[f_index],
1039 ref_index_real(f_index, T_index),
1040 ref_index_imag(f_index, T_index),
1041 precision,
1042 nza,
1043 ndgs,
1044 quiet);
1045 } catch (const std::runtime_error& e) {
1046 //ostringstream os;
1047 //os << "Calculation of SingleScatteringData properties failed for\n"
1048 os << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
1049 << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index]
1050 << "\n\n";
1051 //<< e.what();
1052 //throw std::runtime_error(os.str());
1053 thisfailed = true;
1054 anyfailed = true;
1055 cout << "\n\n";
1056 }
1057
1058 if (!thisfailed) {
1059 mono_pha_mat_data(joker, 0) = f11;
1060 mono_pha_mat_data(joker, 1) = f12;
1061 mono_pha_mat_data(joker, 2) = f22;
1062 mono_pha_mat_data(joker, 3) = f33;
1063 mono_pha_mat_data(joker, 4) = f34;
1064 mono_pha_mat_data(joker, 5) = f44;
1065
1066 mono_pha_mat_data *= csca / 4. / PI;
1067 ssd.pha_mat_data(f_index, T_index, joker, 0, 0, 0, joker) =
1068 mono_pha_mat_data;
1069
1070 ssd.ext_mat_data(f_index, T_index, 0, 0, 0) = cext;
1071 ssd.abs_vec_data(f_index, T_index, 0, 0, 0) = cext - csca;
1072 }
1073 }
1074 if (anyfailed)
1075 if (robust)
1076 cout << os.str();
1077 else
1078 throw std::runtime_error(os.str());
1079 else {
1080 os << "None\n";
1081 }
1082
1083 break;
1084 }
1085 case PTYPE_AZIMUTH_RND: {
1086 const Index nza = ssd.za_grid.nelem();
1087 const Index naa = ssd.aa_grid.nelem();
1088
1089 ssd.pha_mat_data.resize(nf, nT, nza, naa, nza, 1, 16);
1090 ssd.ext_mat_data.resize(nf, nT, nza, 1, 3);
1091 ssd.abs_vec_data.resize(nf, nT, nza, 1, 2);
1092
1093 ssd.ext_mat_data = NAN;
1094 ssd.pha_mat_data = NAN;
1095 ssd.abs_vec_data = NAN;
1096
1097 // Output variables
1098 Numeric cext = NAN;
1099 Numeric csca = NAN;
1100 Index nmax = -1;
1101
1102 Tensor5 csca_data(nf, nT, nza, 1, 2);
1103
1104#pragma omp critical(tmatrix_ssp)
1105 for (Index f_index = 0; f_index < nf; ++f_index) {
1106 const Numeric lam_f = lam[f_index];
1107
1108 for (Index T_index = 0; T_index < nT; ++T_index) {
1109 try {
1111 csca,
1112 nmax,
1113 equiv_radius,
1114 aspect_ratio,
1115 np,
1116 lam_f,
1117 ref_index_real(f_index, T_index),
1118 ref_index_imag(f_index, T_index),
1119 precision);
1120 } catch (const std::runtime_error& e) {
1121 ostringstream os;
1122 os << "Calculation of SingleScatteringData properties failed for\n"
1123 << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
1124 << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index] << "\n"
1125 << e.what();
1126 throw std::runtime_error(os.str());
1127 }
1128
1129 Matrix phamat;
1130 for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index)
1131 for (Index aa_index = 0; aa_index < naa; ++aa_index)
1132 for (Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
1133 if (aspect_ratio < 1.0) {
1134 // Phase matrix for prolate particles
1136 nmax,
1137 lam_f,
1138 ssd.za_grid[za_inc_index],
1139 ssd.za_grid[za_scat_index],
1140 0.0,
1141 ssd.aa_grid[aa_index],
1142 90.0,
1143 0.0,
1144 180.0);
1145 phamat /= 180.;
1146 } else {
1147 // Phase matrix for oblate particles
1148 calc_phamat(phamat,
1149 nmax,
1150 lam_f,
1151 ssd.za_grid[za_inc_index],
1152 ssd.za_grid[za_scat_index],
1153 0.0,
1154 ssd.aa_grid[aa_index],
1155 0.0,
1156 0.0);
1157 }
1158
1159 ssd.pha_mat_data(f_index,
1160 T_index,
1161 za_scat_index,
1162 aa_index,
1163 za_inc_index,
1164 0,
1165 Range(0, 4)) = phamat(0, joker);
1166 ssd.pha_mat_data(f_index,
1167 T_index,
1168 za_scat_index,
1169 aa_index,
1170 za_inc_index,
1171 0,
1172 Range(4, 4)) = phamat(1, joker);
1173 ssd.pha_mat_data(f_index,
1174 T_index,
1175 za_scat_index,
1176 aa_index,
1177 za_inc_index,
1178 0,
1179 Range(8, 4)) = phamat(2, joker);
1180 ssd.pha_mat_data(f_index,
1181 T_index,
1182 za_scat_index,
1183 aa_index,
1184 za_inc_index,
1185 0,
1186 Range(12, 4)) = phamat(3, joker);
1187 }
1188
1189 // Csca integral
1190 for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index) {
1191 Matrix csca_integral;
1192 if (aspect_ratio < 1.0) {
1193 // Csca for prolate particles
1195 nmax,
1196 lam_f,
1197 0,
1198 180,
1199 ssd.za_grid[za_scat_index],
1200 0.,
1201 0.,
1202 180.,
1203 90.,
1204 0.,
1205 180.);
1206 csca_integral /= 180.;
1207 } else {
1208 // Csca for oblate particles
1209 integrate_phamat_theta0_phi10(csca_integral,
1210 nmax,
1211 lam_f,
1212 0,
1213 180,
1214 ssd.za_grid[za_scat_index],
1215 0.,
1216 0.,
1217 180,
1218 0.,
1219 0.);
1220 }
1221 csca_data(f_index, T_index, za_scat_index, 0, joker) =
1222 csca_integral(Range(0, 2), 0);
1223 }
1224
1225 // Extinction matrix
1226 if (aspect_ratio < 1.0) {
1227 // Average T-Matrix for prolate particles
1228 avgtmatrix_(nmax);
1229 }
1230
1231 for (Index za_inc_index = 0; za_inc_index < nza; ++za_inc_index) {
1232 Complex s11;
1233 Complex s12;
1234 Complex s21;
1235 Complex s22;
1236 VectorView K =
1237 ssd.ext_mat_data(f_index, T_index, za_inc_index, 0, joker);
1238
1239 const Numeric beta = 0.;
1240 const Numeric alpha = 0.;
1241 ampl_(nmax,
1242 lam_f,
1243 ssd.za_grid[za_inc_index],
1244 ssd.za_grid[za_inc_index],
1245 0.,
1246 0.,
1247 alpha,
1248 beta,
1249 s11,
1250 s12,
1251 s21,
1252 s22);
1253
1254 K[0] = (Complex(0., -1.) * (s11 + s22)).real();
1255 K[1] = (Complex(0., 1.) * (s22 - s11)).real();
1256 K[2] = (s22 - s11).real();
1257
1258 K *= lam_f;
1259 }
1260 }
1261 }
1262
1263 csca_data *= 2. * PI * PI / 32400.;
1264 ssd.abs_vec_data =
1265 ssd.ext_mat_data(joker, joker, joker, joker, Range(0, 2));
1266 ssd.abs_vec_data -= csca_data;
1267
1268 break;
1269 }
1270 default: {
1271 std::ostringstream os;
1272 os << "Only particle types totally_random and azimuthally_random\n"
1273 << "are currently supported: " << ssd.ptype;
1274 throw std::runtime_error(os.str());
1275 break;
1276 }
1277 }
1278}
1279
1280// Documentation in header file.
1281void tmatrix_ampld_test(const Verbosity& verbosity) {
1283
1284 out0 << "======================================================\n";
1285 out0 << "Test for nonspherical particles in a fixed orientation\n";
1286 out0 << "Output should match 3rdparty/tmatrix/tmatrix_ampld.ref\n";
1287 out0 << "======================================================\n";
1288
1289 // Same inputs as in example included in original ampld.lp.f
1290 Numeric rat = 1.;
1291 Numeric axi = 10.; //[um]
1292 Index np = -1;
1293 Numeric lam = acos(-1.) * 2.; //[um]
1294 Numeric eps = 0.5;
1295 Numeric mrr = 1.5;
1296 Numeric mri = 0.02;
1297 Numeric ddelt = 0.001;
1298
1299 Index quiet = 1;
1300
1301 // Output variables
1302 Index nmax;
1303 Numeric csca;
1304 Numeric cext;
1305 char errmsg[1024] = "";
1306
1307 tmatrix_(
1308 rat, axi, np, lam, eps, mrr, mri, ddelt, quiet, nmax, csca, cext, errmsg);
1309
1310 out0 << "nmax: " << nmax << "\n";
1311 out0 << "csca: " << csca << " um2\n";
1312 out0 << "cext: " << cext << " um2\n";
1313
1314 out0 << "Error message: " << (strlen(errmsg) ? errmsg : "None") << "\n";
1315
1316 // Same inputs as in example included in original ampld.lp.f
1317 Numeric alpha = 145.;
1318 Numeric beta = 52.;
1319 Numeric thet0 = 56.;
1320 Numeric thet = 65.;
1321 Numeric phi0 = 114.;
1322 Numeric phi = 128.;
1323
1324 // Output variables
1325 Complex s11;
1326 Complex s12;
1327 Complex s21;
1328 Complex s22;
1329 ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
1330
1331 out0 << "AMPLITUDE MATRIX (all in [um]): \n";
1332 out0 << "s11: " << s11 << "\n";
1333 out0 << "s12: " << s12 << "\n";
1334 out0 << "s21: " << s21 << "\n";
1335 out0 << "s22: " << s22 << "\n";
1336
1337 Matrix z;
1338 ampmat_to_phamat(z, s11, s12, s21, s22);
1339 //z *= 1e12; // meter^2 to micron^2 for comparison with original results
1340
1341 out0 << "PHASE MATRIX (all un [um2]): \n" << z << "\n";
1342}
1343
1344// Documentation in header file.
1345void tmatrix_tmd_test(const Verbosity& verbosity) {
1347
1348 out0 << "======================================================\n";
1349 out0 << "Test for randomly oriented nonspherical particles\n";
1350 out0 << "Output should match 3rdparty/tmatrix/tmatrix_tmd.ref\n";
1351 out0 << "======================================================\n";
1352
1353 // Same inputs as in example included in original tmd.lp.f
1354 Numeric rat = 0.5;
1355 Index ndistr = 3;
1356 Numeric axmax = 1.; //[um]
1357 Index npnax = 2;
1358 Numeric b = 0.1;
1359 Numeric gam = 0.5;
1360 Index nkmax = 5;
1361 Numeric eps = 2;
1362 Index np = -1;
1363 Numeric lam = 0.5; //[um]
1364 Numeric mrr = 1.53;
1365 Numeric mri = 0.008;
1366 Numeric ddelt = 0.001;
1367 Index npna = 19;
1368 Index ndgs = 2;
1369 Numeric r1rat = 0.89031;
1370 Numeric r2rat = 1.56538;
1371
1372 Index quiet = 1;
1373
1374 // Output variables
1375 Numeric reff;
1376 Numeric veff;
1377 Numeric cext;
1378 Numeric csca;
1379 Numeric walb;
1380 Numeric asymm;
1381 Vector f11(npna, 0.);
1382 Vector f22(npna, 0.);
1383 Vector f33(npna, 0.);
1384 Vector f44(npna, 0.);
1385 Vector f12(npna, 0.);
1386 Vector f34(npna, 0.);
1387 char errmsg[1024] = "";
1388
1389 tmd_(rat,
1390 ndistr,
1391 axmax,
1392 npnax,
1393 b,
1394 gam,
1395 nkmax,
1396 eps,
1397 np,
1398 lam,
1399 mrr,
1400 mri,
1401 ddelt,
1402 npna,
1403 ndgs,
1404 r1rat,
1405 r2rat,
1406 quiet,
1407 reff,
1408 veff,
1409 cext,
1410 csca,
1411 walb,
1412 asymm,
1413 f11.unsafe_data_handle(),
1414 f22.unsafe_data_handle(),
1415 f33.unsafe_data_handle(),
1416 f44.unsafe_data_handle(),
1417 f12.unsafe_data_handle(),
1418 f34.unsafe_data_handle(),
1419 errmsg);
1420
1421 out0 << "reff: " << reff << " um\n";
1422 out0 << "veff: " << veff << "\n";
1423 out0 << "cext: " << cext << " um2\n";
1424 out0 << "csca: " << csca << " um2\n";
1425 out0 << "walb: " << walb << "\n";
1426 out0 << "asymm: " << asymm << "\n";
1427 out0 << "f11: " << f11 << "\n";
1428 out0 << "f22: " << f22 << "\n";
1429 out0 << "f33: " << f33 << "\n";
1430 out0 << "f44: " << f44 << "\n";
1431 out0 << "f12: " << f12 << "\n";
1432 out0 << "f34: " << f34 << "\n";
1433
1434 out0 << "Error message: " << (strlen(errmsg) ? errmsg : "None") << "\n";
1435}
1436
1437// Documentation in header file.
1438void calc_ssp_random_test(const Verbosity& verbosity) {
1440 out0 << "======================================================\n";
1441 out0 << "Test calculation of single scattering data\n";
1442 out0 << "for randomly oriented, oblate particles\n";
1443 out0 << "======================================================\n";
1444
1446
1447 ssd.ptype = PTYPE_TOTAL_RND;
1448 ssd.f_grid = {230e9, 240e9};
1449 ssd.T_grid = {220, 250};
1450 nlinspace(ssd.za_grid, 0, 180, 19);
1451 nlinspace(ssd.aa_grid, 0, 180, 19);
1452
1453 // Refractive index real and imagenary parts
1454 // Dimensions: [nf, nT];
1455 Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
1456 Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
1457
1458 mrr(0, 0) = 1.78031135;
1459 mrr(0, 1) = 1.78150475;
1460 mrr(1, 0) = 1.78037238;
1461 mrr(1, 1) = 1.78147686;
1462
1463 mri(0, 0) = 0.00278706;
1464 mri(0, 1) = 0.00507565;
1465 mri(1, 0) = 0.00287245;
1466 mri(1, 1) = 0.00523012;
1467
1468 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 1.5);
1469
1470 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1471 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1472
1473 out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
1474 out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1475
1476 out0 << "======================================================\n";
1477 out0 << "Test calculation of single scattering data\n";
1478 out0 << "for randomly oriented, prolate particles\n";
1479 out0 << "======================================================\n";
1480
1481 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 0.7);
1482
1483 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1484 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1485
1486 out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
1487 out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1488}
1489
1490// Documentation in header file.
1491void calc_ssp_fixed_test(const Verbosity& verbosity) {
1493 out0 << "======================================================\n";
1494 out0 << "Test calculation of single scattering data\n";
1495 out0 << "for oblate particles with fixed orientation\n";
1496 out0 << "======================================================\n";
1497
1499
1501 ssd.f_grid = {230e9, 240e9};
1502 ssd.T_grid = {220, 250};
1503 nlinspace(ssd.za_grid, 0, 180, 19);
1504 nlinspace(ssd.aa_grid, 0, 180, 19);
1505
1506 // Refractive index real and imagenary parts
1507 // Dimensions: [nf, nT];
1508 Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
1509 Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
1510
1511 mrr(0, 0) = 1.78031135;
1512 mrr(0, 1) = 1.78150475;
1513 mrr(1, 0) = 1.78037238;
1514 mrr(1, 1) = 1.78147686;
1515
1516 mri(0, 0) = 0.00278706;
1517 mri(0, 1) = 0.00507565;
1518 mri(1, 0) = 0.00287245;
1519 mri(1, 1) = 0.00523012;
1520
1521 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 1.5);
1522
1523 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1524 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1525
1526 out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1527 << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
1528
1529 out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1530
1531 out0 << "======================================================\n";
1532 out0 << "Test calculation of single scattering data\n";
1533 out0 << "for prolate particles with fixed orientation\n";
1534 out0 << "======================================================\n";
1535 calcSingleScatteringDataProperties(ssd, mrr, mri, 200.e-6, -1, 0.7);
1536
1537 out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1538 << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1539
1540 out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1541 << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
1542
1543 out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1544}
Constants of physical expressions as constexpr.
constexpr Numeric SPEED_OF_LIGHT
Definition cia.cc:26
constexpr Numeric PI
Definition disort.cc:46
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Declarations having to do with the four output streams.
#define CREATE_OUT0
Definition messages.h:186
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...
Scattering database structure and functions.
@ PTYPE_AZIMUTH_RND
@ PTYPE_TOTAL_RND
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
Definition tmatrix.cc:1491
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:963
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:544
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:925
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:403
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:827
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition tmatrix.cc:1345
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:362
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition tmatrix.cc:1438
void avgtmatrix_(const Index &nmax)
Perform orientation averaging.
Definition tmatrix.cc:396
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition tmatrix.cc:1281
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:380
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:326
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:494
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:434
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:596
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:707
Declarations for the T-Matrix interface.
#define c
#define b