absorption.h

Go to the documentation of this file.

/* Copyright (C) 2000-2008
   Stefan Buehler <sbuehler@ltu.se>
   Axel von Engeln <engeln@uni-bremen.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. */
 
#ifndef absorption_h
#define absorption_h
 
#include <stdexcept>
#include "matpackI.h"
#include "array.h"
#include "mystring.h"
#include "make_array.h"
#include "messages.h"
 
typedef void (*lsf_type)(Vector&,
                         Vector&,
                         Numeric,
                         Numeric,
                         Numeric,
                         VectorView,
                         const Index);
 
class LineshapeRecord{
public:
 
  LineshapeRecord() : mname(),
                      mdescription(),
                      mfunction()
  { /* Nothing to do here. */ }
 
  LineshapeRecord(const String& name,
                  const String& description,
                  lsf_type      function)
    : mname(name),
      mdescription(description),
      mfunction(function)
  { /* Nothing to do here. */ }
  const String&  Name()        const { return mname;        }   
  const String&  Description() const { return mdescription; }
  lsf_type Function() const { return mfunction; }
private:        
  String  mname;        
  String  mdescription; 
  lsf_type mfunction;   
 
};
 
typedef void (*lsnf_type)(Vector&,
                          Numeric,
                          VectorView,
                          const Numeric,
                          const Index);
 
class LineshapeNormRecord{
public:
 
  LineshapeNormRecord() : mname(),
                          mdescription(),
                          mfunction()
  { /* Nothing to do here. */ }
 
  LineshapeNormRecord(const String& name,
                      const String& description,
                      lsnf_type      function)
    : mname(name),
      mdescription(description),
      mfunction(function)
  { /* Nothing to do here. */ }
  const String&  Name()        const { return mname;        }   
  const String&  Description() const { return mdescription; }
  lsnf_type Function() const { return mfunction; }
private:        
  String  mname;        
  String  mdescription; 
  lsnf_type mfunction;  
};
 
class LineshapeSpec{
public:
 
  LineshapeSpec() : mind_ls(-1),
                    mind_lsn(-1),
                    mcutoff(0.)
  { /* Nothing to do here. */ }
 
  LineshapeSpec(const Index&    ind_ls,
                const Index&    ind_lsn,
                const Numeric&   cutoff)
    : mind_ls(ind_ls),
      mind_lsn(ind_lsn),
      mcutoff(cutoff)
  { /* Nothing to do here. */ }
 
  const Index&  Ind_ls()        const { return mind_ls; }   
  void SetInd_ls( Index ind_ls ) { mind_ls = ind_ls; }
 
  const Index&  Ind_lsn()       const { return mind_lsn; }
  void SetInd_lsn( Index ind_lsn ) { mind_lsn = ind_lsn; }
 
  const Numeric& Cutoff() const { return mcutoff; }
  void SetCutoff( Numeric cutoff ) { mcutoff = cutoff; }
private:        
  Index  mind_ls;
  Index  mind_lsn;
  Numeric mcutoff;
};
 
ostream& operator<< (ostream& os, const LineshapeSpec& lsspec);
 
typedef Array<LineshapeSpec> ArrayOfLineshapeSpec;
 
 
 
class IsotopeRecord{
public:
 
  IsotopeRecord() : mname(),
                    mabundance(0.),
                    mmass(0.),
                    mmytrantag(-1),
                    mhitrantag(-1),
                    mjpltags(),
                    mqcoeff()
  { /* Nothing left to do here. */ }
 
  IsotopeRecord(const IsotopeRecord& x) :
    mname(x.mname),
    mabundance(x.mabundance),
    mmass(x.mmass),
    mmytrantag(x.mmytrantag),
    mhitrantag(x.mhitrantag),
    mjpltags(x.mjpltags),
    mqcoeff()
  { /* Nothing left to do here. */ }
 
