ARTS  2.2.66
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 <stdexcept>
28 #include <cstring>
29 #include <cmath>
30 #include "complex.h"
31 #include "messages.h"
32 #include "matpackI.h"
33 #include "make_vector.h"
34 #include "math_funcs.h"
35 #include "optproperties.h"
36 
37 
38 void calc_phamat(Matrix& z,
39  const Index& nmax,
40  const Numeric& lam,
41  const Numeric& thet0,
42  const Numeric& thet,
43  const Numeric& phi0,
44  const Numeric& phi,
45  const Numeric& beta,
46  const Numeric& alpha);
47 
48 void ampmat_to_phamat(Matrix& z,
49  const Complex& s11,
50  const Complex& s12,
51  const Complex& s21,
52  const Complex& s22);
53 
54 void integrate_phamat_alpha10(Matrix& phamat,
55  const Index& nmax,
56  const Numeric& lam,
57  const Numeric& thet0,
58  const Numeric& thet,
59  const Numeric& phi0,
60  const Numeric& phi,
61  const Numeric& beta,
62  const Numeric& alpha_1,
63  const Numeric& alpha_2);
64 
65 void integrate_phamat_alpha6(Matrix& phamat,
66  const Index& nmax,
67  const Numeric& lam,
68  const Numeric& thet0,
69  const Numeric& thet,
70  const Numeric& phi0,
71  const Numeric& phi,
72  const Numeric& beta,
73  const Numeric& alpha_1,
74  const Numeric& alpha_2);
75 
76 void integrate_phamat_theta0_phi10(Matrix& phamat,
77  const Index& nmax,
78  const Numeric& lam,
79  const Numeric& thet0_1,
80  const Numeric& thet0_2,
81  const Numeric& thet,
82  const Numeric& phi0,
83  const Numeric& phi_1,
84  const Numeric& phi_2,
85  const Numeric& beta,
86  const Numeric& alpha);
87 
88 void integrate_phamat_theta0_phi_alpha6(Matrix& phamat,
89  const Index& nmax,
90  const Numeric& lam,
91  const Numeric& thet0_1,
92  const Numeric& thet0_2,
93  const Numeric& thet,
94  const Numeric& phi0,
95  const Numeric& phi_1,
96  const Numeric& phi_2,
97  const Numeric& beta,
98  const Numeric& alpha_1,
99  const Numeric& alpha_2);
100 
101 
102 #ifdef ENABLE_TMATRIX
103 extern "C" {
104 #endif
105 
212  void tmd_(const Numeric& rat,
213  const Index& ndistr,
214  const Numeric& axmax,
215  const Index& npnax,
216  const Numeric& b,
217  const Numeric& gam,
218  const Index& nkmax,
219  const Numeric& eps,
220  const Index& np,
221  const Numeric& lam,
222  const Numeric& mrr,
223  const Numeric& mri,
224  const Numeric& ddelt,
225  const Index& npna,
226  const Index& ndgs,
227  const Numeric& r1rat,
228  const Numeric& r2rat,
229  const Index& quiet,
230  Numeric& reff, // OUT
231  Numeric& veff, // OUT
232  Numeric& cext, // OUT
233  Numeric& csca, // OUT
234  Numeric& walb, // OUT
235  Numeric& asymm, // OUT
236  Numeric* f11, // Array npna
237  Numeric* f22, // Array npna
238  Numeric* f33, // Array npna
239  Numeric* f44, // Array npna
240  Numeric* f12, // Array npna
241  Numeric* f34, // Array npna
242  char* errmsg);
243 
277  void tmatrix_(const Numeric& rat,
278  const Numeric& axi,
279  const Index& np,
280  const Numeric& lam,
281  const Numeric& eps,
282  const Numeric& mrr,
283  const Numeric& mri,
284  const Numeric& ddelt,
285  const Index& quiet,
286  Index& nmax,
287  Numeric& csca,
288  Numeric& cext,
289  char* errmsg);
290 
314  void ampl_(const Index& nmax,
315  const Numeric& lam,
316  const Numeric& thet0,
317  const Numeric& thet,
318  const Numeric& phi0,
319  const Numeric& phi,
320  const Numeric& alpha,
321  const Numeric& beta,
322  Complex& s11,
323  Complex& s12,
324  Complex& s21,
325  Complex& s22);
326 
335  void avgtmatrix_(const Index& nmax);
336 #ifdef ENABLE_TMATRIX
337 }
338 #endif
339 
340 
341 // Define dummy functions that throw runtime errors if ARTS is
342 // compiled without T-Matrix support.
343 #ifndef ENABLE_TMATRIX
344 
345 // T-matrix code for randomly oriented nonspherical particles.
346 void tmd_(const Numeric&,
347  const Index&,
348  const Numeric&,
349  const Index&,
350  const Numeric&,
351  const Numeric&,
352  const Index&,
353  const Numeric&,
354  const Index&,
355  const Numeric&,
356  const Numeric&,
357  const Numeric&,
358  const Numeric&,
359  const Index&,
360  const Index&,
361  const Numeric&,
362  const Numeric&,
363  const Index&,
364  Numeric&,
365  Numeric&,
366  Numeric&,
367  Numeric&,
368  Numeric&,
369  Numeric&,
370  Numeric*,
371  Numeric*,
372  Numeric*,
373  Numeric*,
374  Numeric*,
375  Numeric*,
376  char*)
377 {
378  throw std::runtime_error("This version of ARTS was compiled without T-Matrix support.");
379 }
380 
381 
382 // T-matrix code for nonspherical particles in a fixed orientation
383 void tmatrix_(const Numeric&,
384  const Numeric&,
385  const Index&,
386  const Numeric&,
387  const Numeric&,
388  const Numeric&,
389  const Numeric&,
390  const Numeric&,
391  const Index&,
392  Index&,
393  Numeric&,
394  Numeric&,
395  char*)
396 {
397  throw std::runtime_error("This version of ARTS was compiled without T-Matrix support.");
398 }
399 
400 
401 // T-matrix code for nonspherical particles in a fixed orientation
402 void ampl_(const Index&,
403  const Numeric&,
404  const Numeric&,
405  const Numeric&,
406  const Numeric&,
407  const Numeric&,
408  const Numeric&,
409  const Numeric&,
410  Complex&,
411  Complex&,
412  Complex&,
413  Complex&)
414 {
415  throw std::runtime_error("This version of ARTS was compiled without T-Matrix support.");
416 }
417 
418 void avgtmatrix_(const Index&)
419 {
420  throw std::runtime_error("This version of ARTS was compiled without T-Matrix support.");
421 }
422 
423 #endif
424 
425 
426 void calc_phamat(Matrix& z,
427  const Index& nmax,
428  const Numeric& lam,
429  const Numeric& thet0,
430  const Numeric& thet,
431  const Numeric& phi0,
432  const Numeric& phi,
433  const Numeric& beta,
434  const Numeric& alpha)
435 {
436  Complex s11;
437  Complex s12;
438  Complex s21;
439  Complex s22;
440  ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta, s11, s12, s21, s22);
441 
442  ampmat_to_phamat(z, s11, s12, s21, s22);
443 }
444 
445 
459 void ampmat_to_phamat(Matrix& z,
460  const Complex& s11,
461  const Complex& s12,
462  const Complex& s21,
463  const Complex& s22)
464 {
465  z.resize(4, 4);
466  z(0, 0) = 0.5 * ( s11 * conj(s11) + s12 * conj(s12)
467  + s21 * conj(s21) + s22 * conj(s22)).real();
468  z(0, 1) = 0.5 * ( s11 * conj(s11) - s12 * conj(s12)
469  + s21 * conj(s21) - s22 * conj(s22)).real();
470  z(0, 2) = ( -s11 * conj(s12) - s22 * conj(s21)).real();
471  z(0, 3) = (Complex(0., 1.) * ( s11 * conj(s12) - s22 * conj(s21))).real();
472 
473  z(1, 0) = 0.5 * ( s11 * conj(s11) + s12 * conj(s12)
474  - s21 * conj(s21) - s22 * conj(s22)).real();
475  z(1, 1) = 0.5 * ( s11 * conj(s11) - s12 * conj(s12)
476  - s21 * conj(s21) + s22 * conj(s22)).real();
477  z(1, 2) = ( -s11 * conj(s12) + s22 * conj(s21)).real();
478  z(1, 3) = (Complex(0., 1.) * ( s11 * conj(s12) + s22 * conj(s21))).real();
479 
480  z(2, 0) = ( -s11 * conj(s21) - s22 * conj(s12)).real();
481  z(2, 1) = ( -s11 * conj(s21) + s22 * conj(s12)).real();
482  z(2, 2) = ( s11 * conj(s22) + s12 * conj(s21)).real();
483  z(2, 3) = (Complex(0., -1.) * ( s11 * conj(s22) + s21 * conj(s12))).real();
484 
485  z(3, 0) = (Complex(0., 1.) * ( s21 * conj(s11) + s22 * conj(s12))).real();
486  z(3, 1) = (Complex(0., 1.) * ( s21 * conj(s11) - s22 * conj(s12))).real();
487  z(3, 2) = (Complex(0., -1.) * ( s22 * conj(s11) - s12 * conj(s21))).real();
488  z(3, 3) = ( s22 * conj(s11) - s12 * conj(s21)).real();
489 
490  z *= 1e-12; // micron^2 to meter^2
491 }
492 
493 
494 static const Numeric GaussLeg6[][3] = {
495  {0.23861918, 0.66120939, 0.93246951},
496  {0.46791393, 0.36076157, 0.17132449}
497 };
498 
499 static const Numeric GaussLeg10[][5] = {
500  {0.14887434, 0.43339539, 0.67940957, 0.86506337, 0.97390653},
501  {0.29552422, 0.26926672, 0.21908636, 0.14945135, 0.06667134}
502 };
503 
504 
523 void integrate_phamat_alpha10(Matrix& phamat,
524  const Index& nmax,
525  const Numeric& lam,
526  const Numeric& thet0,
527  const Numeric& thet,
528  const Numeric& phi0,
529  const Numeric& phi,
530  const Numeric& beta,
531  const Numeric& alpha_1,
532  const Numeric& alpha_2)
533 {
534  const Numeric c = 0.5 * (alpha_2+alpha_1);
535  const Numeric m = 0.5 * (alpha_2-alpha_1);
536 
537  phamat.resize(4, 4);
538  phamat = 0.;
539  Matrix z;
540 
541  for (Index i = 0; i < 5; ++i)
542  {
543  const Numeric abscissa = GaussLeg10[0][i];
544  const Numeric weight = GaussLeg10[1][i];
545 
546  calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m*abscissa);
547  z *= weight;
548  phamat += z;
549 
550  calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m*abscissa);
551  z *= weight;
552  phamat += z;
553  }
554  phamat *= m;
555 }
556 
557 
576 void integrate_phamat_alpha6(Matrix& phamat,
577  const Index& nmax,
578  const Numeric& lam,
579  const Numeric& thet0,
580  const Numeric& thet,
581  const Numeric& phi0,
582  const Numeric& phi,
583  const Numeric& beta,
584  const Numeric& alpha_1,
585  const Numeric& alpha_2)
586 {
587  const Numeric c = 0.5 * (alpha_2+alpha_1);
588  const Numeric m = 0.5 * (alpha_2-alpha_1);
589 
590  phamat.resize(4, 4);
591  phamat = 0.;
592  Matrix z;
593 
594  for (Index i = 0; i < 3; ++i)
595  {
596  const Numeric abscissa = GaussLeg6[0][i];
597  const Numeric weight = GaussLeg6[1][i];
598 
599  calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c + m*abscissa);
600  z *= weight;
601  phamat += z;
602 
603  calc_phamat(z, nmax, lam, thet0, thet, phi0, phi, beta, c - m*abscissa);
604  z *= weight;
605  phamat += z;
606  }
607  phamat *= m;
608 }
609 
610 
631 void integrate_phamat_theta0_phi10(Matrix& phamat,
632  const Index& nmax,
633  const Numeric& lam,
634  const Numeric& thet0_1,
635  const Numeric& thet0_2,
636  const Numeric& thet,
637  const Numeric& phi0,
638  const Numeric& phi_1,
639  const Numeric& phi_2,
640  const Numeric& beta,
641  const Numeric& alpha)
642 {
643  extern const Numeric PI;
644 
645  phamat.resize(4, 4);
646  phamat = 0.;
647  Matrix z;
648 
649  const Numeric c_thet0 = 0.5 * (thet0_2+thet0_1);
650  const Numeric m_thet0 = 0.5 * (thet0_2-thet0_1);
651  const Numeric c_phi = 0.5 * (phi_2+phi_1);
652  const Numeric m_phi = 0.5 * (phi_2-phi_1);
653 
654  for (Index t = 0; t < 5; ++t)
655  {
656  const Numeric abscissa_thet0 = GaussLeg10[0][t];
657  const Numeric weight_thet0 = GaussLeg10[1][t];
658 
659  Matrix phamat_phi(4, 4, 0.);
660 
661  for (Index p = 0; p < 5; ++p)
662  {
663  const Numeric abscissa_phi = GaussLeg10[0][p];
664  const Numeric weight_phi = GaussLeg10[1][p];
665 
666  Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
667  calc_phamat(z, nmax, lam, this_thet0, thet, phi0, c_phi + m_phi * abscissa_phi, beta, alpha);
668  z *= weight_phi * sin(this_thet0*PI/180.);
669  phamat_phi += z;
670 
671  calc_phamat(z, nmax, lam, this_thet0, thet, phi0, c_phi - m_phi * abscissa_phi, beta, alpha);
672  z *= weight_phi * sin(this_thet0*PI/180.);
673  phamat_phi += z;
674 
675  this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
676  calc_phamat(z, nmax, lam, this_thet0, thet, phi0, c_phi + m_phi * abscissa_phi, beta, alpha);
677  z *= weight_phi * sin(this_thet0*PI/180.);
678  phamat_phi += z;
679 
680  calc_phamat(z, nmax, lam, this_thet0, thet, phi0, c_phi - m_phi * abscissa_phi, beta, alpha);
681  z *= weight_phi * sin(this_thet0*PI/180.);
682  phamat_phi += z;
683  }
684  phamat_phi *= m_phi * weight_thet0;
685  phamat += phamat_phi;
686  }
687 
688  phamat *= m_thet0;
689 }
690 
691 
714 void integrate_phamat_theta0_phi_alpha6(Matrix& phamat,
715  const Index& nmax,
716  const Numeric& lam,
717  const Numeric& thet0_1,
718  const Numeric& thet0_2,
719  const Numeric& thet,
720  const Numeric& phi0,
721  const Numeric& phi_1,
722  const Numeric& phi_2,
723  const Numeric& beta,
724  const Numeric& alpha_1,
725  const Numeric& alpha_2)
726 {
727  extern const Numeric PI;
728 
729  Matrix z;
730  Matrix phamat_phi(4, 4);
731 
732 
733  const Numeric c_thet0 = 0.5 * (thet0_2+thet0_1);
734  const Numeric m_thet0 = 0.5 * (thet0_2-thet0_1);
735  const Numeric c_phi = 0.5 * (phi_2+phi_1);
736  const Numeric m_phi = 0.5 * (phi_2-phi_1);
737 
738  phamat.resize(4, 4);
739  phamat = 0.;
740 
741  for (Index t = 0; t < 3; ++t)
742  {
743  const Numeric abscissa_thet0 = GaussLeg6[0][t];
744  const Numeric weight_thet0 = GaussLeg6[1][t];
745 
746  phamat_phi = 0.;
747 
748  for (Index p = 0; p < 3; ++p)
749  {
750  const Numeric abscissa_phi = GaussLeg6[0][p];
751  const Numeric weight_phi = GaussLeg6[1][p];
752 
753  Numeric this_thet0 = c_thet0 + m_thet0 * abscissa_thet0;
754  integrate_phamat_alpha6(z, nmax, lam, this_thet0, thet, phi0, c_phi + m_phi * abscissa_phi, beta, alpha_1, alpha_2);
755  z *= weight_phi * sin(this_thet0*PI/180.);
756  phamat_phi += z;
757 
758  integrate_phamat_alpha6(z, nmax, lam, this_thet0, thet, phi0, c_phi - m_phi * abscissa_phi, beta, alpha_1, alpha_2);
759  z *= weight_phi * sin(this_thet0*PI/180.);
760  phamat_phi += z;
761 
762  this_thet0 = c_thet0 - m_thet0 * abscissa_thet0;
763  integrate_phamat_alpha6(z, nmax, lam, this_thet0, thet, phi0, c_phi + m_phi * abscissa_phi, beta, alpha_1, alpha_2);
764  z *= weight_phi * sin(this_thet0*PI/180.);
765  phamat_phi += z;
766 
767  integrate_phamat_alpha6(z, nmax, lam, this_thet0, thet, phi0, c_phi - m_phi * abscissa_phi, beta, alpha_1, alpha_2);
768  z *= weight_phi * sin(this_thet0*PI/180.);
769  phamat_phi += z;
770  }
771  phamat_phi *= m_phi * weight_thet0;
772  phamat += phamat_phi;
773  }
774 
775  phamat *= m_thet0;
776 }
777 
778 
802 void tmatrix_random_orientation(Numeric& cext,
803  Numeric& csca,
804  Vector& f11,
805  Vector& f22,
806  Vector& f33,
807  Vector& f44,
808  Vector& f12,
809  Vector& f34,
810  const Numeric equiv_radius,
811  const Numeric aspect_ratio,
812  const Index np,
813  const Numeric lam,
814  const Numeric ref_index_real,
815  const Numeric ref_index_imag,
816  const Numeric precision,
817  const Index nza,
818  const Index quiet = 1)
819 {
820  Numeric reff;
821  Numeric veff;
822  Numeric walb;
823  Numeric asymm;
824 
825  char errmsg[1024] = "";
826 
827  f11.resize(nza); f11 = NAN;
828  f22.resize(nza); f22 = NAN;
829  f33.resize(nza); f33 = NAN;
830  f44.resize(nza); f44 = NAN;
831  f12.resize(nza); f12 = NAN;
832  f34.resize(nza); f34 = NAN;
833 
834  // It is necessary to make sure that the Fortran code is not
835  // called from different threads at the same time. The common
836  // blocks are not threadsafe.
837 #pragma omp critical(tmatrix_code)
838  tmd_(1.0,
839  4,
840  equiv_radius,
841  1,
842  0.1,
843  1.0,
844  -1,
845  aspect_ratio,
846  np,
847  lam,
848  ref_index_real,
849  ref_index_imag,
850  precision,
851  nza,
852  2,
853  0.9999999,
854  1.0000001,
855  quiet,
856  reff,
857  veff,
858  cext,
859  csca,
860  walb,
861  asymm,
862  f11.get_c_array(),
863  f22.get_c_array(),
864  f33.get_c_array(),
865  f44.get_c_array(),
866  f12.get_c_array(),
867  f34.get_c_array(),
868  errmsg);
869 
870  if (strlen(errmsg))
871  {
872  std::ostringstream os;
873  os << "T-Matrix code failed: " << errmsg;
874  throw std::runtime_error(os.str());
875  }
876 }
877 
878 
896 void tmatrix_fixed_orientation(Numeric& cext,
897  Numeric& csca,
898  Index& nmax,
899  const Numeric equiv_radius,
900  const Numeric aspect_ratio,
901  const Index np,
902  const Numeric lam,
903  const Numeric ref_index_real,
904  const Numeric ref_index_imag,
905  const Numeric precision,
906  const Index quiet = 1)
907 {
908  char errmsg[1024] = "";
909 
910  // It is necessary to make sure that the Fortran code is not
911  // called from different threads at the same time. The common
912  // blocks are not threadsafe.
913 #pragma omp critical(tmatrix_code)
914  tmatrix_(1., equiv_radius, np, lam, aspect_ratio,
915  ref_index_real, ref_index_imag, precision, quiet,
916  nmax, csca, cext, errmsg);
917 
918  if (strlen(errmsg))
919  {
920  std::ostringstream os;
921  os << "T-Matrix code failed: " << errmsg;
922  throw std::runtime_error(os.str());
923  }
924 }
925 
926 
927 void calcSingleScatteringDataProperties(SingleScatteringData& ssd,
928  ConstMatrixView ref_index_real,
929  ConstMatrixView ref_index_imag,
930  const Numeric equiv_radius,
931  const Index np,
932  const Numeric aspect_ratio,
933  const Numeric precision)
934 {
935  const Index nf = ssd.f_grid.nelem();
936  const Index nT = ssd.T_grid.nelem();
937 
938  const Index nza = ssd.za_grid.nelem();
939 
940  extern const Numeric PI;
941  extern const Numeric SPEED_OF_LIGHT;
942 
943  if (ref_index_real.nrows() != nf)
944  throw std::runtime_error("Number of rows of refractive index real part must match ssd.f_grid.");
945  if (ref_index_real.ncols() != nT)
946  throw std::runtime_error("Number of cols of refractive index real part must match ssd.T_grid.");
947  if (ref_index_imag.nrows() != nf)
948  throw std::runtime_error("Number of rows of refractive index imaginary part must match ssd.f_grid.");
949  if (ref_index_imag.ncols() != nT)
950  throw std::runtime_error("Number of cols of refractive index imaginary part must match ssd.T_grid.");
951 
952  // The T-Matrix codes needs wavelength in microns
953  Vector lam(nf, 1e6*SPEED_OF_LIGHT);
954  lam /= ssd.f_grid;
955 
956  switch (ssd.particle_type) {
957  case PARTICLE_TYPE_MACROS_ISO:
958  {
959  ssd.pha_mat_data.resize(nf, nT, nza, 1, 1, 1, 6);
960  ssd.ext_mat_data.resize(nf, nT, 1, 1, 1);
961  ssd.abs_vec_data.resize(nf, nT, 1, 1, 1);
962 
963  ssd.pha_mat_data = NAN;
964  ssd.ext_mat_data = NAN;
965  ssd.abs_vec_data = NAN;
966 
967  // Output variables
968  Numeric cext = NAN;
969  Numeric csca = NAN;
970  Vector f11;
971  Vector f22;
972  Vector f33;
973  Vector f44;
974  Vector f12;
975  Vector f34;
976  Matrix mono_pha_mat_data(nza, 6, NAN);
977 
978 #pragma omp critical(tmatrix_ssp)
979  for (Index f_index = 0; f_index < nf; ++f_index)
980  for (Index T_index = 0; T_index < nT; ++T_index)
981  {
982  try {
983  tmatrix_random_orientation
984  (cext, csca,
985  f11, f22, f33, f44, f12, f34,
986  equiv_radius, aspect_ratio, np, lam[f_index],
987  ref_index_real(f_index, T_index),
988  ref_index_imag(f_index, T_index),
989  precision,
990  nza);
991  } catch (std::runtime_error e) {
992  ostringstream os;
993  os << "Calculation of SingleScatteringData properties failed for\n"
994  << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
995  << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index] << "\n"
996  << e.what();
997  throw std::runtime_error(os.str());
998  }
999 
1000  cext *= 1e-12; // um^2 -> m^2
1001  csca *= 1e-12;
1002 
1003  mono_pha_mat_data(joker, 0) = f11;
1004  mono_pha_mat_data(joker, 1) = f12;
1005  mono_pha_mat_data(joker, 2) = f22;
1006  mono_pha_mat_data(joker, 3) = f33;
1007  mono_pha_mat_data(joker, 4) = f34;
1008  mono_pha_mat_data(joker, 5) = f44;
1009 
1010  mono_pha_mat_data *= csca/4./PI;
1011  ssd.pha_mat_data(f_index, T_index, joker, 0, 0, 0, joker) = mono_pha_mat_data;
1012 
1013  ssd.ext_mat_data(f_index, T_index, 0, 0, 0) = cext;
1014  ssd.abs_vec_data(f_index, T_index, 0, 0, 0) = cext - csca;
1015  }
1016 
1017  break;
1018  }
1019  case PARTICLE_TYPE_HORIZ_AL:
1020  {
1021  Index nza_inc = (nza-1)/2 + 1;
1022  const Index naa = ssd.aa_grid.nelem();
1023 
1024  ssd.ext_mat_data.resize(nf, nT, nza_inc, 1, 3);
1025  ssd.pha_mat_data.resize(nf, nT, nza, naa, nza_inc, 1, 16);
1026  ssd.abs_vec_data.resize(nf, nza_inc, 1, 2, 1);
1027 
1028  ssd.ext_mat_data = NAN;
1029  ssd.pha_mat_data = NAN;
1030  ssd.abs_vec_data = NAN;
1031 
1032  // Output variables
1033  Numeric cext = NAN;
1034  Numeric csca = NAN;
1035  Index nmax = -1;
1036 
1037  Tensor5 csca_data(nf, nT, nza_inc, 1, 2);
1038 
1039 #pragma omp critical(tmatrix_ssp)
1040  for (Index f_index = 0; f_index < nf; ++f_index)
1041  {
1042  const Numeric lam_f = lam[f_index];
1043 
1044  for (Index T_index = 0; T_index < nT; ++T_index)
1045  {
1046  try {
1047  tmatrix_fixed_orientation
1048  (cext, csca, nmax,
1049  equiv_radius, aspect_ratio, np, lam_f,
1050  ref_index_real(f_index, T_index),
1051  ref_index_imag(f_index, T_index),
1052  precision);
1053  } catch (std::runtime_error e) {
1054  ostringstream os;
1055  os << "Calculation of SingleScatteringData properties failed for\n"
1056  << "f_grid[" << f_index << "] = " << ssd.f_grid[f_index] << "\n"
1057  << "T_grid[" << T_index << "] = " << ssd.T_grid[T_index] << "\n"
1058  << e.what();
1059  throw std::runtime_error(os.str());
1060  }
1061 
1062 
1063  Matrix phamat;
1064  for (Index za_scat_index = 0; za_scat_index < nza; ++za_scat_index)
1065  for (Index aa_index = 0; aa_index < naa; ++aa_index)
1066  for (Index za_inc_index = 0; za_inc_index < nza_inc; ++za_inc_index)
1067  {
1068  if (aspect_ratio < 1.0)
1069  {
1070  // Phase matrix for prolate particles
1071  integrate_phamat_alpha10(phamat,
1072  nmax, lam_f,
1073  ssd.za_grid[za_inc_index],
1074  ssd.za_grid[za_scat_index],
1075  0.0,
1076  ssd.aa_grid[aa_index],
1077  90.0, 0.0, 180.0);
1078  phamat /= 180.;
1079  }
1080  else
1081  {
1082  // Phase matrix for oblate particles
1083  calc_phamat(phamat,
1084  nmax, lam_f,
1085  ssd.za_grid[za_inc_index],
1086  ssd.za_grid[za_scat_index],
1087  0.0,
1088  ssd.aa_grid[aa_index],
1089  0.0, 0.0);
1090  }
1091 
1092  ssd.pha_mat_data(f_index, T_index,
1093  za_scat_index, aa_index, za_inc_index,
1094  0, Range(0, 4)) = phamat(0, joker);
1095  ssd.pha_mat_data(f_index, T_index,
1096  za_scat_index, aa_index, za_inc_index,
1097  0, Range(4, 4)) = phamat(1, joker);
1098  ssd.pha_mat_data(f_index, T_index,
1099  za_scat_index, aa_index, za_inc_index,
1100  0, Range(8, 4)) = phamat(2, joker);
1101  ssd.pha_mat_data(f_index, T_index,
1102  za_scat_index, aa_index, za_inc_index,
1103  0, Range(12, 4)) = phamat(3, joker);
1104  }
1105 
1106  // Csca integral
1107  for (Index za_scat_index = 0; za_scat_index < nza_inc; ++za_scat_index)
1108  {
1109  Matrix csca_integral;
1110  if (aspect_ratio < 1.0)
1111  {
1112  // Csca for prolate particles
1113  integrate_phamat_theta0_phi_alpha6(csca_integral, nmax, lam_f,
1114  0, 180,
1115  ssd.za_grid[za_scat_index],
1116  0.,
1117  0., 180.,
1118  90.,
1119  0., 180.);
1120  csca_integral /= 180.;
1121  }
1122  else
1123  {
1124  // Csca for oblate particles
1125  integrate_phamat_theta0_phi10(csca_integral, nmax, lam_f,
1126  0, 180,
1127  ssd.za_grid[za_scat_index],
1128  0.,
1129  0., 180,
1130  0., 0.);
1131  }
1132  csca_data(f_index, T_index, za_scat_index, 0, joker) = csca_integral(Range(0, 2), 0);
1133  }
1134 
1135  // Extinction matrix
1136  if (aspect_ratio < 1.0)
1137  {
1138  // Average T-Matrix for prolate particles
1139  avgtmatrix_(nmax);
1140  }
1141 
1142  for (Index za_inc_index = 0; za_inc_index < nza_inc; ++za_inc_index)
1143  {
1144  Complex s11;
1145  Complex s12;
1146  Complex s21;
1147  Complex s22;
1148  VectorView K = ssd.ext_mat_data(f_index, T_index, za_inc_index, 0, joker);
1149 
1150  const Numeric beta = 0.;
1151  const Numeric alpha = 0.;
1152  ampl_(nmax, lam_f,
1153  ssd.za_grid[za_inc_index],
1154  ssd.za_grid[za_inc_index],
1155  0., 0.,
1156  alpha, beta,
1157  s11, s12, s21, s22);
1158 
1159 
1160  K[0] = (Complex(0., -1.) * (s11+s22)).real();
1161  K[1] = (Complex(0., 1.) * (s22-s11)).real();
1162  K[2] = (s22-s11).real();
1163 
1164  K *= lam_f * 1e-12; // micron^2 to meter^2
1165 
1166  }
1167  }
1168  }
1169 
1170  csca_data *= 2.*PI*PI/32400.;
1171  ssd.abs_vec_data = ssd.ext_mat_data(joker, joker, joker, joker, Range(0, 2));
1172  ssd.abs_vec_data -= csca_data;
1173 
1174  break;
1175  }
1176  default:
1177  {
1178  std::ostringstream os;
1179  os << "Only particle types 20 and 30 are currently supported: "
1180  << ssd.particle_type;
1181  throw std::runtime_error(os.str());
1182  break;
1183  }
1184  }
1185 }
1186 
1187 
1188 // Documentation in header file.
1189 void tmatrix_ampld_test(const Verbosity& verbosity)
1190 {
1191  CREATE_OUT0;
1192 
1193  out0 << "======================================================\n";
1194  out0 << "Test for nonspherical particles in a fixed orientation\n";
1195  out0 << "Output should match 3rdparty/tmatrix/tmatrix_ampld.ref\n";
1196  out0 << "======================================================\n";
1197 
1198  // Same inputs as in example included in original ampld.lp.f
1199  Numeric rat = 1.;
1200  Numeric axi = 10.;
1201  Index np = -1;
1202  Numeric lam = acos(-1.)*2.;
1203  Numeric eps = 0.5;
1204  Numeric mrr = 1.5;
1205  Numeric mri = 0.02;
1206  Numeric ddelt = 0.001;
1207 
1208  Index quiet = 1;
1209 
1210  // Output variables
1211  Index nmax;
1212  Numeric csca;
1213  Numeric cext;
1214  char errmsg[1024] = "";
1215 
1216  tmatrix_(rat, axi, np, lam, eps, mrr, mri, ddelt, quiet,
1217  nmax, csca, cext, errmsg);
1218 
1219  out0 << "nmax: " << nmax << "\n";
1220  out0 << "csca: " << csca << "\n";
1221  out0 << "cext: " << cext << "\n";
1222 
1223  out0 << "Error message: " << (strlen(errmsg)?errmsg:"None") << "\n";
1224 
1225  // Same inputs as in example included in original ampld.lp.f
1226  Numeric alpha = 145.;
1227  Numeric beta = 52.;
1228  Numeric thet0 = 56.;
1229  Numeric thet = 65.;
1230  Numeric phi0 = 114.;
1231  Numeric phi = 128.;
1232 
1233  // Output variables
1234  Complex s11;
1235  Complex s12;
1236  Complex s21;
1237  Complex s22;
1238  ampl_(nmax, lam, thet0, thet, phi0, phi, alpha, beta,
1239  s11, s12, s21, s22);
1240 
1241  out0 << "AMPLITUDE MATRIX: \n";
1242  out0 << "s11: " << s11 << "\n";
1243  out0 << "s12: " << s12 << "\n";
1244  out0 << "s21: " << s21 << "\n";
1245  out0 << "s22: " << s22 << "\n";
1246 
1247  Matrix z;
1248  ampmat_to_phamat(z, s11, s12, s21, s22);
1249  z *= 1e12; // meter^2 to micron^2 for comparison with original results
1250 
1251  out0 << "PHASE MATRIX: \n" << z << "\n";
1252 }
1253 
1254 
1255 // Documentation in header file.
1256 void tmatrix_tmd_test(const Verbosity& verbosity)
1257 {
1258  CREATE_OUT0;
1259 
1260  out0 << "======================================================\n";
1261  out0 << "Test for randomly oriented nonspherical particles\n";
1262  out0 << "Output should match 3rdparty/tmatrix/tmatrix_tmd.ref\n";
1263  out0 << "======================================================\n";
1264 
1265  // Same inputs as in example included in original tmd.lp.f
1266  Numeric rat = 0.5;
1267  Index ndistr = 3;
1268  Numeric axmax = 1.;
1269  Index npnax = 2;
1270  Numeric b = 0.1;
1271  Numeric gam = 0.5;
1272  Index nkmax = 5;
1273  Numeric eps = 2;
1274  Index np = -1;
1275  Numeric lam = 0.5;
1276  Numeric mrr = 1.53;
1277  Numeric mri = 0.008;
1278  Numeric ddelt = 0.001;
1279  Index npna = 19;
1280  Index ndgs = 2;
1281  Numeric r1rat = 0.89031;
1282  Numeric r2rat = 1.56538;
1283 
1284  Index quiet = 1;
1285 
1286  // Output variables
1287  Numeric reff;
1288  Numeric veff;
1289  Numeric cext;
1290  Numeric csca;
1291  Numeric walb;
1292  Numeric asymm;
1293  Vector f11(npna, 0.);
1294  Vector f22(npna, 0.);
1295  Vector f33(npna, 0.);
1296  Vector f44(npna, 0.);
1297  Vector f12(npna, 0.);
1298  Vector f34(npna, 0.);
1299  char errmsg[1024] = "";
1300 
1301  tmd_(rat,
1302  ndistr,
1303  axmax,
1304  npnax,
1305  b,
1306  gam,
1307  nkmax,
1308  eps,
1309  np,
1310  lam,
1311  mrr,
1312  mri,
1313  ddelt,
1314  npna,
1315  ndgs,
1316  r1rat,
1317  r2rat,
1318  quiet,
1319  reff,
1320  veff,
1321  cext,
1322  csca,
1323  walb,
1324  asymm,
1325  f11.get_c_array(),
1326  f22.get_c_array(),
1327  f33.get_c_array(),
1328  f44.get_c_array(),
1329  f12.get_c_array(),
1330  f34.get_c_array(),
1331  errmsg);
1332 
1333  out0 << "reff: " << reff << "\n";
1334  out0 << "veff: " << veff << "\n";
1335  out0 << "cext: " << cext << "\n";
1336  out0 << "csca: " << csca << "\n";
1337  out0 << "walb: " << walb << "\n";
1338  out0 << "asymm: " << asymm << "\n";
1339  out0 << "f11: " << f11 << "\n";
1340  out0 << "f22: " << f22 << "\n";
1341  out0 << "f33: " << f33 << "\n";
1342  out0 << "f44: " << f44 << "\n";
1343  out0 << "f12: " << f12 << "\n";
1344  out0 << "f34: " << f34 << "\n";
1345 
1346  out0 << "Error message: " << (strlen(errmsg)?errmsg:"None") << "\n";
1347 }
1348 
1349 
1350 // Documentation in header file.
1351 void calc_ssp_random_test(const Verbosity& verbosity)
1352 {
1353  CREATE_OUT0;
1354  out0 << "======================================================\n";
1355  out0 << "Test calculation of single scattering data\n";
1356  out0 << "for randomly oriented, oblate particles\n";
1357  out0 << "======================================================\n";
1358 
1359 
1360  SingleScatteringData ssd;
1361 
1362  ssd.particle_type = PARTICLE_TYPE_MACROS_ISO;
1363  ssd.f_grid = MakeVector(230e9, 240e9);
1364  ssd.T_grid = MakeVector(220, 250);
1365  nlinspace(ssd.za_grid, 0, 180, 19);
1366  nlinspace(ssd.aa_grid, 0, 180, 19);
1367 
1368  // Refractive index real and imagenary parts
1369  // Dimensions: [nf, nT];
1370  Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
1371  Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
1372 
1373  mrr(0, 0) = 1.78031135; mrr(0, 1) = 1.78150475;
1374  mrr(1, 0) = 1.78037238; mrr(1, 1) = 1.78147686;
1375 
1376  mri(0, 0) = 0.00278706; mri(0, 1) = 0.00507565;
1377  mri(1, 0) = 0.00287245; mri(1, 1) = 0.00523012;
1378 
1379  calcSingleScatteringDataProperties(ssd, mrr, mri,
1380  200.,
1381  -1,
1382  1.5);
1383 
1384  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1385  << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1386 
1387  out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
1388  out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1389 
1390  out0 << "======================================================\n";
1391  out0 << "Test calculation of single scattering data\n";
1392  out0 << "for randomly oriented, prolate particles\n";
1393  out0 << "======================================================\n";
1394 
1395  calcSingleScatteringDataProperties(ssd, mrr, mri,
1396  200.,
1397  -1,
1398  0.7);
1399 
1400  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1401  << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1402 
1403  out0 << "ssd.ext_mat_data:\n" << ssd.ext_mat_data << "\n\n";
1404  out0 << "ssd.abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1405 }
1406 
1407 
1408 // Documentation in header file.
1409 void calc_ssp_fixed_test(const Verbosity& verbosity)
1410 {
1411  CREATE_OUT0;
1412  out0 << "======================================================\n";
1413  out0 << "Test calculation of single scattering data\n";
1414  out0 << "for oblate particles with fixed orientation\n";
1415  out0 << "======================================================\n";
1416 
1417 
1418  SingleScatteringData ssd;
1419 
1420  ssd.particle_type = PARTICLE_TYPE_HORIZ_AL;
1421  ssd.f_grid = MakeVector(230e9, 240e9);
1422  ssd.T_grid = MakeVector(220, 250);
1423  nlinspace(ssd.za_grid, 0, 180, 19);
1424  nlinspace(ssd.aa_grid, 0, 180, 19);
1425 
1426  // Refractive index real and imagenary parts
1427  // Dimensions: [nf, nT];
1428  Matrix mrr(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 1.78031135);
1429  Matrix mri(ssd.f_grid.nelem(), ssd.T_grid.nelem(), 0.00278706);
1430 
1431  mrr(0, 0) = 1.78031135; mrr(0, 1) = 1.78150475;
1432  mrr(1, 0) = 1.78037238; mrr(1, 1) = 1.78147686;
1433 
1434  mri(0, 0) = 0.00278706; mri(0, 1) = 0.00507565;
1435  mri(1, 0) = 0.00287245; mri(1, 1) = 0.00523012;
1436 
1437  calcSingleScatteringDataProperties(ssd, mrr, mri,
1438  200.,
1439  -1,
1440  1.5);
1441 
1442  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1443  << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1444 
1445  out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1446  << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
1447 
1448  out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1449 
1450  out0 << "======================================================\n";
1451  out0 << "Test calculation of single scattering data\n";
1452  out0 << "for prolate particles with fixed orientation\n";
1453  out0 << "======================================================\n";
1454  calcSingleScatteringDataProperties(ssd, mrr, mri,
1455  200.,
1456  -1,
1457  0.7);
1458 
1459  out0 << "ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker):\n"
1460  << ssd.pha_mat_data(0, 0, joker, 0, 0, joker, joker) << "\n\n";
1461 
1462  out0 << "ssd.ext_mat_data(0, 0, joker, joker, joker):\n"
1463  << ssd.ext_mat_data(0, 0, joker, joker, joker) << "\n\n";
1464 
1465  out0 << "abs_vec_data:\n" << ssd.abs_vec_data << "\n\n";
1466 }
1467 
tmatrix_tmd_test
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1256
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:631
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
precision
#define precision
Definition: logic.cc:45
VectorView
The VectorView class.
Definition: matpackI.h:372
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)
Calculate SingleScatteringData properties.
Definition: tmatrix.cc:927
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:2430
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:346
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:383
Vector
The Vector class.
Definition: matpackI.h:556
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:523
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:426
Range
The range class.
Definition: matpackI.h:148
SingleScatteringData
Structure which describes the single scattering properties of a particle or a particle distribution...
Definition: optproperties.h:84
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
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 quiet=1)
Calculate properties for randomly oriented particles.
Definition: tmatrix.cc:802
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
VectorView::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackI.cc:542
calc_ssp_fixed_test
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
Definition: tmatrix.cc:1409
Complex
std::complex< Numeric > Complex
Definition: complex.h:32
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:459
PI
const Numeric PI
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:402
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:261
calc_ssp_random_test
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition: tmatrix.cc:1351
joker
const Joker joker
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Matrix
The Matrix class.
Definition: matpackI.h:788
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:896
avgtmatrix_
void avgtmatrix_(const Index &nmax)
Perform orientation averaging.
Definition: tmatrix.cc:418
make_vector.h
The class MakeVector is a special kind of Vector that can be initialized explicitly from one or more ...
Vector::resize
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
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:576
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:596
real
float real
Definition: continua.cc:20315
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:211
beta
#define beta
Definition: continua.cc:21194
SingleScatteringData::particle_type
ParticleType particle_type
Definition: optproperties.h:85
tmatrix_ampld_test
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1189
Tensor7::resize
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5379
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:832
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:714
Tensor5
The Tensor5 class.
Definition: matpackV.h:451
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580