ARTS  2.2.66
optproperties.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2012 Claudia Emde <claudia.emde@dlr.de>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  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, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
18 
19 
20 /*===========================================================================
21  === File description
22  ===========================================================================*/
23 
36 /*===========================================================================
37  === External declarations
38  ===========================================================================*/
39 
40 #include <cmath>
41 #include <stdexcept>
42 #include "arts.h"
43 #include "matpackVII.h"
44 #include "array.h"
45 #include "math_funcs.h"
46 #include "messages.h"
47 #include "logic.h"
48 #include "interpolation.h"
49 #include "optproperties.h"
50 #include "xml_io.h"
51 extern const Numeric DEG2RAD;
52 extern const Numeric RAD2DEG;
53 extern const Numeric PI;
54 
55 
56 #define F11 pha_mat_int[0]
57 #define F12 pha_mat_int[1]
58 #define F22 pha_mat_int[2]
59 #define F33 pha_mat_int[3]
60 #define F34 pha_mat_int[4]
61 #define F44 pha_mat_int[5]
62 
64 
85 void abs_vecTransform(//Output and Input
86  VectorView abs_vec_lab,
87  //Input
88  ConstTensor3View abs_vec_data,
89  ConstVectorView za_datagrid,
90  ConstVectorView aa_datagrid _U_,
91  const ParticleType& particle_type,
92  const Numeric& za_sca _U_,
93  const Numeric& aa_sca _U_,
94  const Verbosity& verbosity)
95 {
96  const Index stokes_dim = abs_vec_lab.nelem();
97 
98  if (stokes_dim > 4 || stokes_dim < 1){
99  throw runtime_error("The dimension of the stokes vector \n"
100  "must be 1,2,3 or 4");
101  }
102 
103  switch (particle_type){
104 
105  case PARTICLE_TYPE_GENERAL:
106  {
107  CREATE_OUT0;
108  out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n";
109  break;
110  }
111  case PARTICLE_TYPE_MACROS_ISO:
112  {
113  // The first element of the vector corresponds to the absorption
114  // coefficient which is stored in the database, the others are 0.
115 
116  abs_vec_lab = 0;
117 
118  abs_vec_lab[0] = abs_vec_data(0,0,0);
119  break;
120  }
121 
122  case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis 9/12/03
123  {
124  assert (abs_vec_data.ncols() == 2);
125 
126  // In the case of horizontally oriented particles the absorption
127  // coefficient vector only the first two elements are non-zero.
128  // These values are dependent on the zenith angle of propagation. The
129  // data storage format also makes use of the fact that in this case
130  //K_abs(za_sca)=K_abs(180-za_sca).
131 
132  // 1st interpolate data by za_sca
133  GridPos gp;
134  Vector itw(2);
135 
136  // JM 2010-11-05: reduce za_datagrid to abs_vec_data range,
137  // then gridpos correctly recognizes last grid point and
138  // returns correct index (follows SAB 2010-04-28)
139  ConstVectorView this_za_datagrid = za_datagrid[Range(0,abs_vec_data.npages())];
140 
141  if (za_sca>90)
142  {
143  gridpos(gp,this_za_datagrid,180-za_sca);
144  }
145  else
146  {
147  gridpos(gp,this_za_datagrid,za_sca);
148  }
149  interpweights(itw,gp);
150  abs_vec_lab = 0;
151  abs_vec_lab[0] = interp(itw,abs_vec_data(Range(joker),0,0),gp);
152 
153  if( stokes_dim == 1 ){
154  break;
155  }
156  abs_vec_lab[1] = interp(itw,abs_vec_data(Range(joker),0,1),gp);
157  break;
158  }
159  default:
160  {
161  CREATE_OUT0;
162  out0 << "Not all particle type cases are implemented\n";
163  }
164  }
165 }
166 
167 
169 
190 void ext_matTransform(//Output and Input
191  MatrixView ext_mat_lab,
192  //Input
193  ConstTensor3View ext_mat_data,
194  ConstVectorView za_datagrid,
195  ConstVectorView aa_datagrid _U_,
196  const ParticleType& particle_type,
197  const Numeric& za_sca,
198  const Numeric& aa_sca _U_,
199  const Verbosity& verbosity)
200 {
201  const Index stokes_dim = ext_mat_lab.ncols();
202 
203  if (stokes_dim > 4 || stokes_dim < 1){
204  throw runtime_error("The dimension of the stokes vector \n"
205  "must be 1,2,3 or 4");
206  }
207 
208  switch (particle_type){
209 
210  case PARTICLE_TYPE_GENERAL:
211  {
212  CREATE_OUT0;
213  out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n";
214  break;
215  }
216  case PARTICLE_TYPE_MACROS_ISO:
217  {
218  assert (ext_mat_data.ncols() == 1);
219 
220  // In the case of randomly oriented particles the extinction matrix is
221  // diagonal. The value of each element of the diagonal is the
222  // extinction cross section, which is stored in the database.
223 
224  ext_mat_lab = 0.;
225 
226  ext_mat_lab(0,0) = ext_mat_data(0,0,0);
227 
228 
229  if( stokes_dim == 1 ){
230  break;
231  }
232 
233  ext_mat_lab(1,1) = ext_mat_data(0,0,0);
234 
235  if( stokes_dim == 2 ){
236  break;
237  }
238 
239  ext_mat_lab(2,2) = ext_mat_data(0,0,0);
240 
241  if( stokes_dim == 3 ){
242  break;
243  }
244 
245  ext_mat_lab(3,3) = ext_mat_data(0,0,0);
246  break;
247  }
248 
249  case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis 9/12/03
250  {
251  assert (ext_mat_data.ncols() == 3);
252 
253  // In the case of horizontally oriented particles the extinction matrix
254  // has only 3 independent non-zero elements Kjj, K12=K21, and K34=-K43.
255  // These values are dependent on the zenith angle of propagation. The
256  // data storage format also makes use of the fact that in this case
257  //K(za_sca)=K(180-za_sca).
258 
259  // 1st interpolate data by za_sca
260  GridPos gp;
261  Vector itw(2);
262  Numeric Kjj;
263  Numeric K12;
264  Numeric K34;
265 
266  // JM 2010-11-05: reduce za_datagrid to ext_mat_data range,
267  // then gridpos correctly recognizes last grid point and
268  // returns correct index (follows SAB 2010-04-28)
269  ConstVectorView this_za_datagrid = za_datagrid[Range(0,ext_mat_data.npages())];
270 
271  if (za_sca>90)
272  {
273  gridpos(gp,this_za_datagrid,180-za_sca);
274  }
275  else
276  {
277  gridpos(gp,this_za_datagrid,za_sca);
278  }
279 
280  interpweights(itw,gp);
281 
282  ext_mat_lab=0.0;
283  Kjj=interp(itw,ext_mat_data(Range(joker),0,0),gp);
284  ext_mat_lab(0,0)=Kjj;
285 
286  if( stokes_dim == 1 ){
287  break;
288  }
289 
290  K12=interp(itw,ext_mat_data(Range(joker),0,1),gp);
291  ext_mat_lab(1,1)=Kjj;
292  ext_mat_lab(0,1)=K12;
293  ext_mat_lab(1,0)=K12;
294 
295  if( stokes_dim == 2 ){
296  break;
297  }
298 
299  ext_mat_lab(2,2)=Kjj;
300 
301  if( stokes_dim == 3 ){
302  break;
303  }
304 
305  K34=interp(itw,ext_mat_data(Range(joker),0,2),gp);
306  ext_mat_lab(2,3)=K34;
307  ext_mat_lab(3,2)=-K34;
308  ext_mat_lab(3,3)=Kjj;
309  break;
310 
311  }
312  default:
313  {
314  CREATE_OUT0;
315  out0 << "Not all particle type cases are implemented\n";
316  }
317  }
318 }
319 
320 
322 
345 void pha_matTransform(//Output
346  MatrixView pha_mat_lab,
347  //Input
348  ConstTensor5View pha_mat_data,
349  ConstVectorView za_datagrid,
350  ConstVectorView aa_datagrid,
351  const ParticleType& particle_type,
352  const Index& za_sca_idx,
353  const Index& aa_sca_idx,
354  const Index& za_inc_idx,
355  const Index& aa_inc_idx,
356  ConstVectorView scat_za_grid,
357  ConstVectorView scat_aa_grid,
358  const Verbosity& verbosity)
359 {
360  const Index stokes_dim = pha_mat_lab.ncols();
361 
362  Numeric za_sca = scat_za_grid[za_sca_idx];
363  Numeric aa_sca = scat_aa_grid[aa_sca_idx];
364  Numeric za_inc = scat_za_grid[za_inc_idx];
365  Numeric aa_inc = scat_aa_grid[aa_inc_idx];
366 
367  if (stokes_dim > 4 || stokes_dim < 1){
368  throw runtime_error("The dimension of the stokes vector \n"
369  "must be 1,2,3 or 4");
370  }
371 
372  switch (particle_type){
373 
374  case PARTICLE_TYPE_GENERAL:
375  {
376  CREATE_OUT0;
377  out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n";
378  break;
379  }
380  case PARTICLE_TYPE_MACROS_ISO:
381  {
382  // Calculate the scattering and interpolate the data on the scattering
383  // angle:
384 
385  Vector pha_mat_int(6);
386  Numeric theta_rad;
387 
388  // Interpolation of the data on the scattering angle:
389  interpolate_scat_angle(pha_mat_int, theta_rad, pha_mat_data,
390  za_datagrid, za_sca, aa_sca,
391  za_inc, aa_inc);
392 
393  // Caclulate the phase matrix in the laboratory frame:
394  pha_mat_labCalc(pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc,
395  aa_inc, theta_rad);
396 
397  break;
398  }
399 
400  case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis
401  //Data is already stored in the laboratory frame, but it is compressed
402  //a little. Details elsewhere
403  {
404  // SAB 2010-04-28: For the incoming angle, not the whole of za_datagrid
405  // is used in this case, only the range [0,90] degrees.
406  // (Inclusive interval at both ends.)
407  // How much is needed can be seen from pha_mat_data.npages().
408  ConstVectorView this_za_datagrid = za_datagrid[Range(0,pha_mat_data.npages())];
409 
410  assert (pha_mat_data.ncols()==16);
411  Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
412  (aa_sca-aa_inc>180)*360;//delta_aa corresponds to the "books"
413  //dimension of pha_mat_data
414  GridPos za_sca_gp;
415  GridPos delta_aa_gp;
416  GridPos za_inc_gp;
417  Vector itw(8);
418 
419  gridpos(delta_aa_gp,aa_datagrid,abs(delta_aa));
420  if (za_inc>90)
421  {
422  gridpos(za_inc_gp,this_za_datagrid,180-za_inc);
423  gridpos(za_sca_gp,za_datagrid,180-za_sca);
424  }
425  else
426  {
427  gridpos(za_inc_gp,this_za_datagrid,za_inc);
428  gridpos(za_sca_gp,za_datagrid,za_sca);
429  }
430 
431  interpweights(itw,za_sca_gp,delta_aa_gp,za_inc_gp);
432 
433  pha_mat_lab(0,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
434  Range(joker),0,0),
435  za_sca_gp,delta_aa_gp,za_inc_gp);
436  if( stokes_dim == 1 ){
437  break;
438  }
439  pha_mat_lab(0,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
440  Range(joker),0,1),
441  za_sca_gp,delta_aa_gp,za_inc_gp);
442  pha_mat_lab(1,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
443  Range(joker),0,4),
444  za_sca_gp,delta_aa_gp,za_inc_gp);
445  pha_mat_lab(1,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
446  Range(joker),0,5),
447  za_sca_gp,delta_aa_gp,za_inc_gp);
448  if( stokes_dim == 2 ){
449  break;
450  }
451  if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
452  {
453  pha_mat_lab(0,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
454  Range(joker),0,2),
455  za_sca_gp,delta_aa_gp,za_inc_gp);
456  pha_mat_lab(1,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
457  Range(joker),0,6),
458  za_sca_gp,delta_aa_gp,za_inc_gp);
459  pha_mat_lab(2,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
460  Range(joker),0,8),
461  za_sca_gp,delta_aa_gp,za_inc_gp);
462  pha_mat_lab(2,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
463  Range(joker),0,9),
464  za_sca_gp,delta_aa_gp,za_inc_gp);
465  }
466  else
467  {
468  pha_mat_lab(0,2)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
469  Range(joker),0,2),
470  za_sca_gp,delta_aa_gp,za_inc_gp);
471  pha_mat_lab(1,2)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
472  Range(joker),0,6),
473  za_sca_gp,delta_aa_gp,za_inc_gp);
474  pha_mat_lab(2,0)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
475  Range(joker),0,8),
476  za_sca_gp,delta_aa_gp,za_inc_gp);
477  pha_mat_lab(2,1)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
478  Range(joker),0,9),
479  za_sca_gp,delta_aa_gp,za_inc_gp);
480  }
481  pha_mat_lab(2,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
482  Range(joker),0,10),
483  za_sca_gp,delta_aa_gp,za_inc_gp);
484  if( stokes_dim == 3 ){
485  break;
486  }
487  if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
488  {
489  pha_mat_lab(0,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
490  Range(joker),0,3),
491  za_sca_gp,delta_aa_gp,za_inc_gp);
492  pha_mat_lab(1,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
493  Range(joker),0,7),
494  za_sca_gp,delta_aa_gp,za_inc_gp);
495  pha_mat_lab(3,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
496  Range(joker),0,12),
497  za_sca_gp,delta_aa_gp,za_inc_gp);
498  pha_mat_lab(3,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
499  Range(joker),0,13),
500  za_sca_gp,delta_aa_gp,za_inc_gp);
501  }
502  else
503  {
504  pha_mat_lab(0,3)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
505  Range(joker),0,3),
506  za_sca_gp,delta_aa_gp,za_inc_gp);
507  pha_mat_lab(1,3)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
508  Range(joker),0,7),
509  za_sca_gp,delta_aa_gp,za_inc_gp);
510  pha_mat_lab(3,0)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
511  Range(joker),0,12),
512  za_sca_gp,delta_aa_gp,za_inc_gp);
513  pha_mat_lab(3,1)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
514  Range(joker),0,13),
515  za_sca_gp,delta_aa_gp,za_inc_gp);
516  }
517  pha_mat_lab(2,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
518  Range(joker),0,11),
519  za_sca_gp,delta_aa_gp,za_inc_gp);
520  pha_mat_lab(3,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
521  Range(joker),0,14),
522  za_sca_gp,delta_aa_gp,za_inc_gp);
523  pha_mat_lab(3,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
524  Range(joker),0,15),
525  za_sca_gp,delta_aa_gp,za_inc_gp);
526  break;
527 
528  }
529  default:
530  {
531  CREATE_OUT0;
532  out0 << "Not all particle type cases are implemented\n";
533  }
534  }
535 }
536 
537 
538 
540 
568 void ext_matFromabs_vec(//Output:
569  MatrixView ext_mat,
570  //Input:
571  ConstVectorView abs_vec,
572  const Index& stokes_dim)
573 {
574  assert( stokes_dim>=1 && stokes_dim<=4 );
575  assert( ext_mat.nrows() == stokes_dim );
576  assert( ext_mat.ncols() == stokes_dim );
577  assert( abs_vec.nelem() == stokes_dim );
578 
579  // first: diagonal elements
580  for (Index is=0; is<stokes_dim; is++)
581  {
582  ext_mat(is,is) += abs_vec[0];
583  }
584  // second: off-diagonal elements, namely first row and column
585  for (Index is=1; is<stokes_dim; is++)
586  {
587  ext_mat(0,is) += abs_vec[is];
588  ext_mat(is,0) += abs_vec[is];
589  }
590 }
591 
592 
593 
595 
619 void interpolate_scat_angleDOIT(//Output:
620  VectorView pha_mat_int,
621  //Input:
622  ConstTensor5View pha_mat_data,
623  const Index& za_sca_idx,
624  const Index& aa_sca_idx,
625  const Index& za_inc_idx,
626  const Index& aa_inc_idx,
627  const ArrayOfArrayOfArrayOfArrayOfGridPos&
628  scat_theta_gps,
629  ConstTensor5View scat_theta_itws
630  )
631 {
632 
633  ConstVectorView itw = scat_theta_itws(za_sca_idx, aa_sca_idx, za_inc_idx, aa_inc_idx, joker);
634 
635  for (Index i = 0; i < 6; i++)
636  {
637  pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i),
638  scat_theta_gps[za_sca_idx][aa_sca_idx][za_inc_idx][aa_inc_idx]);
639  }
640 
641 }
642 
643 
645 
668 void interpolate_scat_angle(//Output:
669  VectorView pha_mat_int,
670  Numeric& theta_rad,
671  //Input:
672  ConstTensor5View pha_mat_data,
673  ConstVectorView za_datagrid,
674  const Numeric& za_sca,
675  const Numeric& aa_sca,
676  const Numeric& za_inc,
677  const Numeric& aa_inc
678  )
679 {
680  Numeric ANG_TOL=1e-7;
681 
682  //Calculate scattering angle from incident and scattered directions.
683  //The two special cases are implemented here to avoid NaNs that can
684  //sometimes occur in in the acos... formula in forward and backscatter
685  //cases. CPD 5/10/03.
686  //
687  // Consider not only aa_sca-aa_inc ~= 0, but also aa_sca-aa_inc ~= 360.
688  // GH 2011-05-31
689 
690  if( (abs(aa_sca-aa_inc)<ANG_TOL) || (abs(abs(aa_sca-aa_inc)-360) < ANG_TOL) )
691  {
692  theta_rad=DEG2RAD*abs(za_sca-za_inc);
693  }
694  else if (abs(abs(aa_sca-aa_inc)-180)<ANG_TOL)
695  {
696  theta_rad=DEG2RAD*(za_sca+za_inc);
697  if (theta_rad>PI){theta_rad=2*PI-theta_rad;}
698  }
699  else
700  {
701  const Numeric za_sca_rad = za_sca * DEG2RAD;
702  const Numeric za_inc_rad = za_inc * DEG2RAD;
703  const Numeric aa_sca_rad = aa_sca * DEG2RAD;
704  const Numeric aa_inc_rad = aa_inc * DEG2RAD;
705 
706  // cout << "Interpolation on scattering angle" << endl;
707  assert (pha_mat_data.ncols() == 6);
708  // Calculation of the scattering angle:
709  theta_rad = acos(cos(za_sca_rad) * cos(za_inc_rad) +
710  sin(za_sca_rad) * sin(za_inc_rad) *
711  cos(aa_sca_rad - aa_inc_rad));
712  }
713  const Numeric theta = RAD2DEG * theta_rad;
714 
715  // Interpolation of the data on the scattering angle:
716 
717  GridPos thet_gp;
718  gridpos(thet_gp, za_datagrid, theta);
719 
720  Vector itw(2);
721  interpweights(itw, thet_gp);
722 
723  for (Index i = 0; i < 6; i++)
724  {
725  pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i),
726  thet_gp);
727  }
728 }
729 
730 
731 
732 
734 
759 void pha_mat_labCalc(//Output:
760  MatrixView pha_mat_lab,
761  //Input:
762  ConstVectorView pha_mat_int,
763  const Numeric& za_sca,
764  const Numeric& aa_sca,
765  const Numeric& za_inc,
766  const Numeric& aa_inc,
767  const Numeric& theta_rad)
768 {
769  Numeric za_sca_rad = za_sca * DEG2RAD;
770  Numeric za_inc_rad = za_inc * DEG2RAD;
771  Numeric aa_sca_rad = aa_sca * DEG2RAD;
772  Numeric aa_inc_rad = aa_inc * DEG2RAD;
773 
774  const Numeric theta = RAD2DEG * theta_rad;
775  const Index stokes_dim = pha_mat_lab.ncols();
776 
777  // cout << "Transformation of phase matrix:" <<endl;
778 
779 
780 
781  // For stokes_dim = 1, we only need Z11=F11:
782  pha_mat_lab(0,0) = F11;
783 
784  if( stokes_dim > 1 ){
785  //
786  // Several cases have to be considered:
787  //
788  const Numeric ANGTOL = 1e-6; //CPD: this constant is used to adjust zenith angles
789  //close to 0 and PI. This is also used to avoid
790  //float == float statements.
791 
792  if(
793  // Forward scattering
794  ((theta > -.01) && (theta < .01) ) ||
795  // Backward scattering
796  ((theta > 179.99) && (theta < 180.01)) ||
797  // "Grosskreis" through poles: no rotation required
798  ((aa_sca == aa_inc) || (aa_sca == 360-aa_inc) || (aa_inc == 360-aa_sca) ||
799  (aa_sca == 180-aa_inc) || (aa_inc == 180-aa_sca) )
800  )
801  {
802  pha_mat_lab(0,1) = F12;
803  pha_mat_lab(1,0) = F12;
804  pha_mat_lab(1,1) = F22;
805 
806  if( stokes_dim > 2 ){
807  pha_mat_lab(0,2) = 0;
808  pha_mat_lab(1,2) = 0;
809  pha_mat_lab(2,0) = 0;
810  pha_mat_lab(2,1) = 0;
811  pha_mat_lab(2,2) = F33;
812 
813  if( stokes_dim > 3 ){
814  pha_mat_lab(0,3) = 0;
815  pha_mat_lab(1,3) = 0;
816  pha_mat_lab(2,3) = F34;
817  pha_mat_lab(3,0) = 0;
818  pha_mat_lab(3,1) = 0;
819  pha_mat_lab(3,2) = -F34;
820  pha_mat_lab(3,3) = F44;
821  }
822  }
823  }
824 
825  else
826  {
827  Numeric sigma1;
828  Numeric sigma2;
829 
830  Numeric s1, s2;
831 
832  // In these cases we have to take limiting values.
833 
834  if (za_inc_rad < ANGTOL)
835  {
836  sigma1 = PI + aa_sca_rad - aa_inc_rad;
837  sigma2 = 0;
838  }
839  else if (za_inc_rad > PI-ANGTOL)
840  {
841  sigma1 = aa_sca_rad - aa_inc_rad;
842  sigma2 = PI;
843  }
844  else if (za_sca_rad < ANGTOL)
845  {
846  sigma1 = 0;
847  sigma2 = PI + aa_sca_rad - aa_inc_rad;
848  }
849  else if (za_sca_rad > PI - ANGTOL)
850  {
851  sigma1 = PI;
852  sigma2 = aa_sca_rad - aa_inc_rad;
853  }
854  else
855  {
856  s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad))
857  /(sin(za_inc_rad)*sin(theta_rad));
858  s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos (theta_rad))/
859  (sin(za_sca_rad)*sin(theta_rad));
860 
861  sigma1 = acos(s1);
862  sigma2 = acos(s2);
863 
864  // Arccos is only defined in the range from -1 ... 1
865  // Numerical problems can appear for values close to 1 or -1
866  if ( isnan(sigma1) || isnan(sigma2) )
867  {
868  if ( abs(s1 - 1) < ANGTOL)
869  sigma1 = 0;
870  if ( abs(s1 + 1) < ANGTOL)
871  sigma1 = PI;
872  if ( abs(s2 - 1) < ANGTOL)
873  sigma2 = 0;
874  if ( abs(s2 + 1) < ANGTOL)
875  sigma2 = PI;
876  }
877  }
878 
879  const Numeric C1 = cos(2*sigma1);
880  const Numeric C2 = cos(2*sigma2);
881 
882  const Numeric S1 = sin(2*sigma1);
883  const Numeric S2 = sin(2*sigma2);
884 
885  pha_mat_lab(0,1) = C1 * F12;
886  pha_mat_lab(1,0) = C2 * F12;
887  pha_mat_lab(1,1) = C1 * C2 * F22 - S1 * S2 * F33;
888 
889  assert(!isnan(pha_mat_lab(0,1)));
890  assert(!isnan(pha_mat_lab(1,0)));
891  assert(!isnan(pha_mat_lab(1,1)));
892 
893  if( stokes_dim > 2 ){
894  /*CPD: For skokes_dim > 2 some of the transformation formula
895  for each element have a different sign depending on whether or
896  not 0<aa_scat-aa_inc<180. For details see pages 94 and 95 of
897  Mishchenkos chapter in :
898  Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic
899  scattering by nonspherical particles. In Exploring the Atmosphere
900  by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag,
901  Berlin, pp. 77-127.
902  This is available at http://www.giss.nasa.gov/~crmim/publications/ */
903  Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
904  (aa_sca-aa_inc>180)*360;
905  if(delta_aa>=0)
906  {
907  pha_mat_lab(0,2) = S1 * F12;
908  pha_mat_lab(1,2) = S1 * C2 * F22 + C1 * S2 * F33;
909  pha_mat_lab(2,0) = -S2 * F12;
910  pha_mat_lab(2,1) = -C1 * S2 * F22 - S1 * C2 * F33;
911  }
912  else
913  {
914  pha_mat_lab(0,2) = -S1 * F12;
915  pha_mat_lab(1,2) = -S1 * C2 * F22 - C1 * S2 * F33;
916  pha_mat_lab(2,0) = S2 * F12;
917  pha_mat_lab(2,1) = C1 * S2 * F22 + S1 * C2 * F33;
918  }
919  pha_mat_lab(2,2) = -S1 * S2 * F22 + C1 * C2 * F33;
920 
921  if( stokes_dim > 3 ){
922  if(delta_aa>=0)
923  {
924  pha_mat_lab(1,3) = S2 * F34;
925  pha_mat_lab(3,1) = S1 * F34;
926  }
927  else
928  {
929  pha_mat_lab(1,3) = -S2 * F34;
930  pha_mat_lab(3,1) = -S1 * F34;
931  }
932  pha_mat_lab(0,3) = 0;
933  pha_mat_lab(2,3) = C2 * F34;
934  pha_mat_lab(3,0) = 0;
935  pha_mat_lab(3,2) = -C1 * F34;
936  pha_mat_lab(3,3) = F44;
937  }
938  }
939  }
940  }
941 }
942 
943 
944 ostream& operator<< (ostream& os, const SingleScatteringData& /*ssd*/)
945 {
946  os << "SingleScatteringData: Output operator not implemented";
947  return os;
948 }
949 
950 
951 ostream& operator<< (ostream& os, const ArrayOfSingleScatteringData& /*assd*/)
952 {
953  os << "ArrayOfSingleScatteringData: Output operator not implemented";
954  return os;
955 }
956 
957 ostream& operator<< (ostream& os, const ScatteringMetaData& /*ssd*/)
958 {
959  os << "ScatteringMetaData: Output operator not implemented";
960  return os;
961 }
962 
963 
964 ostream& operator<< (ostream& os, const ArrayOfScatteringMetaData& /*assd*/)
965 {
966  os << "ArrayOfScatteringMetaData: Output operator not implemented";
967  return os;
968 }
969 
970 
972 
987 void opt_prop_sum_propmat_clearsky(//Output:
988  Tensor3& ext_mat,
989  Matrix& abs_vec,
990  //Input:
991  const Tensor4 propmat_clearsky)
992 {
993 
994  Index stokes_dim = propmat_clearsky.ncols();
995 
996  Index freq_dim = propmat_clearsky.npages();
997 
998  // old abs_vecInit
999  abs_vec.resize( freq_dim, stokes_dim );
1000  abs_vec = 0; // Initialize to zero!
1001 
1002  // old ext_matInit
1003  ext_mat.resize( freq_dim, stokes_dim, stokes_dim );
1004  ext_mat = 0; // Initialize to zero!
1005 
1006  // old ext_matAddGas and abs_vecAddGas for 0 vector and matrix
1007  for ( Index iv=0; iv<freq_dim; ++iv )
1008  for ( Index is1=0; is1<stokes_dim; ++is1 )
1009  {
1010  abs_vec(iv,is1) += propmat_clearsky(joker,iv,is1,0).sum();
1011  for ( Index is2=0; is2<stokes_dim; ++is2 )
1012  ext_mat(iv,is1,is2) += propmat_clearsky(joker,iv,is1,is2).sum();
1013  }
1014 }
1015 
1016 
1018 
1026 ParticleType ParticleTypeFromString(const String& particle_type_string)
1027 {
1028  ParticleType particle_type;
1029  if (particle_type_string == "general")
1030  particle_type = PARTICLE_TYPE_GENERAL;
1031  else if (particle_type_string == "macroscopically_isotropic")
1032  particle_type = PARTICLE_TYPE_MACROS_ISO;
1033  else if (particle_type_string == "horizontally_aligned")
1034  particle_type = PARTICLE_TYPE_HORIZ_AL;
1035  else if (particle_type_string == "spherical")
1036  particle_type = PARTICLE_TYPE_SPHERICAL;
1037  else
1038  {
1039  ostringstream os;
1040  os << "Unknown particle type: " << particle_type_string << endl
1041  << "Valid types are: general, macroscopically_isotropic, "
1042  << "horizontally_aligned and spherical.";
1043  throw std::runtime_error(os.str());
1044  }
1045 
1046  return particle_type;
1047 }
1048 
1049 
1051 
1059 String ParticleTypeToString(const ParticleType& particle_type)
1060 {
1061  String particle_type_string;
1062 
1063  switch (particle_type) {
1064  case PARTICLE_TYPE_GENERAL:
1065  particle_type_string = "general";
1066  break;
1067  case PARTICLE_TYPE_MACROS_ISO:
1068  particle_type_string = "macroscopically_isotropic";
1069  break;
1070  case PARTICLE_TYPE_HORIZ_AL:
1071  particle_type_string = "horizontally_aligned";
1072  break;
1073  case PARTICLE_TYPE_SPHERICAL:
1074  particle_type_string = "spherical";
1075  break;
1076  default:
1077  ostringstream os;
1078  os << "Internal error: Cannot map ParticleType enum value "
1079  << particle_type << " to String.";
1080  throw std::runtime_error(os.str());
1081  break;
1082  }
1083 
1084  return particle_type_string;
1085 }
1086 
1087 
1089 
1097 ParticleSSDMethod ParticleSSDMethodFromString(const String& particle_ssdmethod_string)
1098 {
1099  ParticleSSDMethod particle_ssdmethod;
1100  if (particle_ssdmethod_string == "tmatrix")
1101  particle_ssdmethod = PARTICLE_SSDMETHOD_TMATRIX;
1102  else
1103  {
1104  ostringstream os;
1105  os << "Unknown particle SSD method: " << particle_ssdmethod_string << endl
1106  << "Valid methods: tmatrix";
1107  throw std::runtime_error(os.str());
1108  }
1109 
1110  return particle_ssdmethod;
1111 }
1112 
1113 
1115 
1123 String ParticleTypeToString(const ParticleSSDMethod& particle_ssdmethod)
1124 {
1125  String particle_ssdmethod_string;
1126 
1127  switch (particle_ssdmethod) {
1128  case PARTICLE_SSDMETHOD_TMATRIX:
1129  particle_ssdmethod_string = "tmatrix";
1130  break;
1131  default:
1132  ostringstream os;
1133  os << "Internal error: Cannot map ParticleSSDMethod enum value "
1134  << particle_ssdmethod << " to String.";
1135  throw std::runtime_error(os.str());
1136  break;
1137  }
1138 
1139  return particle_ssdmethod_string;
1140 }
1141 
Matrix
The Matrix class.
Definition: matpackI.h:788
RAD2DEG
const Numeric RAD2DEG
ANGTOL
const Numeric ANGTOL
Definition: ppath.h:103
F11
#define F11
Definition: optproperties.cc:56
MatrixView
The MatrixView class.
Definition: matpackI.h:679
ParticleTypeToString
String ParticleTypeToString(const ParticleType &particle_type)
Convert particle type enum value to String.
Definition: optproperties.cc:1059
PARTICLE_TYPE_GENERAL
@ PARTICLE_TYPE_GENERAL
Definition: optproperties.h:56
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:167
Tensor3
The Tensor3 class.
Definition: matpackIII.h:348
F22
#define F22
Definition: optproperties.cc:58
ConstTensor5View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackV.cc:59
PARTICLE_SSDMETHOD_TMATRIX
@ PARTICLE_SSDMETHOD_TMATRIX
Definition: optproperties.h:70
joker
const Joker joker
interpolation.h
Header file for interpolation.cc.
interpolate_scat_angle
void interpolate_scat_angle(VectorView pha_mat_int, Numeric &theta_rad, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc)
Interpolate data on the scattering angle.
Definition: optproperties.cc:668
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:1046
ParticleType
ParticleType
An attribute to classify the particle type in a SingleScatteringData structure.
Definition: optproperties.h:55
F34
#define F34
Definition: optproperties.cc:60
PARTICLE_TYPE_MACROS_ISO
@ PARTICLE_TYPE_MACROS_ISO
Definition: optproperties.h:57
PARTICLE_TYPE_SPHERICAL
@ PARTICLE_TYPE_SPHERICAL
Definition: optproperties.h:59
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
Tensor4
The Tensor4 class.
Definition: matpackIV.h:383
array.h
This file contains the definition of Array.
Tensor3::resize
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
ConstTensor5View::npages
Index npages() const
Returns the number of pages.
Definition: matpackV.cc:47
SingleScatteringData
Structure which describes the single scattering properties of a particle or a particle distribution.
Definition: optproperties.h:84
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
DEG2RAD
const Numeric DEG2RAD
ext_matTransform
void ext_matTransform(MatrixView ext_mat_lab, ConstTensor3View ext_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const ParticleType &particle_type, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of extinction matrix.
Definition: optproperties.cc:190
Array
This can be used to make arrays out of anything.
Definition: array.h:107
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:211
messages.h
Declarations having to do with the four output streams.
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
_U_
#define _U_
Definition: config.h:167
my_basic_string
The implementation for String, the ARTS string class.
Definition: mystring.h:64
abs
#define abs(x)
Definition: continua.cc:20458
ScatteringMetaData
Definition: optproperties.h:109
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
VectorView
The VectorView class.
Definition: matpackI.h:372
ConstTensor4View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
optproperties.h
Scattering database structure and functions.
pha_matTransform
void pha_matTransform(MatrixView pha_mat_lab, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const ParticleType &particle_type, const Index &za_sca_idx, const Index &aa_sca_idx, const Index &za_inc_idx, const Index &aa_inc_idx, ConstVectorView scat_za_grid, ConstVectorView scat_aa_grid, const Verbosity &verbosity)
Transformation of phase matrix.
Definition: optproperties.cc:345
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Verbosity
Definition: messages.h:50
ConstTensor4View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
math_funcs.h
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
abs_vecTransform
void abs_vecTransform(VectorView abs_vec_lab, ConstTensor3View abs_vec_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const ParticleType &particle_type, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of absorption vector.
Definition: optproperties.cc:85
interpolate_scat_angleDOIT
void interpolate_scat_angleDOIT(VectorView pha_mat_int, ConstTensor5View pha_mat_data, const Index &za_sca_idx, const Index &aa_sca_idx, const Index &za_inc_idx, const Index &aa_inc_idx, const ArrayOfArrayOfArrayOfArrayOfGridPos &scat_theta_gps, ConstTensor5View scat_theta_itws)
Interpolate data on the scattering angle.
Definition: optproperties.cc:619
ParticleSSDMethod
ParticleSSDMethod
An attribute to classify the method to be used for SingleScatteringData.
Definition: optproperties.h:68
ext_matFromabs_vec
void ext_matFromabs_vec(MatrixView ext_mat, ConstVectorView abs_vec, const Index &stokes_dim)
Derive extinction matrix from absorption vector.
Definition: optproperties.cc:568
Range
The range class.
Definition: matpackI.h:148
GridPos
Structure to store a grid position.
Definition: interpolation.h:74
ParticleSSDMethodFromString
ParticleSSDMethod ParticleSSDMethodFromString(const String &particle_ssdmethod_string)
Convert particle ssd method name to enum value.
Definition: optproperties.cc:1097
logic.h
Header file for logic.cc.
pha_mat_labCalc
void pha_mat_labCalc(MatrixView pha_mat_lab, ConstVectorView pha_mat_int, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc, const Numeric &theta_rad)
Calculate phase matrix in laboratory coordinate system.
Definition: optproperties.cc:759
F33
#define F33
Definition: optproperties.cc:59
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:139
PARTICLE_TYPE_HORIZ_AL
@ PARTICLE_TYPE_HORIZ_AL
Definition: optproperties.h:58
operator<<
ostream & operator<<(ostream &os, const SingleScatteringData &)
Definition: optproperties.cc:944
ParticleTypeFromString
ParticleType ParticleTypeFromString(const String &particle_type_string)
Convert particle name to enum value.
Definition: optproperties.cc:1026
F12
#define F12
Definition: optproperties.cc:57
F44
#define F44
Definition: optproperties.cc:61
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
PI
const Numeric PI
opt_prop_sum_propmat_clearsky
void opt_prop_sum_propmat_clearsky(Tensor3 &ext_mat, Matrix &abs_vec, const Tensor4 propmat_clearsky)
Get optical properties from propmat_clearsky.
Definition: optproperties.cc:987
Vector
The Vector class.
Definition: matpackI.h:556
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:292
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:802
ConstTensor5View
A constant view of a Tensor5.
Definition: matpackV.h:152
matpackVII.h
arts.h
The global header file for ARTS.
xml_io.h
This file contains basic functions to handle XML data files.