optproperties.cc

Go to the documentation of this file.

/* Copyright (C) 2003-2012 Claudia Emde <claudia.emde@dlr.de>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */



/*===========================================================================
  === File description
  ===========================================================================*/

/*===========================================================================
  === External declarations
  ===========================================================================*/

#include <cmath>
#include <stdexcept>
#include "arts.h"
#include "matpackVII.h"
#include "array.h"
#include "math_funcs.h"
#include "messages.h"
#include "logic.h"
#include "interpolation.h"
#include "optproperties.h"
#include "xml_io.h"
extern const Numeric DEG2RAD;
extern const Numeric RAD2DEG;
extern const Numeric PI;


#define F11 pha_mat_int[0]
#define F12 pha_mat_int[1]
#define F22 pha_mat_int[2]
#define F33 pha_mat_int[3]
#define F34 pha_mat_int[4]
#define F44 pha_mat_int[5]


void abs_vecTransform(//Output and Input
                      VectorView abs_vec_lab,
                      //Input
                      ConstTensor3View abs_vec_data,
                      ConstVectorView za_datagrid,
                      ConstVectorView aa_datagrid _U_,
                      const ParticleType& particle_type,
                      const Numeric& za_sca _U_,
                      const Numeric& aa_sca _U_,
                      const Verbosity& verbosity)
{
  const Index stokes_dim = abs_vec_lab.nelem();
  
  if (stokes_dim > 4 || stokes_dim < 1){
    throw runtime_error("The dimension of the stokes vector \n"
                        "must be 1,2,3 or 4");
  }
  
  switch (particle_type){
      
    case PARTICLE_TYPE_GENERAL:
    {
      CREATE_OUT0;
      out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n"; 
      break;
    }
    case PARTICLE_TYPE_MACROS_ISO:
    {
      // The first element of the vector corresponds to the absorption 
      // coefficient which is stored in the database, the others are 0.
      
      abs_vec_lab = 0;
      
      abs_vec_lab[0] = abs_vec_data(0,0,0);
      break;
    }
      
    case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis 9/12/03
    {
      assert (abs_vec_data.ncols() == 2);
      
      // In the case of horizontally oriented particles the absorption
      // coefficient vector only the first two elements are non-zero.
      // These values are dependent on the zenith angle of propagation. The 
      // data storage format also makes use of the fact that in this case
      //K_abs(za_sca)=K_abs(180-za_sca). 
      
      // 1st interpolate data by za_sca
      GridPos gp;
      Vector itw(2);
      
      // JM 2010-11-05: reduce za_datagrid to abs_vec_data range,
      // then gridpos correctly recognizes last grid point and
      // returns correct index (follows SAB 2010-04-28)
      ConstVectorView this_za_datagrid = za_datagrid[Range(0,abs_vec_data.npages())];
      
      if (za_sca>90)
      {
        gridpos(gp,this_za_datagrid,180-za_sca);
      }
      else
      {
        gridpos(gp,this_za_datagrid,za_sca);
      }
      interpweights(itw,gp);
      abs_vec_lab = 0;
      abs_vec_lab[0] = interp(itw,abs_vec_data(Range(joker),0,0),gp);
      
      if( stokes_dim == 1 ){
        break;
      }
      abs_vec_lab[1] = interp(itw,abs_vec_data(Range(joker),0,1),gp);
      break;
    }
    default:
    {
      CREATE_OUT0;
      out0 << "Not all particle type cases are implemented\n";
    }
  }  
}



