cia.cc

Go to the documentation of this file.

/* Copyright (C) 2012 Stefan Buehler <sbuehler@ltu.se>
 
 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. */

#include <cmath>
#include "cia.h"
#include "interpolation_poly.h"
#include "absorption.h"
#include "abs_species_tags.h"
#include "file.h"

extern const Numeric SPEED_OF_LIGHT;

void cia_interpolation(VectorView result,
                       ConstVectorView f_grid,
                       const Numeric& temperature,
                       const GriddedField2& cia_data,
                       const Numeric& T_extrapolfac,
                       const Index& robust,
                       const Verbosity& verbosity)
{
    CREATE_OUTS;

    const Index nf = f_grid.nelem();
    
    // Assert that result vector has right size:
    assert(result.nelem()==nf);
    
    // Get data grids:
    ConstVectorView data_f_grid = cia_data.get_numeric_grid(0);
    ConstVectorView data_T_grid = cia_data.get_numeric_grid(1);

    if (out3.sufficient_priority())
    {
        // Some detailed information to the most verbose output stream:
        ostringstream os;
        os   << "    f_grid:      " << f_grid[0] << " - "
        << f_grid[nf-1] << " Hz\n"
        << "    data_f_grid: " << data_f_grid[0] << " - "
        << data_f_grid[data_f_grid.nelem()-1] << " Hz\n"
        << "    temperature: " << temperature << " K\n"
        << "    data_T_grid: " << data_T_grid[0] << " - "
        << data_T_grid[data_T_grid.nelem()-1] << " K\n";
        out3 << os.str();
    }
    
    // Initialize result to zero (important for those frequencies outside the data grid).
    result = 0;
    
    // We want to return result zero for all f_grid points that are outside the
    // data_f_grid, because some CIA datasets are defined only where the absorption
    // is not zero. So, we have to find out which part of f_grid is inside
    // data_f_grid.
    Index i_fstart, i_fstop;
    
    for (i_fstart=0; i_fstart<nf; ++i_fstart)
        if (f_grid[i_fstart] >= data_f_grid[0]) break;
    
    // Return directly if all frequencies are below data_f_grid:
    if (i_fstart==nf) return;
    
    for (i_fstop=nf-1; i_fstop>=0; --i_fstop)
        if (f_grid[i_fstop] <= data_f_grid[data_f_grid.nelem()-1]) break;
    
    // Return directly if all frequencies are above data_f_grid:
    if (i_fstop==-1) return;
    
    // Extent for active frequency vector:
    const Index f_extent = i_fstop-i_fstart+1;

    if (out3.sufficient_priority())
    {
        ostringstream os;
        os << "    " << f_extent << " frequency extraction points starting at "
           << "frequency index " << i_fstart << ".\n";
        out3 << os.str();
    }
    
    // If f_extent is less than one, then the entire data_f_grid is between two
    // grid points of f_grid. (So that we do not have any f_grid points inside
    // data_f_grid.) Return also in this case.
    if (f_extent < 1) return;
    
    
    // This is the part of f_grid for which we have to do the interpolation.
    ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
    
    // We have to create a matching view on the result vector:
    VectorView result_active = result[Range(i_fstart, f_extent)];
    
    
    // Decide on interpolation orders:
    const Index f_order = 3;
    
    // The frequency grid has to have enough points for this interpolation
    // order, otherwise throw a runtime error.
    if ( data_f_grid.nelem() < f_order+1 )
      {
        ostringstream os;
        os << "Not enough frequency grid points in CIA data.\n"
           << "You have only " << data_f_grid.nelem() << " grid points.\n"
           << "But need at least " << f_order+1 << ".";
        throw runtime_error(os.str());
      }


    // For T we have to be adaptive, since sometimes there is only one T in
    // the data
    Index T_order;
    switch (data_T_grid.nelem()) {
        case 1:
            T_order = 0;
            break;
        case 2:
            T_order = 1;
            break;
        case 3:
            T_order = 2;
            break;
        default:
            T_order = 3;
            break;
    }
    
    // Check if frequency is inside the range covered by the data:
    chk_interpolation_grids("Frequency interpolation for CIA continuum",
                                data_f_grid,
                                f_grid_active,
                                f_order);
    

    // Check if temperature is inside the range covered by the data:
    if (T_order > 0) {
        try {
            chk_interpolation_grids("Temperature interpolation for CIA continuum",
                                    data_T_grid,
                                    temperature,
                                    T_order,
                                    T_extrapolfac);
        } catch (runtime_error e) {
            //            cout << "Gotcha!\n";
            if (robust) {
                // Just return NANs, but continue.
                result_active = NAN;
                return;
            } else {
                // Re-throw the exception.
                throw runtime_error(e.what());
            }
        }
    }
    
    // Find frequency grid positions:
    ArrayOfGridPosPoly f_gp(f_grid_active.nelem()), T_gp(1);
    gridpos_poly(f_gp, data_f_grid, f_grid_active, f_order);
 
    // Do the rest of the interpolation.
    if (T_order == 0)
      {
        // No temperature interpolation in this case, just a frequency interpolation.
        
        Matrix itw(f_gp.nelem(),f_order+1);
        interpweights(itw,f_gp);
        interp(result_active, itw, cia_data.data(joker,0), f_gp);
      }
    else
      {
        // Temperature and frequency interpolation.
        
        // Find temperature grid position:
        gridpos_poly(T_gp, data_T_grid, temperature, T_order, T_extrapolfac);
    
        // Calculate combined interpolation weights:
        Tensor3 itw(f_gp.nelem(),T_gp.nelem(),(f_order+1)*(T_order+1));
        interpweights(itw, f_gp, T_gp);
        
        // Make a matrix view of the results vector:
        MatrixView result_matrix(result_active);
      
//        cout << "result_matrix r/c: " << result_matrix.nrows() << " / "
//             << result_matrix.ncols() << endl;
//        cout << "result_matrix before: " << result_matrix << endl;
        
        //        cout << "cia_data: " << cia_data.data << endl;
        
        // Actual interpolation:
        interp(result_matrix, itw, cia_data.data, f_gp, T_gp);

//        cout << "result_matrix after: " << result_matrix << endl;
      }
    
//    cout << "result_active before: " << result_active << endl;

    // Set negative values to zero. (These could happen due to overshooting
    // of the higher order interpolation.)
    for (Index i=0; i<result_active.nelem(); ++i)
        if (result_active[i]<0)
            result_active[i] = 0;

//    cout << "result_active after: " << result_active << endl;
    
}


