ARTS  2.0.49
optproperties.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2008 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& ptype,
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 (ptype){
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& ptype,
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 (ptype){
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& ptype,
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 (ptype){
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 
539 
541 
565 void interpolate_scat_angleDOIT(//Output:
566  VectorView pha_mat_int,
567  //Input:
568  ConstTensor5View pha_mat_data,
569  const Index& za_sca_idx,
570  const Index& aa_sca_idx,
571  const Index& za_inc_idx,
572  const Index& aa_inc_idx,
573  const ArrayOfArrayOfArrayOfArrayOfGridPos&
574  scat_theta_gps,
575  ConstTensor5View scat_theta_itws
576  )
577 {
578 
579  ConstVectorView itw = scat_theta_itws(za_sca_idx, aa_sca_idx, za_inc_idx, aa_inc_idx, joker);
580 
581  for (Index i = 0; i < 6; i++)
582  {
583  pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i),
584  scat_theta_gps[za_sca_idx][aa_sca_idx][za_inc_idx][aa_inc_idx]);
585  }
586 
587 }
588 
589 
591 
614 void interpolate_scat_angle(//Output:
615  VectorView pha_mat_int,
616  Numeric& theta_rad,
617  //Input:
618  ConstTensor5View pha_mat_data,
619  ConstVectorView za_datagrid,
620  const Numeric& za_sca,
621  const Numeric& aa_sca,
622  const Numeric& za_inc,
623  const Numeric& aa_inc
624  )
625 {
626  Numeric ANG_TOL=1e-7;
627 
628  //Calculate scattering angle from incident and scattered directions.
629  //The two special cases are implemented here to avoid NaNs that can
630  //sometimes occur in in the acos... formula in forward and backscatter
631  //cases. CPD 5/10/03.
632  //
633  // Consider not only aa_sca-aa_inc ~= 0, but also aa_sca-aa_inc ~= 360.
634  // GH 2011-05-31
635 
636  if( (abs(aa_sca-aa_inc)<ANG_TOL) || (abs(abs(aa_sca-aa_inc)-360) < ANG_TOL) )
637  {
638  theta_rad=DEG2RAD*abs(za_sca-za_inc);
639  }
640  else if (abs(abs(aa_sca-aa_inc)-180)<ANG_TOL)
641  {
642  theta_rad=DEG2RAD*(za_sca+za_inc);
643  if (theta_rad>PI){theta_rad=2*PI-theta_rad;}
644  }
645  else
646  {
647  const Numeric za_sca_rad = za_sca * DEG2RAD;
648  const Numeric za_inc_rad = za_inc * DEG2RAD;
649  const Numeric aa_sca_rad = aa_sca * DEG2RAD;
650  const Numeric aa_inc_rad = aa_inc * DEG2RAD;
651 
652  // cout << "Interpolation on scattering angle" << endl;
653  assert (pha_mat_data.ncols() == 6);
654  // Calculation of the scattering angle:
655  theta_rad = acos(cos(za_sca_rad) * cos(za_inc_rad) +
656  sin(za_sca_rad) * sin(za_inc_rad) *
657  cos(aa_sca_rad - aa_inc_rad));
658  }
659  const Numeric theta = RAD2DEG * theta_rad;
660 
661  // Interpolation of the data on the scattering angle:
662 
663  GridPos thet_gp;
664  gridpos(thet_gp, za_datagrid, theta);
665 
666  Vector itw(2);
667  interpweights(itw, thet_gp);
668 
669  for (Index i = 0; i < 6; i++)
670  {
671  pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i),
672  thet_gp);
673  }
674 }
675 
676 
677 
678 
680 
705 void pha_mat_labCalc(//Output:
706  MatrixView pha_mat_lab,
707  //Input:
708  ConstVectorView pha_mat_int,
709  const Numeric& za_sca,
710  const Numeric& aa_sca,
711  const Numeric& za_inc,
712  const Numeric& aa_inc,
713  const Numeric& theta_rad)
714 {
715  Numeric za_sca_rad = za_sca * DEG2RAD;
716  Numeric za_inc_rad = za_inc * DEG2RAD;
717  Numeric aa_sca_rad = aa_sca * DEG2RAD;
718  Numeric aa_inc_rad = aa_inc * DEG2RAD;
719 
720  const Numeric theta = RAD2DEG * theta_rad;
721  const Index stokes_dim = pha_mat_lab.ncols();
722 
723  // cout << "Transformation of phase matrix:" <<endl;
724 
725 
726 
727  // For stokes_dim = 1, we only need Z11=F11:
728  pha_mat_lab(0,0) = F11;
729 
730  if( stokes_dim > 1 ){
731  //
732  // Several cases have to be considered:
733  //
734  const Numeric ANGTOL = 1e-6; //CPD: this constant is used to adjust zenith angles
735  //close to 0 and PI. This is also used to avoid
736  //float == float statements.
737 
738  if(
739  // Forward scattering
740  ((theta > -.01) && (theta < .01) ) ||
741  // Backward scattering
742  ((theta > 179.99) && (theta < 180.01)) ||
743  // "Grosskreis" through poles: no rotation required
744  ((aa_sca == aa_inc) || (aa_sca == 360-aa_inc) || (aa_inc == 360-aa_sca) ||
745  (aa_sca == 180-aa_inc) || (aa_inc == 180-aa_sca) )
746  )
747  {
748  pha_mat_lab(0,1) = F12;
749  pha_mat_lab(1,0) = F12;
750  pha_mat_lab(1,1) = F22;
751 
752  if( stokes_dim > 2 ){
753  pha_mat_lab(0,2) = 0;
754  pha_mat_lab(1,2) = 0;
755  pha_mat_lab(2,0) = 0;
756  pha_mat_lab(2,1) = 0;
757  pha_mat_lab(2,2) = F33;
758 
759  if( stokes_dim > 3 ){
760  pha_mat_lab(0,3) = 0;
761  pha_mat_lab(1,3) = 0;
762  pha_mat_lab(2,3) = F34;
763  pha_mat_lab(3,0) = 0;
764  pha_mat_lab(3,1) = 0;
765  pha_mat_lab(3,2) = -F34;
766  pha_mat_lab(3,3) = F44;
767  }
768  }
769  }
770 
771  else
772  {
773  Numeric sigma1;
774  Numeric sigma2;
775 
776  Numeric s1, s2;
777 
778  // In these cases we have to take limiting values.
779 
780  if (za_inc_rad < ANGTOL)
781  {
782  sigma1 = PI + aa_sca_rad - aa_inc_rad;
783  sigma2 = 0;
784  }
785  else if (za_inc_rad > PI-ANGTOL)
786  {
787  sigma1 = aa_sca_rad - aa_inc_rad;
788  sigma2 = PI;
789  }
790  else if (za_sca_rad < ANGTOL)
791  {
792  sigma1 = 0;
793  sigma2 = PI + aa_sca_rad - aa_inc_rad;
794  }
795  else if (za_sca_rad > PI - ANGTOL)
796  {
797  sigma1 = PI;
798  sigma2 = aa_sca_rad - aa_inc_rad;
799  }
800  else
801  {
802  s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad))
803  /(sin(za_inc_rad)*sin(theta_rad));
804  s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos (theta_rad))/
805  (sin(za_sca_rad)*sin(theta_rad));
806 
807  sigma1 = acos(s1);
808  sigma2 = acos(s2);
809 
810  // Arccos is only defined in the range from -1 ... 1
811  // Numerical problems can appear for values close to 1 or -1
812  if ( isnan(sigma1) || isnan(sigma2) )
813  {
814  if ( abs(s1 - 1) < ANGTOL)
815  sigma1 = 0;
816  if ( abs(s1 + 1) < ANGTOL)
817  sigma1 = PI;
818  if ( abs(s2 - 1) < ANGTOL)
819  sigma2 = 0;
820  if ( abs(s2 + 1) < ANGTOL)
821  sigma2 = PI;
822  }
823  }
824 
825  const Numeric C1 = cos(2*sigma1);
826  const Numeric C2 = cos(2*sigma2);
827 
828  const Numeric S1 = sin(2*sigma1);
829  const Numeric S2 = sin(2*sigma2);
830 
831  pha_mat_lab(0,1) = C1 * F12;
832  pha_mat_lab(1,0) = C2 * F12;
833  pha_mat_lab(1,1) = C1 * C2 * F22 - S1 * S2 * F33;
834 
835  assert(!isnan(pha_mat_lab(0,1)));
836  assert(!isnan(pha_mat_lab(1,0)));
837  assert(!isnan(pha_mat_lab(1,1)));
838 
839  if( stokes_dim > 2 ){
840  /*CPD: For skokes_dim > 2 some of the transformation formula
841  for each element have a different sign depending on whether or
842  not 0<aa_scat-aa_inc<180. For details see pages 94 and 95 of
843  Mishchenkos chapter in :
844  Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic
845  scattering by nonspherical particles. In Exploring the Atmosphere
846  by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag,
847  Berlin, pp. 77-127.
848  This is available at http://www.giss.nasa.gov/~crmim/publications/ */
849  Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
850  (aa_sca-aa_inc>180)*360;
851  if(delta_aa>=0)
852  {
853  pha_mat_lab(0,2) = S1 * F12;
854  pha_mat_lab(1,2) = S1 * C2 * F22 + C1 * S2 * F33;
855  pha_mat_lab(2,0) = -S2 * F12;
856  pha_mat_lab(2,1) = -C1 * S2 * F22 - S1 * C2 * F33;
857  }
858  else
859  {
860  pha_mat_lab(0,2) = -S1 * F12;
861  pha_mat_lab(1,2) = -S1 * C2 * F22 - C1 * S2 * F33;
862  pha_mat_lab(2,0) = S2 * F12;
863  pha_mat_lab(2,1) = C1 * S2 * F22 + S1 * C2 * F33;
864  }
865  pha_mat_lab(2,2) = -S1 * S2 * F22 + C1 * C2 * F33;
866 
867  if( stokes_dim > 3 ){
868  if(delta_aa>=0)
869  {
870  pha_mat_lab(1,3) = S2 * F34;
871  pha_mat_lab(3,1) = S1 * F34;
872  }
873  else
874  {
875  pha_mat_lab(1,3) = -S2 * F34;
876  pha_mat_lab(3,1) = -S1 * F34;
877  }
878  pha_mat_lab(0,3) = 0;
879  pha_mat_lab(2,3) = C2 * F34;
880  pha_mat_lab(3,0) = 0;
881  pha_mat_lab(3,2) = -C1 * F34;
882  pha_mat_lab(3,3) = F44;
883  }
884  }
885  }
886  }
887 }
888 
889 
890 ostream& operator<< (ostream& os, const SingleScatteringData& /*ssd*/)
891 {
892  os << "SingleScatteringData: Output operator not implemented";
893  return os;
894 }
895 
896 
897 ostream& operator<< (ostream& os, const ArrayOfSingleScatteringData& /*assd*/)
898 {
899  os << "ArrayOfSingleScatteringData: Output operator not implemented";
900  return os;
901 }
902 
903 ostream& operator<< (ostream& os, const ScatteringMetaData& /*ssd*/)
904 {
905  os << "ScatteringMetaData: Output operator not implemented";
906  return os;
907 }
908 
909 
910 ostream& operator<< (ostream& os, const ArrayOfScatteringMetaData& /*assd*/)
911 {
912  os << "ArrayOfScatteringMetaData: Output operator not implemented";
913  return os;
914 }
RAD2DEG
const Numeric RAD2DEG
F11
#define F11
Definition: optproperties.cc:56
MatrixView
The MatrixView class.
Definition: matpackI.h:668
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
F22
#define F22
Definition: optproperties.cc:58
ConstTensor5View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackV.cc:56
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:614
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:1000
ext_matTransform
void ext_matTransform(MatrixView ext_mat_lab, ConstTensor3View ext_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const ParticleType &ptype, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of extinction matrix.
Definition: optproperties.cc:190
ParticleType
ParticleType
An attribute to classify the particle type in a SingleScatteringData.
Definition: optproperties.h:55
F34
#define F34
Definition: optproperties.cc:60
PARTICLE_TYPE_MACROS_ISO
@ PARTICLE_TYPE_MACROS_ISO
Definition: optproperties.h:57
abs_vecTransform
void abs_vecTransform(VectorView abs_vec_lab, ConstTensor3View abs_vec_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const ParticleType &ptype, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of absorption vector.
Definition: optproperties.cc:85
array.h
This file contains the definition of Array.
ConstTensor5View::npages
Index npages() const
Returns the number of pages.
Definition: matpackV.cc:44
SingleScatteringData
Structure which describes the single scattering properties of a.
Definition: optproperties.h:74
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:151
DEG2RAD
const Numeric DEG2RAD
Array
This can be used to make arrays out of anything.
Definition: array.h:103
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:205
pha_matTransform
void pha_matTransform(MatrixView pha_mat_lab, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const ParticleType &ptype, 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
ANGTOL
const double ANGTOL
Definition: ppath.cc:98
messages.h
Declarations having to do with the four output streams.
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:802
_U_
#define _U_
Definition: config.h:158
abs
#define abs(x)
Definition: continua.cc:13094
ScatteringMetaData
Definition: optproperties.h:99
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:175
VectorView
The VectorView class.
Definition: matpackI.h:373
optproperties.h
Scattering database structure and functions.
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:50
math_funcs.h
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:565
Range
The range class.
Definition: matpackI.h:148
GridPos
Structure to store a grid position.
Definition: interpolation.h:74
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:705
F33
#define F33
Definition: optproperties.cc:59
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:157
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:147
PARTICLE_TYPE_HORIZ_AL
@ PARTICLE_TYPE_HORIZ_AL
Definition: optproperties.h:58
operator<<
ostream & operator<<(ostream &os, const SingleScatteringData &)
Definition: optproperties.cc:890
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:39
PI
const Numeric PI
Vector
The Vector class.
Definition: matpackI.h:555
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:300
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:756
ConstTensor5View
A constant view of a Tensor5.
Definition: matpackV.h:160
matpackVII.h
arts.h
The global header file for ARTS.
xml_io.h
This file contains basic functions to handle XML data files.