void ext_matTransform(//Output and Input
                      MatrixView ext_mat_lab,
                      //Input
                      ConstTensor3View ext_mat_data,
                      ConstVectorView za_datagrid,
                      ConstVectorView aa_datagrid _U_,
                      const ParticleType& particle_type,
                      const Numeric& za_sca,
                      const Numeric& aa_sca _U_,
                      const Verbosity& verbosity)
{
  const Index stokes_dim = ext_mat_lab.ncols();

  if (stokes_dim > 4 || stokes_dim < 1){
    throw runtime_error("The dimension of the stokes vector \n"
                        "must be 1,2,3 or 4");
  }
  
  switch (particle_type){
      
    case PARTICLE_TYPE_GENERAL:
    {
      CREATE_OUT0;
      out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n"; 
      break;
    }
    case PARTICLE_TYPE_MACROS_ISO:
    {
      assert (ext_mat_data.ncols() == 1);
      
      // In the case of randomly oriented particles the extinction matrix is 
      // diagonal. The value of each element of the diagonal is the
      // extinction cross section, which is stored in the database.
      
      ext_mat_lab = 0.;
      
      ext_mat_lab(0,0) = ext_mat_data(0,0,0);
      
      
      if( stokes_dim == 1 ){
        break;
      }
      
      ext_mat_lab(1,1) = ext_mat_data(0,0,0);
      
      if( stokes_dim == 2 ){
        break;
      }
      
      ext_mat_lab(2,2) = ext_mat_data(0,0,0);
      
      if( stokes_dim == 3 ){
        break;
      }
      
      ext_mat_lab(3,3) = ext_mat_data(0,0,0);
      break;
    }
      
    case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis 9/12/03
    {
      assert (ext_mat_data.ncols() == 3);
      
      // In the case of horizontally oriented particles the extinction matrix
      // has only 3 independent non-zero elements Kjj, K12=K21, and K34=-K43.
      // These values are dependent on the zenith angle of propagation. The 
      // data storage format also makes use of the fact that in this case
      //K(za_sca)=K(180-za_sca). 
      
      // 1st interpolate data by za_sca
      GridPos gp;
      Vector itw(2);
      Numeric Kjj;
      Numeric K12;
      Numeric K34;
      
      // JM 2010-11-05: reduce za_datagrid to ext_mat_data range,
      // then gridpos correctly recognizes last grid point and
      // returns correct index (follows SAB 2010-04-28)
      ConstVectorView this_za_datagrid = za_datagrid[Range(0,ext_mat_data.npages())];
      
      if (za_sca>90)
      {
        gridpos(gp,this_za_datagrid,180-za_sca);
      }
      else
      {
        gridpos(gp,this_za_datagrid,za_sca);
      }
      
      interpweights(itw,gp);
      
      ext_mat_lab=0.0;
      Kjj=interp(itw,ext_mat_data(Range(joker),0,0),gp);
      ext_mat_lab(0,0)=Kjj;
      
      if( stokes_dim == 1 ){
        break;
      }
      
      K12=interp(itw,ext_mat_data(Range(joker),0,1),gp);
      ext_mat_lab(1,1)=Kjj;
      ext_mat_lab(0,1)=K12;
      ext_mat_lab(1,0)=K12;
      
      if( stokes_dim == 2 ){
        break;
      }
      
      ext_mat_lab(2,2)=Kjj;
      
      if( stokes_dim == 3 ){
        break;
      }
      
      K34=interp(itw,ext_mat_data(Range(joker),0,2),gp);
      ext_mat_lab(2,3)=K34;
      ext_mat_lab(3,2)=-K34;
      ext_mat_lab(3,3)=Kjj;
      break;
      
    }
    default:
    {
      CREATE_OUT0;
      out0 << "Not all particle type cases are implemented\n";
    }
  }
}  
 