Index cia_get_index(const ArrayOfCIARecord& cia_data,
                    const Index sp1, const Index sp2)
{
    for (Index i = 0; i < cia_data.nelem(); i++)
        if ((cia_data[i].Species(0) == sp1 && cia_data[i].Species(1) ==  sp2)
            || (cia_data[i].Species(0) == sp2 && cia_data[i].Species(1) ==  sp1))
            return i;

    return -1;
}


// Documentation in header file.
void CIARecord::Extract(VectorView      result,
                        ConstVectorView f_grid,
                        const Numeric&  temperature,
                        const Index&    dataset,
                        const Numeric&  T_extrapolfac,
                        const Index&    robust,
                        const Verbosity& verbosity) const
{
    // If there is more than one dataset available for this species pair,
    // we have to decide on which one to use. The rest is done by helper function
    // cia_interpolate.
    
    // Make sure dataset index exists
    if ( dataset >= mdata.nelem() )
      {
        ostringstream os;
        os << "There are only " << mdata.nelem() << " datasets in this CIA file.\n"
           << "But you are trying to use dataset " << dataset << ". (Zero-based indexing.)";
        throw runtime_error(os.str());
      }


    // Get a handle on this dataset:
    const GriddedField2& this_cia = mdata[dataset];

    cia_interpolation(result,
                      f_grid,
                      temperature,
                      this_cia,
                      T_extrapolfac,
                      robust,
                      verbosity);
}


// Documentation in header file.
String CIARecord::MoleculeName(const Index i) const
{
    // Assert that i is 0 or 1:
    assert(i>=0);
    assert(i<=1);

    // The function species_name_from_species_index internally does an assertion
    // that the species with this index really exists.
    return species_name_from_species_index( mspecies[i] );
}

// Documentation in header file.
void CIARecord::SetMoleculeName(const Index i,
                                const String& name)
{
    // Assert that i is 0 or 1:
    assert(i>=0);
    assert(i<=1);
    
    // Find out the species index for name:
    Index spec_ind = species_index_from_species_name(name);

    // Function species_index_from_species_name returns -1 if the species does
    // not exist. Check this:
    if ( spec_ind < 0 )
      {
        ostringstream os;
        os << "Species does not exist in ARTS: " << name;
        throw runtime_error(os.str());
      }

    // Assign species:
    mspecies[i] = spec_ind;
}



