absorption.h

Go to the documentation of this file.

/* Copyright (C) 2000, 2001 Stefan Buehler <sbuehler@uni-bremen.de>
                            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. */

// This file contains 

#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(){};

  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(){};

  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(){};

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

typedef Array<LineshapeSpec> ArrayOfLineshapeSpec;



class IsotopeRecord{
public:

  IsotopeRecord() { /* Nothing 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)
  { /* 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)
  {
    // 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 << "Negative value of partition function was obtained for "
           << mname << ".\nProbably caused by low or negative temperature.";
        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(){}  ;
  
  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>& 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        )
  {
    // 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.
    extern Array<SpeciesRecord> species_data;
    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 {
    // The species lookup data:
    extern const Array<SpeciesRecord> species_data;
    const SpeciesRecord& sr = species_data[mspecies];
    return sr.Name() + "-" + sr.Isotope()[misotope].Name();
  }

  const SpeciesRecord& SpeciesData() const {
    // The species lookup data:
    extern const Array<SpeciesRecord> species_data;
    return species_data[mspecies];
  }

  const IsotopeRecord& IsotopeData() const {
    // The species lookup data:
    extern const Array<SpeciesRecord> species_data;
    return species_data[mspecies].Isotope()[misotope];
  }

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



//------------------------------< Tag Group Stuff >------------------------------

class OneTag {
public:
  OneTag() { /* Nothing to be done here. */ }

  OneTag(String def); 

  String Name() const;
    
  Index Species() const { return mspecies; }

  Index Isotope() const { return misotope; }

  Numeric Lf() const { return mlf; }

  Numeric Uf() const { return muf; }

private:
  // Molecular species index: 
  Index mspecies;
  // Isotopic species index.
  // If this is equal to the number of isotopes (one more than
  // allowed) it means all isotopes of this species.
  Index misotope;
  // The lower line center frequency in Hz.
  // If this is <0 it means no lower limit. 
  Numeric mlf;
  // The upper line center frequency in Hz:
  // If this is <0 it means no upper limit. 
  Numeric muf;
};


ostream& operator << (ostream& os, const OneTag& ot);


typedef  Array< Array<OneTag> > TagGroups;


void get_tagindex_for_Strings( 
              ArrayOfIndex&   tags1_index, 
        const TagGroups&      tags1, 
        const ArrayOfString&  tags2_Strings );

void get_tag_group_index_for_tag_group( 
              Index&         tags1_index, 
        const TagGroups&      tags1, 
        const Array<OneTag>&  tags2 );

String get_tag_group_name( const Array<OneTag>& tg );

// Doc header in absorption.cc
void write_lines_to_stream(ostream& os,
                           const ArrayOfLineRecord& lines);


void xsec_species( MatrixView              xsec,
                   ConstVectorView         f_mono,
                   ConstVectorView         p_abs,
                   ConstVectorView         t_abs,           
                   ConstVectorView         h2o_abs,           
                   ConstVectorView            vmr,
                   const ArrayOfLineRecord& 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   p_abs,
              ConstVectorView   t_abs );

void refr_index_Boudouris (
                    Vector&     refr_index,
              ConstVectorView   p_abs,
              ConstVectorView   t_abs,
              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

---