void pha_matTransform(//Output
                      MatrixView pha_mat_lab,
                      //Input
                      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)
{
  const Index stokes_dim = pha_mat_lab.ncols();
  
  Numeric za_sca = scat_za_grid[za_sca_idx]; 
  Numeric aa_sca = scat_aa_grid[aa_sca_idx];
  Numeric za_inc = scat_za_grid[za_inc_idx]; 
  Numeric aa_inc = scat_aa_grid[aa_inc_idx];
  
  if (stokes_dim > 4 || stokes_dim < 1){
    throw runtime_error("The dimension of the stokes vector \n"
                        "must be 1,2,3 or 4");
  }
  
  switch (particle_type){
      
    case PARTICLE_TYPE_GENERAL:
    {
      CREATE_OUT0;
      out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n"; 
      break;
    }
    case PARTICLE_TYPE_MACROS_ISO:
    {
      // Calculate the scattering and interpolate the data on the scattering
      // angle:
      
      Vector pha_mat_int(6);
      Numeric theta_rad;
      
      // Interpolation of the data on the scattering angle:
      interpolate_scat_angle(pha_mat_int, theta_rad, pha_mat_data,
                             za_datagrid, za_sca, aa_sca,
                             za_inc, aa_inc);
      
      // Caclulate the phase matrix in the laboratory frame:
      pha_mat_labCalc(pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, 
                      aa_inc, theta_rad);
      
      break;
    }
      
    case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis
                                //Data is already stored in the laboratory frame, but it is compressed
                                //a little.  Details elsewhere
    {
      // SAB 2010-04-28: For the incoming angle, not the whole of za_datagrid 
      // is used in this case, only the range [0,90] degrees. 
      // (Inclusive interval at both ends.)
      // How much is needed can be seen from pha_mat_data.npages().
      ConstVectorView this_za_datagrid = za_datagrid[Range(0,pha_mat_data.npages())];
      
      assert (pha_mat_data.ncols()==16);
      Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
      (aa_sca-aa_inc>180)*360;//delta_aa corresponds to the "books" 
                              //dimension of pha_mat_data
      GridPos za_sca_gp;
      GridPos delta_aa_gp;
      GridPos za_inc_gp;
      Vector itw(8);
      
      gridpos(delta_aa_gp,aa_datagrid,abs(delta_aa));
      if (za_inc>90)
      {
        gridpos(za_inc_gp,this_za_datagrid,180-za_inc);
        gridpos(za_sca_gp,za_datagrid,180-za_sca);
      }
      else
      {
        gridpos(za_inc_gp,this_za_datagrid,za_inc);
        gridpos(za_sca_gp,za_datagrid,za_sca);
      }
      
      interpweights(itw,za_sca_gp,delta_aa_gp,za_inc_gp);
      
      pha_mat_lab(0,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,0),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      if( stokes_dim == 1 ){
        break;
      }
      pha_mat_lab(0,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,1),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      pha_mat_lab(1,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,4),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      pha_mat_lab(1,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,5),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      if( stokes_dim == 2 ){
        break;
      }
      if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
      {
        pha_mat_lab(0,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,2),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(1,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,6),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(2,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,8),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(2,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,9),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
      }
      else
      {
        pha_mat_lab(0,2)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,2),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(1,2)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,6),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(2,0)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,8),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(2,1)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,9),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
      }                             
      pha_mat_lab(2,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,10),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      if( stokes_dim == 3 ){
        break;
      }
      if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
      {
        pha_mat_lab(0,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,3),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(1,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,7),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(3,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,12),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(3,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                 Range(joker),0,13),
                                za_sca_gp,delta_aa_gp,za_inc_gp);
      }
      else
      {
        pha_mat_lab(0,3)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,3),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(1,3)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,7),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(3,0)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,12),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
        pha_mat_lab(3,1)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
                                                  Range(joker),0,13),
                                 za_sca_gp,delta_aa_gp,za_inc_gp);
      }
      pha_mat_lab(2,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,11),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      pha_mat_lab(3,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,14),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      pha_mat_lab(3,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
                                               Range(joker),0,15),
                              za_sca_gp,delta_aa_gp,za_inc_gp);
      break;
      
    }  
    default:
    {
      CREATE_OUT0;
      out0 << "Not all particle type cases are implemented\n";
    }
  }
}




void ext_matFromabs_vec(//Output:
                        MatrixView ext_mat,
                        //Input:
                        ConstVectorView abs_vec,
                        const Index& stokes_dim)
{
  assert( stokes_dim>=1  &&  stokes_dim<=4 );
  assert( ext_mat.nrows() == stokes_dim );
  assert( ext_mat.ncols() == stokes_dim );
  assert( abs_vec.nelem() == stokes_dim );

  // first: diagonal elements
  for (Index is=0; is<stokes_dim; is++)
    {
      ext_mat(is,is) += abs_vec[0];
    }
  // second: off-diagonal elements, namely first row and column
  for (Index is=1; is<stokes_dim; is++)
    {
      ext_mat(0,is) += abs_vec[is];
      ext_mat(is,0) += abs_vec[is];
    }
}




void interpolate_scat_angleDOIT(//Output:
                            VectorView pha_mat_int,
                            //Input:
                            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
                            )
{
  
  ConstVectorView itw = scat_theta_itws(za_sca_idx, aa_sca_idx, za_inc_idx, aa_inc_idx, joker);
  
  for (Index i = 0; i < 6; i++)
    {
      pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), 
                             scat_theta_gps[za_sca_idx][aa_sca_idx][za_inc_idx][aa_inc_idx]);
    }
  
} 