void CIARecord::ReadFromCIA(const String& filename, const Verbosity& verbosity)
{
    CREATE_OUT2;

    ifstream is;

    out2 << "  Reading file: " << filename << "\n";
    open_input_file(is, filename);

    // Number of points for spectral range in current dataset
    Index npoints = -1;

    // Min/max wave numbers for this dataset
    Numeric wave_min = -1.;
    Numeric wave_max = -1.;

    // Current dataset index
    Index ndataset = -1;

    // Current set in dataset
    Index nset = -1;

    // Frequency, temp and cia values for current dataset
    Vector freq;
    ArrayOfNumeric temp;
    ArrayOfVector cia;

    // Keep track of current line in file
    Index nline = 0;

    mdata.resize(0);
    istringstream istr;

    while (is)
    {
        String line;

        // Extract needed information from dataset header line
        nline++;
        getline(is,line);
        if (is.eof()) continue;

        if (line.nelem() < 100)
        {
            ostringstream os;
            os << "Error in line " << nline
            << " reading CIA catalog file " << filename << endl
            << "Header line unexpectedly short: " << endl << line;

            throw runtime_error(os.str());
        }

        if (is.bad())
        {
            ostringstream os;
            os << "Error in line " << nline
            << " reading CIA catalog file " << filename << endl;

            throw runtime_error(os.str());
        }

        line.erase(0, 20);

        // Data for current set
        Index set_npoints;
        Numeric set_temp = NAN;
        Numeric set_wave_min = NAN;
        Numeric set_wave_max = NAN;

        istr.str(line);
        istr.clear();
        istr >> set_wave_min >> set_wave_max >> set_npoints >> set_temp;

        if (!istr || isnan(set_temp) || isnan(set_wave_min) || isnan(set_wave_max))
        {
            ostringstream os;
            os << "Error in line " << nline
            << " reading CIA catalog file " << filename << endl;

            throw runtime_error(os.str());
        }

        // If the min/max wave numbers of this set are different from the
        // previous one, a new dataset starts
        if (npoints == -1 || wave_min != set_wave_min || wave_max != set_wave_max)
        {
            if (ndataset != -1)
                AppendDataset(freq, temp, cia);

            npoints = set_npoints;
            ndataset++;
            wave_min = set_wave_min;
            wave_max = set_wave_max;
            nset = 0;
            freq.resize(set_npoints);
            temp.resize(0);
            cia.resize(0);
        }
        if (npoints != set_npoints)
        {
            ostringstream os;
            os << "Error in line " << nline
            << " reading CIA catalog file " << filename << endl
            << "Inconsistent number of data points. Expected " << npoints
            << ", got " << set_npoints;

            throw runtime_error(os.str());
        }

        temp.push_back(set_temp);
        cia.push_back(Vector(npoints));

        // Read dataset
        for (Index i = 0; i < npoints; i++)
        {
            Numeric w, c;

            nline++;
            getline(is,line);

            w = c = NAN;
            istr.str(line);
            istr.clear();
            istr >> w >> c;

            if (isnan(w) || isnan(c) || is.bad() || istr.bad())
            {
                ostringstream os;
                os << "Error in line " << nline
                << " reading CIA catalog file " << filename << ":" << endl
                << line;

                throw runtime_error(os.str());
            }

            // Convert wavenumbers to Herz:
            freq[i] = 100. * w * SPEED_OF_LIGHT;
          
            // Convert binary absorption cross-sections from
            // cm^5 molec^(-2) to m^5 molec^(-2):
            cia[nset][i] = c / 1e10;
        }

        nset++;
    }

    if (is.bad())
    {
        ostringstream os;
        os << "Error in line " << nline
        << " reading CIA catalog file " << filename << endl;

        throw runtime_error(os.str());
    }

    AppendDataset(freq, temp, cia);

//    // For debugging
//    for(Index i = 0; i < mdata.nelem(); i++)
//    {
//        cout << i << " ";
//        cout << mdata[i].get_numeric_grid(0).nelem() << " ";
//        cout << mdata[i].get_numeric_grid(1).nelem() << endl;
//    }
}


void CIARecord::AppendDataset(const Vector& freq,
                              const ArrayOfNumeric& temp,
                              const ArrayOfVector& cia)
{
    GriddedField2 dataset;
    dataset.resize(freq.nelem(), temp.nelem());
    dataset.set_grid(0, freq);
    dataset.set_grid_name(0, "Frequency");

    Vector temp_t;
    temp_t = temp;
    dataset.set_grid(1, temp_t);
    dataset.set_grid_name(1, "Temperature");

    for (Index t = 0; t < temp.nelem(); t++)
        dataset.data(joker, t) = cia[t];
    mdata.push_back(dataset);
}



ostream& operator<<(ostream& os, const CIARecord& /* cr */)
{
    os << "CIARecord output operator not yet implemented." << endl;
    return os;
}