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"

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()
    : 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 */)
    : 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 Version() const
  {
    ostringstream os;
    os << "ARTSCAT-" << mversion;
    return os.str();
  }

  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; }


  bool ReadFromHitranStream(istream& is);



  bool ReadFromHitran2004Stream(istream& is);




  bool ReadFromMytran2Stream(istream& is);


  bool ReadFromJplStream(istream& is);

  bool ReadFromArtsStream(istream& is);


private:
  // Versin number:
  static const Index mversion = 3;
  // 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;
};

// 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);


// 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

---