void interpolate_scat_angle(//Output:
                            VectorView pha_mat_int,
                            Numeric& theta_rad,
                            //Input:
                            ConstTensor5View pha_mat_data,
                            ConstVectorView za_datagrid,
                            const Numeric& za_sca,
                            const Numeric& aa_sca,
                            const Numeric& za_inc,
                            const Numeric& aa_inc
                            )
{
  Numeric ANG_TOL=1e-7;

  //Calculate scattering angle from incident and scattered directions.
  //The two special cases are implemented here to avoid NaNs that can 
  //sometimes occur in in the acos... formula in forward and backscatter
  //cases. CPD 5/10/03.
  //
  // Consider not only aa_sca-aa_inc ~= 0, but also aa_sca-aa_inc ~= 360.
  // GH 2011-05-31
  
  if( (abs(aa_sca-aa_inc)<ANG_TOL) || (abs(abs(aa_sca-aa_inc)-360) < ANG_TOL) )
    {
      theta_rad=DEG2RAD*abs(za_sca-za_inc);
    }
  else if (abs(abs(aa_sca-aa_inc)-180)<ANG_TOL)
    {
      theta_rad=DEG2RAD*(za_sca+za_inc);
      if (theta_rad>PI){theta_rad=2*PI-theta_rad;}
    }
  else
    {
      const Numeric za_sca_rad = za_sca * DEG2RAD;
      const Numeric za_inc_rad = za_inc * DEG2RAD;
      const Numeric aa_sca_rad = aa_sca * DEG2RAD;
      const Numeric aa_inc_rad = aa_inc * DEG2RAD;
      
      // cout << "Interpolation on scattering angle" << endl;
      assert (pha_mat_data.ncols() == 6);
      // Calculation of the scattering angle:
      theta_rad = acos(cos(za_sca_rad) * cos(za_inc_rad) + 
                       sin(za_sca_rad) * sin(za_inc_rad) * 
                       cos(aa_sca_rad - aa_inc_rad));
   }
      const Numeric theta = RAD2DEG * theta_rad;
      
  // Interpolation of the data on the scattering angle:
 
      GridPos thet_gp;
      gridpos(thet_gp, za_datagrid, theta);
      
      Vector itw(2);
      interpweights(itw, thet_gp);
      
      for (Index i = 0; i < 6; i++)
        {
          pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), 
                              thet_gp);
    }
} 