  IsotopeRecord(const String&           name,
                const Numeric&          abundance,
                const Numeric&          mass,
                const Index&            mytrantag,
                const Index&            hitrantag,
                const MakeArray<Index>& jpltags) :
    mname(name),
    mabundance(abundance),
    mmass(mass),
    mmytrantag(mytrantag),
    mhitrantag(hitrantag),
    mjpltags(jpltags),
    mqcoeff()
  {
    // With Matpack, initialization of mjpltags from jpltags should now work correctly.
 
    // Some consistency checks whether the given data makes sense.
#ifndef NDEBUG
      {
        /* 1. All the tags must be positive or -1 */
        assert( (0<mmytrantag) || (-1==mmytrantag) );
        assert( (0<mhitrantag) || (-1==mhitrantag) );
        for ( Index i=0; i<mjpltags.nelem(); ++i )
          assert( (0<mjpltags[i]) || (-1==mjpltags[i]) );
      }
#endif // ifndef NDEBUG
  }
 
  const String&       Name()         const { return mname;  }
  const Numeric&      Abundance()    const { return mabundance; }
  const Numeric&      Mass()         const { return mmass;    }
  const Index&          MytranTag()    const { return mmytrantag;    }
  const Index&          HitranTag()    const { return mhitrantag;    }
  const ArrayOfIndex&   JplTags()      const { return mjpltags;      }
 
  void SetPartitionFctCoeff( const ArrayOfNumeric& qcoeff )
  {
    mqcoeff = qcoeff;
  }
 
 
  Numeric CalculatePartitionFctRatio( Numeric reference_temperature,
                                      Numeric actual_temperature ) const
  {
    Numeric qcoeff_at_t_ref =
      CalculatePartitionFctAtTemp( reference_temperature );
    Numeric qtemp =
      CalculatePartitionFctAtTemp( actual_temperature    );
 
    if ( qtemp < 0. ) 
      {
        ostringstream os;
        os << "Partition function of "
           << "Isotope = " << mname
           << "is unknown.";
        throw runtime_error(os.str());
      }
    return qcoeff_at_t_ref / qtemp;
  }
 
private:
 
  // calculate the partition fct at a certain temperature
  // this is only the prototyping
  Numeric CalculatePartitionFctAtTemp( Numeric temperature ) const;
 
  String mname;
  Numeric mabundance;
  Numeric mmass;
  Index mmytrantag;
  Index mhitrantag;
  ArrayOfIndex mjpltags;
  ArrayOfNumeric mqcoeff;
};
 
 
class SpeciesRecord{
public:
 
  SpeciesRecord() : mname(),
                    mdegfr(-1),
                    misotope() { /* Nothing to do here */ }
  
  SpeciesRecord(const char name[],
                const Index degfr,
                const MakeArray<IsotopeRecord>& isotope)
    : mname(name),
      mdegfr(degfr),
      misotope(isotope)
  {
 
    // Thanks to Matpack, initialization of misotope with isotope
    // should now work correctly.  
 
#ifndef NDEBUG
      {
        /* Check that the isotopes are correctly sorted. */
        for ( Index i=0; i<misotope.nelem()-1; ++i )
          {
            assert( misotope[i].Abundance() >= misotope[i+1].Abundance() );
          }
 
        /* Check that the Mytran tags are correctly sorted. */
        for ( Index i=0; i<misotope.nelem()-1; ++i )
          {
            if ( (0<misotope[i].MytranTag()) && (0<misotope[i+1].MytranTag()) )
              {
                assert( misotope[i].MytranTag() < misotope[i+1].MytranTag() );
            
                // Also check that the tags have the same base number:
                assert( misotope[i].MytranTag()/10 == misotope[i].MytranTag()/10 );
              }
          }
 
        /* Check that the Hitran tags are correctly sorted. */
        for ( Index i=0; i<misotope.nelem()-1; ++i )
          {
            if ( (0<misotope[i].HitranTag()) && (0<misotope[i+1].HitranTag()) )
              {
                assert( misotope[i].HitranTag() < misotope[i+1].HitranTag() );
            
                // Also check that the tags have the same base number:
                assert( misotope[i].HitranTag()/10 == misotope[i+1].HitranTag()/10 );
              }
          }
      }
#endif // #ifndef NDEBUG
  }
 
  const String&               Name()     const { return mname;     }   
  Index                         Degfr()    const { return mdegfr;    }
  const Array<IsotopeRecord>& Isotope()  const { return misotope;  }
  Array<IsotopeRecord>&       Isotope()        { return misotope;  }
  
private:
  String mname;
  Index mdegfr;
  Array<IsotopeRecord> misotope;
};
 
