ARTS: optproperties.cc Source File

00001 /* Copyright (C) 2003-2008 Claudia Emde <claudia.emde@dlr.de>
00002 
00003    This program is free software; you can redistribute it and/or modify it
00004    under the terms of the GNU General Public License as published by the
00005    Free Software Foundation; either version 2, or (at your option) any
00006    later version.
00007 
00008    This program is distributed in the hope that it will be useful,
00009    but WITHOUT ANY WARRANTY; without even the implied warranty of
00010    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011    GNU General Public License for more details.
00012 
00013    You should have received a copy of the GNU General Public License
00014    along with this program; if not, write to the Free Software
00015    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
00016    USA. */
00017 
00018 
00019 
00020 /*===========================================================================
00021   === File description
00022   ===========================================================================*/
00023 
00036 /*===========================================================================
00037   === External declarations
00038   ===========================================================================*/
00039 
00040 #include <cmath>
00041 #include <stdexcept>
00042 #include "arts.h"
00043 #include "matpackVII.h"
00044 #include "array.h"
00045 #include "math_funcs.h"
00046 #include "messages.h"
00047 #include "logic.h"
00048 #include "interpolation.h"
00049 #include "optproperties.h"
00050 #include "xml_io.h"
00051 extern const Numeric DEG2RAD;
00052 extern const Numeric RAD2DEG;
00053 extern const Numeric PI;
00054 
00055 
00056 #define F11 pha_mat_int[0]
00057 #define F12 pha_mat_int[1]
00058 #define F22 pha_mat_int[2]
00059 #define F33 pha_mat_int[3]
00060 #define F34 pha_mat_int[4]
00061 #define F44 pha_mat_int[5]
00062 
00064 
00085 void abs_vecTransform(//Output and Input
00086                       VectorView abs_vec_lab,
00087                       //Input
00088                       ConstTensor3View abs_vec_data,
00089                       ConstVectorView za_datagrid,
00090                       ConstVectorView aa_datagrid _U_,
00091                       const PType& ptype,
00092                       const Numeric& za_sca _U_,
00093                       const Numeric& aa_sca _U_)
00094 {
00095  const Index stokes_dim = abs_vec_lab.nelem();
00096     
00097   if (stokes_dim > 4 || stokes_dim < 1){
00098     throw runtime_error("The dimension of the stokes vector \n"
00099                          "must be 1,2,3 or 4");
00100   }
00101 
00102   switch (ptype){
00103 
00104   case PTYPE_GENERAL:
00105 
00106     out0 << "Case PTYPE_GENERAL not yet implemented. \n"; 
00107     break;
00108     
00109   case PTYPE_MACROS_ISO:
00110     {
00111       // The first element of the vector corresponds to the absorption 
00112       // coefficient which is stored in the database, the others are 0.
00113       
00114       abs_vec_lab = 0;
00115 
00116       abs_vec_lab[0] = abs_vec_data(0,0,0);
00117       break;
00118     }
00119 
00120 case PTYPE_HORIZ_AL://Added by Cory Davis 9/12/03
00121     {
00122       assert (abs_vec_data.ncols() == 2);
00123       
00124       // In the case of horizontally oriented particles the absorption
00125       // coefficient vector only the first two elements are non-zero.
00126       // These values are dependent on the zenith angle of propagation. The 
00127       // data storage format also makes use of the fact that in this case
00128       //K_abs(za_sca)=K_abs(180-za_sca). 
00129 
00130       // 1st interpolate data by za_sca
00131       GridPos gp;
00132       Vector itw(2);
00133       if (za_sca>90)
00134         {
00135           gridpos(gp,za_datagrid,180-za_sca);
00136         }
00137       else
00138         {
00139           gridpos(gp,za_datagrid,za_sca);
00140         }
00141       interpweights(itw,gp);
00142       abs_vec_lab = 0;
00143       abs_vec_lab[0] = interp(itw,abs_vec_data(Range(joker),0,0),gp);
00144 
00145       if( stokes_dim == 1 ){
00146         break;
00147       }
00148       abs_vec_lab[1] = interp(itw,abs_vec_data(Range(joker),0,1),gp);
00149       break;
00150     }
00151 
00152    default:
00153     out0 << "Not all particle type cases are implemented\n";
00154     
00155   }  
00156     
00157 }
00158 
00159 
00161 
00182 void ext_matTransform(//Output and Input
00183                       MatrixView ext_mat_lab,
00184                       //Input
00185                       ConstTensor3View ext_mat_data,
00186                       ConstVectorView za_datagrid,
00187                       ConstVectorView aa_datagrid _U_,
00188                       const PType& ptype,
00189                       const Numeric& za_sca,
00190                       const Numeric& aa_sca _U_)
00191 {
00192  const Index stokes_dim = ext_mat_lab.ncols();
00193     
00194   if (stokes_dim > 4 || stokes_dim < 1){
00195     throw runtime_error("The dimension of the stokes vector \n"
00196                          "must be 1,2,3 or 4");
00197   }
00198 
00199   switch (ptype){
00200 
00201   case PTYPE_GENERAL:
00202 
00203     out0 << "Case PTYPE_GENERAL not yet implemented. \n"; 
00204     break;
00205     
00206   case PTYPE_MACROS_ISO:
00207     {
00208       assert (ext_mat_data.ncols() == 1);
00209       
00210       // In the case of randomly oriented particles the extinction matrix is 
00211       // diagonal. The value of each element of the diagonal is the
00212       // extinction cross section, which is stored in the database.
00213      
00214       ext_mat_lab = 0.;
00215       
00216       ext_mat_lab(0,0) = ext_mat_data(0,0,0);
00217       
00218       
00219       if( stokes_dim == 1 ){
00220         break;
00221       }
00222       
00223       ext_mat_lab(1,1) = ext_mat_data(0,0,0);
00224       
00225       if( stokes_dim == 2 ){
00226         break;
00227       }
00228       
00229       ext_mat_lab(2,2) = ext_mat_data(0,0,0);
00230       
00231       if( stokes_dim == 3 ){
00232         break;
00233       }
00234       
00235       ext_mat_lab(3,3) = ext_mat_data(0,0,0);
00236       break;
00237     }
00238 
00239   case PTYPE_HORIZ_AL://Added by Cory Davis 9/12/03
00240     {
00241       assert (ext_mat_data.ncols() == 3);
00242       
00243       // In the case of horizontally oriented particles the extinction matrix
00244       // has only 3 independent non-zero elements Kjj, K12=K21, and K34=-K43.
00245       // These values are dependent on the zenith angle of propagation. The 
00246       // data storage format also makes use of the fact that in this case
00247       //K(za_sca)=K(180-za_sca). 
00248 
00249       // 1st interpolate data by za_sca
00250       GridPos gp;
00251       Vector itw(2);
00252       Numeric Kjj;
00253       Numeric K12;
00254       Numeric K34;
00255       
00256       if (za_sca>90)
00257         {
00258           gridpos(gp,za_datagrid,180-za_sca);
00259         }
00260       else
00261         {
00262           gridpos(gp,za_datagrid,za_sca);
00263         }
00264 
00265       interpweights(itw,gp);
00266       
00267       ext_mat_lab=0.0;
00268       Kjj=interp(itw,ext_mat_data(Range(joker),0,0),gp);
00269       ext_mat_lab(0,0)=Kjj;
00270 
00271       if( stokes_dim == 1 ){
00272         break;
00273       }
00274       
00275       K12=interp(itw,ext_mat_data(Range(joker),0,1),gp);
00276       ext_mat_lab(1,1)=Kjj;
00277       ext_mat_lab(0,1)=K12;
00278       ext_mat_lab(1,0)=K12;
00279 
00280       if( stokes_dim == 2 ){
00281         break;
00282       }
00283       
00284       ext_mat_lab(2,2)=Kjj;
00285       
00286       if( stokes_dim == 3 ){
00287         break;
00288       }
00289       
00290       K34=interp(itw,ext_mat_data(Range(joker),0,2),gp);
00291       ext_mat_lab(2,3)=K34;
00292       ext_mat_lab(3,2)=-K34;
00293       ext_mat_lab(3,3)=Kjj;
00294       break;
00295 
00296     }
00297   default:
00298     out0 << "Not all particle type cases are implemented\n";
00299     
00300   }
00301 }  
00302  
00303 
00305 
00328 void pha_matTransform(//Output
00329                       MatrixView pha_mat_lab,
00330                       //Input
00331                       ConstTensor5View pha_mat_data,
00332                       ConstVectorView za_datagrid,
00333                       ConstVectorView aa_datagrid,
00334                       const PType& ptype,
00335                       const Index& za_sca_idx,
00336                       const Index& aa_sca_idx,
00337                       const Index& za_inc_idx,
00338                       const Index& aa_inc_idx,
00339                       ConstVectorView scat_za_grid,
00340                       ConstVectorView scat_aa_grid
00341                       )
00342 {
00343 
00344   
00345   const Index stokes_dim = pha_mat_lab.ncols();
00346   
00347   Numeric za_sca = scat_za_grid[za_sca_idx]; 
00348   Numeric aa_sca = scat_aa_grid[aa_sca_idx];
00349   Numeric za_inc = scat_za_grid[za_inc_idx]; 
00350   Numeric aa_inc = scat_aa_grid[aa_inc_idx];
00351                
00352   if (stokes_dim > 4 || stokes_dim < 1){
00353     throw runtime_error("The dimension of the stokes vector \n"
00354                          "must be 1,2,3 or 4");
00355   }
00356 
00357   switch (ptype){
00358 
00359   case PTYPE_GENERAL:
00360 
00361     out0 << "Case PTYPE_GENERAL not yet implemented. \n"; 
00362     break;
00363     
00364   case PTYPE_MACROS_ISO:
00365     {
00366       // Calculate the scattering and interpolate the data on the scattering
00367       // angle:
00368       
00369       Vector pha_mat_int(6);
00370       Numeric theta_rad;
00371             
00372       // Interpolation of the data on the scattering angle:
00373       interpolate_scat_angle(pha_mat_int, theta_rad, pha_mat_data,
00374                              za_datagrid, za_sca, aa_sca,
00375                              za_inc, aa_inc);
00376       
00377       // Caclulate the phase matrix in the laboratory frame:
00378       pha_mat_labCalc(pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, 
00379                       aa_inc, theta_rad);
00380       
00381       break;
00382     }
00383 
00384   case PTYPE_HORIZ_AL://Added by Cory Davis
00385     //Data is already stored in the laboratory frame, but it is compressed
00386     //a little.  Details elsewhere
00387     {
00388       assert (pha_mat_data.ncols()==16);
00389       Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
00390         (aa_sca-aa_inc>180)*360;//delta_aa corresponds to the "pages" 
00391                                 //dimension of pha_mat_data
00392       GridPos za_sca_gp;
00393       GridPos delta_aa_gp;
00394       GridPos za_inc_gp;
00395       Vector itw(8);
00396 
00397       gridpos(delta_aa_gp,aa_datagrid,abs(delta_aa));
00398       if (za_inc>90)
00399         {
00400           gridpos(za_inc_gp,za_datagrid,180-za_inc);
00401           gridpos(za_sca_gp,za_datagrid,180-za_sca);
00402         }
00403       else
00404         {
00405           gridpos(za_inc_gp,za_datagrid,za_inc);
00406           gridpos(za_sca_gp,za_datagrid,za_sca);
00407         }
00408 
00409       interpweights(itw,za_sca_gp,delta_aa_gp,za_inc_gp);
00410 
00411       pha_mat_lab(0,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00412                                                Range(joker),0,0),
00413                               za_sca_gp,delta_aa_gp,za_inc_gp);
00414       if( stokes_dim == 1 ){
00415         break;
00416       }
00417       pha_mat_lab(0,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00418                                                Range(joker),0,1),
00419                               za_sca_gp,delta_aa_gp,za_inc_gp);
00420       pha_mat_lab(1,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00421                                                Range(joker),0,4),
00422                               za_sca_gp,delta_aa_gp,za_inc_gp);
00423       pha_mat_lab(1,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00424                                                Range(joker),0,5),
00425                               za_sca_gp,delta_aa_gp,za_inc_gp);
00426       if( stokes_dim == 2 ){
00427         break;
00428       }
00429       if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
00430         {
00431           pha_mat_lab(0,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00432                                                    Range(joker),0,2),
00433                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00434           pha_mat_lab(1,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00435                                                    Range(joker),0,6),
00436                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00437           pha_mat_lab(2,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00438                                                    Range(joker),0,8),
00439                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00440           pha_mat_lab(2,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00441                                                    Range(joker),0,9),
00442                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00443         }
00444       else
00445         {
00446           pha_mat_lab(0,2)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00447                                                    Range(joker),0,2),
00448                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00449           pha_mat_lab(1,2)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00450                                                    Range(joker),0,6),
00451                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00452           pha_mat_lab(2,0)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00453                                                    Range(joker),0,8),
00454                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00455           pha_mat_lab(2,1)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00456                                                    Range(joker),0,9),
00457                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00458         }                             
00459       pha_mat_lab(2,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00460                                                Range(joker),0,10),
00461                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00462       if( stokes_dim == 3 ){
00463         break;
00464       }
00465       if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
00466         {
00467           pha_mat_lab(0,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00468                                                    Range(joker),0,3),
00469                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00470           pha_mat_lab(1,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00471                                                    Range(joker),0,7),
00472                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00473           pha_mat_lab(3,0)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00474                                                    Range(joker),0,12),
00475                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00476           pha_mat_lab(3,1)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00477                                                    Range(joker),0,13),
00478                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00479         }
00480       else
00481         {
00482           pha_mat_lab(0,3)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00483                                                    Range(joker),0,3),
00484                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00485           pha_mat_lab(1,3)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00486                                                    Range(joker),0,7),
00487                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00488           pha_mat_lab(3,0)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00489                                                    Range(joker),0,12),
00490                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00491           pha_mat_lab(3,1)=-interp(itw,pha_mat_data(Range(joker),Range(joker),
00492                                                    Range(joker),0,13),
00493                                   za_sca_gp,delta_aa_gp,za_inc_gp);
00494         }
00495       pha_mat_lab(2,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00496                                                Range(joker),0,11),
00497                               za_sca_gp,delta_aa_gp,za_inc_gp);
00498       pha_mat_lab(3,2)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00499                                                Range(joker),0,14),
00500                               za_sca_gp,delta_aa_gp,za_inc_gp);
00501       pha_mat_lab(3,3)=interp(itw,pha_mat_data(Range(joker),Range(joker),
00502                                                Range(joker),0,15),
00503                               za_sca_gp,delta_aa_gp,za_inc_gp);
00504       break;
00505       
00506     }  
00507   default:
00508     out0 << "Not all particle type cases are implemented\n";
00509     
00510   }
00511 }
00512 
00513 
00514 
00515 
00517 
00541 void interpolate_scat_angleDOIT(//Output:
00542                             VectorView pha_mat_int,
00543                             //Input:
00544                             ConstTensor5View pha_mat_data,
00545                             const Index& za_sca_idx,
00546                             const Index& aa_sca_idx,
00547                             const Index& za_inc_idx,
00548                             const Index& aa_inc_idx,
00549                             const ArrayOfArrayOfArrayOfArrayOfGridPos&
00550                                 scat_theta_gps,
00551                             ConstTensor5View scat_theta_itws
00552                             )
00553 {
00554   
00555   ConstVectorView itw = scat_theta_itws(za_sca_idx, aa_sca_idx, za_inc_idx, aa_inc_idx, joker);
00556   
00557   for (Index i = 0; i < 6; i++)
00558     {
00559       pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), 
00560                              scat_theta_gps[za_sca_idx][aa_sca_idx][za_inc_idx][aa_inc_idx]);
00561     }
00562   
00563 } 
00564 
00565 
00567 
00590 void interpolate_scat_angle(//Output:
00591                             VectorView pha_mat_int,
00592                             Numeric& theta_rad,
00593                             //Input:
00594                             ConstTensor5View pha_mat_data,
00595                             ConstVectorView za_datagrid,
00596                             const Numeric& za_sca,
00597                             const Numeric& aa_sca,
00598                             const Numeric& za_inc,
00599                             const Numeric& aa_inc
00600                             )
00601 {
00602   Numeric ANG_TOL=1e-7;
00603 
00604   //Calculate scattering angle from incident and scattered directions.
00605   //The two special cases are implemented here to avoid NaNs that can 
00606   //sometimes occur in in the acos... formula in forward and backscatter
00607   //cases. CPD 5/10/03.
00608   
00609   if(abs(aa_sca-aa_inc)<ANG_TOL)
00610     {
00611       theta_rad=DEG2RAD*abs(za_sca-za_inc);
00612     }
00613   else if (abs(abs(aa_sca-aa_inc)-180)<ANG_TOL)
00614     {
00615       theta_rad=DEG2RAD*(za_sca+za_inc);
00616       if (theta_rad>PI){theta_rad=2*PI-theta_rad;}
00617     }
00618   else
00619     {
00620       const Numeric za_sca_rad = za_sca * DEG2RAD;
00621       const Numeric za_inc_rad = za_inc * DEG2RAD;
00622       const Numeric aa_sca_rad = aa_sca * DEG2RAD;
00623       const Numeric aa_inc_rad = aa_inc * DEG2RAD;
00624       
00625       // cout << "Interpolation on scattering angle" << endl;
00626       assert (pha_mat_data.ncols() == 6);
00627       // Calculation of the scattering angle:
00628       theta_rad = acos(cos(za_sca_rad) * cos(za_inc_rad) + 
00629                        sin(za_sca_rad) * sin(za_inc_rad) * 
00630                        cos(aa_sca_rad - aa_inc_rad));
00631    }
00632       const Numeric theta = RAD2DEG * theta_rad;
00633       
00634   // Interpolation of the data on the scattering angle:
00635  
00636       GridPos thet_gp;
00637       gridpos(thet_gp, za_datagrid, theta);
00638       
00639       Vector itw(2);
00640       interpweights(itw, thet_gp);
00641       
00642       for (Index i = 0; i < 6; i++)
00643         {
00644           pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), 
00645                               thet_gp);
00646     }
00647 } 
00648 
00649 
00650 
00651 
00653 
00678 void pha_mat_labCalc(//Output:
00679                       MatrixView pha_mat_lab,
00680                       //Input:
00681                       ConstVectorView pha_mat_int,
00682                       const Numeric& za_sca,
00683                       const Numeric& aa_sca,
00684                       const Numeric& za_inc,
00685                       const Numeric& aa_inc,
00686                       const Numeric& theta_rad)
00687 {
00688   Numeric za_sca_rad = za_sca * DEG2RAD;
00689   Numeric za_inc_rad = za_inc * DEG2RAD;
00690   Numeric aa_sca_rad = aa_sca * DEG2RAD;
00691   Numeric aa_inc_rad = aa_inc * DEG2RAD;
00692 
00693   const Numeric theta = RAD2DEG * theta_rad;
00694   const Index stokes_dim = pha_mat_lab.ncols();
00695 
00696   // cout << "Transformation of phase matrix:" <<endl; 
00697   
00698  
00699   
00700   // For stokes_dim = 1, we only need Z11=F11:
00701   pha_mat_lab(0,0) = F11;
00702   
00703   if( stokes_dim > 1 ){
00704     //
00705     // Several cases have to be considered:
00706     //
00707     const Numeric ANGTOL = 1e-6; //CPD: this constant is used to adjust zenith angles 
00708                                //close to 0 and PI.  This is also used to avoid
00709                                //float == float statements.  
00710 
00711    if(
00712         // Forward scattering
00713         ((theta > -.01) && (theta < .01) ) ||
00714         // Backward scattering
00715         ((theta > 179.99) && (theta < 180.01)) ||
00716         // "Grosskreis" through poles: no rotation required
00717         ((aa_sca == aa_inc) || (aa_sca == 360-aa_inc) || (aa_inc == 360-aa_sca) ||
00718          (aa_sca == 180-aa_inc) || (aa_inc == 180-aa_sca) )  
00719         )
00720       {
00721         pha_mat_lab(0,1) = F12;
00722         pha_mat_lab(1,0) = F12;
00723         pha_mat_lab(1,1) = F22;
00724         
00725         if( stokes_dim > 2 ){
00726           pha_mat_lab(0,2) = 0;
00727           pha_mat_lab(1,2) = 0;
00728           pha_mat_lab(2,0) = 0;
00729           pha_mat_lab(2,1) = 0;
00730           pha_mat_lab(2,2) = F33;
00731           
00732           if( stokes_dim > 3 ){
00733             pha_mat_lab(0,3) = 0;
00734             pha_mat_lab(1,3) = 0;
00735             pha_mat_lab(2,3) = F34;
00736             pha_mat_lab(3,0) = 0;
00737             pha_mat_lab(3,1) = 0;
00738             pha_mat_lab(3,2) = -F34;
00739             pha_mat_lab(3,3) = F44;
00740           }
00741         }
00742       }
00743    
00744    else 
00745      {
00746        Numeric sigma1;
00747        Numeric sigma2;
00748 
00749        Numeric s1, s2;
00750 
00751        // In these cases we have to take limiting values.
00752  
00753        if (za_inc_rad < ANGTOL)
00754          {
00755            sigma1 = PI + aa_sca_rad - aa_inc_rad;
00756            sigma2 = 0;
00757          }
00758        else if (za_inc_rad > PI-ANGTOL)
00759          {
00760            sigma1 = aa_sca_rad - aa_inc_rad;
00761            sigma2 = PI; 
00762          }
00763        else if (za_sca_rad < ANGTOL)
00764          {
00765            sigma1 = 0;
00766            sigma2 = PI + aa_sca_rad - aa_inc_rad;
00767          }
00768        else if (za_sca_rad > PI - ANGTOL)
00769          {
00770            sigma1 = PI;
00771            sigma2 = aa_sca_rad - aa_inc_rad; 
00772          }
00773        else
00774          {
00775            s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad))
00776               /(sin(za_inc_rad)*sin(theta_rad));
00777            s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos (theta_rad))/
00778              (sin(za_sca_rad)*sin(theta_rad)); 
00779        
00780            sigma1 =  acos(s1);
00781            sigma2 =  acos(s2);
00782            
00783            // Arccos is only defined in the range from -1 ... 1
00784            // Numerical problems can appear for values close to 1 or -1      
00785            if ( isnan(sigma1) || isnan(sigma2) )
00786              {
00787                if ( abs(s1 - 1) < ANGTOL)
00788                  sigma1 = 0;
00789                if ( abs(s1 + 1) < ANGTOL)
00790                  sigma1 = PI;
00791                if ( abs(s2 - 1) < ANGTOL)
00792                  sigma2 = 0;
00793                if ( abs(s2 + 1) < ANGTOL)
00794                  sigma2 = PI;
00795              }
00796          }
00797       
00798        const Numeric C1 = cos(2*sigma1);
00799        const Numeric C2 = cos(2*sigma2);
00800         
00801        const Numeric S1 = sin(2*sigma1); 
00802        const Numeric S2 = sin(2*sigma2);
00803         
00804         pha_mat_lab(0,1) = C1 * F12;
00805         pha_mat_lab(1,0) = C2 * F12;
00806         pha_mat_lab(1,1) = C1 * C2 * F22 - S1 * S2 * F33;
00807         
00808         assert(!isnan(pha_mat_lab(0,1)));        
00809         assert(!isnan(pha_mat_lab(1,0)));
00810         assert(!isnan(pha_mat_lab(1,1)));
00811 
00812         if( stokes_dim > 2 ){
00813           /*CPD: For skokes_dim > 2 some of the transformation formula 
00814             for each element have a different sign depending on whether or
00815             not 0<aa_scat-aa_inc<180.  For details see pages 94 and 95 of 
00816             Mishchenkos chapter in : 
00817             Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic 
00818             scattering by nonspherical particles. In Exploring the Atmosphere 
00819             by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag, 
00820             Berlin, pp. 77-127. 
00821             This is available at http://www.giss.nasa.gov/~crmim/publications/ */
00822           Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
00823             (aa_sca-aa_inc>180)*360;
00824           if(delta_aa>=0)
00825             {
00826               pha_mat_lab(0,2) = S1 * F12;
00827               pha_mat_lab(1,2) = S1 * C2 * F22 + C1 * S2 * F33;
00828               pha_mat_lab(2,0) = -S2 * F12;
00829               pha_mat_lab(2,1) = -C1 * S2 * F22 - S1 * C2 * F33;
00830             }
00831           else
00832             {
00833               pha_mat_lab(0,2) = -S1 * F12;
00834               pha_mat_lab(1,2) = -S1 * C2 * F22 - C1 * S2 * F33;
00835               pha_mat_lab(2,0) = S2 * F12;
00836               pha_mat_lab(2,1) = C1 * S2 * F22 + S1 * C2 * F33;
00837             }
00838           pha_mat_lab(2,2) = -S1 * S2 * F22 + C1 * C2 * F33;
00839 
00840           if( stokes_dim > 3 ){
00841             if(delta_aa>=0)
00842               {
00843                 pha_mat_lab(1,3) = S2 * F34;
00844                 pha_mat_lab(3,1) = S1 * F34;
00845               }
00846             else
00847               {
00848                 pha_mat_lab(1,3) = -S2 * F34;
00849                 pha_mat_lab(3,1) = -S1 * F34;
00850               }
00851             pha_mat_lab(0,3) = 0;
00852             pha_mat_lab(2,3) = C2 * F34;
00853             pha_mat_lab(3,0) = 0;
00854             pha_mat_lab(3,2) = -C1 * F34;
00855             pha_mat_lab(3,3) = F44;
00856           }
00857         }     
00858      }
00859    }
00860 }
00861      
00862 
00863 ostream& operator<< (ostream& os, const SingleScatteringData& /*ssd*/)
00864 {
00865   os << "SingleScatteringData: Output operator not implemented";
00866   return os;
00867 }
00868 
00869 
00870 ostream& operator<< (ostream& os, const ArrayOfSingleScatteringData& /*assd*/)
00871 {
00872   os << "ArrayOfSingleScatteringData: Output operator not implemented";
00873   return os;
00874 }
00875