void pha_mat_labCalc(//Output:
                      MatrixView pha_mat_lab,
                      //Input:
                      ConstVectorView pha_mat_int,
                      const Numeric& za_sca,
                      const Numeric& aa_sca,
                      const Numeric& za_inc,
                      const Numeric& aa_inc,
                      const Numeric& theta_rad)
{
  Numeric za_sca_rad = za_sca * DEG2RAD;
  Numeric za_inc_rad = za_inc * DEG2RAD;
  Numeric aa_sca_rad = aa_sca * DEG2RAD;
  Numeric aa_inc_rad = aa_inc * DEG2RAD;

  const Numeric theta = RAD2DEG * theta_rad;
  const Index stokes_dim = pha_mat_lab.ncols();

  // cout << "Transformation of phase matrix:" <<endl; 
  
 
  
  // For stokes_dim = 1, we only need Z11=F11:
  pha_mat_lab(0,0) = F11;
  
  if( stokes_dim > 1 ){
    //
    // Several cases have to be considered:
    //
    const Numeric ANGTOL = 1e-6; //CPD: this constant is used to adjust zenith angles 
                               //close to 0 and PI.  This is also used to avoid
                               //float == float statements.  

   if(
        // Forward scattering
        ((theta > -.01) && (theta < .01) ) ||
        // Backward scattering
        ((theta > 179.99) && (theta < 180.01)) ||
        // "Grosskreis" through poles: no rotation required
        ((aa_sca == aa_inc) || (aa_sca == 360-aa_inc) || (aa_inc == 360-aa_sca) ||
         (aa_sca == 180-aa_inc) || (aa_inc == 180-aa_sca) )  
        )
      {
        pha_mat_lab(0,1) = F12;
        pha_mat_lab(1,0) = F12;
        pha_mat_lab(1,1) = F22;
        
        if( stokes_dim > 2 ){
          pha_mat_lab(0,2) = 0;
          pha_mat_lab(1,2) = 0;
          pha_mat_lab(2,0) = 0;
          pha_mat_lab(2,1) = 0;
          pha_mat_lab(2,2) = F33;
          
          if( stokes_dim > 3 ){
            pha_mat_lab(0,3) = 0;
            pha_mat_lab(1,3) = 0;
            pha_mat_lab(2,3) = F34;
            pha_mat_lab(3,0) = 0;
            pha_mat_lab(3,1) = 0;
            pha_mat_lab(3,2) = -F34;
            pha_mat_lab(3,3) = F44;
          }
        }
      }
   
   else 
     {
       Numeric sigma1;
       Numeric sigma2;

       Numeric s1, s2;

       // In these cases we have to take limiting values.
 
       if (za_inc_rad < ANGTOL)
         {
           sigma1 = PI + aa_sca_rad - aa_inc_rad;
           sigma2 = 0;
         }
       else if (za_inc_rad > PI-ANGTOL)
         {
           sigma1 = aa_sca_rad - aa_inc_rad;
           sigma2 = PI; 
         }
       else if (za_sca_rad < ANGTOL)
         {
           sigma1 = 0;
           sigma2 = PI + aa_sca_rad - aa_inc_rad;
         }
       else if (za_sca_rad > PI - ANGTOL)
         {
           sigma1 = PI;
           sigma2 = aa_sca_rad - aa_inc_rad; 
         }
       else
         {
           s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad))
              /(sin(za_inc_rad)*sin(theta_rad));
           s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos (theta_rad))/
             (sin(za_sca_rad)*sin(theta_rad)); 
       
           sigma1 =  acos(s1);
           sigma2 =  acos(s2);
           
           // Arccos is only defined in the range from -1 ... 1
           // Numerical problems can appear for values close to 1 or -1      
           if ( isnan(sigma1) || isnan(sigma2) )
             {
               if ( abs(s1 - 1) < ANGTOL)
                 sigma1 = 0;
               if ( abs(s1 + 1) < ANGTOL)
                 sigma1 = PI;
               if ( abs(s2 - 1) < ANGTOL)
                 sigma2 = 0;
               if ( abs(s2 + 1) < ANGTOL)
                 sigma2 = PI;
             }
         }
      
       const Numeric C1 = cos(2*sigma1);
       const Numeric C2 = cos(2*sigma2);
        
       const Numeric S1 = sin(2*sigma1); 
       const Numeric S2 = sin(2*sigma2);
        
        pha_mat_lab(0,1) = C1 * F12;
        pha_mat_lab(1,0) = C2 * F12;
        pha_mat_lab(1,1) = C1 * C2 * F22 - S1 * S2 * F33;
        
        assert(!isnan(pha_mat_lab(0,1)));        
        assert(!isnan(pha_mat_lab(1,0)));
        assert(!isnan(pha_mat_lab(1,1)));

        if( stokes_dim > 2 ){
          /*CPD: For skokes_dim > 2 some of the transformation formula 
            for each element have a different sign depending on whether or
            not 0<aa_scat-aa_inc<180.  For details see pages 94 and 95 of 
            Mishchenkos chapter in : 
            Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic 
            scattering by nonspherical particles. In Exploring the Atmosphere 
            by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag, 
            Berlin, pp. 77-127. 
            This is available at http://www.giss.nasa.gov/~crmim/publications/ */
          Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
            (aa_sca-aa_inc>180)*360;
          if(delta_aa>=0)
            {
              pha_mat_lab(0,2) = S1 * F12;
              pha_mat_lab(1,2) = S1 * C2 * F22 + C1 * S2 * F33;
              pha_mat_lab(2,0) = -S2 * F12;
              pha_mat_lab(2,1) = -C1 * S2 * F22 - S1 * C2 * F33;
            }
          else
            {
              pha_mat_lab(0,2) = -S1 * F12;
              pha_mat_lab(1,2) = -S1 * C2 * F22 - C1 * S2 * F33;
              pha_mat_lab(2,0) = S2 * F12;
              pha_mat_lab(2,1) = C1 * S2 * F22 + S1 * C2 * F33;
            }
          pha_mat_lab(2,2) = -S1 * S2 * F22 + C1 * C2 * F33;

          if( stokes_dim > 3 ){
            if(delta_aa>=0)
              {
                pha_mat_lab(1,3) = S2 * F34;
                pha_mat_lab(3,1) = S1 * F34;
              }
            else
              {
                pha_mat_lab(1,3) = -S2 * F34;
                pha_mat_lab(3,1) = -S1 * F34;
              }
            pha_mat_lab(0,3) = 0;
            pha_mat_lab(2,3) = C2 * F34;
            pha_mat_lab(3,0) = 0;
            pha_mat_lab(3,2) = -C1 * F34;
            pha_mat_lab(3,3) = F44;
          }
        }     
     }
   }
}
     

ostream& operator<< (ostream& os, const SingleScatteringData& /*ssd*/)
{
  os << "SingleScatteringData: Output operator not implemented";
  return os;
}