class LineRecord {
public:
 
  friend void linesElowToJoule(Array<LineRecord>& abs_lines);
 
  LineRecord()
    : mversion (3),
      mspecies (1000000),
      misotope (1000000),
      mf       (0.     ),
      mpsf     (0.     ),
      mi0      (0.     ),
      mti0     (0.     ),
      melow    (0.     ),
      magam    (0.     ),
      msgam    (0.     ),
      mnair    (0.     ),
      mnself   (0.     ),
      mtgam    (0.     ),
      maux     (       ),
      mdf      (-1.    ),
      mdi0     (-1.    ),
      mdagam   (-1.    ),
      mdsgam   (-1.    ),
      mdnair   (-1.    ),
      mdnself  (-1.    ),
      mdpsf    (-1.    )
 { /* Nothing to do here. */ }
 
  LineRecord( Index                 species,
              Index                 isotope,
              Numeric               f,
              Numeric               psf,
              Numeric               i0,
              Numeric               ti0,
              Numeric               elow,
              Numeric               agam,
              Numeric               sgam,
              Numeric               nair,
              Numeric               nself,
              Numeric               tgam,
              const ArrayOfNumeric& aux,
              Numeric               /* df */,
              Numeric               /* di0 */,
              Numeric               /* dagam */,
              Numeric               /* dsgam */,
              Numeric               /* dnair */,
              Numeric               /* dnself */,
              Numeric               /* dpsf */)
    : mversion (3),
      mspecies (species    ),
      misotope (isotope    ),
      mf       (f          ),
      mpsf     (psf        ),
      mi0      (i0         ),
      mti0     (ti0        ),
      melow    (elow       ),
      magam    (agam       ),
      msgam    (sgam       ),
      mnair    (nair       ),
      mnself   (nself      ),
      mtgam    (tgam       ), 
      maux     (aux        ),
      mdf      (-1.    ),
      mdi0     (-1.    ),
      mdagam   (-1.    ),
      mdsgam   (-1.    ),
      mdnair   (-1.    ),
      mdnself  (-1.    ),
      mdpsf    (-1.    )
  {
    // Thanks to Matpack, initialization of misotope with isotope
    // should now work correctly.  
 
    // Check if this species is legal, i.e., if species and isotope
    // data exists.
    //assert( mspecies < species_data.nelem() );
    //assert( misotope < species_data[mspecies].Isotope().nelem() );
  }
 
  String VersionString() const
  {
    ostringstream os;
    os << "ARTSCAT-" << mversion;
    return os.str();
  }
  
  Index Version() const { return mversion; }
  
  Index Species() const { return mspecies; }
 
  Index Isotope() const { return misotope; }
 
  String Name() const;
 
  const SpeciesRecord& SpeciesData() const;
 
  const IsotopeRecord& IsotopeData() const;
 
  Numeric F() const     { return mf; }
 
  void setF( Numeric new_mf ) { mf = new_mf; }
 
  Numeric Psf() const   { return mpsf; }
 
  void setPsf( Numeric new_mpsf ) { mpsf = new_mpsf; }
 
  Numeric I0() const    { return mi0; }
 
  void setI0( Numeric new_mi0 ) { mi0 = new_mi0; }
 
  Numeric Ti0() const   { return mti0; }
 
  Numeric Elow() const  { return melow; }
 
  Numeric Agam() const  { return magam; }
 
  void setAgam( Numeric new_agam ) { magam = new_agam; }
 
  Numeric Sgam() const  { return msgam; }
 
  void setSgam( Numeric new_sgam ) { msgam = new_sgam; }
 
  Numeric Nair() const  { return mnair; }
 
  void setNair( Numeric new_mnair ) { mnair = new_mnair; }
 
  Numeric Nself() const { return mnself; }
 
  void setNself( Numeric new_mnself ) { mnself = new_mnself; }
 
  Numeric Tgam() const  { return mtgam; }
 
  Index Naux() const   { return maux.nelem(); }
 
  const ArrayOfNumeric& Aux() const { return maux; }
  //
 
  Numeric dF() const  { return mdf; }
 
  Numeric dI0() const  { return mdi0; }
 
