absorption.h

Go to the documentation of this file.

/* Copyright (C) 2000-2012
   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 <cmath>
#include "matpackI.h"
#include "array.h"
#include "mystring.h"
#include "make_array.h"
#include "messages.h"
#include "abs_species_tags.h"
#include "linerecord.h"
#include "linemixingrecord.h"


typedef void (*lsf_type)(Vector&,
                         Vector&,
                         Vector&,
                         const Numeric,
                         const Numeric,
                         const Numeric,
                         ConstVectorView);

class LineshapeRecord{
public:

  LineshapeRecord() : mname(),
                      mdescription(),
                      mphase(),
                      mfunction()
  { /* Nothing to do here. */ }

  LineshapeRecord(const String& name,
                  const String& description,
                  lsf_type      function,
                  const bool    phase)
    : mname(name),
      mdescription(description),
      mphase(phase),
      mfunction(function)
  { /* Nothing to do here. */ }
  const String&  Name()        const { return mname;        }   
  const String&  Description() const { return mdescription; }
  lsf_type Function() const { return mfunction; }
  bool Phase() const { return mphase; }
private:        
  String  mname;        
  String  mdescription; 
  bool    mphase;       
  lsf_type mfunction;   

};

typedef void (*lsnf_type)(Vector&,
                          const Numeric,
                          ConstVectorView,
                          const Numeric);

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 IsotopologueRecord{
public:

  IsotopologueRecord() : mname(),
                    mabundance(0.),
                    mmass(0.),
                    mmytrantag(-1),
                    mhitrantag(-1),
                    mjpltags(),
                    mqcoeff(),
                    mqcoefftype(-1),
                    mqcoeffgrid(),
                    mqcoeffinterporder(-1)
  { /* Nothing left to do here. */ }

  IsotopologueRecord(const IsotopologueRecord& x) :
    mname(x.mname),
    mabundance(x.mabundance),
    mmass(x.mmass),
    mmytrantag(x.mmytrantag),
    mhitrantag(x.mhitrantag),
    mjpltags(x.mjpltags),
    mqcoeff(),
    mqcoefftype(),
    mqcoeffgrid(),
    mqcoeffinterporder()
  { /* Nothing left to do here. */ }

  IsotopologueRecord(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(),
    mqcoefftype(PF_NOTHING),
    mqcoeffgrid(),
    mqcoeffinterporder()
  {
    // 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;      }


  bool isContinuum() const { return mname.length() && !isdigit(mname[0]); }

  void SetPartitionFctCoeff( const ArrayOfNumeric& qcoeff, const Index& qcoefftype )
  {
    mqcoeff = qcoeff;
    mqcoefftype = qcoefftype;
  }


  Numeric CalculatePartitionFctRatio( Numeric reference_temperature,
                                      Numeric actual_temperature ) const
  {
      
      Numeric qcoeff_at_t_ref, qtemp;
      
      switch(mqcoefftype)
        {
          case PF_FROMCOEFF:
            qcoeff_at_t_ref =
            CalculatePartitionFctAtTempFromCoeff( reference_temperature );
            qtemp =
            CalculatePartitionFctAtTempFromCoeff( actual_temperature    );
            break;
          case PF_FROMTEMP:
              qcoeff_at_t_ref =
              CalculatePartitionFctAtTempFromData( actual_temperature    );
              qtemp =
              CalculatePartitionFctAtTempFromData( actual_temperature    );
              break;
          default:
              throw runtime_error("The partition functions are incorrect.\n");
              break;
        }
/*        cout << "ref_t: " << reference_temperature << ", act_t:" <<
          actual_temperature << "\n";
        cout << "ref_q: " << qcoeff_at_t_ref << ", act_q:" <<
          qtemp << "\n";
*/
        if ( qtemp > 0. ) 
            return qcoeff_at_t_ref / qtemp;
        else
          {
            ostringstream os;
            os << "Partition function of "
               << "Isotopologue " << mname
//               << " is unknown.";
               << " at T=" << actual_temperature << "K is zero or negative.";
            throw runtime_error(os.str());
          }
  }
  
  
  enum {
      PF_FROMCOEFF,   // Partition function will be from coefficients
      PF_FROMTEMP,    // Partition function will be from temperature field
      PF_NOTHING      // This will be the designated starter value
  };

private:

  // calculate the partition fct at a certain temperature
  // this is only the prototyping
  Numeric CalculatePartitionFctAtTempFromCoeff( Numeric temperature ) const;
  Numeric CalculatePartitionFctAtTempFromData( Numeric temperature ) const;

  String mname;
  Numeric mabundance;
  Numeric mmass;
  Index mmytrantag;
  Index mhitrantag;
  ArrayOfIndex mjpltags;
  Vector mqcoeff;
  Index mqcoefftype;
  Vector mqcoeffgrid;
  Index mqcoeffinterporder;
};


class SpeciesRecord{
public:

  SpeciesRecord() : mname(),
                    mdegfr(-1),
                    misotopologue() { /* Nothing to do here */ }
  
  SpeciesRecord(const char name[],
                const Index degfr,
                const MakeArray<IsotopologueRecord>& isotopologue)
    : mname(name),
      mdegfr(degfr),
      misotopologue(isotopologue)
  {

    // Thanks to Matpack, initialization of misotopologue with isotopologue
    // should now work correctly.  

#ifndef NDEBUG
      {
        /* Check that the isotopologues are correctly sorted. */
        for ( Index i=0; i<misotopologue.nelem()-1; ++i )
          {
              assert(isnan(misotopologue[i].Abundance()) || isnan(misotopologue[i+1].Abundance())
                     || misotopologue[i].Abundance() >= misotopologue[i+1].Abundance());
          }

        /* Check that the Mytran tags are correctly sorted. */
        for ( Index i=0; i<misotopologue.nelem()-1; ++i )
          {
            if ( (0<misotopologue[i].MytranTag()) && (0<misotopologue[i+1].MytranTag()) )
              {
                assert( misotopologue[i].MytranTag() < misotopologue[i+1].MytranTag() );
            
                // Also check that the tags have the same base number:
                assert( misotopologue[i].MytranTag()/10 == misotopologue[i].MytranTag()/10 );
              }
          }

        /* Check that the Hitran tags are correctly sorted. */
        for ( Index i=0; i<misotopologue.nelem()-1; ++i )
          {
            if ( (0<misotopologue[i].HitranTag()) && (0<misotopologue[i+1].HitranTag()) )
              {
//                assert( misotopologue[i].HitranTag() < misotopologue[i+1].HitranTag() );
            
                // Also check that the tags have the same base number:
                assert( misotopologue[i].HitranTag()/10 == misotopologue[i+1].HitranTag()/10 );
              }
          }
      }
#endif // #ifndef NDEBUG
  }

  const String&               Name()     const { return mname;     }   
  Index                         Degfr()    const { return mdegfr;    }
  const Array<IsotopologueRecord>& Isotopologue()  const { return misotopologue;  }
  Array<IsotopologueRecord>&       Isotopologue()        { return misotopologue;  }
  
private:
  String mname;
  Index mdegfr;
  Array<IsotopologueRecord> misotopologue;
};


class SpeciesAuxData
{
public:
    SpeciesAuxData() : mparams() { }

    void initParams(Index nparams);

    Numeric getParam(Index species, Index isotopologue, Index col) const
    {
        return mparams[species](isotopologue, col);
    }

    void setParam(Index species, Index isotopologue, Index col, Numeric v)
    {
        mparams[species](isotopologue, col) = v;
    }

    const ArrayOfMatrix& getParams() const { return mparams; }

    bool ReadFromStream(String& artsid, istream& is, Index nparams,
                        const Verbosity& verbosity);

private:
    ArrayOfMatrix mparams;
};


void checkIsotopologueRatios(const ArrayOfArrayOfSpeciesTag& abs_species,
                             const SpeciesAuxData& sad);

void fillSpeciesAuxDataWithIsotopologueRatiosFromSpeciesData(SpeciesAuxData& sad);


// is needed to map jpl tags/arts identifier to the species/isotopologue data within arts
class SpecIsoMap{
public:
  SpecIsoMap():mspeciesindex(0), misotopologueindex(0){}
  SpecIsoMap(const Index& speciesindex,
                const Index& isotopologueindex)
    : mspeciesindex(speciesindex),
      misotopologueindex(isotopologueindex) 
  {}

  // Return the index to the species 
  const Index& Speciesindex() const { return mspeciesindex; }
  // Return the index to the isotopologue
  const Index& Isotopologueindex() const { return misotopologueindex; }

private:
  Index mspeciesindex;
  Index misotopologueindex;
};



ostream& operator<< (ostream& os, const SpeciesRecord& sr);

ostream& operator<< (ostream& os, const SpeciesAuxData& sad);



void define_species_map();


void xsec_species(MatrixView               xsec_attenuation,
                  MatrixView               xsec_phase,
                  ConstVectorView          f_grid,
                  ConstVectorView          abs_p,
                  ConstVectorView          abs_t,
                  ConstMatrixView          all_vmrs,
                  const ArrayOfArrayOfSpeciesTag& abs_species,
                  const Index              this_species,
                  const ArrayOfLineRecord& abs_lines,
                  const Index              ind_ls,
                  const Index              ind_lsn,
                  const Numeric            cutoff,
                  const SpeciesAuxData&    isotopologue_ratios,
                  const Verbosity&         verbosity );


void xsec_species_line_mixing_wrapper(  MatrixView               xsec_attenuation,
                                        MatrixView               xsec_phase,
                                        const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
                                        const ArrayOfArrayOfIndex& line_mixing_data_lut,
                                        ConstVectorView          f_grid,
                                        ConstVectorView          abs_p,
                                        ConstVectorView          abs_t,
                                        ConstMatrixView          all_vmrs,
                                        const ArrayOfArrayOfSpeciesTag& abs_species,
                                        const Index              this_species,
                                        const ArrayOfLineRecord& abs_lines,
                                        const Index              ind_ls,
                                        const Index              ind_lsn,
                                        const Numeric            cutoff,
                                        const SpeciesAuxData&    isotopologue_ratios,
                                        const Verbosity&         verbosity );


void xsec_species_line_mixing_2nd_order(    MatrixView               xsec_attenuation,
                                            MatrixView               xsec_phase,
                                            const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
                                            const ArrayOfArrayOfIndex& line_mixing_data_lut,
                                            ConstVectorView          f_grid,
                                            ConstVectorView          abs_p,
                                            ConstVectorView          abs_t,
                                            ConstMatrixView          all_vmrs,
                                            const ArrayOfArrayOfSpeciesTag& abs_species,
                                            const Index              this_species,
                                            const ArrayOfLineRecord& abs_lines,
                                            const Index              ind_ls,
                                            const Index              ind_lsn,
                                            const Numeric            cutoff,
                                            const SpeciesAuxData&    isotopologue_ratios,
                                            const Verbosity&         verbosity );


// A helper function for energy conversion:
Numeric wavenumber_to_joule(Numeric e);


//======================================================================
//             Functions related to species
//======================================================================

Index species_index_from_species_name( String name );

String species_name_from_species_index( const Index spec_ind );


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


// Functions to set abs_n2 and abs_h2o:

void abs_n2Set(Vector&            abs_n2,
               const ArrayOfArrayOfSpeciesTag& abs_species,
               const Matrix&    abs_vmrs,
               const Verbosity&);

void abs_h2oSet(Vector&          abs_h2o,
                const ArrayOfArrayOfSpeciesTag& abs_species,
                const Matrix&    abs_vmrs,
                const Verbosity&);

#endif // absorption_h