ARTS  2.2.66
absorption.cc
Go to the documentation of this file.
1 /* Copyright (C) 2000-2012
2  Stefan Buehler <sbuehler@ltu.se>
3  Axel von Engeln <engeln@uni-bremen.de>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
33 #include "arts.h"
34 #include <map>
35 #include <cfloat>
36 #include <algorithm>
37 #include <cmath>
38 #include <cstdlib>
39 #include "file.h"
40 #include "absorption.h"
41 #include "math_funcs.h"
42 #include "messages.h"
43 #include "logic.h"
44 #include "interpolation_poly.h"
45 
46 #include "global_data.h"
47 
48 
50 std::map<String, Index> SpeciesMap;
51 
52 
53 // member fct of isotopologuerecord, calculates the partition fct at the
54 // given temperature from the partition fct coefficients (3rd order
55 // polynomial in temperature)
57  temperature ) const
58 {
59  Numeric result = 0.;
60  Numeric exponent = 1.;
61 
63 
64 // cout << "T: " << temperature << "\n";
65  for (it=mqcoeff.begin(); it != mqcoeff.end(); ++it)
66  {
67  result += *it * exponent;
68  exponent *= temperature;
69 // cout << "it: " << it << "\n";
70 // cout << "res: " << result << ", exp: " << exponent << "\n";
71  }
72  return result;
73 }
74 
75 
76 // member fct of isotopologuerecord, calculates the partition fct at the
77 // given temperature from the partition fct coefficients (3rd order
78 // polynomial in temperature)
80 temperature ) const
81 {
82  GridPosPoly gp;
83  gridpos_poly( gp, mqcoeffgrid, temperature, mqcoeffinterporder);
84  Vector itw;
85  interpweights(itw, gp );
86  return interp( itw, mqcoeff, gp );
87 }
88 
89 
91 {
93 
94  mparams.resize(species_data.nelem());
95 
96  for (Index isp = 0; isp < species_data.nelem(); isp++)
97  {
98  mparams[isp].resize(species_data[isp].Isotopologue().nelem(), nparams);
99  mparams[isp] = NAN;
100  }
101 }
102 
103 
104 bool SpeciesAuxData::ReadFromStream(String& artsid, istream& is, Index nparams, const Verbosity& verbosity)
105 {
106  CREATE_OUT3;
107 
108  // Global species lookup data:
110 
111  // We need a species index sorted by Arts identifier. Keep this in a
112  // static variable, so that we have to do this only once. The ARTS
113  // species index is ArtsMap[<Arts String>].
114  static map<String, SpecIsoMap> ArtsMap;
115 
116  // Remember if this stuff has already been initialized:
117  static bool hinit = false;
118 
119  if ( !hinit )
120  {
121 
122  out3 << " ARTS index table:\n";
123 
124  for ( Index i=0; i<species_data.nelem(); ++i )
125  {
126  const SpeciesRecord& sr = species_data[i];
127 
128 
129  for ( Index j=0; j<sr.Isotopologue().nelem(); ++j)
130  {
131 
132  SpecIsoMap indicies(i,j);
133  String buf = sr.Name()+"-"+sr.Isotopologue()[j].Name();
134 
135  ArtsMap[buf] = indicies;
136 
137  // Print the generated data structures (for debugging):
138  // The explicit conversion of Name to a c-String is
139  // necessary, because setw does not work correctly for
140  // stl Strings.
141  const Index& i1 = ArtsMap[buf].Speciesindex();
142  const Index& i2 = ArtsMap[buf].Isotopologueindex();
143 
144  out3 << " Arts Identifier = " << buf << " Species = "
145  << setw(10) << setiosflags(ios::left)
146  << species_data[i1].Name().c_str()
147  << "iso = "
148  << species_data[i1].Isotopologue()[i2].Name().c_str()
149  << "\n";
150 
151  }
152  }
153  hinit = true;
154  }
155 
156 
157  // This always contains the rest of the line to parse. At the
158  // beginning the entire line. Line gets shorter and shorter as we
159  // continue to extract stuff from the beginning.
160  String line;
161 
162  // Look for more comments?
163  bool comment = true;
164 
165  while (comment)
166  {
167  // Return true if eof is reached:
168  if (is.eof()) return true;
169 
170  // Throw runtime_error if stream is bad:
171  if (!is) throw runtime_error ("Stream bad.");
172 
173  // Read line from file into linebuffer:
174  getline(is,line);
175 
176  // It is possible that we were exactly at the end of the file before
177  // calling getline. In that case the previous eof() was still false
178  // because eof() evaluates only to true if one tries to read after the
179  // end of the file. The following check catches this.
180  if (line.nelem() == 0 && is.eof()) return true;
181 
182  // @ as first character marks catalogue entry
183  char c;
184  extract(c,line,1);
185 
186  // check for empty line
187  if (c == '@')
188  {
189  comment = false;
190  }
191  }
192 
193 
194  // read the arts identifier String
195  istringstream icecream(line);
196 
197  icecream >> artsid;
198 
199  if (artsid.length() != 0)
200  {
201  Index mspecies;
202  Index misotopologue;
203 
204  // ok, now for the cool index map:
205  // is this arts identifier valid?
206  const map<String, SpecIsoMap>::const_iterator i = ArtsMap.find(artsid);
207  if ( i == ArtsMap.end() )
208  {
209  ostringstream os;
210  os << "ARTS Tag: " << artsid << " is unknown.";
211  throw runtime_error(os.str());
212  }
213 
214  SpecIsoMap id = i->second;
215 
216 
217  // Set mspecies:
218  mspecies = id.Speciesindex();
219 
220  // Set misotopologue:
221  misotopologue = id.Isotopologueindex();
222 
223  Matrix& params = mparams[mspecies];
224  // Extract accuracies:
225  try
226  {
227  Numeric p;
228  for (Index ip = 0; ip < nparams; ip++)
229  {
230  icecream >> double_imanip() >> p;
231  params(misotopologue, ip) = p;
232  }
233  }
234  catch (runtime_error)
235  {
236  throw runtime_error("Error reading SpeciesAuxData.");
237  }
238  }
239 
240  // That's it!
241  return false;
242 }
243 
244 
246  const SpeciesAuxData& isoratios)
247 {
249 
250  // Check total number of species:
251  if (species_data.nelem() != isoratios.getParams().nelem())
252  {
253  ostringstream os;
254  os << "Number of species in SpeciesAuxData (" << isoratios.getParams().nelem()
255  << "does not fit builtin species data (" << species_data.nelem() << ").";
256  throw runtime_error(os.str());
257  }
258 
259  // For the selected species, we check all isotopes by looping over the
260  // species data. (Trying to check only the isotopes actually used gets
261  // quite complicated, actually, so we do the simple thing here.)
262 
263  // Loop over the absorption species:
264  for (Index i=0; i<abs_species.nelem(); i++)
265  {
266  // sp is the index of this species in the internal lookup table
267  const Index sp = abs_species[i][0].Species();
268 
269  // Get handle on species data for this species:
270  const SpeciesRecord& this_sd = species_data[sp];
271 
272  // Check number of isotopologues:
273  if (this_sd.Isotopologue().nelem() != isoratios.getParams()[sp].nrows())
274  {
275  ostringstream os;
276  os << "Incorrect number of isotopologues in isotopologue data.\n"
277  << "Species: " << this_sd.Name() << ".\n"
278  << "Number of isotopes in SpeciesAuxData ("
279  << isoratios.getParams()[i].nrows() << ") "
280  << "does not fit builtin species data (" << this_sd.Isotopologue().nelem() << ").";
281  throw runtime_error(os.str());
282  }
283 
284  for (Index iso = 0; iso < this_sd.Isotopologue().nelem(); ++iso)
285  {
286  // For "real" species (not representing continau) the isotopologue
287  // ratio must not be NAN or below zero.
288  if (!this_sd.Isotopologue()[iso].isContinuum()) {
289  if (isnan(isoratios.getParam(sp, iso, 0)) ||
290  isoratios.getParam(sp, iso, 0) < 0.) {
291 
292  ostringstream os;
293  os << "Invalid isotopologue ratio.\n"
294  << "Species: " << this_sd.Name() << "-"
295  << this_sd.Isotopologue()[iso].Name() << "\n"
296  << "Ratio: " << isoratios.getParam(sp, iso, 0);
297  throw runtime_error(os.str());
298  }
299  }
300  }
301  }
302 }
303 
304 
306 {
308 
309  sad.initParams(1);
310 
311  for (Index isp = 0; isp < species_data.nelem(); isp++)
312  for (Index iiso = 0; iiso < species_data[isp].Isotopologue().nelem(); iiso++)
313  {
314  sad.setParam(isp, iiso, 0, species_data[isp].Isotopologue()[iiso].Abundance());
315  }
316 }
317 
318 
323 {
325 
326  for ( Index i=0 ; i<species_data.nelem() ; ++i)
327  {
328  SpeciesMap[species_data[i].Name()] = i;
329  }
330 }
331 
332 
333 ostream& operator<< (ostream& os, const SpeciesRecord& sr)
334 {
335  for ( Index j=0; j<sr.Isotopologue().nelem(); ++j)
336  {
337  os << sr.Name() << "-" << sr.Isotopologue()[j].Name() << "\n";
338  }
339  return os;
340 }
341 
342 
343 ostream& operator<< (ostream& os, const SpeciesAuxData& sad)
344 {
345  os << sad.getParams();
346  return os;
347 }
348 
349 
350 
372 void find_broad_spec_locations(ArrayOfIndex& broad_spec_locations,
373  const ArrayOfArrayOfSpeciesTag& abs_species,
374  const Index this_species)
375 {
376  // Make sure this_species points to somewhere inside abs_species:
377  assert(this_species>=0);
378  assert(this_species<abs_species.nelem());
379 
380  // Number of broadening species:
381  const Index nbs = LineRecord::NBroadSpec();
382 
383  // Resize output array:
384  broad_spec_locations.resize(nbs);
385 
386 // // Set broadening species names, using the enums defined in absorption.h.
387 // // This is hardwired here and quite primitive, but should do the job.
388 // ArrayOfString broad_spec_names(nbs);
389 // broad_spec_names[LineRecord::SPEC_POS_N2] = "N2";
390 // broad_spec_names[LineRecord::SPEC_POS_O2] = "O2";
391 // broad_spec_names[LineRecord::SPEC_POS_H2O] = "H2O";
392 // broad_spec_names[LineRecord::SPEC_POS_CO2] = "CO2";
393 // broad_spec_names[LineRecord::SPEC_POS_H2] = "H2";
394 // broad_spec_names[LineRecord::SPEC_POS_He] = "He";
395 
396  // Loop over all broadening species and see if we can find them in abs_species.
397  for (Index i=0; i<nbs; ++i) {
398  // Find associated internal species index (we do the lookup by index, not by name).
399  const Index isi = LineRecord::BroadSpecSpecIndex(i);
400 
401  // First check if this broadening species is the same as this_species
402  if ( isi == abs_species[this_species][0].Species() )
403  broad_spec_locations[i] = -2;
404  else
405  {
406  // Find position of broadening species isi in abs_species. The called
407  // function returns -1 if not found, which is already the correct
408  // treatment for this case that we also want here.
409  broad_spec_locations[i] = find_first_species_tg(abs_species,isi);
410  }
411  }
412 }
413 
433  Numeric& deltaf,
434  const Numeric p,
435  const Numeric t,
436  ConstVectorView vmrs,
437  const Index this_species,
438  const ArrayOfIndex& broad_spec_locations,
439  const LineRecord& l_l,
440  const Verbosity& verbosity)
441 {
442  CREATE_OUT2;
443 
444  // Number of broadening species:
445  const Index nbs = LineRecord::NBroadSpec();
446  assert(nbs==broad_spec_locations.nelem());
447 
448  // Theta is reference temperature divided by local temperature. Used in
449  // several place by the broadening and shift formula.
450  const Numeric theta = l_l.Ti0() / t;
451 
452  // Split total pressure in self and foreign part:
453  const Numeric p_self = vmrs[this_species] * p;
454  const Numeric p_foreign = p-p_self;
455 
456  // Calculate sum of VMRs of all available foreign broadening species (we need this
457  // for normalization). The species "Self" will not be included in the sum!
458  Numeric broad_spec_vmr_sum = 0;
459 
460  // Gamma is the line width. We first initialize gamma with the self width
461  gamma = l_l.Sgam() * pow(theta, l_l.Nself()) * p_self;
462 
463  // and treat foreign width separately:
464  Numeric gamma_foreign = 0;
465 
466  // There is no self shift parameter (or rather, we do not have it), so
467  // we do not need separate treatment of self and foreign for the shift:
468  deltaf = 0;
469 
470  // Add up foreign broadening species, where available:
471  for (Index i=0; i<nbs; ++i) {
472  if ( broad_spec_locations[i] < -1 ) {
473  // -2 means that this broadening species is identical to Self.
474  // Throw runtime errors if the parameters are not identical.
475  if (l_l.Gamma_foreign(i)!=l_l.Sgam() ||
476  l_l.N_foreign(i)!=l_l.Nself())
477  {
478  ostringstream os;
479  os << "Inconsistency in LineRecord, self broadening and line "
480  << "broadening for " << LineRecord::BroadSpecName(i) << "\n"
481  << "should be identical.\n"
482  << "LineRecord:\n"
483  << l_l;
484  throw runtime_error(os.str());
485  }
486  } else if ( broad_spec_locations[i] >= 0 ) {
487 
488  // Add to VMR sum:
489  broad_spec_vmr_sum += vmrs[broad_spec_locations[i]];
490 
491  // foreign broadening:
492  gamma_foreign += l_l.Gamma_foreign(i) * pow(theta, l_l.N_foreign(i))
493  * vmrs[broad_spec_locations[i]];
494 
495  // Pressure shift:
496  // The T dependence is connected to that of the corresponding
497  // broadening parameter by:
498  // n_shift = .25 + 1.5 * n_gamma
499  deltaf += l_l.Delta_foreign(i)
500  * pow( theta , (Numeric).25 + (Numeric)1.5*l_l.N_foreign(i) )
501  * vmrs[broad_spec_locations[i]];
502  }
503  }
504 
505  // Check that sum of self and all foreign VMRs is not too far from 1:
506  if ( abs(vmrs[this_species]+broad_spec_vmr_sum-1) > 0.1
507  && out2.sufficient_priority() )
508  {
509  ostringstream os;
510  os << "Warning: The total VMR of all your defined broadening\n"
511  << "species (including \"self\") is "
512  << vmrs[this_species]+broad_spec_vmr_sum
513  << ", more than 10% " << "different from 1.\n";
514  out2 << os.str();
515  }
516 
517  // Normalize foreign gamma and deltaf with the foreign VMR sum (but only if
518  // we have any foreign broadening species):
519  if (broad_spec_vmr_sum != 0.)
520  {
521  gamma_foreign /= broad_spec_vmr_sum;
522  deltaf /= broad_spec_vmr_sum;
523  }
524  else if (p_self > 0.)
525  // If there are no foreign broadening species present, the best assumption
526  // we can make is to use gamma_self in place of gamma_foreign. for deltaf
527  // there is no equivalent solution, as we don't have a Delta_self and don't
528  // know which other Delta we should apply (in this case delta_f gets 0,
529  // which should be okayish):
530  {
531  gamma_foreign = gamma/p_self;
532  }
533  // It can happen that broad_spec_vmr_sum==0 AND p_self==0 (e.g., when p_grid
534  // exceeds the given atmosphere and zero-padding is applied). In this case,
535  // both gamma_foreign and deltaf are 0 and we leave it like that.
536 
537  // Multiply by pressure. For the width we take only the foreign pressure.
538  // This is consistent with that we have scaled with the sum of all foreign
539  // broadening VMRs. In this way we make sure that the total foreign broadening
540  // scales with the total foreign pressure.
541  gamma_foreign *= p_foreign;
542  // For the shift we simply take the total pressure, since there is no self part.
543  deltaf *= p;
544 
545  // For the width, add foreign parts:
546  gamma += gamma_foreign;
547 
548  // That's it, we're done.
549 }
550 
594 void xsec_species( MatrixView xsec_attenuation,
595  MatrixView xsec_phase,
596  ConstVectorView f_grid,
597  ConstVectorView abs_p,
598  ConstVectorView abs_t,
599  ConstMatrixView all_vmrs,
600  const ArrayOfArrayOfSpeciesTag& abs_species,
601  const Index this_species,
602  const ArrayOfLineRecord& abs_lines,
603  const Index ind_ls,
604  const Index ind_lsn,
605  const Numeric cutoff,
606  const SpeciesAuxData& isotopologue_ratios,
607  const Verbosity& verbosity )
608 {
609  // Make lineshape and species lookup data visible:
612 
613  // speed of light constant
614  extern const Numeric SPEED_OF_LIGHT;
615 
616  // Boltzmann constant
617  extern const Numeric BOLTZMAN_CONST;
618 
619  // Avogadros constant
620  extern const Numeric AVOGADROS_NUMB;
621 
622  // Planck constant
623  extern const Numeric PLANCK_CONST;
624 
625  // sqrt(ln(2))
626  // extern const Numeric SQRT_NAT_LOG_2;
627 
628  // Constant within the Doppler Broadening calculation:
629  const Numeric doppler_const = sqrt( 2.0 * BOLTZMAN_CONST *
630  AVOGADROS_NUMB) / SPEED_OF_LIGHT;
631 
632  // dimension of f_grid, abs_lines
633  const Index nf = f_grid.nelem();
634  const Index nl = abs_lines.nelem();
635 
636  // number of pressure levels:
637  const Index np = abs_p.nelem();
638 
639  // Define the vector for the line shape function and the
640  // normalization factor of the lineshape here, so that we don't need
641  // so many free store allocations. the last element is used to
642  // calculate the value at the cutoff frequency
643  Vector ls_attenuation(nf+1);
644  Vector ls_phase(nf+1);
645  Vector fac(nf+1);
646 
647  const bool cut = (cutoff != -1) ? true : false;
648 
649  const bool calc_phase = lineshape_data[ind_ls].Phase();
650 
651  // Check that the frequency grid is sorted in the case of lineshape
652  // with cutoff. Duplicate frequency values are allowed.
653  if (cut)
654  {
655  if ( ! is_sorted( f_grid ) )
656  {
657  ostringstream os;
658  os << "If you use a lineshape function with cutoff, your\n"
659  << "frequency grid *f_grid* must be sorted.\n"
660  << "(Duplicate values are allowed.)";
661  throw runtime_error(os.str());
662  }
663  }
664 
665  // Check that all temperatures are non-negative
666  bool negative = false;
667 
668  for (Index i = 0; !negative && i < abs_t.nelem (); i++)
669  {
670  if (abs_t[i] < 0.)
671  negative = true;
672  }
673 
674  if (negative)
675  {
676  ostringstream os;
677  os << "abs_t contains at least one negative temperature value.\n"
678  << "This is not allowed.";
679  throw runtime_error(os.str());
680  }
681 
682  // We need a local copy of f_grid which is 1 element longer, because
683  // we append a possible cutoff to it.
684  // The initialization of this has to be inside the line loop!
685  Vector f_local( nf + 1 );
686 
687  // Voigt generally needs a different frequency grid. If we allocate
688  // that in the outer loop, instead of in voigt, we don't have the
689  // free store allocation at each lineshape call. Calculation is
690  // still done in the voigt routine itself, this is just an auxillary
691  // parameter, passed to lineshape. For selected lineshapes (e.g.,
692  // Rosenkranz) it is used additionally to pass parameters needed in
693  // the lineshape (e.g., overlap, ...). Consequently we have to
694  // assure that aux has a dimension not less then the number of
695  // parameters passed.
696  Index ii = (nf+1 < 10) ? 10 : nf+1;
697  Vector aux(ii);
698 
699  // Check that abs_p, abs_t, and abs_vmrs have consistent
700  // dimensions. This could be a user error, so we throw a
701  // runtime_error.
702 
703  if ( abs_t.nelem() != np )
704  {
705  ostringstream os;
706  os << "Variable abs_t must have the same dimension as abs_p.\n"
707  << "abs_t.nelem() = " << abs_t.nelem() << '\n'
708  << "abs_p.nelem() = " << np;
709  throw runtime_error(os.str());
710  }
711 
712  // all_vmrs should have dimensions [nspecies, np]:
713 
714  if ( all_vmrs.ncols() != np )
715  {
716  ostringstream os;
717  os << "Number of columns of all_vmrs must match abs_p.\n"
718  << "all_vmrs.ncols() = " << all_vmrs.ncols() << '\n'
719  << "abs_p.nelem() = " << np;
720  throw runtime_error(os.str());
721  }
722 
723  const Index nspecies = abs_species.nelem();
724 
725  if ( all_vmrs.nrows() != nspecies)
726  {
727  ostringstream os;
728  os << "Number of rows of all_vmrs must match abs_species.\n"
729  << "all_vmrs.nrows() = " << all_vmrs.nrows() << '\n'
730  << "abs_species.nelem() = " << nspecies;
731  throw runtime_error(os.str());
732  }
733 
734  // With abs_h2o it is different. We do not really need this in most
735  // cases, only the Rosenkranz lineshape for oxygen uses it. There is
736  // a global (scalar) default value of -1, that we are expanding to a
737  // vector here if we find it. The Rosenkranz lineshape does a check
738  // to make sure that the value is actually set, and not the default
739  // value.
740 // Vector abs_h2o(np);
741 // if ( abs_h2o_orig.nelem() == np )
742 // {
743 // abs_h2o = abs_h2o_orig;
744 // }
745 // else
746 // {
747 // if ( ( 1 == abs_h2o_orig.nelem()) &&
748 // ( -.99 > abs_h2o_orig[0]) )
749 // {
750 // // We have found the global default value. Expand this to a
751 // // vector with the right length, by copying -1 to all
752 // // elements of abs_h2o.
753 // abs_h2o = -1;
754 // }
755 // else
756 // {
757 // ostringstream os;
758 // os << "Variable abs_h2o must have default value -1 or the\n"
759 // << "same dimension as abs_p.\n"
760 // << "abs_h2o.nelem() = " << abs_h2o.nelem() << '\n'
761 // << "abs_p.nelem() = " << np;
762 // throw runtime_error(os.str());
763 // }
764 // }
765 
766  // Check that the dimension of xsec is indeed [f_grid.nelem(),
767  // abs_p.nelem()]:
768  if ( xsec_attenuation.nrows() != nf || xsec_attenuation.ncols() != np )
769  {
770  ostringstream os;
771  os << "Variable xsec must have dimensions [f_grid.nelem(),abs_p.nelem()].\n"
772  << "[xsec_attenuation.nrows(),xsec_attenuation.ncols()] = [" << xsec_attenuation.nrows()
773  << ", " << xsec_attenuation.ncols() << "]\n"
774  << "f_grid.nelem() = " << nf << '\n'
775  << "abs_p.nelem() = " << np;
776  throw runtime_error(os.str());
777  }
778  if ( xsec_phase.nrows() != nf || xsec_phase.ncols() != np )
779  {
780  ostringstream os;
781  os << "Variable xsec must have dimensions [f_grid.nelem(),abs_p.nelem()].\n"
782  << "[xsec_phase.nrows(),xsec_phase.ncols()] = [" << xsec_phase.nrows()
783  << ", " << xsec_phase.ncols() << "]\n"
784  << "f_grid.nelem() = " << nf << '\n'
785  << "abs_p.nelem() = " << np;
786  throw runtime_error(os.str());
787  }
788 
789  // Find the location of all broadening species in abs_species. Set to -1 if
790  // not found. The length of array broad_spec_locations is the number of allowed
791  // broadening species (in ARTSCAT-4 Self, N2, O2, H2O, CO2, H2, He). The value
792  // means:
793  // -1 = not in abs_species
794  // -2 = in abs_species, but should be ignored because it is identical to Self
795  // N = species is number N in abs_species
796  ArrayOfIndex broad_spec_locations;
797  find_broad_spec_locations(broad_spec_locations,
798  abs_species,
799  this_species);
800 
801  String fail_msg;
802  bool failed = false;
803 
804  // Loop all pressures:
805  if (np)
806 #pragma omp parallel for \
807  if (!arts_omp_in_parallel() \
808  && np >= arts_omp_get_max_threads()) \
809  firstprivate(ls_attenuation, ls_phase, fac, f_local, aux)
810  for ( Index i=0; i<np; ++i )
811  {
812  if (failed) continue;
813 
814  // Store input profile variables, this is perhaps slightly faster.
815  const Numeric p_i = abs_p[i];
816  const Numeric t_i = abs_t[i];
817  const Numeric vmr_i = all_vmrs(this_species,i);
818 
819  //out3 << " p = " << p_i << " Pa\n";
820 
821  // Calculate total number density from pressure and temperature.
822  // n = n0*T0/p0 * p/T or n = p/kB/t, ideal gas law
823  // const Numeric n = p_i / BOLTZMAN_CONST / t_i;
824  // This is not needed anymore, since we now calculate true cross
825  // sections, which do not contain the n.
826 
827  // For the pressure broadening, we also need the partial pressure:
828  const Numeric p_partial = p_i * vmr_i;
829 
830  // Get handle on xsec for this pressure level i.
831  // Watch out! This is output, we have to be careful not to
832  // introduce race conditions when writing to it.
833  VectorView xsec_i_attenuation = xsec_attenuation(Range(joker),i);
834  VectorView xsec_i_phase = xsec_phase(Range(joker),i);
835 
836 
837 // if (omp_in_parallel())
838 // cout << "omp_in_parallel: true\n";
839 // else
840 // cout << "omp_in_parallel: false\n";
841 
842 
843  // Prepare a variable that can be used by the individual LBL
844  // threads to add up absorption:
845  Index n_lbl_threads;
846  if (arts_omp_in_parallel())
847  {
848  // If we already are running parallel, then the LBL loop
849  // will not be parallelized.
850  n_lbl_threads = 1;
851  }
852  else
853  {
854  n_lbl_threads = arts_omp_get_max_threads();
855  }
856  Matrix xsec_accum_attenuation(n_lbl_threads, xsec_i_attenuation.nelem(), 0);
857  Matrix xsec_accum_phase(n_lbl_threads, xsec_i_phase.nelem(), 0);
858 
859 
860  // Loop all lines:
861  if (nl)
862 #pragma omp parallel for \
863  if (!arts_omp_in_parallel() \
864  && nl >= arts_omp_get_max_threads()) \
865  firstprivate(ls_attenuation, ls_phase, fac, f_local, aux)
866  for ( Index l=0; l< nl; ++l )
867  {
868  // Skip remaining iterations if an error occurred
869  if (failed) continue;
870 
871 // if (omp_in_parallel())
872 // cout << "LBL: omp_in_parallel: true\n";
873 // else
874 // cout << "LBL: omp_in_parallel: false\n";
875 
876 
877  // The try block here is necessary to correctly handle
878  // exceptions inside the parallel region.
879  try
880  {
881  // Copy f_grid to the beginning of f_local. There is one
882  // element left at the end of f_local.
883  // THIS HAS TO BE INSIDE THE LINE LOOP, BECAUSE THE CUTOFF
884  // FREQUENCY IS ALWAYS PUT IN A DIFFERENT PLACE!
885  f_local[Range(0,nf)] = f_grid;
886 
887  // This will hold the actual number of frequencies to add to
888  // xsec later on:
889  Index nfl = nf;
890 
891  // This will hold the actual number of frequencies for the
892  // call to the lineshape functions later on:
893  Index nfls = nf;
894 
895  // abs_lines[l] is used several times, this construct should be
896  // faster (Oliver Lemke)
897  const LineRecord& l_l = abs_lines[l]; // which line are we dealing with
898 
899  // Make sure that catalogue version is either 3 or 4 (no other
900  // versions are supported yet):
901  if ( 3!=l_l.Version() && 4!=l_l.Version() )
902  {
903  ostringstream os;
904  os << "Unknown spectral line catalogue version (artscat-"
905  << l_l.Version() << ").\n"
906  << "Allowed are artscat-3 and artscat-4.";
907  throw runtime_error(os.str());
908  }
909 
910  // Center frequency in vacuum:
911  Numeric F0 = l_l.F();
912 
913  // Intensity is already in the right units (Hz*m^2). It also
914  // includes already the isotopologue ratio. Needs only to be
915  // coverted to the actual temperature and multiplied by total
916  // number density and lineshape.
917  Numeric intensity = l_l.I0();
918 
919  // Lower state energy is already in the right unit (Joule).
920  Numeric e_lower = l_l.Elow();
921 
922  // Upper state energy:
923  Numeric e_upper = e_lower + F0 * PLANCK_CONST;
924 
925  // Get the ratio of the partition function.
926  // This will throw a runtime error if no data exists.
927  // Important: This function needs both the reference
928  // temperature and the actual temperature, because the
929  // reference temperature can be different for each line,
930  // even of the same species.
931  Numeric part_fct_ratio =
933  t_i );
934 
935  // Boltzmann factors
936  Numeric nom = exp(- e_lower / ( BOLTZMAN_CONST * t_i ) ) -
937  exp(- e_upper / ( BOLTZMAN_CONST * t_i ) );
938 
939  Numeric denom = exp(- e_lower / ( BOLTZMAN_CONST * l_l.Ti0() ) ) -
940  exp(- e_upper / ( BOLTZMAN_CONST * l_l.Ti0() ) );
941 
942 
943  // intensity at temperature
944  // (calculate the line intensity according to the standard
945  // expression given in HITRAN)
946  intensity *= part_fct_ratio * nom / denom;
947 
948  if (lineshape_norm_data[ind_lsn].Name() == "quadratic")
949  {
950  // in case of the quadratic normalization factor use the
951  // so called 'microwave approximation' of the line intensity
952  // given by
953  // P. W. Rosenkranz, Chapter 2, Eq.(2.16), in M. A. Janssen,
954  // Atmospheric Remote Sensing by Microwave Radiometry,
955  // John Wiley & Sons, Inc., 1993
956  Numeric mafac = (PLANCK_CONST * F0) / (2.000e0 * BOLTZMAN_CONST
957  * t_i);
958  intensity = intensity * mafac / sinh(mafac);
959  }
960 
961  // 2. Calculate the pressure broadened line width and the pressure
962  // shifted center frequency.
963  //
964  // Here there is a difference betweeen catalogue version 4
965  // (from ESA planetary study) and earlier catalogue versions.
966  Numeric gamma; // The line width.
967  Numeric deltaf; // Pressure shift.
968  if (l_l.Version() == 4)
969  {
971  deltaf,
972  p_i,
973  t_i,
974  all_vmrs(joker,i),
975  this_species,
976  broad_spec_locations,
977  l_l,
978  verbosity);
979  }
980  else if (l_l.Version() == 3)
981  {
982  // Numeric Tgam;
983  // Numeric Nair;
984  // Numeric Agam;
985  // Numeric Sgam;
986  // Numeric Nself;
987  // Numeric Psf;
988 
989  // if (l_l.Version() == 3)
990  // {
991  const Numeric Tgam = l_l.Tgam();
992  const Numeric Agam = l_l.Agam();
993  const Numeric Nair = l_l.Nair();
994  const Numeric Sgam = l_l.Sgam();
995  const Numeric Nself = l_l.Nself();
996  const Numeric Psf = l_l.Psf();
997  // }
998  // else
999  // {
1000  // static bool warn = false;
1001  // if (!warn)
1002  // {
1003  // CREATE_OUT0;
1004  // warn = true;
1005  // out0 << " WARNING: Using artscat version 4 for calculations is currently\n"
1006  // << " just a hack and results are most likely wrong!!!\n";
1007  // }
1008  // Tgam = l_l.Tgam();
1009  // // Use hardcoded mixing ratios for air
1010  // Agam = l_l.Gamma_N2() * 0.79 + l_l.Gamma_O2() * 0.21;
1011  // Nair = l_l.Gam_N_N2() * 0.79 + l_l.Gam_N_O2() * 0.21;
1012  // Sgam = l_l.Gamma_self();
1013  // Nself = l_l.Gam_N_self();
1014  // Psf = l_l.Delta_N2() * 0.79 + l_l.Delta_O2() * 0.21;
1015  // }
1016 
1017  // Get pressure broadened line width:
1018  // (Agam is in Hz/Pa, abs_p is in Pa, gamma is in Hz)
1019  const Numeric theta = Tgam / t_i;
1020  const Numeric theta_Nair = pow(theta, Nair);
1021 
1022  gamma = Agam * theta_Nair * (p_i - p_partial)
1023  + Sgam * pow(theta, Nself) * p_partial;
1024 
1025 
1026  // Pressure shift:
1027  // The T dependence is connected to that of agam by:
1028  // n_shift = .25 + 1.5 * n_agam
1029  // Theta has been initialized above.
1030  deltaf = Psf * p_i *
1031  std::pow( theta , (Numeric).25 + (Numeric)1.5*Nair );
1032 
1033  }
1034  else
1035  {
1036  // There is a runtime error check for allowed artscat versions
1037  // further up.
1038  assert(false);
1039  }
1040 
1041  // Apply pressure shift:
1042  F0 += deltaf;
1043 
1044  // 3. Doppler broadening without the sqrt(ln(2)) factor, which
1045  // seems to be redundant.
1046  const Numeric sigma = F0 * doppler_const *
1047  sqrt( t_i / l_l.IsotopologueData().Mass());
1048 
1049 // cout << l_l.IsotopologueData().Name() << " " << l_l.F() << " "
1050 // << Nair << " " << Agam << " " << Sgam << " " << Nself << endl;
1051 
1052 
1053 
1054  // Indices pointing at begin/end frequencies of f_grid or at
1055  // the elements that have to be calculated in case of cutoff
1056  Index i_f_min = 0;
1057  Index i_f_max = nf-1;
1058 
1059  // cutoff ?
1060  if ( cut )
1061  {
1062  // Check whether we have elements in ls that can be
1063  // ignored at lower frequencies of f_grid.
1064  //
1065  // Loop through all frequencies, finding min value and
1066  // set all values to zero on that way.
1067  while ( i_f_min < nf && (F0 - cutoff) > f_grid[i_f_min] )
1068  {
1069  // ls[i_f_min] = 0;
1070  ++i_f_min;
1071  }
1072 
1073 
1074  // Check whether we have elements in ls that can be
1075  // ignored at higher frequencies of f_grid.
1076  //
1077  // Loop through all frequencies, finding max value and
1078  // set all values to zero on that way.
1079  while ( i_f_max >= 0 && (F0 + cutoff) < f_grid[i_f_max] )
1080  {
1081  // ls[i_f_max] = 0;
1082  --i_f_max;
1083  }
1084 
1085  // Append the cutoff frequency to f_local:
1086  ++i_f_max;
1087  f_local[i_f_max] = F0 + cutoff;
1088 
1089  // Number of frequencies to calculate:
1090  nfls = i_f_max - i_f_min + 1; // Add one because indices
1091  // are pointing to first and
1092  // last valid element. This
1093  // is for the lineshape
1094  // calls.
1095  nfl = nfls -1; // This is for xsec.
1096  }
1097  else
1098  {
1099  // Nothing to do here. Note that nfl and nfls are both still set to nf.
1100  }
1101 
1102  // cout << "nf, nfl, nfls = " << nf << ", " << nfl << ", " << nfls << ".\n";
1103 
1104  // Maybe there are no frequencies left to compute? Note that
1105  // the number that counts here is nfl, since only these are
1106  // the `real' frequencies, for which xsec is changed. nfls
1107  // will always be at least one, because it contains the cutoff.
1108  if ( nfl > 0 )
1109  {
1110  // cout << ls << endl
1111  // << "aux / F0 / gamma / sigma" << aux << "/" << F0 << "/" << gamma << "/" << sigma << endl
1112  // << f_local[Range(i_f_min,nfls)] << endl
1113  // << nfls << endl;
1114 
1115  // Calculate the line shape:
1116  lineshape_data[ind_ls].Function()(ls_attenuation,ls_phase,
1117  aux,F0,gamma,sigma,
1118  f_local[Range(i_f_min,nfls)]);
1119 
1120  // Calculate the chosen normalization factor:
1121  lineshape_norm_data[ind_lsn].Function()(fac,F0,
1122  f_local[Range(i_f_min,nfls)],
1123  t_i);
1124 
1125  // Get a handle on the range of xsec that we want to change.
1126  // We use nfl here, which could be one less than nfls.
1127  VectorView this_xsec_attenuation = xsec_accum_attenuation(arts_omp_get_thread_num(), Range(i_f_min,nfl));
1128  VectorView this_xsec_phase = xsec_accum_phase(arts_omp_get_thread_num(), Range(i_f_min,nfl));
1129 
1130  // Get handles on the range of ls and fac that we need.
1131  VectorView this_ls_attenuation = ls_attenuation[Range(0,nfl)];
1132  VectorView this_ls_phase = ls_phase[Range(0,nfl)];
1133  VectorView this_fac = fac[Range(0,nfl)];
1134 
1135  // cutoff ?
1136  if ( cut )
1137  {
1138  // Subtract baseline for cutoff frequency
1139  // The index nfls-1 should be exactly the index pointing
1140  // to the value at the cutoff frequency.
1141  // Subtract baseline from xsec.
1142  // this_xsec -= base;
1143  this_ls_attenuation -= ls_attenuation[nfls-1];
1144  if (calc_phase) this_ls_phase -= ls_phase[nfls-1];
1145  }
1146 
1147  // Add line to xsec.
1148  {
1149  // To make the loop a bit faster, precompute all constant
1150  // factors. These are:
1151  // 1. Total number density of the air. --> Not
1152  // anymore, we now to real cross-sections
1153  // 2. Line intensity.
1154  // 3. Isotopologue ratio.
1155  //
1156  // The isotopologue ratio must be applied here, since we are
1157  // summing up lines belonging to different isotopologues.
1158 
1159  // const Numeric factors = n * intensity * l_l.IsotopologueData().Abundance();
1160 // const Numeric factors = intensity * l_l.IsotopologueData().Abundance();
1161  const Numeric factors = intensity
1162  * isotopologue_ratios.getParam(l_l.Species(), l_l.Isotopologue(), 0);
1163 
1164  // We have to do:
1165  // xsec(j,i) += factors * ls[j] * fac[j];
1166  //
1167  // We use ls as a dummy to compute the product, then add it
1168  // to this_xsec.
1169 
1170  this_ls_attenuation *= this_fac;
1171  this_ls_attenuation *= factors;
1172  this_xsec_attenuation += this_ls_attenuation;
1173 
1174  if (calc_phase)
1175  {
1176  this_ls_phase *= this_fac;
1177  this_ls_phase *= factors;
1178  this_xsec_phase += this_ls_phase;
1179  }
1180  }
1181  }
1182 
1183  } // end of try block
1184  catch (runtime_error e)
1185  {
1186 #pragma omp critical (xsec_species_fail)
1187  { fail_msg = e.what(); failed = true; }
1188  }
1189 
1190  } // end of parallel LBL loop
1191 
1192  // Bail out if an error occurred in the LBL loop
1193  if (failed) continue;
1194 
1195  // Now we just have to add up all the rows of xsec_accum:
1196  for (Index j=0; j<xsec_accum_attenuation.nrows(); ++j)
1197  {
1198  xsec_i_attenuation += xsec_accum_attenuation(j, Range(joker));
1199  }
1200 
1201  if (calc_phase)
1202  for (Index j=0; j<xsec_accum_phase.nrows(); ++j)
1203  {
1204  xsec_i_phase += xsec_accum_phase(j, Range(joker));
1205  }
1206  } // end of parallel pressure loop
1207 
1208  if (failed) throw runtime_error("Run-time error in function: xsec_species\n" + fail_msg);
1209 
1210 }
1211 
1212 
1213 
1214 
1215 
1228 {
1229  // Planck constant [Js]
1230  extern const Numeric PLANCK_CONST;
1231 
1232  // Speed of light [m/s]
1233  extern const Numeric SPEED_OF_LIGHT;
1234 
1235  // Constant to convert lower state energy from cm^-1 to J
1236  const Numeric lower_energy_const = PLANCK_CONST * SPEED_OF_LIGHT * 1E2;
1237 
1238  return e*lower_energy_const;
1239 }
1240 
1241 
1242 //======================================================================
1243 // Functions related to species
1244 //======================================================================
1245 
1247 
1261 {
1262  // For the return value:
1263  Index mspecies;
1264 
1265  // Trim whitespace
1266  name.trim();
1267 
1268  // cout << "name / def = " << name << " / " << def << endl;
1269 
1270  // Look for species name in species map:
1271  map<String, Index>::const_iterator mi = SpeciesMap.find(name);
1272  if ( mi != SpeciesMap.end() )
1273  {
1274  // Ok, we've found the species. Set mspecies.
1275  mspecies = mi->second;
1276  }
1277  else
1278  {
1279  // The species does not exist!
1280  mspecies = -1;
1281  }
1282 
1283  return mspecies;
1284 }
1285 
1286 
1288 
1302 {
1303  // Species lookup data:
1305 
1306  // Assert that spec_ind is inside species data. (This is an assertion,
1307  // because species indices should never be user input, but set by the
1308  // program automatically, based on species names.)
1309  assert( spec_ind>=0 );
1310  assert( spec_ind<species_data.nelem() );
1311 
1312  // A reference to the relevant record of the species data:
1313  const SpeciesRecord& spr = species_data[spec_ind];
1314 
1315  return spr.Name();
1316 }
1317 
1318 
1319 //======================================================================
1320 // Functions to convert the accuracy index to ARTS units
1321 //======================================================================
1322 
1323 // ********* for HITRAN database *************
1324 //convert HITRAN index for line position accuracy to ARTS
1325 //units (Hz).
1326 
1328  Numeric& mdf,
1329  const Index& df
1330  )
1331 {
1332  switch ( df )
1333  {
1334  case 0:
1335  {
1336  mdf = -1;
1337  break;
1338  }
1339  case 1:
1340  {
1341  mdf = 30*1E9;
1342  break;
1343  }
1344  case 2:
1345  {
1346  mdf = 3*1E9;
1347  break;
1348  }
1349  case 3:
1350  {
1351  mdf = 300*1E6;
1352  break;
1353  }
1354  case 4:
1355  {
1356  mdf = 30*1E6;
1357  break;
1358  }
1359  case 5:
1360  {
1361  mdf = 3*1E6;
1362  break;
1363  }
1364  case 6:
1365  {
1366  mdf = 0.3*1E6;
1367  break;
1368  }
1369  }
1370 }
1371 
1372 //convert HITRAN index for intensity and halfwidth accuracy to ARTS
1373 //units (relative difference).
1375  Numeric& mdh,
1376  const Index& dh
1377  )
1378 {
1379  switch ( dh )
1380  {
1381  case 0:
1382  {
1383  mdh = -1;
1384  break;
1385  }
1386  case 1:
1387  {
1388  mdh = -1;
1389  break;
1390  }
1391  case 2:
1392  {
1393  mdh = -1;
1394  break;
1395  }
1396  case 3:
1397  {
1398  mdh = 30;
1399  break;
1400  }
1401  case 4:
1402  {
1403  mdh = 20;
1404  break;
1405  }
1406  case 5:
1407  {
1408  mdh = 10;
1409  break;
1410  }
1411  case 6:
1412  {
1413  mdh =5;
1414  break;
1415  }
1416  case 7:
1417  {
1418  mdh =2;
1419  break;
1420  }
1421  case 8:
1422  {
1423  mdh =1;
1424  break;
1425  }
1426  }
1427  mdh=mdh/100;
1428 }
1429 
1430 // ********* for MYTRAN database *************
1431 //convert MYTRAN index for intensity and halfwidth accuracy to ARTS
1432 //units (relative difference).
1434  Numeric& mdh,
1435  const Index & dh
1436  )
1437 {
1438  switch ( dh )
1439  {
1440  case 0:
1441  {
1442  mdh = 200;
1443  break;
1444  }
1445  case 1:
1446  {
1447  mdh = 100;
1448  break;
1449  }
1450  case 2:
1451  {
1452  mdh = 50;
1453  break;
1454  }
1455  case 3:
1456  {
1457  mdh = 30;
1458  break;
1459  }
1460  case 4:
1461  {
1462  mdh = 20;
1463  break;
1464  }
1465  case 5:
1466  {
1467  mdh = 10;
1468  break;
1469  }
1470  case 6:
1471  {
1472  mdh =5;
1473  break;
1474  }
1475  case 7:
1476  {
1477  mdh = 2;
1478  break;
1479  }
1480  case 8:
1481  {
1482  mdh = 1;
1483  break;
1484  }
1485  case 9:
1486  {
1487  mdh = 0.5;
1488  break;
1489  }
1490  }
1491  mdh=mdh/100;
1492 }
1493 
1494 ostream& operator<< (ostream &os, const LineshapeSpec& lsspec)
1495 {
1496  os << "LineshapeSpec Index: " << lsspec.Ind_ls() << ", NormIndex: " << lsspec.Ind_lsn()
1497  << ", Cutoff: " << lsspec.Cutoff() << endl;
1498 
1499  return os;
1500 }
1501 
1502 
1531  MatrixView xsec_phase,
1532  const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
1533  const ArrayOfArrayOfIndex& line_mixing_data_lut,
1534  ConstVectorView f_grid,
1535  ConstVectorView abs_p,
1536  ConstVectorView abs_t,
1537  ConstMatrixView all_vmrs,
1538  const ArrayOfArrayOfSpeciesTag& abs_species,
1539  const Index this_species,
1540  const ArrayOfLineRecord& abs_lines,
1541  const Index ind_ls,
1542  const Index ind_lsn,
1543  const Numeric cutoff,
1544  const SpeciesAuxData& isotopologue_ratios,
1545  const Verbosity& verbosity )
1546 {
1547  if(abs_species[this_species][0].LineMixingType() == SpeciesTag::LINE_MIXING_TYPE_NONE) // If no linemixing data exists
1548  xsec_species(xsec_attenuation, xsec_phase, f_grid, abs_p, abs_t, all_vmrs,
1549  abs_species, this_species, abs_lines, ind_ls, ind_lsn, cutoff,
1550  isotopologue_ratios, verbosity);
1551  else // If linemixing data exists
1552  {
1553  switch(abs_species[this_species][0].LineMixingType())
1554  {
1556  xsec_species_line_mixing_2nd_order(xsec_attenuation, xsec_phase, line_mixing_data,
1557  line_mixing_data_lut, f_grid, abs_p, abs_t, all_vmrs,
1558  abs_species, this_species, abs_lines, ind_ls, ind_lsn,
1559  cutoff, isotopologue_ratios, verbosity);
1560  break;
1561  // default:
1562  // //ERROR
1563  }
1564  }
1565 
1566 }
1567 
1568 
1599  MatrixView xsec_phase,
1600  const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
1601  const ArrayOfArrayOfIndex& line_mixing_data_lut,
1602  ConstVectorView f_grid,
1603  ConstVectorView abs_p,
1604  ConstVectorView abs_t,
1605  ConstMatrixView all_vmrs,
1606  const ArrayOfArrayOfSpeciesTag& abs_species,
1607  const Index this_species,
1608  const ArrayOfLineRecord& abs_lines,
1609  const Index ind_ls,
1610  const Index ind_lsn,
1611  const Numeric cutoff,
1612  const SpeciesAuxData& isotopologue_ratios,
1613  const Verbosity& verbosity )
1614 {
1615  // Must have the phase
1617  if (! lineshape_data[ind_ls].Phase())
1618  {
1619  ostringstream os;
1620  os << "This is an error message. You are using " << lineshape_data[ind_ls].Name() <<
1621  ".\n"<<"This line shape does not include phase in its calculations and\nis therefore invalid for " <<
1622  "second order line mixing.\n\n";
1623  throw runtime_error(os.str());
1624  }
1625  const ArrayOfLineMixingRecord& data = line_mixing_data[this_species];
1626  const ArrayOfIndex& lut = line_mixing_data_lut[this_species];
1627 
1628  // Helper variables
1629  Matrix attenuation(f_grid.nelem(),1), phase(f_grid.nelem(),1);
1630 
1631  //Since we need to be able to modify the line.
1633  ll.resize(1);
1634  // Variable containing modifying elements
1635  Vector a(2);
1636 
1637  for(Index jj=0;jj<abs_p.nelem(); jj++)
1638  {
1639  const Numeric p = abs_p[jj], t = abs_t[jj];
1640 
1641  for(Index ii = 0; ii < lut.nelem(); ii++)
1642  {
1643  // Since we need cross section per line to do the mixing.
1644  ll[0] = abs_lines[ii];//Possible speed-up by separating all non-split lines and doing those as a batch?
1645 
1646  // Since we need line mixing only for selected lines the rest should ignore by doing nothing.
1647  a=0;
1648 
1649  if(lut[ii]!=-1)
1650  {
1651  const Vector& dat = data[lut[ii]].Data();
1652  if( dat.nelem() != 10 )
1653  {
1654  ostringstream os;
1655  os << "This is an error message. You have not defined the linemixing data vector " <<
1656  "properly. It should be of length 10, yet your version is of length: " << dat.nelem() << ".\n" <<
1657  "The structure of the line mixing data Vector should be as follows:\n" <<
1658  "data[0] is the first order zeroth phase correction\n" <<
1659  "data[1] is the first order first phase correction\n" <<
1660  "data[2] is the second order zeroth phase correction\n" <<
1661  "data[3] is the second order first phase correction \n" <<
1662  "data[4] is the second order zeroth frequency correction\n" <<
1663  "data[5] is the second order first frequency correction \n" <<
1664  "data[6] is the reference temperature\n" <<
1665  "data[7] is the first order phase correction temperature dependence exponent\n" <<
1666  "data[8] is the second order attenuation correction temperature dependence exponent\n" <<
1667  "data[9] is the second order frequency correction temperature dependence exponent";
1668  throw std::runtime_error(os.str());
1669  }
1670 
1671  // First order phase correction
1672  a[0] = p
1673  * ( ( dat[0] + dat[1] * ( dat[6]/t-1. ) )
1674  * pow( dat[6]/t, dat[7] ) );
1675 
1676  // Second order attenuation correction
1677  a[1] = p * p
1678  * ( ( dat[2] + dat[3] * ( dat[6]/t-1. ) )
1679  * pow( dat[6]/t, dat[8] ) );
1680 
1681  // Second order line-center correction
1682  ll[0].setF( p * p
1683  * ( ( dat[4] + dat[5] * ( dat[6]/t-1. ) )
1684  * pow( dat[6]/t, dat[9] ) ) + ll[0].F());
1685  }
1686 
1687  attenuation=0;
1688  phase=0;
1689 
1690  Vector P(1,p),T(1,t);
1691  Matrix V(all_vmrs.nrows(),1);V(joker,0)=all_vmrs(joker,jj);
1692 
1693  xsec_species(attenuation, phase, f_grid, P, T, V,
1694  abs_species, this_species, ll, ind_ls, ind_lsn, cutoff,
1695  isotopologue_ratios, verbosity);
1696 
1697  // Do the actual line mixing and add this to xsec_attenuation.
1698  xsec_phase(joker,jj) += phase(joker,0);
1699  phase *= a[0];
1700  xsec_attenuation(joker,jj) += phase(joker,0); // First order phase correction
1701  xsec_attenuation(joker,jj) += attenuation(joker,0); // Zeroth order attenuation
1702  attenuation *= a[1];
1703  xsec_attenuation(joker,jj) += attenuation(joker,0); // Second order attenuation correction
1704 
1705  }
1706  }
1707 }
Numeric Ti0() const
Reference temperature for I0 in K:
Definition: linerecord.h:368
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
const Numeric PLANCK_CONST
Global constant, the Planck constant [Js].
The VectorView class.
Definition: matpackI.h:372
Numeric N_foreign(const Index i) const
ARTSCAT-4 foreign temperature exponents (dimensionless):
Definition: linerecord.h:445
Index nelem() const
Number of elements.
Definition: array.h:176
const Numeric & Mass() const
Mass of the isotopologue.
Definition: absorption.h:260
Declarations having to do with the four output streams.
const Numeric BOLTZMAN_CONST
Global constant, the Boltzmann constant [J/K].
Numeric Tgam() const
Reference temperature for AGAM and SGAM in K:
Definition: linerecord.h:398
Lineshape related specification like which lineshape to use, the normalizationfactor, and the cutoff.
Definition: absorption.h:141
bool arts_omp_in_parallel()
Wrapper for omp_in_parallel.
Definition: arts_omp.cc:66
The Vector class.
Definition: matpackI.h:556
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)
This will work as a wrapper for linemixing when abs_species contain relevant data.
Definition: absorption.cc:1530
Numeric Agam() const
Air broadened width in Hz/Pa:
Definition: linerecord.h:374
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:68
int arts_omp_get_thread_num()
Wrapper for omp_get_thread_num.
Definition: arts_omp.cc:83
const IsotopologueRecord & IsotopologueData() const
The matching IsotopologueRecord from species_data.
Definition: linerecord.cc:57
The range class.
Definition: matpackI.h:148
bool ReadFromStream(String &artsid, istream &is, Index nparams, const Verbosity &verbosity)
Read parameters from input stream.
Definition: absorption.cc:104
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:47
ConstIterator1D begin() const
Return const iterator to first element.
Definition: matpackI.cc:346
ConstIterator1D end() const
Return const iterator behind last element.
Definition: matpackI.cc:354
This file contains basic functions to handle ASCII files.
void iso(Array< IsotopologueRecord >::iterator &ii, String name, const ArrayOfNumeric &coeff, const Index &coefftype)
ostream & operator<<(ostream &os, const SpeciesRecord &sr)
Output operator for SpeciesRecord.
Definition: absorption.cc:333
Numeric Elow() const
Lower state energy in cm^-1:
Definition: linerecord.h:371
const Index & Ind_lsn() const
Return the index of the normalization factor.
Definition: absorption.h:165
void setParam(Index species, Index isotopologue, Index col, Numeric v)
Set parameter.
Definition: absorption.h:455
Numeric Psf() const
The pressure shift parameter in Hz/Pa.
Definition: linerecord.h:345
void convMytranIER(Numeric &mdh, const Index &dh)
Definition: absorption.cc:1433
void find_broad_spec_locations(ArrayOfIndex &broad_spec_locations, const ArrayOfArrayOfSpeciesTag &abs_species, const Index this_species)
Find the location of all broadening species in abs_species.
Definition: absorption.cc:372
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
const Array< LineshapeRecord > lineshape_data
The lookup data for the different lineshapes.
Definition: lineshapes.cc:1925
const Array< SpeciesRecord > species_data
Species Data.
String species_name_from_species_index(const Index spec_ind)
Return species name for given species index.
Definition: absorption.cc:1301
const Array< IsotopologueRecord > & Isotopologue() const
Definition: absorption.h:425
void calc_gamma_and_deltaf_artscat4(Numeric &gamma, Numeric &deltaf, const Numeric p, const Numeric t, ConstVectorView vmrs, const Index this_species, const ArrayOfIndex &broad_spec_locations, const LineRecord &l_l, const Verbosity &verbosity)
Calculate line width and pressure shift for artscat4.
Definition: absorption.cc:432
Numeric Nself() const
SGAM temperature exponent (dimensionless):
Definition: linerecord.h:392
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The global header file for ARTS.
Index Version() const
Return the version number.
Definition: linerecord.h:298
bool is_sorted(ConstVectorView x)
Checks if a vector is sorted in ascending order.
Definition: logic.cc:253
void convHitranIERSH(Numeric &mdh, const Index &dh)
Definition: absorption.cc:1374
Index mqcoeffinterporder
Definition: absorption.h:360
static String BroadSpecName(const Index i)
Return the name of an artscat-4 broadening species, as function of its broadening species index...
Definition: linerecord.h:559
Numeric Delta_foreign(const Index i) const
ARTSCAT-4 pressure shift parameters in Hz/Pa :
Definition: linerecord.h:448
Contains the lookup data for one species.
Definition: absorption.h:367
const ArrayOfMatrix & getParams() const
Return a constant reference to the parameters.
Definition: absorption.h:461
Index Isotopologue() const
The index of the isotopologue species that this line belongs to.
Definition: linerecord.h:307
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
void extract(T &x, String &line, Index n)
Extract something from the beginning of a string.
Definition: mystring.h:333
#define abs(x)
Definition: continua.cc:20458
const Joker joker
const Numeric & Cutoff() const
Return the cutoff frequency (in Hz).
Definition: absorption.h:172
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
std::map< String, Index > SpeciesMap
The map associated with species_data.
Definition: absorption.cc:50
The Matrix class.
Definition: matpackI.h:788
const Index & Ind_ls() const
Return the index of this lineshape.
Definition: absorption.h:160
Index nelem() const
Number of elements.
Definition: mystring.h:278
const String & Name() const
Definition: absorption.h:423
Declarations required for the calculation of absorption coefficients.
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
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)
Calculate line absorption cross sections for one tag group.
Definition: absorption.cc:594
const Array< LineshapeNormRecord > lineshape_norm_data
The lookup data for the different normalization factors to the lineshapes.
Definition: lineshapes.cc:2030
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:1260
void trim()
Trim leading and trailing whitespace.
Definition: mystring.h:253
const Numeric AVOGADROS_NUMB
Global constant, the Avogadro&#39;s number [molec/kg].
Header file for logic.cc.
Numeric Sgam() const
Self broadened width in Hz/Pa:
Definition: linerecord.h:380
Header file for interpolation_poly.cc.
Numeric CalculatePartitionFctAtTempFromCoeff(Numeric temperature) const
Definition: absorption.cc:56
void checkIsotopologueRatios(const ArrayOfArrayOfSpeciesTag &abs_species, const SpeciesAuxData &isoratios)
Check that isotopologue ratios for the given species are correctly defined.
Definition: absorption.cc:245
void initParams(Index nparams)
Resize according to builtin isotopologues in species data.
Definition: absorption.cc:90
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)
This is the second order line mixing correction as presented by: Makarov, D.S, M.Yu.
Definition: absorption.cc:1598
A constant view of a Vector.
Definition: matpackI.h:292
A constant view of a Matrix.
Definition: matpackI.h:596
static Index BroadSpecSpecIndex(const Index i)
Return the internal species index (index in species_data) of an artscat-4 broadening species...
Definition: linerecord.cc:64
const Index & Speciesindex() const
Definition: absorption.h:491
Numeric getParam(Index species, Index isotopologue, Index col) const
Get single parameter value.
Definition: absorption.h:449
void convHitranIERF(Numeric &mdf, const Index &df)
Definition: absorption.cc:1327
Structure to store a grid position for higher order interpolation.
The constant iterator class for sub vectors.
Definition: matpackI.h:249
void fillSpeciesAuxDataWithIsotopologueRatiosFromSpeciesData(SpeciesAuxData &sad)
Fill SpeciesAuxData with default isotopologue ratios from species data.
Definition: absorption.cc:305
#define CREATE_OUT3
Definition: messages.h:214
const Numeric & Abundance() const
Normal abundance ( = isotopologue ratio).
Definition: absorption.h:257
Auxiliary data for isotopologues.
Definition: absorption.h:439
Index Species() const
The index of the molecular species that this line belongs to.
Definition: linerecord.h:302
#define ll
Definition: continua.cc:21705
Numeric I0() const
The line intensity in m^2*Hz at the reference temperature Ti0.
Definition: linerecord.h:362
static Index NBroadSpec()
Return the number of artscat-4 foreign broadening species (6).
Definition: linerecord.h:554
const String & Name() const
Isotopologue name.
Definition: absorption.h:255
Numeric CalculatePartitionFctRatio(Numeric reference_temperature, Numeric actual_temperature) const
Calculate partition function ratio.
Definition: absorption.h:295
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define CREATE_OUT2
Definition: messages.h:213
void define_species_map()
Define the species data map.
Definition: absorption.cc:322
const Numeric SPEED_OF_LIGHT
Numeric Gamma_foreign(const Index i) const
ARTSCAT-4 foreign broadening parameters in Hz/Pa :
Definition: linerecord.h:442
Numeric F() const
The line center frequency in Hz.
Definition: linerecord.h:339
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
Numeric CalculatePartitionFctAtTempFromData(Numeric temperature) const
Definition: absorption.cc:79
Input manipulator class for doubles to enable nan and inf parsing.
Definition: file.h:124
Spectral line catalog data.
Definition: linerecord.h:196
Numeric wavenumber_to_joule(Numeric e)
A little helper function to convert energy from units of wavenumber (cm^-1) to Joule (J)...
Definition: absorption.cc:1227
Numeric Nair() const
AGAM temperature exponent (dimensionless):
Definition: linerecord.h:386