  Numeric dAgam() const  { return mdagam; }
 
  Numeric dSgam() const  { return mdsgam; }
 
  Numeric dNair() const  { return mdnair; }
 
  Numeric dNself() const { return mdnself; }
 
  Numeric dPsf() const { return mdpsf; }
 
  Numeric A() const { return ma; }
  
  Numeric G_upper() const { return mgupper; }
  
  Numeric G_lower() const { return mglower; }
  
  Numeric Gamma_self() const { return mgamma_self; }
  
  Numeric Gamma_N2() const { return mgamma_n2; }
  
  Numeric Gamma_O2() const { return mgamma_o2; }
  
  Numeric Gamma_H2O() const { return mgamma_h2o; }
  
  Numeric Gamma_CO2() const { return mgamma_co2; }
  
  Numeric Gamma_H2() const { return mgamma_h2; }
  
  Numeric Gamma_He() const { return mgamma_he; }
  
  Numeric Gam_N_self() const { return mn_self; }
  
  Numeric Gam_N_N2() const { return mn_n2; }
  
  Numeric Gam_N_O2() const { return mn_o2; }
  
  Numeric Gam_N_H2O() const { return mn_h2o; }
  
  Numeric Gam_N_CO2() const { return mn_co2; }
  
  Numeric Gam_N_H2() const { return mn_h2; }
  
  Numeric Gam_N_He() const { return mn_he; }
  
  Numeric Delta_N2() const { return mdelta_n2; }
  
  Numeric Delta_O2() const { return mdelta_o2; }
  
  Numeric Delta_H2O() const { return mdelta_h2o; }
  
  Numeric Delta_CO2() const { return mdelta_co2; }
  
  Numeric Delta_H2() const { return mdelta_h2; }
  
  Numeric Delta_He() const { return mdelta_he; }
  
  bool ReadFromHitranStream(istream& is, const Verbosity& verbosity);
 
 
 
  bool ReadFromHitran2004Stream(istream& is, const Verbosity& verbosity);
 
 
 
 
  bool ReadFromMytran2Stream(istream& is, const Verbosity& verbosity);
 
 
  bool ReadFromJplStream(istream& is, const Verbosity& verbosity);
 
  bool ReadFromArtscat3Stream(istream& is, const Verbosity& verbosity);
 
  bool ReadFromArtscat4Stream(istream& is, const Verbosity& verbosity);
  
 
private:
  // Version number:
  Index mversion;
  // Molecular species index: 
  Index mspecies;
  // Isotopic species index:
  Index misotope;
  // The line center frequency in Hz:
  Numeric mf;
  // The pressure shift parameter in Hz/Pa:
  Numeric mpsf;
  // The line intensity in m^2/Hz:
  Numeric mi0;
  // Reference temperature for I0 in K:
  Numeric mti0;
  // Lower state energy in cm^-1:
  Numeric melow;
  // Air broadened width in Hz/Pa:
  Numeric magam;
  // Self broadened width in Hz/Pa:
  Numeric msgam;
  // AGAM temperature exponent (dimensionless):
  Numeric mnair;
  // SGAM temperature exponent (dimensionless):
  Numeric mnself;
  // Reference temperature for AGAM and SGAM in K:
  Numeric mtgam;
  // Array to hold auxiliary parameters:
  ArrayOfNumeric maux;
  //
  // Fields for the spectroscopic parameters accuracies
  //
  // Accuracy for line center frequency in Hz:
  Numeric mdf;
  // Accuracy for line intensity in %:
  Numeric mdi0;
  // Accuracy for air broadened width in %:
  Numeric mdagam;
  // Accuracy for self broadened width in %:
  Numeric mdsgam;
  // Accuracy for AGAM temperature exponent in %:
  Numeric mdnair;
  //  Accuracy for SGAM temperature exponent in %:
  Numeric mdnself;
 //  Accuracy for pressure shift in %:
  Numeric mdpsf;
  
  
  // Einstein A-coefficient in 1/s:
  Numeric ma;
  // Upper state stat. weight:
  Numeric mgupper;
  // Lower state stat. weight:
  Numeric mglower;
  
