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