ostream& operator<< (ostream& os, const ArrayOfSingleScatteringData& /*assd*/)
{
  os << "ArrayOfSingleScatteringData: Output operator not implemented";
  return os;
}

ostream& operator<< (ostream& os, const ScatteringMetaData& /*ssd*/)
{
  os << "ScatteringMetaData: Output operator not implemented";
  return os;
}


ostream& operator<< (ostream& os, const ArrayOfScatteringMetaData& /*assd*/)
{
  os << "ArrayOfScatteringMetaData: Output operator not implemented";
  return os;
}



void opt_prop_sum_propmat_clearsky(//Output:
                                      Tensor3&         ext_mat,
                                      Matrix&          abs_vec,
                                      //Input:
                                      const Tensor4    propmat_clearsky)
{

    Index stokes_dim = propmat_clearsky.ncols();
    
    Index freq_dim = propmat_clearsky.npages();

    // old abs_vecInit
    abs_vec.resize( freq_dim, stokes_dim );
    abs_vec = 0;                  // Initialize to zero!

    // old ext_matInit
    ext_mat.resize( freq_dim, stokes_dim, stokes_dim );
    ext_mat = 0;                  // Initialize to zero!

   // old ext_matAddGas and abs_vecAddGas for 0 vector and matrix
    for ( Index iv=0; iv<freq_dim; ++iv )
        for ( Index is1=0; is1<stokes_dim; ++is1 )
        {
            abs_vec(iv,is1) += propmat_clearsky(joker,iv,is1,0).sum();
            for ( Index is2=0; is2<stokes_dim; ++is2 )
                ext_mat(iv,is1,is2) += propmat_clearsky(joker,iv,is1,is2).sum();
        }
}



ParticleType ParticleTypeFromString(const String& particle_type_string)
{
    ParticleType particle_type;
    if (particle_type_string == "general")
        particle_type = PARTICLE_TYPE_GENERAL;
    else if (particle_type_string == "macroscopically_isotropic")
        particle_type = PARTICLE_TYPE_MACROS_ISO;
    else if (particle_type_string == "horizontally_aligned")
        particle_type = PARTICLE_TYPE_HORIZ_AL;
    else if (particle_type_string == "spherical")
        particle_type = PARTICLE_TYPE_SPHERICAL;
    else
    {
        ostringstream os;
        os << "Unknown particle type: " << particle_type_string << endl
           << "Valid types are: general, macroscopically_isotropic, "
           << "horizontally_aligned and spherical.";
        throw std::runtime_error(os.str());
    }

    return particle_type;
}



String ParticleTypeToString(const ParticleType& particle_type)
{
    String particle_type_string;

    switch (particle_type) {
        case PARTICLE_TYPE_GENERAL:
            particle_type_string = "general";
            break;
        case PARTICLE_TYPE_MACROS_ISO:
            particle_type_string = "macroscopically_isotropic";
            break;
        case PARTICLE_TYPE_HORIZ_AL:
            particle_type_string = "horizontally_aligned";
            break;
        case PARTICLE_TYPE_SPHERICAL:
            particle_type_string = "spherical";
            break;
        default:
            ostringstream os;
            os << "Internal error: Cannot map ParticleType enum value "
            << particle_type << " to String.";
            throw std::runtime_error(os.str());
            break;
    }

    return particle_type_string;
}



ParticleSSDMethod ParticleSSDMethodFromString(const String& particle_ssdmethod_string)
{
    ParticleSSDMethod particle_ssdmethod;
    if (particle_ssdmethod_string == "tmatrix")
        particle_ssdmethod = PARTICLE_SSDMETHOD_TMATRIX;
    else
    {
        ostringstream os;
        os << "Unknown particle SSD method: " << particle_ssdmethod_string << endl
        << "Valid methods: tmatrix";
        throw std::runtime_error(os.str());
    }

    return particle_ssdmethod;
}



String ParticleTypeToString(const ParticleSSDMethod& particle_ssdmethod)
{
    String particle_ssdmethod_string;

    switch (particle_ssdmethod) {
        case PARTICLE_SSDMETHOD_TMATRIX:
            particle_ssdmethod_string = "tmatrix";
            break;
        default:
            ostringstream os;
            os << "Internal error: Cannot map ParticleSSDMethod enum value "
            << particle_ssdmethod << " to String.";
            throw std::runtime_error(os.str());
            break;
    }

    return particle_ssdmethod_string;
}