  // Broadening parameter self in Hz/Pa:
  Numeric mgamma_self;
  // Broadening parameter N2 in Hz/Pa:
  Numeric mgamma_n2;
  // Broadening parameter O2 in Hz/Pa:
  Numeric mgamma_o2;
  // Broadening parameter H2O in Hz/Pa:
  Numeric mgamma_h2o;
  // Broadening parameter CO2 in Hz/Pa:
  Numeric mgamma_co2;
  // Broadening parameter H2 in Hz/Pa:
  Numeric mgamma_h2;
  // Broadening parameter He in Hz/Pa:
  Numeric mgamma_he;
 
  // GAM temp. exponent self:
  Numeric mn_self;
  // GAM temp. exponent N2:
  Numeric mn_n2;
  // GAM temp. exponent O2:
  Numeric mn_o2;
  // GAM temp. exponent H2O:
  Numeric mn_h2o;
  // GAM temp. exponent CO2:
  Numeric mn_co2;
  // GAM temp. exponent H2:
  Numeric mn_h2;
  // GAM temp. exponent He:
  Numeric mn_he;
  
  // F Pressure shift N2 in Hz/Pa:
  Numeric mdelta_n2;
  // F Pressure shift O2 in Hz/Pa:
  Numeric mdelta_o2;
  // F Pressure shift H2O in Hz/Pa:
  Numeric mdelta_h2o;
  // F Pressure shift CO2 in Hz/Pa:
  Numeric mdelta_co2;
  // F Pressure shift H2 in Hz/Pa:
  Numeric mdelta_h2;
  // F Pressure shift He in Hz/Pa:
  Numeric mdelta_he;
  
};
 
// is needed to map jpl tags/arts identifier to the species/isotope data within arts
class SpecIsoMap{
public:
  SpecIsoMap():mspeciesindex(0), misotopeindex(0){}
  SpecIsoMap(const Index& speciesindex,
                const Index& isotopeindex)
    : mspeciesindex(speciesindex),
      misotopeindex(isotopeindex) 
  {}
 
  // Return the index to the species 
  const Index& Speciesindex() const { return mspeciesindex; }
  // Return the index to the isotope
  const Index& Isotopeindex() const { return misotopeindex; }
 
private:
  Index mspeciesindex;
  Index misotopeindex;
};
 
 
 
typedef Array<LineRecord> ArrayOfLineRecord;
 
typedef Array< Array<LineRecord> > ArrayOfArrayOfLineRecord;
 
 
 
ostream& operator<< (ostream& os, const LineRecord& lr);
 
 
 
void define_species_map();
 
 
// Doc header in absorption.cc
void write_lines_to_stream(ostream& os,
                           const ArrayOfLineRecord& abs_lines);
 
 
void xsec_species( MatrixView               xsec,
                   ConstVectorView          f_mono,
                   ConstVectorView          abs_p,
                   ConstVectorView          abs_t,           
                   ConstVectorView          h2o_abs,           
                   ConstVectorView          vmr,
                   const ArrayOfLineRecord& abs_lines,
                   const Index              ind_ls,
                   const Index              ind_lsn,
                   const Numeric            cutoff,
                   const Verbosity&         verbosity );
 
 
// A helper function for energy conversion:
Numeric wavenumber_to_joule(Numeric e);
 
 
//======================================================================
//             Functions related to refraction
//======================================================================
 
void refr_index_BoudourisDryAir (
                    Vector&     refr_index,
              ConstVectorView   abs_p,
              ConstVectorView   abs_t );
 
void refr_index_Boudouris (
                    Vector&     refr_index,
              ConstVectorView   abs_p,
              ConstVectorView   abs_t,
              ConstVectorView   h2o_abs );
 
//======================================================================
//             Functions to convert the accuracy index
//======================================================================
 
// ********* for HITRAN database *************
// convert index for the frequency accuracy.
void convHitranIERF(     
                    Numeric&     mdf,
              const Index&       df 
                    );
 
// convert to percents index for intensity and halfwidth accuracy.
 
void convHitranIERSH(     
                    Numeric&     mdh,
              const Index&       dh 
                    );
 
// ********* for MYTRAN database *************
// convert index for the halfwidth accuracy.
void convMytranIER(     
                    Numeric&     mdh,
              const Index  &      dh 
                    );
 
#endif // absorption_h