ARTS  2.2.66
m_abs.cc
Go to the documentation of this file.
1 
2 /* Copyright (C) 2000-2012
3  Stefan Buehler <sbuehler@ltu.se>
4  Patrick Eriksson <patrick.eriksson@chalmers.se>
5  Axel von Engeln <engeln@uni-bremen.de>
6  Thomas Kuhn <tkuhn@uni-bremen.de>
7 
8  This program is free software; you can redistribute it and/or modify it
9  under the terms of the GNU General Public License as published by the
10  Free Software Foundation; either version 2, or (at your option) any
11  later version.
12 
13  This program is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
21  USA. */
22 
23 //
24 
33 #include <cmath>
34 #include <algorithm>
35 #include "arts.h"
36 #include "matpackI.h"
37 #include "array.h"
38 #include "messages.h"
39 #include "file.h"
40 #include "auto_md.h"
41 #include "math_funcs.h"
42 #include "make_array.h"
43 #include "make_vector.h"
44 #include "global_data.h"
45 #include "physics_funcs.h"
46 #include "absorption.h"
47 #include "continua.h"
48 #include "check_input.h"
49 #include "montecarlo.h"
50 #include "m_xml.h"
51 #include "optproperties.h"
52 #include "parameters.h"
53 #include "rte.h"
54 #include "xml_io.h"
55 
56 #ifdef ENABLE_NETCDF
57 #include <netcdf.h>
58 #include "nc_io.h"
59 #endif
60 
61 
62 extern const Numeric ELECTRON_CHARGE;
63 extern const Numeric ELECTRON_MASS;
64 extern const Numeric PI;
65 extern const Numeric SPEED_OF_LIGHT;
66 extern const Numeric VACUUM_PERMITTIVITY;
67 
68 
69 
70 
71 
72 /* Workspace method: Doxygen documentation will be auto-generated */
73 void AbsInputFromRteScalars(// WS Output:
74  Vector& abs_p,
75  Vector& abs_t,
76  Matrix& abs_vmrs,
77  // WS Input:
78  const Numeric& rtp_pressure,
79  const Numeric& rtp_temperature,
80  const Vector& rtp_vmr,
81  const Verbosity&)
82 {
83  // Prepare abs_p:
84  abs_p.resize(1);
85  abs_p = rtp_pressure;
86 
87  // Prepare abs_t:
88  abs_t.resize(1);
89  abs_t = rtp_temperature;
90 
91  // Prepare abs_vmrs:
92  abs_vmrs.resize(rtp_vmr.nelem(),1);
93  abs_vmrs = rtp_vmr;
94 }
95 
96 
97 /* Workspace method: Doxygen documentation will be auto-generated */
99  ArrayOfArrayOfLineRecord& abs_lines_per_species,
100  // WS Input:
101  const ArrayOfArrayOfSpeciesTag& tgs,
102  const Verbosity&)
103 {
104  // Make abs_lines_per_species the right size:
105  abs_lines_per_species.resize( tgs.nelem() );
106 
107  for (Index i=0; i<tgs.nelem(); ++i)
108  {
109  abs_lines_per_species[i].resize(0);
110  }
111 }
112 
113 /* Workspace method: Doxygen documentation will be auto-generated */
115  ArrayOfLineRecord& abs_lines,
116  // Verbosity object:
117  const Verbosity&)
118 {
119  String fail_msg;
120  bool failed = false;
121 
122  // Loop over all lines, use member function to do conversion.
123 #pragma omp parallel for \
124  if (!arts_omp_in_parallel() \
125  && abs_lines.nelem() >= arts_omp_get_max_threads())
126  for ( Index i=0; i<abs_lines.nelem(); ++i )
127  {
128  try
129  {
130  abs_lines[i].ARTSCAT4FromARTSCAT3();
131  }
132  catch (runtime_error e)
133  {
134 #pragma omp critical (abs_linesArtscat4FromArtscat3_fail)
135  { fail_msg = e.what(); failed = true; }
136  }
137  }
138 
139  if (failed) throw runtime_error(fail_msg);
140 }
141 
142 
143 /* Workspace method: Doxygen documentation will be auto-generated */
145  // WS Output:
146  ArrayOfLineRecord& abs_lines,
147  // Control Parameters:
148  const String& filename,
149  const Numeric& fmin,
150  const Numeric& fmax,
151  const Verbosity& verbosity)
152 {
153  CREATE_OUT2;
154 
155  ifstream is;
156 
157  // Reset lines in case it already existed:
158  abs_lines.resize(0);
159 
160  out2 << " Reading file: " << filename << "\n";
161  open_input_file(is, filename);
162 
163  bool go_on = true;
164  while ( go_on )
165  {
166  LineRecord lr;
167  if ( lr.ReadFromHitran2001Stream(is, verbosity) )
168  {
169  // If we are here the read function has reached eof and has
170  // returned no data.
171  go_on = false;
172  }
173  else
174  {
175  if ( fmin <= lr.F() )
176  {
177  if ( lr.F() <= fmax )
178  abs_lines.push_back(lr);
179  else
180  go_on = false;
181  }
182  }
183  }
184  out2 << " Read " << abs_lines.nelem() << " lines.\n";
185 }
186 
187 
188 /* Workspace method: Doxygen documentation will be auto-generated */
189 void abs_linesReadFromHitran(// WS Output:
190  ArrayOfLineRecord& abs_lines,
191  // Control Parameters:
192  const String& filename,
193  const Numeric& fmin,
194  const Numeric& fmax,
195  const Verbosity& verbosity)
196 {
197  CREATE_OUT2;
198 
199  ifstream is;
200 
201  // Reset lines in case it already existed:
202  abs_lines.resize(0);
203 
204  String filename_lower = filename;
205  filename_lower.tolower();
206  ArrayOfString splitted_fname;
207  filename_lower.split(splitted_fname, "/");
208  if (splitted_fname.nelem())
209  {
210  filename_lower = splitted_fname[splitted_fname.nelem()-1];
211  }
212  else
213  {
214  throw std::runtime_error("Catalog filename is empty");
215  }
216 #ifdef USE_HITRAN2008
217  if (filename_lower.nelem() < 8
218  || (filename_lower.substr(0, 8) != "hitran08"
219  && filename_lower.substr(0, 8) != "hitran04"))
220  {
221  ostringstream os;
222  os << "'" << filename << "'\n"
223  << "does not appear to be a HITRAN 2008 catalogue. The catalog filename\n"
224  << "name must start with HITRAN08. This version of arts was compiled with\n"
225  << "support only for HITRAN 2008. To switch back to the latest HITRAN\n"
226  << "run 'cmake -DWITH_HITRAN2008=0 ..' and recompile arts";
227  throw std::runtime_error(os.str());
228  }
229 #else
230  if (filename_lower.nelem() < 10
231  || (filename_lower.substr(0, 10) != "hitran2012"
232  && filename_lower.substr(0, 10) != "hitran2016"))
233  {
234  ostringstream os;
235  os << "'" << filename << "'\n"
236  << "does not appear to be a HITRAN 2012 catalogue. The catalog filename\n"
237  << "must start with HITRAN2012. If you intend to use a HITRAN 2008 catalog\n"
238  << "run 'cmake -DWITH_HITRAN2008 ..' and recompile arts";
239  throw std::runtime_error(os.str());
240  }
241 #endif
242 
243  out2 << " Reading file: " << filename << "\n";
244  open_input_file(is, filename);
245 
246  bool go_on = true;
247  Index linenr = 0;
248  while ( go_on )
249  {
250  linenr++;
251  try {
252  LineRecord lr;
253  if ( lr.ReadFromHitran2004Stream(is, verbosity, fmin) )
254  {
255  // If we are here the read function has reached eof and has
256  // returned no data.
257  go_on = false;
258  }
259  else
260  {
261  if ( fmin <= lr.F() )
262  {
263  if ( lr.F() <= fmax )
264  abs_lines.push_back(lr);
265  else
266  go_on = false;
267  }
268  }
269  }
270  catch (runtime_error e)
271  {
272  ostringstream os;
273  os << e.what() << "\n";
274  os << "Error parsing line " << linenr << " from catalog.\n";
275  throw runtime_error(os.str());
276  }
277  }
278  out2 << " Read " << abs_lines.nelem() << " lines.\n";
279 }
280 
281 
282 /* Workspace method: Doxygen documentation will be auto-generated */
283 void abs_linesReadFromMytran2(// WS Output:
284  ArrayOfLineRecord& abs_lines,
285  // Control Parameters:
286  const String& filename,
287  const Numeric& fmin,
288  const Numeric& fmax,
289  const Verbosity& verbosity)
290 {
291  CREATE_OUT2;
292 
293  ifstream is;
294 
295  // Reset lines in case it already existed:
296  abs_lines.resize(0);
297 
298  out2 << " Reading file: " << filename << "\n";
299  open_input_file(is, filename);
300 
301  bool go_on = true;
302  while ( go_on )
303  {
304  LineRecord lr;
305  if ( lr.ReadFromMytran2Stream(is, verbosity) )
306  {
307  // If we are here the read function has reached eof and has
308  // returned no data.
309  go_on = false;
310  }
311  else
312  {
313  // lines are not necessarily frequency sorted
314  if ( fmin <= lr.F() )
315  if ( lr.F() <= fmax )
316  abs_lines.push_back(lr);
317  }
318  }
319  out2 << " Read " << abs_lines.nelem() << " lines.\n";
320 }
321 
322 
323 /* Workspace method: Doxygen documentation will be auto-generated */
324 void abs_linesReadFromJpl(// WS Output:
325  ArrayOfLineRecord& abs_lines,
326  // Control Parameters:
327  const String& filename,
328  const Numeric& fmin,
329  const Numeric& fmax,
330  const Verbosity& verbosity)
331 {
332  CREATE_OUT2;
333 
334  ifstream is;
335 
336  // Reset lines in case it already existed:
337  abs_lines.resize(0);
338 
339  out2 << " Reading file: " << filename << "\n";
340  open_input_file(is, filename);
341 
342  bool go_on = true;
343  while ( go_on )
344  {
345  LineRecord lr;
346  if ( lr.ReadFromJplStream(is, verbosity) )
347  {
348  // If we are here the read function has reached eof and has
349  // returned no data.
350  go_on = false;
351  }
352  else
353  {
354  // we expect lines to be sorted
355  if ( fmin <= lr.F() )
356  {
357  if ( lr.F() <= fmax )
358  abs_lines.push_back(lr);
359  else
360  go_on = false;
361  }
362  }
363  }
364  out2 << " Read " << abs_lines.nelem() << " lines.\n";
365 }
366 
367 
368 /* Workspace method: Doxygen documentation will be auto-generated */
369 void abs_linesReadFromArts(// WS Output:
370  ArrayOfLineRecord& abs_lines,
371  // Control Parameters:
372  const String& filename,
373  const Numeric& fmin,
374  const Numeric& fmax,
375  const Verbosity& verbosity)
376 {
377  xml_read_arts_catalogue_from_file (filename, abs_lines, fmin, fmax, verbosity);
378 }
379 
380 
381 /* Workspace method: Doxygen documentation will be auto-generated */
383  const String& output_file_format,
384  const ArrayOfArrayOfLineRecord& abs_lines_per_species,
385  // Control Parameters:
386  const String& basename,
387  const Verbosity& verbosity)
388 {
390 
391  for (ArrayOfArrayOfLineRecord::const_iterator it = abs_lines_per_species.begin();
392  it != abs_lines_per_species.end();
393  it++)
394  {
395  if (it->nelem())
396  {
397  String species_filename = basename;
398  if (basename.length() && basename[basename.length()-1] != '/')
399  species_filename += ".";
400 
401  species_filename += species_data[(*it)[0].Species()].Name() + ".xml";
402  WriteXML(output_file_format, *it, species_filename, 0,
403  "", "", "", verbosity);
404  }
405  }
406 }
407 
408 
409 /* Workspace method: Doxygen documentation will be auto-generated */
411  ArrayOfLineRecord& abs_lines,
412  const ArrayOfArrayOfSpeciesTag& abs_species,
413  // Control Parameters:
414  const String& basename,
415  const Numeric& fmin,
416  const Numeric& fmax,
417  const Verbosity& verbosity)
418 {
419  CREATE_OUT2;
421 
422  // Build a set of species indices. Duplicates are ignored.
423  set<Index> unique_species;
424  for (ArrayOfArrayOfSpeciesTag::const_iterator asp = abs_species.begin();
425  asp != abs_species.end(); asp++)
426  for (ArrayOfSpeciesTag::const_iterator sp = asp->begin();
427  sp != asp->end(); sp++)
428  {
429  // Of the four different types of SpeciesTags, only PLAIN and ZEEMAN
430  // actually require explicit spectral lines.
431  if (sp->Type()==SpeciesTag::TYPE_PLAIN ||
432  sp->Type()==SpeciesTag::TYPE_ZEEMAN) {
433  unique_species.insert(sp->Species());
434  }
435 
436 // if (sp->Isotopologue() == species_data[sp->Species()].Isotopologue().nelem()
437 // || !species_data[sp->Species()].Isotopologue()[sp->Isotopologue()].isContinuum())
438 // {
439 // unique_species.insert(sp->Species());
440 // }
441  }
442 
443  // Read catalogs for each identified species and put them all into
444  // abs_lines.
445  abs_lines.resize(0);
446  for (set<Index>::const_iterator it = unique_species.begin();
447  it != unique_species.end(); it++)
448  {
449  ArrayOfLineRecord more_abs_lines;
450  String tmpbasename = basename;
451  if (basename.length() && basename[basename.length()-1] != '/')
452  {
453  tmpbasename += '.';
454  }
455 
456  abs_linesReadFromArts(more_abs_lines, tmpbasename + (species_data[*it].Name()) + ".xml",
457  fmin, fmax, verbosity);
458  abs_lines.insert(abs_lines.end(), more_abs_lines.begin(), more_abs_lines.end());
459  }
460 
461  out2 << " Read " << abs_lines.nelem() << " lines in total.\n";
462 }
463 
464 
466 
477  ArrayOfLineRecord& abs_lines,
478  // Control Parameters:
479  const String& filename,
480  const Numeric& fmin,
481  const Numeric& fmax,
482  const Verbosity& verbosity)
483 {
484  CREATE_OUT2;
485 
486  // The input stream:
487  ifstream is;
488 
489  // We will use this line record to read in the line records in the
490  // file one after another:
491  LineRecord lr;
492 
493  // Reset lines in case it already existed:
494  abs_lines.resize(0);
495 
496  out2 << " Reading file: " << filename << "\n";
497  open_input_file(is, filename);
498 
499  // Get version tag and check that it corresponds to the current version.
500  {
501  String v;
502  is >> v;
503  if ( v!=lr.VersionString() )
504  {
505  ostringstream os;
506 
507  // If what we read is the version String, it should have at elast 9 characters.
508  if ( 9 <= v.nelem() )
509  {
510  if ( "ARTSCAT" == v.substr(0,7) )
511  {
512  os << "The ARTS line file you are trying contains a version tag "
513  << "different from the current version.\n"
514  << "Tag in file: " << v << "\n"
515  << "Current version: " << lr.Version();
516  throw runtime_error(os.str());
517  }
518  }
519 
520  os << "The ARTS line file you are trying to read does not contain a valid version tag.\n"
521  << "Probably it was created with an older version of ARTS that used different units.";
522  throw runtime_error(os.str());
523  }
524  }
525 
526  bool go_on = true;
527  while ( go_on )
528  {
529  if ( lr.ReadFromArtscat3Stream(is, verbosity) )
530  {
531  // If we are here the read function has reached eof and has
532  // returned no data.
533  go_on = false;
534  }
535  else
536  {
537  // lines are not necessarily frequency sorted
538  if ( fmin <= lr.F() && lr.F() <= fmax )
539  {
540  abs_lines.push_back(lr);
541  }
542  }
543  }
544  out2 << " Read " << abs_lines.nelem() << " lines.\n";
545 }
546 
547 
548 /* Workspace method: Doxygen documentation will be auto-generated */
550  ArrayOfArrayOfLineRecord& abs_lines_per_species,
551  // WS Input:
552  const ArrayOfArrayOfSpeciesTag& tgs,
553  // Control Parameters:
554  const ArrayOfString& filenames,
555  const ArrayOfString& formats,
556  const Vector& fmin,
557  const Vector& fmax,
558  const Verbosity& verbosity)
559 {
560  CREATE_OUT3;
561 
562  const Index n_tg = tgs.nelem(); // # tag groups
563  const Index n_cat = filenames.nelem(); // # tag Catalogues
564 
565  // Check that dimensions of the keyword parameters are consistent
566  // (must all be the same).
567 
568  if ( n_cat != formats.nelem() ||
569  n_cat != fmin.nelem() ||
570  n_cat != fmax.nelem() )
571  {
572  ostringstream os;
573  os << "abs_lines_per_speciesReadFromCatalogues: All keyword\n"
574  << "parameters must get the same number of arguments.\n"
575  << "You specified:\n"
576  << "filenames: " << n_cat << "\n"
577  << "formats: " << formats.nelem() << "\n"
578  << "fmin: " << fmin.nelem() << "\n"
579  << "fmax: " << fmax.nelem();
580  throw runtime_error(os.str());
581  }
582 
583  // Furthermore, the dimension must be
584  // smaller than or equal to the number of tag groups.
585 
586  if ( n_cat > n_tg )
587  {
588  ostringstream os;
589  os << "abs_lines_per_speciesReadFromCatalogues: You specified more\n"
590  << "catalugues than you have tag groups.\n"
591  << "You specified:\n"
592  << "Catalogues: " << n_cat << "\n"
593  << "tgs: " << n_tg;
594  throw runtime_error(os.str());
595  }
596 
597  // There must be at least one tag group and at least one catalogue:
598 
599  if ( n_cat < 1 ||
600  n_tg < 1 )
601  {
602  ostringstream os;
603  os << "abs_lines_per_speciesReadFromCatalogues: You must have at\n"
604  << "least one catalogue and at least one tag group.\n"
605  << "You specified:\n"
606  << "Catalogues: " << n_cat << "\n"
607  << "tgs: " << n_tg;
608  throw runtime_error(os.str());
609  }
610 
611  // There can be repetitions in the keyword paramters. We want to read
612  // and process each catalogue only once, so we'll compile a set of
613  // real catalogues, along with a data structure that tells us which
614  // tag groups should use this catalogue.
615 
616  MakeArray<String> real_filenames ( filenames[0] );
617  MakeArray<String> real_formats ( formats[0] );
618  MakeArray<Numeric> real_fmin ( fmin[0] );
619  MakeArray<Numeric> real_fmax ( fmax[0] );
620 
621  Array< ArrayOfIndex > real_tgs(1);
622  real_tgs[0].resize(1);
623  real_tgs[0][0] = 0;
624 
625  // The last specified catalogue, to which we should assign all
626  // remaining lines. Index of this one in real_ arrays.
627  Index last_cat = 0;
628 
629  for ( Index i=1; i<n_tg; ++i )
630  {
631  // Is there a catalogue specified?
632  if ( n_cat > i )
633  {
634  // Yes, there is a catalogue.
635 
636  // Has this been specified before?
637  // We use the STL find algorithm to look for the catalogue
638  // name in the real_catalogues. Find returns an iterator, so
639  // to get an index we have to take the difference to
640  // .begin().
641  const Index that_cat = find( real_filenames.begin(),
642  real_filenames.end(),
643  filenames[i] ) - real_filenames.begin();
644  if ( that_cat < real_filenames.nelem() )
645  {
646  // Yes, it has been specified before
647  // ==> Assign to that catalogue
648  real_tgs[that_cat].push_back(i);
649 
650  // Verify, that format, fmin, and fmax are consistent:
651  if ( formats[i] != real_formats[that_cat] ||
652  fmin[i] != real_fmin[that_cat] ||
653  fmax[i] != real_fmax[that_cat] )
654  {
655  ostringstream os;
656  os << "abs_lines_per_speciesReadFromCatalogues: If you specify the\n"
657  << "same catalogue repeatedly, format, fmin, and fmax must be\n"
658  << "consistent. There is an inconsistency between\n"
659  << "catalogue " << that_cat << " and " << i << ".";
660  throw runtime_error(os.str());
661  }
662  }
663  else
664  {
665  // No, it has not been specified before.
666  // ==> Add an entry to real_tgs and the other real_ variables:
667  real_tgs.push_back( MakeArray<Index>(i) );
668 
669  real_filenames.push_back( filenames[i] );
670  real_formats.push_back ( formats[i] );
671  real_fmin.push_back ( fmin[i] );
672  real_fmax.push_back ( fmax[i] );
673 
674  last_cat = i; // assign remainder of lines to this
675  // catalogue, if there is no other catalogue.
676  }
677  }
678  else
679  {
680  // No, there is no catalogue.
681  // ==> Assign to the last catalogue
682  real_tgs[last_cat].push_back(i);
683  }
684  }
685 
686  Index n_real_cat = real_filenames.nelem(); // # real catalogues to read
687 
688  // Some output to low priority stream:
689  out3 << " Catalogues to read and tgs for which these will be used:\n";
690  for ( Index i=0; i<n_real_cat; ++i )
691  {
692  out3 << " " << real_filenames[i] << ":";
693  for ( Index s=0; s<real_tgs[i].nelem(); ++s )
694  out3 << " " << real_tgs[i][s];
695  out3 << "\n";
696  }
697 
698  // Make abs_lines_per_species the right size:
699  abs_lines_per_species.resize( tgs.nelem() );
700 
701  // Loop through the catalogues to read:
702  for ( Index i=0; i<n_real_cat; ++i )
703  {
704  ArrayOfLineRecord abs_lines;
705 
706  // Read catalogue:
707 
708  if ( "HITRAN96"==real_formats[i] )
709  {
710  abs_linesReadFromHitranPre2004( abs_lines, real_filenames[i], real_fmin[i], real_fmax[i], verbosity );
711  }
712  else if ( "HITRAN04"==real_formats[i] )
713  {
714  abs_linesReadFromHitran( abs_lines, real_filenames[i], real_fmin[i], real_fmax[i], verbosity );
715  }
716  else if ( "MYTRAN2"==real_formats[i] )
717  {
718  abs_linesReadFromMytran2( abs_lines, real_filenames[i], real_fmin[i], real_fmax[i], verbosity );
719  }
720  else if ( "JPL"==real_formats[i] )
721  {
722  abs_linesReadFromJpl( abs_lines, real_filenames[i], real_fmin[i], real_fmax[i], verbosity );
723  }
724  else if ( "ARTS"==real_formats[i] )
725  {
726  abs_linesReadFromArts( abs_lines, real_filenames[i], real_fmin[i], real_fmax[i], verbosity );
727  }
728  else
729  {
730  ostringstream os;
731  os << "abs_lines_per_speciesReadFromCatalogues: You specified the\n"
732  << "format `" << real_formats[i] << "', which is unknown.\n"
733  << "Allowd formats are: HITRAN96, HITRAN04, MYTRAN2, JPL, ARTS.";
734  throw runtime_error(os.str());
735  }
736 
737  // We need to make subset tgs for the groups that should
738  // be read from this catalogue.
739  ArrayOfArrayOfSpeciesTag these_tgs(real_tgs[i].nelem());
740  for ( Index s=0; s<real_tgs[i].nelem(); ++s )
741  {
742  these_tgs[s] = tgs[real_tgs[i][s]];
743  }
744 
745  // Create these_abs_lines_per_species:
746  ArrayOfArrayOfLineRecord these_abs_lines_per_species;
747  abs_lines_per_speciesCreateFromLines( these_abs_lines_per_species, abs_lines, these_tgs, verbosity );
748 
749  // Put these lines in the right place in abs_lines_per_species:
750  for ( Index s=0; s<real_tgs[i].nelem(); ++s )
751  {
752  abs_lines_per_species[real_tgs[i][s]] = these_abs_lines_per_species[s];
753  }
754  }
755 }
756 
757 
758 /* Workspace method: Doxygen documentation will be auto-generated */
760  ArrayOfArrayOfLineRecord& abs_lines_per_species,
761  // WS Input:
762  const ArrayOfLineRecord& abs_lines,
763  const ArrayOfArrayOfSpeciesTag& tgs,
764  const Verbosity& verbosity)
765 {
766  CREATE_OUT3;
767 
768  // The species lookup data:
770 
771  // As a safety feature, we will watch out for the case that a
772  // species is included in the calculation, but not all lines are
773  // used. For this we need an array to flag the used species:
774 
775  // For some weird reason, Arrays of bool do not work, although all
776  // other types seem to work fine. So in this single case, I'll use
777  // the stl vector directly. The other place where this is done is in
778  // the function executor in main.cc.
779  // FIXME: Fix this when Array<bool> works.
780  std::vector<bool> species_used (species_data.nelem(),false);
781 
782  // Make abs_lines_per_species the right size:
783  abs_lines_per_species = ArrayOfArrayOfLineRecord(tgs.nelem());
784 
785  // Unfortunately, MTL conatains a bug that leads to all elements of
786  // the outer Array of an Array<Array>> pointing to the same data
787  // after creation. So we need to fix this explicitly:
788  for ( Index i=0; i<abs_lines_per_species.nelem(); ++i )
789  abs_lines_per_species[i] = ArrayOfLineRecord();
790 
791  // Loop all lines in the input line list:
792  for ( Index i=0; i<abs_lines.nelem(); ++i )
793  {
794  // Get a convenient reference to the current line:
795  const LineRecord& this_line = abs_lines[i];
796 
797  // We have to test each tag group in turn (in the order in which
798  // they appear in the controlfile). The line is assigned to the
799  // first tag group that fits.
800 
801  // The flag found is used to break the for loops when the right
802  bool found = false;
803 
804  // We need to define j here, since we need the value outside the
805  // for loop:
806  Index j;
807 
808  // Loop the tag groups:
809  for ( j=0; j<tgs.nelem() && !found ; ++j )
810  {
811  // A tag group can contain several tags:
812  for ( Index k=0; k<tgs[j].nelem() && !found; ++k )
813  {
814  // Get a reference to the current tag (not really
815  // necessary, but makes for nicer notation below):
816  const SpeciesTag& this_tag = tgs[j][k];
817 
818  // Now we will test different attributes of the line
819  // against matching attributes of the tag. If any
820  // attribute does not match, we continue with the next tag
821  // in the tag group. (Exception: Species, see just below.)
822 
823  // Test species. If this attribute does not match we don`t
824  // have to test the other tags in this group, since all
825  // tags must belong to the same species.
826  if ( this_tag.Species() != this_line.Species() ) break;
827 
828  // Test isotopologue. The isotopologue can either match directly, or
829  // the Isotopologue of the tag can be one larger than the
830  // number of isotopologues, which means `all'. Test the second
831  // condition first, since this will probably be more often
832  // used.
833  if ( this_tag.Isotopologue() != this_line.SpeciesData().Isotopologue().nelem() )
834  if ( this_tag.Isotopologue() != this_line.Isotopologue() )
835  continue;
836 
837  // Test frequncy range we take both the lower (Lf) and the
838  // upper (Uf) border frequency to include the `equal' case.
839  // Both Lf and Uf can also be negative, which means `no limit'
840 
841  // Take the lower limit first:
842  if ( this_tag.Lf() >= 0 )
843  if ( this_tag.Lf() > this_line.F() )
844  continue;
845 
846  // Then the upper limit:
847  if ( this_tag.Uf() >= 0 )
848  if ( this_tag.Uf() < this_line.F() )
849  continue;
850 
851  // When we get here, this_tag has survived all tests. That
852  // means it matches the line perfectly!
853  found = true;
854  }
855  }
856 
857  // If a matching tag was found, this line can be used in the
858  // calculation. Add it to the line list for this tag group.
859  if (found)
860  {
861  // We have to use j-1 here, since j was still increased by
862  // one after the matching tag has been found.
863  abs_lines_per_species[j-1].push_back(this_line);
864 
865  // Flag this species as used, if not already done:
866  if ( !species_used[this_line.Species()] )
867  species_used[this_line.Species()] = true;
868  }
869  else
870  {
871  // Safety feature: Issue a warning messages if the lines for a
872  // species are only partly covered by tags.
873  if ( species_used[this_line.Species()] )
874  {
875  CREATE_OUT2;
876  out2 << " Your tags include other lines of species "
877  << this_line.SpeciesData().Name()
878  << ",\n"
879  << " Why do you not include line "
880  << i
881  << " (at "
882  << this_line.F()
883  << " Hz)?\n";
884  }
885  }
886  }
887 
888  // Write some information to the lowest priority output stream.
889  for (Index i=0; i<tgs.nelem(); ++i)
890  {
891  out3 << " " << i << ":";
892 
893  for (Index s=0; s<tgs[i].nelem(); ++s)
894  out3 << " " << tgs[i][s].Name();
895 
896  out3 << ": " << abs_lines_per_species[i].nelem() << " lines\n";
897  }
898 
899 }
900 
902  ArrayOfArrayOfLineRecord& abs_lines_per_species,
903  const Numeric& max_f,
904  const Verbosity&)
905 {
906  // We will simply append the mirror lines after the original
907  // lines. This way we don't have to make a backup copy of the
908  // original lines.
909 
910  for ( Index i=0; i<abs_lines_per_species.nelem(); ++i )
911  {
912  // Get a reference to the current list of lines to save typing:
913  ArrayOfLineRecord& ll = abs_lines_per_species[i];
914 
915  // It is important that we determine the size of ll *before*
916  // we start the loop. After all, we are adding elements. And
917  // we cerainly don't want to continue looping the newly added
918  // elements, we want to loop only the original elements.
919  //
920  const Index n = ll.nelem();
921 
922  // For increased efficiency, reserve the necessary space:
923  //
924  Index nnew = n;
925  //
926  if( max_f >= 0 )
927  {
928  nnew = 0;
929  for ( Index j=0; j<n; ++j )
930  {
931  if( ll[j].F() <= max_f )
932  { nnew += 1; }
933  }
934  }
935  //
936  ll.reserve( n + nnew );
937 
938  // Loop through all lines of this tag group:
939  for ( Index j=0; j<n; ++j )
940  {
941  if( max_f < 0 || ll[j].F() <= max_f )
942  {
943  LineRecord dummy = ll[j];
944  dummy.setF( -dummy.F() );
945  ll.push_back(dummy);
946  }
947  }
948  }
949 }
950 
951 
952 /* Workspace method: Doxygen documentation will be auto-generated */
954  ArrayOfArrayOfLineRecord& abs_lines_per_species,
955  // WS Input:
956  const ArrayOfLineshapeSpec& abs_lineshape,
957  const Vector& f_grid,
958  const Verbosity&)
959 {
960  // Make sure abs_lines_per_species and abs_lineshape have the same
961  // dimension. As a special case, abs_lineshape can have length 1, in
962  // which case it is valid for all species.
963  if ( (abs_lines_per_species.nelem() != abs_lineshape.nelem()) &&
964  (abs_lineshape.nelem() != 1) )
965  {
966  ostringstream os;
967  os << "Dimension of abs_lines_per_species does\n"
968  << "not match that of abs_lineshape.\n"
969  << "(Dimension of abs_lineshape must be 1 or\n"
970  << "equal to abs_lines_per_species.)";
971  throw runtime_error(os.str());
972  }
973 
974  // Make sure that the frequency grid is properly sorted:
975  for ( Index s=0; s<f_grid.nelem()-1; ++s )
976  {
977  if ( f_grid[s+1] <= f_grid[s] )
978  {
979  ostringstream os;
980  os << "The frequency grid f_grid is not properly sorted.\n"
981  << "f_grid[" << s << "] = " << f_grid[s] << "\n"
982  << "f_grid[" << s+1 << "] = " << f_grid[s+1];
983  throw runtime_error(os.str());
984  }
985  }
986 
987  // Cycle through all tag groups:
988  for ( Index i=0; i<abs_lines_per_species.nelem(); ++i )
989  {
990  // Get cutoff frequency of this tag group.
991  // We have to also handle the special case that there is only
992  // one lineshape for all species.
993  Numeric cutoff;
994  if (1==abs_lineshape.nelem())
995  cutoff = abs_lineshape[0].Cutoff();
996  else
997  cutoff = abs_lineshape[i].Cutoff();
998 
999  // Check whether cutoff is defined:
1000  if ( cutoff != -1)
1001  {
1002  // Get a reference to the current list of lines to save typing:
1003  ArrayOfLineRecord& ll = abs_lines_per_species[i];
1004 
1005  // Calculate the borders:
1006  Numeric upp = f_grid[f_grid.nelem()-1] + cutoff;
1007  Numeric low = f_grid[0] - cutoff;
1008 
1009  // Cycle through all lines within this tag group.
1011  for ( ArrayOfLineRecord::iterator j=ll.begin(); j<ll.end(); ++j )
1012  {
1013  // Center frequency:
1014  const Numeric F0 = j->F();
1015 
1016  if ( ( F0 >= low) && ( F0 <= upp) )
1017  {
1018  // Build list of elements which should be kept
1019  keep.push_back (j);
1020  // The original implementation just erased the non-wanted
1021  // elements from the array. Problem: On every erase the
1022  // whole array will be copied which actually kills
1023  // performance.
1024  }
1025  }
1026 
1027  // If next comparison is false, some elements have to be removed
1028  if (keep.nelem () != ll.nelem ())
1029  {
1030  ArrayOfLineRecord nll;
1031  // Copy all elements that should be kept to a new Array
1033  = keep.begin(); j != keep.end(); j++)
1034  {
1035  nll.push_back (**j);
1036  }
1037  // Overwrite the old array with the new one
1038  ll.resize (nll.nelem ());
1039  ll = nll;
1040  }
1041  }
1042  }
1043 }
1044 
1045 
1046 /* Workspace method: Doxygen documentation will be auto-generated */
1049  Index& propmat_clearsky_agenda_checked,
1050  Index& abs_xsec_agenda_checked,
1051  // Control Parameters:
1052  const String& basename,
1053  const Verbosity& verbosity)
1054 {
1055  CREATE_OUT2;
1056 
1057  // Invalidate agenda check flags
1058  propmat_clearsky_agenda_checked = false;
1059  abs_xsec_agenda_checked = false;
1060 
1061  // Species lookup data:
1063 
1064  // We want to make lists of included and excluded species:
1065  ArrayOfString included(0), excluded(0);
1066 
1067  tgs.resize(0);
1068 
1069  for ( Index i=0; i<species_data.nelem(); ++i )
1070  {
1071  const String specname = species_data[i].Name();
1072 
1073  String filename = basename;
1074  if (basename.length() && basename[basename.length()-1] != '/')
1075  filename += ".";
1076  filename += specname;
1077 
1078  try {
1079  find_xml_file(filename, verbosity);
1080  // Add to included list:
1081  included.push_back(specname);
1082 
1083  // Convert name of species to a SpeciesTag object:
1084  SpeciesTag this_tag(specname);
1085 
1086  // Create Array of SpeciesTags with length 1
1087  // (our tag group has only one tag):
1088  ArrayOfSpeciesTag this_group(1);
1089  this_group[0] = this_tag;
1090 
1091  // Add this tag group to tgs:
1092  tgs.push_back(this_group);
1093  }
1094  catch (runtime_error e)
1095  {
1096  // The file for the species could not be found.
1097  excluded.push_back(specname);
1098  }
1099  }
1100 
1101  // Some nice output:
1102  out2 << " Included Species (" << included.nelem() << "):\n";
1103  for ( Index i=0; i<included.nelem(); ++i )
1104  out2 << " " << included[i] << "\n";
1105 
1106  out2 << " Excluded Species (" << excluded.nelem() << "):\n";
1107  for ( Index i=0; i<excluded.nelem(); ++i )
1108  out2 << " " << excluded[i] << "\n";
1109 }
1110 
1111 
1112 /* Workspace method: Doxygen documentation will be auto-generated */
1113 void abs_lineshapeDefine(// WS Output:
1114  ArrayOfLineshapeSpec& abs_lineshape,
1115  // WS Input:
1116  const String& shape,
1117  const String& normalizationfactor,
1118  const Numeric& cutoff,
1119  const Verbosity& verbosity)
1120 {
1121  CREATE_OUT2;
1122 
1123  // Make lineshape and normalization factor data visible:
1126 
1127 
1128  // generate the right number of elements
1129  // Index tag_sz = tgs.nelem();
1130  // We generate only 1 copy of the lineshape settings. Absorption
1131  // routines check for this case and use it for all species.
1132  Index tag_sz = 1;
1133  abs_lineshape.resize(tag_sz);
1134 
1135  // Is this lineshape available?
1136  Index found0=-1;
1137  for ( Index i=0; i<lineshape_data.nelem() && (found0 == -1) ; ++i )
1138  {
1139  const String& str = lineshape_data[i].Name();
1140  if (str == shape)
1141  {
1142  out2 << " Selected lineshape: " << str << "\n";
1143  found0=i;
1144  }
1145  }
1146 
1147  // Is this normalization to the lineshape available?
1148  Index found1=-1;
1149  for ( Index i=0; i<lineshape_norm_data.nelem() && (found1 == -1); ++i )
1150  {
1151  const String& str = lineshape_norm_data[i].Name();
1152  if (str == normalizationfactor)
1153  {
1154  out2 << " Selected normalization factor : " << normalizationfactor << "\n";
1155 
1156  if ( (cutoff != -1) && (cutoff < 0.0) )
1157  throw runtime_error(" Cutoff must be -1 or positive.");
1158  out2 << " Selected cutoff frequency [Hz] : " << cutoff << "\n";
1159  found1=i;
1160  }
1161  }
1162 
1163 
1164  // did we find the lineshape and normalization factor?
1165  if (found0 == -1)
1166  throw runtime_error("Selected lineshape not available.");
1167  if (found1 == -1)
1168  throw runtime_error("Selected normalization to lineshape not available.");
1169 
1170 
1171  // now set the lineshape
1172  for (Index i=0; i<tag_sz; i++)
1173  {
1174  abs_lineshape[i].SetInd_ls( found0 );
1175  abs_lineshape[i].SetInd_lsn( found1 );
1176  abs_lineshape[i].SetCutoff( cutoff );
1177  }
1178 }
1179 
1180 
1181 /* Workspace method: Doxygen documentation will be auto-generated */
1182 void abs_lineshape_per_tgDefine(// WS Output:
1183  ArrayOfLineshapeSpec& abs_lineshape,
1184  // WS Input:
1185  const ArrayOfArrayOfSpeciesTag& tgs,
1186  const ArrayOfString& shape,
1187  const ArrayOfString& normalizationfactor,
1188  const Vector& cutoff,
1189  const Verbosity& verbosity)
1190 {
1191  CREATE_OUT2;
1192 
1193  // Make lineshape and normalization factor data visible:
1196 
1197  // check that the number of elements are equal
1198  Index tg_sz = tgs.nelem();
1199  if ( (tg_sz != shape.nelem()) ||
1200  (tg_sz != normalizationfactor.nelem()) ||
1201  (tg_sz != cutoff.nelem()) )
1202  {
1203  ostringstream os;
1204  os << "abs_lineshape_per_tgDefine: number of elements does\n"
1205  << "not match the number of tag groups defined.";
1206  throw runtime_error(os.str());
1207  }
1208 
1209 
1210  // generate the right number of elements
1211  abs_lineshape.resize(tg_sz);
1212 
1213  // Is this lineshape available?
1214  for (Index k=0; k<tg_sz; ++k)
1215  {
1216  Index found0=-1;
1217  for ( Index i=0; i<lineshape_data.nelem() && (found0 == -1); ++i )
1218  {
1219  const String& str = lineshape_data[i].Name();
1220  if (str == shape[k])
1221  {
1222  out2 << " Tag Group: [";
1223  for (Index s=0; s<tgs[k].nelem()-1; ++s)
1224  out2 << tgs[k][s].Name() << ", ";
1225  out2 << tgs[k][tgs[k].nelem()-1].Name() << "]\n";
1226  out2 << " Selected lineshape: " << str << "\n";
1227  found0=i;
1228  }
1229  }
1230 
1231  // Is this normalization to the lineshape available?
1232  Index found1=-1;
1233  for ( Index i=0; i<lineshape_norm_data.nelem() && (found1 == -1); ++i )
1234  {
1235  const String& str = lineshape_norm_data[i].Name();
1236  if (str == normalizationfactor[k])
1237  {
1238  out2 << " Selected normalization factor: " << normalizationfactor[k] << "\n";
1239  if ( (cutoff[k] != -1) && (cutoff[k] < 0.0) )
1240  throw runtime_error(" Cutoff must be -1 or positive.");
1241  out2 << " Selected cutoff frequency : " << cutoff[k] << "\n";
1242  found1=i;
1243  }
1244  }
1245 
1246 
1247  // did we find the lineshape and normalization factor?
1248  if (found0 == -1)
1249  {
1250  ostringstream os;
1251  os << "Selected lineshape not available: "<< shape[k] <<"\n";
1252  throw runtime_error(os.str());
1253  }
1254  if (found1 == -1)
1255  {
1256  ostringstream os;
1257  os << "Selected normalization to lineshape not available: "<<
1258  normalizationfactor[k] <<"\n";
1259  throw runtime_error(os.str());
1260  }
1261 
1262  // now set the lineshape variables
1263  abs_lineshape[k].SetInd_ls( found0 );
1264  abs_lineshape[k].SetInd_lsn( found1 );
1265  abs_lineshape[k].SetCutoff( cutoff[k] );
1266  }
1267 }
1268 
1269 
1271 
1287 void abs_h2oSet(Vector& abs_h2o,
1288  const ArrayOfArrayOfSpeciesTag& abs_species,
1289  const Matrix& abs_vmrs,
1290  const Verbosity&)
1291 {
1292  const Index h2o_index
1293  = find_first_species_tg( abs_species,
1295 
1296  abs_h2o.resize( abs_vmrs.ncols() );
1297  if ( h2o_index < 0 )
1298  abs_h2o = -99;
1299  else
1300  abs_h2o = abs_vmrs(h2o_index,Range(joker));
1301 }
1302 
1303 
1305 
1315 void abs_n2Set(Vector& abs_n2,
1316  const ArrayOfArrayOfSpeciesTag& abs_species,
1317  const Matrix& abs_vmrs,
1318  const Verbosity&)
1319 {
1320  const Index n2_index
1321  = find_first_species_tg( abs_species,
1323 
1324  abs_n2.resize( abs_vmrs.ncols() );
1325  if ( n2_index < 0 )
1326  abs_n2 = -99;
1327  else
1328  abs_n2 = abs_vmrs(n2_index,Range(joker));
1329 }
1330 
1332 
1342 void abs_o2Set(Vector& abs_o2,
1343  const ArrayOfArrayOfSpeciesTag& abs_species,
1344  const Matrix& abs_vmrs,
1345  const Verbosity&)
1346 {
1347  const Index o2_index
1348  = find_first_species_tg( abs_species,
1350 
1351  abs_o2.resize( abs_vmrs.ncols() );
1352  if ( o2_index < 0 )
1353  abs_o2 = -99;
1354  else
1355  abs_o2 = abs_vmrs(o2_index,Range(joker));
1356 }
1357 
1358 
1359 /* Workspace method: Doxygen documentation will be auto-generated */
1360 void AbsInputFromAtmFields (// WS Output:
1361  Vector& abs_p,
1362  Vector& abs_t,
1363  Matrix& abs_vmrs,
1364  // WS Input:
1365  const Index& atmosphere_dim,
1366  const Vector& p_grid,
1367  const Tensor3& t_field,
1368  const Tensor4& vmr_field,
1369  const Verbosity&
1370  )
1371 {
1372  // First, make sure that we really have a 1D atmosphere:
1373  if ( 1 != atmosphere_dim )
1374  {
1375  ostringstream os;
1376  os << "Atmospheric dimension must be 1D, but atmosphere_dim is "
1377  << atmosphere_dim << ".";
1378  throw runtime_error(os.str());
1379  }
1380 
1381  abs_p = p_grid;
1382  abs_t = t_field (joker, 0, 0);
1383  abs_vmrs = vmr_field (joker, joker, 0, 0);
1384 }
1385 
1386 
1387 /* Workspace method: Doxygen documentation will be auto-generated */
1388 void abs_coefCalcFromXsec(// WS Output:
1389  Matrix& abs_coef,
1390  ArrayOfMatrix& abs_coef_per_species,
1391  // WS Input:
1392  const ArrayOfMatrix& abs_xsec_per_species,
1393  const Matrix& abs_vmrs,
1394  const Vector& abs_p,
1395  const Vector& abs_t,
1396  const Verbosity& verbosity)
1397 {
1398  CREATE_OUT3;
1399 
1400  // Check that abs_vmrs and abs_xsec_per_species really have compatible
1401  // dimensions. In abs_vmrs there should be one row for each tg:
1402  if ( abs_vmrs.nrows() != abs_xsec_per_species.nelem() )
1403  {
1404  ostringstream os;
1405  os << "Variable abs_vmrs must have compatible dimension to abs_xsec_per_species.\n"
1406  << "abs_vmrs.nrows() = " << abs_vmrs.nrows() << "\n"
1407  << "abs_xsec_per_species.nelem() = " << abs_xsec_per_species.nelem();
1408  throw runtime_error(os.str());
1409  }
1410 
1411  // Check that number of altitudes are compatible. We only check the
1412  // first element, this is possilble because within arts all elements
1413  // are on the same altitude grid.
1414  if ( abs_vmrs.ncols() != abs_xsec_per_species[0].ncols() )
1415  {
1416  ostringstream os;
1417  os << "Variable abs_vmrs must have same numbers of altitudes as abs_xsec_per_species.\n"
1418  << "abs_vmrs.ncols() = " << abs_vmrs.ncols() << "\n"
1419  << "abs_xsec_per_species[0].ncols() = " << abs_xsec_per_species[0].ncols();
1420  throw runtime_error(os.str());
1421  }
1422 
1423  // Check dimensions of abs_p and abs_t:
1424  chk_size("abs_p", abs_p, abs_vmrs.ncols());
1425  chk_size("abs_t", abs_t, abs_vmrs.ncols());
1426 
1427 
1428  // Initialize abs_coef and abs_coef_per_species. The array dimension of abs_coef_per_species
1429  // is the same as that of abs_xsec_per_species. The dimension of abs_coef should
1430  // be equal to one of the abs_xsec_per_species enries.
1431  abs_coef.resize( abs_xsec_per_species[0].nrows(), abs_xsec_per_species[0].ncols() );
1432  abs_coef = 0; // Matpack can set all elements like this.
1433 
1434  abs_coef_per_species.resize( abs_xsec_per_species.nelem() );
1435 
1436  out3 << " Computing abs_coef and abs_coef_per_species from abs_xsec_per_species.\n";
1437  // Loop through all tag groups
1438  for ( Index i=0; i<abs_xsec_per_species.nelem(); ++i )
1439  {
1440  out3 << " Tag group " << i << "\n";
1441 
1442  // Make this element of abs_xsec_per_species the right size:
1443  abs_coef_per_species[i].resize( abs_xsec_per_species[i].nrows(), abs_xsec_per_species[i].ncols() );
1444  abs_coef_per_species[i] = 0; // Initialize all elements to 0.
1445 
1446  // Loop through all altitudes
1447  for ( Index j=0; j<abs_xsec_per_species[i].ncols(); j++)
1448  {
1449  // Calculate total number density from pressure and temperature.
1450  const Numeric n = number_density(abs_p[j],abs_t[j]);
1451 
1452  // Loop through all frequencies
1453  for ( Index k=0; k<abs_xsec_per_species[i].nrows(); k++)
1454  {
1455  abs_coef_per_species[i](k,j) = abs_xsec_per_species[i](k,j) * n * abs_vmrs(i,j);
1456  }
1457  }
1458 
1459  // Add up to the total absorption:
1460  abs_coef += abs_coef_per_species[i]; // In Matpack you can use the +=
1461  // operator to do elementwise addition.
1462  }
1463 }
1464 
1465 
1466 
1467 /* Workspace method: Doxygen documentation will be auto-generated */
1468 void abs_xsec_per_speciesInit(// WS Output:
1469  ArrayOfMatrix& abs_xsec_per_species,
1470  // WS Input:
1471  const ArrayOfArrayOfSpeciesTag& tgs,
1472  const ArrayOfIndex& abs_species_active,
1473  const Vector& f_grid,
1474  const Vector& abs_p,
1475  const Index& abs_xsec_agenda_checked,
1476  const Verbosity& verbosity
1477  )
1478 {
1479  CREATE_OUT3;
1480 
1481  if (!abs_xsec_agenda_checked)
1482  throw runtime_error("You must call *abs_xsec_agenda_checkedCalc* before calling this method.");
1483 
1484  // We need to check that abs_species_active doesn't have more elements than
1485  // abs_species (abs_xsec_agenda_checkedCalc doesn't know abs_species_active.
1486  // Usually we come here through an agenda call, where abs_species_active has
1487  // been properly created somewhere internally. But we might get here by
1488  // direct call, and then need to be safe!).
1489  if ( tgs.nelem() < abs_species_active.nelem() )
1490  {
1491  ostringstream os;
1492  os << "abs_species_active (n=" << abs_species_active.nelem()
1493  << ") not allowed to have more elements than abs_species (n="
1494  << tgs.nelem() << ")!\n";
1495  throw runtime_error(os.str());
1496  }
1497 
1498  // Initialize abs_xsec_per_species. The array dimension of abs_xsec_per_species
1499  // is the same as that of abs_lines_per_species.
1500  abs_xsec_per_species.resize( tgs.nelem() );
1501 
1502  // Loop abs_xsec_per_species and make each matrix the right size,
1503  // initializing to zero.
1504  // But skip inactive species, loop only over the active ones.
1505  for ( Index ii=0; ii<abs_species_active.nelem(); ++ii )
1506  {
1507  const Index i = abs_species_active[ii];
1508  // Check that abs_species_active index is not higher than the number
1509  // of species
1510  if (i >= tgs.nelem())
1511  {
1512  ostringstream os;
1513  os << "*abs_species_active* contains an invalid species index.\n"
1514  << "Species index must be between 0 and " << tgs.nelem()-1;
1515  throw std::runtime_error(os.str());
1516  }
1517  // Make this element of abs_xsec_per_species the right size:
1518  abs_xsec_per_species[i].resize( f_grid.nelem(), abs_p.nelem() );
1519  abs_xsec_per_species[i] = 0; // Matpack can set all elements like this.
1520  }
1521 
1522  ostringstream os;
1523  os << " Initialized abs_xsec_per_species.\n"
1524  << " Number of frequencies : " << f_grid.nelem() << "\n"
1525  << " Number of pressure levels : " << abs_p.nelem() << "\n";
1526  out3 << os.str();
1527 }
1528 
1529 
1530 /* Workspace method: Doxygen documentation will be auto-generated */
1532  ArrayOfMatrix& abs_xsec_per_species,
1533  // WS Input:
1534  const ArrayOfArrayOfSpeciesTag& tgs,
1535  const ArrayOfIndex& abs_species_active,
1536  const Vector& f_grid,
1537  const Vector& abs_p,
1538  const Vector& abs_t,
1539  const Matrix& abs_vmrs,
1540  const ArrayOfArrayOfLineRecord& abs_lines_per_species,
1541  const ArrayOfLineshapeSpec& abs_lineshape,
1542  const SpeciesAuxData& isotopologue_ratios,
1543  const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
1544  const ArrayOfArrayOfIndex& line_mixing_data_lut,
1545  const Verbosity& verbosity)
1546 {
1547  CREATE_OUT3;
1548 
1549  // Check that correct isotopologue ratios are defined for the species
1550  // we want to calculate
1551  checkIsotopologueRatios(tgs, isotopologue_ratios);
1552 
1553  // Check that all temperatures are at least 0 K. (Negative Kelvin
1554  // temperatures are unphysical.)
1555  if ( min(abs_t) < 0 )
1556  {
1557  ostringstream os;
1558  os << "Temperature must be at least 0 K. But you request an absorption\n"
1559  << "calculation at " << min(abs_t) << " K!";
1560  throw runtime_error(os.str());
1561  }
1562 
1563  // Check that all parameters that should have the number of tag
1564  // groups as a dimension are consistent.
1565  {
1566  const Index n_tgs = tgs.nelem();
1567  const Index n_xsec = abs_xsec_per_species.nelem();
1568  const Index n_vmrs = abs_vmrs.nrows();
1569  const Index n_lines = abs_lines_per_species.nelem();
1570  const Index n_shapes = abs_lineshape.nelem();
1571 
1572  if ( n_tgs != n_xsec ||
1573  n_tgs != n_vmrs ||
1574  n_tgs != n_lines ||
1575  ( n_tgs != n_shapes &&
1576  1 != n_shapes ) )
1577  {
1578  ostringstream os;
1579  os << "The following variables must all have the same dimension:\n"
1580  << "tgs: " << tgs.nelem() << "\n"
1581  << "abs_xsec_per_species: " << abs_xsec_per_species.nelem() << "\n"
1582  << "abs_vmrs: " << abs_vmrs.nrows() << "\n"
1583  << "abs_lines_per_species: " << abs_lines_per_species.nelem() << "\n"
1584  << "abs_lineshape: " << abs_lineshape.nelem() << "\n"
1585  << "(As a special case, abs_lineshape is allowed to have only one element.)";
1586  throw runtime_error(os.str());
1587  }
1588  }
1589 
1590  // Print information:
1591  //
1592  out3 << " Calculating line spectra.\n";
1593  //
1594  // Uncomment the part below if you temporarily want detailed info about
1595  // transitions to be done
1596 
1597  // The variables defined here (in particular the frequency
1598  // conversion) are just to make the output nice. They are not used
1599  // in subsequent calculations.
1600  // cout << " Transitions to do: \n";
1601  // Index nlines = 0;
1602  // String funit;
1603  // Numeric ffac;
1604  // if ( f_grid[0] < 3e12 )
1605  // {
1606  // funit = "GHz"; ffac = 1e9;
1607  // }
1608  // else
1609  // {
1610  // extern const Numeric SPEED_OF_LIGHT;
1611  // funit = "cm-1"; ffac = SPEED_OF_LIGHT*100;
1612  // }
1613  // for ( Index i=0; i<abs_lines_per_species.nelem(); ++i )
1614  // {
1615  // for ( Index l=0; l<abs_lines_per_species[i].nelem(); ++l )
1616  // {
1617  // cout << " " << abs_lines_per_species[i][l].Name() << " @ "
1618  // << abs_lines_per_species[i][l].F()/ffac << " " << funit << " ("
1619  // << abs_lines_per_species[i][l].I0() << " "
1620  // << abs_lines_per_species[i][l].Agam() << ")\n";
1621  // nlines++;
1622  // }
1623  // }
1624  // out2 << " Total number of transistions : " << nlines << "\n";
1625 
1626  // Call xsec_species for each tag group.
1627  for ( Index ii=0; ii<abs_species_active.nelem(); ++ii )
1628  {
1629  const Index i = abs_species_active[ii];
1630 
1631  // Get a pointer to the line list for the current species. This
1632  // is just so that we don't have to type abs_lines_per_species[i] over
1633  // and over again.
1634  const ArrayOfLineRecord& ll = abs_lines_per_species[i];
1635 
1636  // Also get a pointer to the lineshape specification. This
1637  // requires special treatment: If there is only 1 lineshape
1638  // given, the same line shape should be used for all species.
1639  LineshapeSpec ls;
1640  if (1==abs_lineshape.nelem())
1641  ls = abs_lineshape[0];
1642  else
1643  ls = abs_lineshape[i];
1644 
1645 
1646  // We do the LBL calculation only if:
1647  // - The line list is not empty, and
1648  // - The species is not a Zeeman species.
1649  if ( 0 < ll.nelem() && tgs[i].nelem() && !is_zeeman(tgs[i]) )
1650  {
1651  // As a safety check, check that the species of the first
1652  // line matches the species we should have according to
1653  // tgs. (This in case the order in tgs has been changed and
1654  // abs_lines_per_species has not been changed consistently.)
1655  if (ll[0].Species() != tgs[i][0].Species() )
1656  {
1657  ostringstream os;
1658  os << "The species in the line list does not match the species\n"
1659  << "for which you want to calculate absorption:\n"
1660  << "abs_species: " << get_tag_group_name(tgs[i]) << "\n"
1661  << "abs_lines_per_species: " << ll[0].Name();
1662  throw runtime_error(os.str());
1663  }
1664 
1665  // Get the name of the species. The member function name of a
1666  // LineRecord returns the full name (species + isotopologue). So
1667  // it is for us more convenient to get the species index
1668  // from the LineRecord member function Species(), and then
1669  // use this to look up the species name in species_data.
1671  String species_name = species_data[ll[0].Species()].Name();
1672 
1673  // Get the name of the lineshape. For that we use the member
1674  // function Ind_ls() to the lineshape data ls, which returns
1675  // an index. With that index we can go into lineshape_data
1676  // to get the name.
1677  // We will need this for safety checks later on.
1679  String lineshape_name = lineshape_data[ls.Ind_ls()].Name();
1680 
1681 
1682  // The use of overlap parameters for oxygen lines only makes
1683  // sense if the special Rosenkranz lineshape is used
1684  // (Rosenkranz_Voigt_Drayson or Rosenkranz_Voigt_Kuntz6).
1685  // Likewise, the use of these lineshapes only makes sense if
1686  // overlap parameters are available. We will test both these
1687  // conditions.
1688 
1689  // First of all, let's find out if the species we are
1690  // dealing with is oxygen.
1691  if ( "O2" == species_name )
1692  {
1693  // Do we have overlap parameters in the aux fields of
1694  // the LineRecord?
1695  if ( 0 != ll[0].Naux() )
1696  {
1697  // Yes. So let's make sure that the lineshape is one
1698  // that can use these parameters.
1699  if ( "Rosenkranz_Voigt_Drayson" != lineshape_name &&
1700  "Rosenkranz_Voigt_Kuntz6" != lineshape_name )
1701  {
1702  ostringstream os;
1703  os
1704  << "You are using a line catalogue that contains auxiliary parameters to\n"
1705  << "take care of overlap for oxygen lines. But you are not using a\n"
1706  << "lineshape that uses these parameters. Use Rosenkranz_Voigt_Drayson or\n"
1707  << "Rosenkranz_Voigt_Kuntz6.";
1708  throw runtime_error(os.str());
1709  }
1710  }
1711  }
1712 
1713  // Now we go the opposite way. Let's see if the Rosenkranz
1714  // lineshapes are used.
1715  if ( "Rosenkranz_Voigt_Drayson" == lineshape_name ||
1716  "Rosenkranz_Voigt_Kuntz6" == lineshape_name )
1717  {
1718  // Is the species oxygen, as it should be?
1719  if ( "O2" != species_name )
1720  {
1721  ostringstream os;
1722  os
1723  << "You are trying to use one of the special Rosenkranz lineshapes with\n"
1724  << "overlap correction for a species other than oxygen. Your species is\n"
1725  << species_name << ".\n"
1726  << "Select another lineshape for this species.";
1727  throw runtime_error(os.str());
1728  }
1729  else
1730  {
1731  // Do we have overlap parameters in the aux fields of
1732  // the LineRecord?
1733  if ( 0 == ll[0].Naux() )
1734  {
1735  ostringstream os;
1736  os
1737  << "You are trying to use one of the special Rosenkranz lineshapes with\n"
1738  << "overlap correction. But your line file does not contain aux\n"
1739  << "parameters. (I've checked only the first LineRecord). Use a line file\n"
1740  << "with overlap parameters.";
1741  throw runtime_error(os.str());
1742  }
1743  }
1744  }
1745  if (tgs[i][0].LineMixingType() == SpeciesTag::LINE_MIXING_TYPE_NONE)
1746  {
1747  Matrix dummy_phase(abs_xsec_per_species[i].nrows(),
1748  abs_xsec_per_species[i].ncols(),
1749  0.);
1750  xsec_species(abs_xsec_per_species[i],
1751  dummy_phase,
1752  f_grid,
1753  abs_p,
1754  abs_t,
1755  abs_vmrs,
1756  tgs,
1757  i,
1758  ll,
1759  ls.Ind_ls(),
1760  ls.Ind_lsn(),
1761  ls.Cutoff(),
1762  isotopologue_ratios,
1763  verbosity );
1764  }
1765  else
1766  {
1767  Matrix dummy_phase(abs_xsec_per_species[i].nrows(),
1768  abs_xsec_per_species[i].ncols(),
1769  0.);
1770  xsec_species_line_mixing_wrapper(abs_xsec_per_species[i],
1771  dummy_phase,
1772  line_mixing_data,
1773  line_mixing_data_lut,
1774  f_grid,
1775  abs_p,
1776  abs_t,
1777  abs_vmrs,
1778  tgs,
1779  i,
1780  ll,
1781  ls.Ind_ls(),
1782  ls.Ind_lsn(),
1783  ls.Cutoff(),
1784  isotopologue_ratios,
1785  verbosity);
1786  }
1787  // Note that we call xsec_species with a row of abs_vmrs,
1788  // selected by the above Matpack expression. This is
1789  // possible, because xsec_species is using Views.
1790 
1791  }
1792 
1793  if (out3.sufficient_priority())
1794  {
1795  ostringstream os;
1796  os << " Tag group " << i
1797  << " (" << get_tag_group_name(tgs[i]) << "): "
1798  << ll.nelem() << " transitions\n";
1799  out3 << os.str();
1800  }
1801 
1802  } // End of species for loop.
1803 }
1804 
1805 /* Workspace method: Doxygen documentation will be auto-generated */
1807  ArrayOfMatrix& abs_xsec_per_species,
1808  // WS Input:
1809  const ArrayOfArrayOfSpeciesTag& tgs,
1810  const ArrayOfIndex& abs_species_active,
1811  const Vector& f_grid,
1812  const Vector& abs_p,
1813  const Vector& abs_t,
1814  const Matrix& abs_vmrs,
1815  const ArrayOfString& abs_cont_names,
1816  const ArrayOfVector& abs_cont_parameters,
1817  const ArrayOfString& abs_cont_models,
1818  const Verbosity& verbosity)
1819 {
1820  CREATE_OUT3;
1821 
1822  // Needed for some continua, and set here from abs_vmrs:
1823  Vector abs_h2o, abs_n2, abs_o2;
1824 
1825  // Check that all paramters that should have the number of tag
1826  // groups as a dimension are consistent.
1827  {
1828  const Index n_tgs = tgs.nelem();
1829  const Index n_xsec = abs_xsec_per_species.nelem();
1830  const Index n_vmrs = abs_vmrs.nrows();
1831 
1832  if ( n_tgs != n_xsec || n_tgs != n_vmrs )
1833  {
1834  ostringstream os;
1835  os << "The following variables must all have the same dimension:\n"
1836  << "tgs: " << tgs.nelem() << "\n"
1837  << "abs_xsec_per_species: " << abs_xsec_per_species.nelem() << "\n"
1838  << "abs_vmrs.nrows(): " << abs_vmrs.nrows();
1839  throw runtime_error(os.str());
1840  }
1841  }
1842 
1843  // Check, that dimensions of abs_cont_names and
1844  // abs_cont_parameters are consistent...
1845  if ( abs_cont_names.nelem() !=
1846  abs_cont_parameters.nelem() )
1847  {
1848  ostringstream os;
1849 
1850  for (Index i=0; i < abs_cont_names.nelem(); ++i)
1851  os << "abs_xsec_per_speciesAddConts: " << i << " name : " << abs_cont_names[i] << "\n";
1852 
1853  for (Index i=0; i < abs_cont_parameters.nelem(); ++i)
1854  os << "abs_xsec_per_speciesAddConts: " << i << " param: " << abs_cont_parameters[i] << "\n";
1855 
1856  for (Index i=0; i < abs_cont_models.nelem(); ++i)
1857  os << "abs_xsec_per_speciesAddConts: " << i << " option: " << abs_cont_models[i] << "\n";
1858 
1859  os << "The following variables must have the same dimension:\n"
1860  << "abs_cont_names: " << abs_cont_names.nelem() << "\n"
1861  << "abs_cont_parameters: " << abs_cont_parameters.nelem();
1862 
1863  throw runtime_error(os.str());
1864  }
1865 
1866  // Check that abs_p, abs_t, and abs_vmrs have the same
1867  // dimension. This could be a user error, so we throw a
1868  // runtime_error.
1869 
1870  if ( abs_t.nelem() != abs_p.nelem() )
1871  {
1872  ostringstream os;
1873  os << "Variable abs_t must have the same dimension as abs_p.\n"
1874  << "abs_t.nelem() = " << abs_t.nelem() << '\n'
1875  << "abs_p.nelem() = " << abs_p.nelem();
1876  throw runtime_error(os.str());
1877  }
1878 
1879  if ( abs_vmrs.ncols() != abs_p.nelem() )
1880  {
1881  ostringstream os;
1882  os << "Variable dimension abs_vmrs.ncols() must\n"
1883  << "be the same as abs_p.nelem().\n"
1884  << "abs_vmrs.ncols() = " << abs_vmrs.ncols() << '\n'
1885  << "abs_p.nelem() = " << abs_p.nelem();
1886  throw runtime_error(os.str());
1887  }
1888 
1889  // We set abs_h2o, abs_n2, and abs_o2 later, because we only want to
1890  // do it if the parameters are really needed.
1891 
1892 
1893  out3 << " Calculating continuum spectra.\n";
1894 
1895  // Loop tag groups:
1896  for ( Index ii=0; ii<abs_species_active.nelem(); ++ii )
1897  {
1898  const Index i = abs_species_active[ii];
1899 
1901 
1902  // Go through the tags in the current tag group to see if they
1903  // are continuum tags:
1904  for ( Index s=0; s<tgs[i].nelem(); ++s )
1905  {
1906  // Continuum tags in the sense that we talk about here
1907  // (including complete absorption models) are marked by a special type.
1908  if (tgs[i][s].Type() == SpeciesTag::TYPE_PREDEF)
1909  {
1910  // We have identified a continuum tag!
1911 
1912  // Get only the continuum name. The full tag name is something like:
1913  // H2O-HITRAN96Self-*-*. We want only the `H2O-HITRAN96Self' part:
1914  const String name =
1915  species_data[tgs[i][s].Species()].Name() + "-"
1916  + species_data[tgs[i][s].Species()].Isotopologue()[tgs[i][s].Isotopologue()].Name();
1917 
1918  // Check, if we have parameters for this model. For
1919  // this, the model name must be listed in
1920  // abs_cont_names.
1921  const Index n =
1922  find( abs_cont_names.begin(),
1923  abs_cont_names.end(),
1924  name ) - abs_cont_names.begin();
1925 
1926  // n==abs_cont_names.nelem() indicates that
1927  // the name was not found.
1928  if ( n==abs_cont_names.nelem() )
1929  {
1930  ostringstream os;
1931  os << "Cannot find model " << name
1932  << " in abs_cont_names.";
1933  throw runtime_error(os.str());
1934  }
1935 
1936  // Ok, the tag specifies a valid continuum model and
1937  // we have continuum parameters.
1938 
1939  if (out3.sufficient_priority())
1940  {
1941  ostringstream os;
1942  os << " Adding " << name
1943  << " to tag group " << i << ".\n";
1944  out3 << os.str();
1945  }
1946 
1947  // find the options for this continuum tag from the input array
1948  // of options. The actual field of the array is n:
1949  const String ContOption = abs_cont_models[n];
1950 
1951  // Set abs_h2o, abs_n2, and abs_o2 from the first matching species.
1952  abs_h2oSet(abs_h2o, tgs, abs_vmrs, verbosity);
1953  abs_n2Set(abs_n2, tgs, abs_vmrs, verbosity);
1954  abs_o2Set(abs_o2, tgs, abs_vmrs, verbosity);
1955 
1956  // Add the continuum for this tag. The parameters in
1957  // this call should be clear. The vmr is in
1958  // abs_vmrs(i,Range(joker)). The other vmr variables,
1959  // abs_h2o, abs_n2, and abs_o2 contains the real vmr of H2O,
1960  // N2, nad O2, which are needed as additional information for
1961  // certain continua:
1962  // abs_h2o for
1963  // O2-PWR88, O2-PWR93, O2-PWR98,
1964  // O2-MPM85, O2-MPM87, O2-MPM89, O2-MPM92, O2-MPM93,
1965  // O2-TRE05,
1966  // O2-SelfContStandardType, O2-SelfContMPM93, O2-SelfContPWR93,
1967  // N2-SelfContMPM93, N2-DryContATM01,
1968  // N2-CIArotCKDMT250, N2-CIAfunCKDMT250
1969  // abs_n2 for
1970  // H2O-SelfContCKD24, H2O-ForeignContCKD24,
1971  // O2-v0v0CKDMT100,
1972  // CO2-ForeignContPWR93, CO2-ForeignContHo66
1973  // abs_o2 for
1974  // N2-CIArotCKDMT250, N2-CIAfunCKDMT250
1975  xsec_continuum_tag( abs_xsec_per_species[i],
1976  name,
1977  abs_cont_parameters[n],
1978  abs_cont_models[n],
1979  f_grid,
1980  abs_p,
1981  abs_t,
1982  abs_n2,
1983  abs_h2o,
1984  abs_o2,
1985  abs_vmrs(i,Range(joker)),
1986  verbosity );
1987  // Calling this function with a row of Matrix abs_vmrs
1988  // is possible because it uses Views.
1989  }
1990  }
1991  }
1992 
1993 
1994 }
1995 
1996 
1997 
1998 //======================================================================
1999 // Methods related to continua
2000 //======================================================================
2001 
2002 /* Workspace method: Doxygen documentation will be auto-generated */
2003 void abs_cont_descriptionInit(// WS Output:
2004  ArrayOfString& abs_cont_names,
2005  ArrayOfString& abs_cont_options,
2006  ArrayOfVector& abs_cont_parameters,
2007  const Verbosity& verbosity)
2008 {
2009  CREATE_OUT2;
2010 
2011  abs_cont_names.resize(0);
2012  abs_cont_options.resize(0);
2013  abs_cont_parameters.resize(0);
2014  out2 << " Initialized abs_cont_names \n"
2015  " abs_cont_models\n"
2016  " abs_cont_parameters.\n";
2017 }
2018 
2019 
2020 /* Workspace method: Doxygen documentation will be auto-generated */
2021 void abs_cont_descriptionAppend(// WS Output:
2022  ArrayOfString& abs_cont_names,
2023  ArrayOfString& abs_cont_models,
2024  ArrayOfVector& abs_cont_parameters,
2025  // Control Parameters:
2026  const String& tagname,
2027  const String& model,
2028  const Vector& userparameters,
2029  const Verbosity&)
2030 {
2031  // First we have to check that name matches a continuum species tag.
2032  check_continuum_model(tagname);
2033 
2034  //cout << " + tagname: " << tagname << "\n";
2035  //cout << " + model: " << model << "\n";
2036  //cout << " + parameters: " << userparameters << "\n";
2037 
2038  // Add name and parameters to the apropriate variables:
2039  abs_cont_names.push_back(tagname);
2040  abs_cont_models.push_back(model);
2041  abs_cont_parameters.push_back(userparameters);
2042 }
2043 
2044 
2045 /* Workspace method: Doxygen documentation will be auto-generated */
2047  Tensor4& propmat_clearsky,
2048  // WS Input:
2049  const ArrayOfMatrix& abs_coef_per_species,
2050  const Verbosity&)
2051 {
2052  // propmat_clearsky has format
2053  // [ abs_species, f_grid, stokes_dim, stokes_dim ].
2054  // abs_coef_per_species has format ArrayOfMatrix (over species),
2055  // where for each species the matrix has format [f_grid, abs_p].
2056 
2057 
2058  // Set stokes_dim (and check that the last two dimensions of
2059  // propmat_clearsky really are equal).
2060  Index nr, nc, stokes_dim;
2061  // In the two stokes dimensions, propmat_clearsky should be a
2062  // square matrix of stokes_dim*stokes_dim. Check this, and determine
2063  // stokes_dim:
2064  nr = propmat_clearsky.nrows();
2065  nc = propmat_clearsky.ncols();
2066  if ( nr!=nc )
2067  {
2068  ostringstream os;
2069  os << "The last two dimensions of propmat_clearsky must be equal (stokes_dim).\n"
2070  << "But here they are " << nr << " and " << nc << ".";
2071  throw runtime_error( os.str() );
2072  }
2073  stokes_dim = nr; // Could be nc here too, since they are the same.
2074 
2075 
2076  Index n_species = abs_coef_per_species.nelem(); // # species
2077 
2078  if (0==n_species)
2079  {
2080  ostringstream os;
2081  os << "Must have at least one species.";
2082  throw runtime_error(os.str());
2083  }
2084 
2085  Index n_f = abs_coef_per_species[0].nrows(); // # frequencies
2086 
2087  // # pressures must be 1:
2088  if (1!=abs_coef_per_species[0].ncols())
2089  {
2090  ostringstream os;
2091  os << "Must have exactly one pressure.";
2092  throw runtime_error(os.str());
2093  }
2094 
2095  // Check species dimension of propmat_clearsky
2096  if ( propmat_clearsky.nbooks()!=n_species )
2097  {
2098  ostringstream os;
2099  os << "Species dimension of propmat_clearsky does not\n"
2100  << "match abs_coef_per_species.";
2101  throw runtime_error( os.str() );
2102  }
2103 
2104  // Check frequency dimension of propmat_clearsky
2105  if ( propmat_clearsky.npages()!=n_f )
2106  {
2107  ostringstream os;
2108  os << "Frequency dimension of propmat_clearsky does not\n"
2109  << "match abs_coef_per_species.";
2110  throw runtime_error( os.str() );
2111  }
2112 
2113  // Loop species and stokes dimensions, and add to propmat_clearsky:
2114  for ( Index si=0; si<n_species; ++si )
2115  for ( Index ii=0; ii<stokes_dim; ++ii )
2116  propmat_clearsky(si,joker,ii, ii) += abs_coef_per_species[si](joker,0);
2117 
2118 }
2119 
2120 
2121 /* Workspace method: Doxygen documentation will be auto-generated */
2122 void propmat_clearskyInit(//WS Output
2123  Tensor4& propmat_clearsky,
2124  //WS Input
2125  const ArrayOfArrayOfSpeciesTag& abs_species,
2126  const Vector& f_grid,
2127  const Index& stokes_dim,
2128  const Index& propmat_clearsky_agenda_checked,
2129  const Verbosity&
2130  )
2131 {
2132  if (!propmat_clearsky_agenda_checked)
2133  throw runtime_error("You must call *propmat_clearsky_agenda_checkedCalc* before calling this method.");
2134 
2135  Index nf = f_grid.nelem();
2136 
2137  if(abs_species.nelem() > 0 )
2138  {
2139  if(nf > 0)
2140  {
2141  if(stokes_dim > 0)
2142  {
2143  propmat_clearsky.resize(abs_species.nelem(),nf, stokes_dim, stokes_dim);
2144  propmat_clearsky = 0;
2145  }
2146  else throw runtime_error("stokes_dim = 0");
2147  }
2148  else throw runtime_error("nf = 0");
2149  }
2150  else throw runtime_error("abs_species.nelem() = 0");
2151 }
2152 
2153 
2154 /* Workspace method: Doxygen documentation will be auto-generated */
2156  Tensor4& propmat_clearsky,
2157  const Index& stokes_dim,
2158  const Index& atmosphere_dim,
2159  const Vector& f_grid,
2160  const ArrayOfArrayOfSpeciesTag& abs_species,
2161  const Vector& rtp_vmr,
2162  const Vector& rtp_los,
2163  const Vector& rtp_mag,
2164  const Verbosity& )
2165 {
2166  // All the physical constants joined into one static constant:
2167  // (abs as e defined as negative)
2168  static const Numeric FRconst = abs(
2170  ( 8 * PI * PI * SPEED_OF_LIGHT * VACUUM_PERMITTIVITY *
2172 
2173  if( stokes_dim < 3 )
2174  throw runtime_error(
2175  "To include Faraday rotation, stokes_dim >= 3 is required." );
2176 
2177  if( atmosphere_dim==1 && rtp_los.nelem() < 1 )
2178  {
2179  ostringstream os;
2180  os << "For applying propmat_clearskyAddFaraday, los needs to be specified\n"
2181  << "(at least zenith angle component for atmosphere_dim==1),\n"
2182  << "but it is not.\n";
2183  throw runtime_error( os.str() );
2184  }
2185  else if( atmosphere_dim>1 && rtp_los.nelem() < 2 )
2186  {
2187  ostringstream os;
2188  os << "For applying propmat_clearskyAddFaraday, los needs to be specified\n"
2189  << "(both zenith and azimuth angle components for atmosphere_dim>1),\n"
2190  << "but it is not.\n";
2191  throw runtime_error( os.str() );
2192  }
2193 
2194  Index ife = -1;
2195  for( Index sp = 0; sp < abs_species.nelem() && ife < 0; sp++ )
2196  {
2197  if (abs_species[sp][0].Type() == SpeciesTag::TYPE_FREE_ELECTRONS)
2198  { ife = sp; }
2199  }
2200 
2201  if( ife < 0 )
2202  {
2203  throw runtime_error( "Free electrons not found in *abs_species* and "
2204  "Faraday rotation can not be calculated." );
2205  }
2206  else
2207  {
2208  const Numeric ne = rtp_vmr[ife];
2209 
2210  if( ne!=0 && ( rtp_mag[0]!=0 || rtp_mag[1]!=0 || rtp_mag[2]!=0 ) )
2211  {
2212  // Include remaining terms, beside /f^2
2213  const Numeric c1 = 2 * FRconst * ne * dotprod_with_los(
2214  rtp_los, rtp_mag[0], rtp_mag[1], rtp_mag[2], atmosphere_dim );
2215 
2216  for( Index iv=0; iv<f_grid.nelem(); iv++ )
2217  {
2218  const Numeric r = c1 / ( f_grid[iv] * f_grid[iv] );
2219  propmat_clearsky(ife,iv,1,2) = r;
2220  propmat_clearsky(ife,iv,2,1) = -r;
2221  }
2222  }
2223  }
2224 }
2225 
2226 
2227 /* Workspace method: Doxygen documentation will be auto-generated */
2229  // WS Output:
2230  Tensor4& propmat_clearsky,
2231  // WS Input:
2232  const Index& stokes_dim,
2233  const Index& atmosphere_dim,
2234  const Vector& f_grid,
2235  const ArrayOfArrayOfSpeciesTag& abs_species,
2236  const Vector& rtp_vmr,
2237  const Vector& rtp_los,
2238  const Numeric& rtp_temperature,
2239  const ArrayOfSingleScatteringData& scat_data_array,
2240  // Verbosity object:
2241  const Verbosity& verbosity)
2242 {
2243  const Index ns = scat_data_array.nelem();
2244  Index np = 0;
2245  for( Index sp = 0; sp < abs_species.nelem(); sp++ )
2246  {
2247  if (abs_species[sp][0].Type() == SpeciesTag::TYPE_PARTICLES)
2248  {
2249  np++;
2250  }
2251  }
2252 
2253  if( np == 0 )
2254  {
2255  ostringstream os;
2256  os << "For applying propmat_clearskyAddParticles, abs_species needs to"
2257  << "contain species 'particles', but it does not.\n";
2258  throw runtime_error( os.str() );
2259  }
2260 
2261  if ( ns != np )
2262  {
2263  ostringstream os;
2264  os << "Number of 'particles' entries in abs_species and of elements in\n"
2265  << "scat_data_array needs to be identical. But you have " << np
2266  << " 'particles' entries\n"
2267  << "and " << ns << " scat_data_array elements.\n";
2268  throw runtime_error( os.str() );
2269  }
2270 
2271  if( atmosphere_dim==1 && rtp_los.nelem() < 1 )
2272  {
2273  ostringstream os;
2274  os << "For applying propmat_clearskyAddParticles, los needs to be specified\n"
2275  << "(at least zenith angle component for atmosphere_dim==1),\n"
2276  << "but it is not.\n";
2277  throw runtime_error( os.str() );
2278  }
2279  else if( atmosphere_dim>1 && rtp_los.nelem() < 2 )
2280  {
2281  ostringstream os;
2282  os << "For applying propmat_clearskyAddParticles, los needs to be specified\n"
2283  << "(both zenith and azimuth angle components for atmosphere_dim>1),\n"
2284  << "but it is not.\n";
2285  throw runtime_error( os.str() );
2286  }
2287 
2288  const Index nf = f_grid.nelem();
2289  Vector rtp_los_back;
2290  mirror_los( rtp_los_back, rtp_los, atmosphere_dim );
2291  ArrayOfSingleScatteringData scat_data_array_mono;
2292  Matrix pnd_ext_mat(stokes_dim,stokes_dim);
2293  Vector pnd_abs_vec(stokes_dim);
2294 
2295  for( Index iv=0; iv<nf; iv++ )
2296  {
2297  // first, get the scat_data at the required frequency. we can do that for
2298  // all particle types at once.
2299  scat_data_array_monoCalc( scat_data_array_mono, scat_data_array, f_grid, iv,
2300  verbosity );
2301 
2302  // now we again need to loop over the abs_species entries. we need to do
2303  // this as there is no guarantee that all the 'particles' are
2304  // one-after-the-other. and we need separate output per abs_species entry,
2305  // i.e., per particle type.
2306  np = -1;
2307  for( Index sp = 0; sp < abs_species.nelem(); sp++ )
2308  {
2309  if (abs_species[sp][0].Type() == SpeciesTag::TYPE_PARTICLES)
2310  {
2311  np++;
2312 
2313  if ( rtp_vmr[sp] > 0. )
2314  {
2315  // second, get extinction matrix and absorption vector at required
2316  // temperature and direction for the individual particle type and
2317  // multiply with their occurence.
2318  opt_propExtract(pnd_ext_mat, pnd_abs_vec, scat_data_array_mono[np],
2319  rtp_los_back[0], rtp_los_back[1],
2320  rtp_temperature, stokes_dim, verbosity);
2321  //pnd_ext_mat *= rtp_vmr[sp];
2322  pnd_abs_vec *= rtp_vmr[sp];
2323 
2324  // last, sort the extracted absorption vector data into propmat_clearsky,
2325  // which is of extinction matrix type
2326 
2327  ext_matFromabs_vec(propmat_clearsky(sp,iv,joker,joker),
2328  pnd_abs_vec, stokes_dim);
2329  }
2330  }
2331  }
2332  }
2333 }
2334 
2335 
2336 /* Workspace method: Doxygen documentation will be auto-generated */
2337 void propmat_clearskyAddOnTheFly(// Workspace reference:
2338  Workspace& ws,
2339  // WS Output:
2340  Tensor4& propmat_clearsky,
2341  // WS Input:
2342  const Vector& f_grid,
2343  const ArrayOfArrayOfSpeciesTag& abs_species,
2344  const Numeric& rtp_pressure,
2345  const Numeric& rtp_temperature,
2346  const Vector& rtp_vmr,
2347  const Agenda& abs_xsec_agenda,
2348  // Verbosity object:
2349  const Verbosity& verbosity)
2350 {
2351  CREATE_OUT3;
2352 
2353  // Define communication variables for the actual absorption calculation:
2354 
2355  // Output of AbsInputFromRteScalars:
2356  Vector abs_p;
2357  Vector abs_t;
2358  Matrix abs_vmrs;
2359  // Output of abs_h2oSet:
2360  Vector abs_h2o;
2361  // Output of abs_coefCalc:
2362  Matrix abs_coef;
2363  ArrayOfMatrix abs_coef_per_species;
2364 
2365  AbsInputFromRteScalars(abs_p,
2366  abs_t,
2367  abs_vmrs,
2368  rtp_pressure,
2369  rtp_temperature,
2370  rtp_vmr,
2371  verbosity);
2372 
2373  // Absorption cross sections per tag group.
2374  ArrayOfMatrix abs_xsec_per_species;
2375 
2376  // Make all species active:
2377  ArrayOfIndex abs_species_active(abs_species.nelem());
2378  for (Index i=0; i<abs_species.nelem(); ++i)
2379  abs_species_active[i] = i;
2380 
2381 
2382  // Call agenda to calculate absorption:
2384  abs_xsec_per_species,
2385  abs_species,
2386  abs_species_active,
2387  f_grid,
2388  abs_p,
2389  abs_t,
2390  abs_vmrs,
2391  abs_xsec_agenda);
2392 
2393  // Calculate absorption coefficients from cross sections:
2394  abs_coefCalcFromXsec(abs_coef, abs_coef_per_species,
2395  abs_xsec_per_species,
2396  abs_vmrs, abs_p, abs_t, verbosity);
2397 
2398 
2399  // Now add abs_coef_per_species to propmat_clearsky:
2401  abs_coef_per_species,
2402  verbosity);
2403 }
2404 
2405 
2406 /* Workspace method: Doxygen documentation will be auto-generated */
2408  Tensor4& propmat_clearsky,
2409  const Vector& f_grid,
2410  const Index& stokes_dim,
2411  const Verbosity& )
2412 {
2413  propmat_clearsky.resize( 1, f_grid.nelem(), stokes_dim, stokes_dim );
2414  propmat_clearsky = 0;
2415 }
2416 
2417 
2418 
2419 
2420 /* Workspace method: Doxygen documentation will be auto-generated */
2422  const Verbosity&)
2423 {
2425 }
2426 
2427 
2428 #ifdef ENABLE_NETCDF
2429 /* Workspace method: Doxygen documentation will be auto-generated */
2430 /* Included by Claudia Emde, 20100707 */
2431 void WriteMolTau(//WS Input
2432  const Vector& f_grid,
2433  const Tensor3& z_field,
2434  const Tensor7& propmat_clearsky_field,
2435  const Index& atmosphere_dim,
2436  //Keyword
2437  const String& filename,
2438  const Verbosity&)
2439 {
2440 
2441  int retval, ncid;
2442  int nlev_dimid, nlyr_dimid, nwvl_dimid, stokes_dimid, none_dimid;
2443  int dimids[4];
2444  int wvlmin_varid, wvlmax_varid, z_varid, wvl_varid, tau_varid;
2445 
2446  if (atmosphere_dim != 1)
2447  throw runtime_error("WriteMolTau can only be used for atmsophere_dim=1");
2448 
2449 #pragma omp critical(netcdf__critical_region)
2450  {
2451  // Open file
2452  if ((retval = nc_create(filename.c_str(), NC_CLOBBER, &ncid)))
2453  nca_error (retval, "nc_create");
2454 
2455  // Define dimensions
2456  if ((retval = nc_def_dim(ncid, "nlev", (int) z_field.npages(), &nlev_dimid)))
2457  nca_error (retval, "nc_def_dim");
2458 
2459  if ((retval = nc_def_dim(ncid, "nlyr", (int) z_field.npages() - 1, &nlyr_dimid)))
2460  nca_error (retval, "nc_def_dim");
2461 
2462  if ((retval = nc_def_dim(ncid, "nwvl", (int) f_grid.nelem(), &nwvl_dimid)))
2463  nca_error (retval, "nc_def_dim");
2464 
2465  if ((retval = nc_def_dim(ncid, "none", 1, &none_dimid)))
2466  nca_error (retval, "nc_def_dim");
2467 
2468  if ((retval = nc_def_dim(ncid, "nstk", (int) propmat_clearsky_field.nbooks(), &stokes_dimid)))
2469  nca_error (retval, "nc_def_dim");
2470 
2471  // Define variables
2472  if ((retval = nc_def_var(ncid, "wvlmin", NC_DOUBLE, 1, &none_dimid, &wvlmin_varid)))
2473  nca_error (retval, "nc_def_var wvlmin");
2474 
2475  if ((retval = nc_def_var(ncid, "wvlmax", NC_DOUBLE, 1, &none_dimid, &wvlmax_varid)))
2476  nca_error (retval, "nc_def_var wvlmax");
2477 
2478  if ((retval = nc_def_var(ncid, "z", NC_DOUBLE, 1, &nlev_dimid, &z_varid)))
2479  nca_error (retval, "nc_def_var z");
2480 
2481  if ((retval = nc_def_var(ncid, "wvl", NC_DOUBLE, 1, &nwvl_dimid, &wvl_varid)))
2482  nca_error (retval, "nc_def_var wvl");
2483 
2484  dimids[0]=nlyr_dimid;
2485  dimids[1]=nwvl_dimid;
2486  dimids[2]=stokes_dimid;
2487  dimids[3]=stokes_dimid;
2488 
2489  if ((retval = nc_def_var(ncid, "tau", NC_DOUBLE, 4, &dimids[0], &tau_varid)))
2490  nca_error (retval, "nc_def_var tau");
2491 
2492  // Units
2493  if ((retval = nc_put_att_text(ncid, wvlmin_varid, "units", 2, "nm")))
2494  nca_error (retval, "nc_put_att_text");
2495 
2496  if ((retval = nc_put_att_text(ncid, wvlmax_varid, "units", 2, "nm")))
2497  nca_error (retval, "nc_put_att_text");
2498 
2499  if ((retval = nc_put_att_text(ncid, z_varid, "units", 2, "km")))
2500  nca_error (retval, "nc_put_att_text");
2501 
2502  if ((retval = nc_put_att_text(ncid, wvl_varid, "units", 2, "nm")))
2503  nca_error (retval, "nc_put_att_text");
2504 
2505  if ((retval = nc_put_att_text(ncid, tau_varid, "units", 1, "-")))
2506  nca_error (retval, "nc_put_att_text");
2507 
2508  // End define mode. This tells netCDF we are done defining
2509  // metadata.
2510  if ((retval = nc_enddef(ncid)))
2511  nca_error (retval, "nc_enddef");
2512 
2513  // Assign data
2514  double wvlmin[1];
2515  wvlmin[0]= SPEED_OF_LIGHT/f_grid[f_grid.nelem()-1]*1e9;
2516  if ((retval = nc_put_var_double (ncid, wvlmin_varid, &wvlmin[0])))
2517  nca_error (retval, "nc_put_var");
2518 
2519  double wvlmax[1];
2520  wvlmax[0]= SPEED_OF_LIGHT/f_grid[0]*1e9;
2521  if ((retval = nc_put_var_double (ncid, wvlmax_varid, &wvlmax[0])))
2522  nca_error (retval, "nc_put_var");
2523 
2524  double z[z_field.npages()];
2525  for (int iz=0; iz<z_field.npages(); iz++)
2526  z[iz]=z_field(z_field.npages()-1-iz, 0, 0)*1e-3;
2527 
2528  if ((retval = nc_put_var_double (ncid, z_varid, &z[0])))
2529  nca_error (retval, "nc_put_var");
2530 
2531  double wvl[f_grid.nelem()];
2532  for (int iv=0; iv<f_grid.nelem(); iv++)
2533  wvl[iv]=SPEED_OF_LIGHT/f_grid[f_grid.nelem()-1-iv]*1e9;
2534 
2535  if ((retval = nc_put_var_double (ncid, wvl_varid, &wvl[0])))
2536  nca_error (retval, "nc_put_var");
2537 
2538  const Index zfnp = z_field.npages()-1;
2539  const Index fgne = f_grid.nelem();
2540  const Index amfnb = propmat_clearsky_field.nbooks();
2541 
2542  Tensor4 tau(zfnp, fgne, amfnb, amfnb, 0.);
2543 
2544  // Calculate average tau for layers
2545  for (int is=0; is<propmat_clearsky_field.nlibraries(); is++)
2546  for (int iz=0; iz<zfnp; iz++)
2547  for (int iv=0; iv<fgne; iv++)
2548  for (int is1=0; is1<amfnb; is1++)
2549  for (int is2=0; is2<amfnb; is2++)
2550  // sum up all species
2551  tau(iz, iv, is1, is2) += 0.5 * (propmat_clearsky_field(is,f_grid.nelem()-1-iv,is1,is2,z_field.npages()-1-iz,0,0)+
2552  propmat_clearsky_field(is,f_grid.nelem()-1-iv,is1,is2,z_field.npages()-2-iz,0,0))
2553  *(z_field(z_field.npages()-1-iz,0,0)-z_field(z_field.npages()-2-iz,0,0));
2554 
2555 
2556  if ((retval = nc_put_var_double (ncid, tau_varid, tau.get_c_array())))
2557  nca_error (retval, "nc_put_var");
2558 
2559  // Close the file
2560  if ((retval = nc_close(ncid)))
2561  nca_error (retval, "nc_close");
2562  }
2563 }
2564 
2565 #else
2566 
2567 void WriteMolTau(//WS Input
2568  const Vector& f_grid _U_,
2569  const Tensor3& z_field _U_,
2570  const Tensor7& propmat_clearsky_field _U_,
2571  const Index& atmosphere_dim _U_,
2572  //Keyword
2573  const String& filename _U_,
2574  const Verbosity&)
2575 {
2576  throw runtime_error("The workspace method WriteMolTau is not available"
2577  "because ARTS was compiled without NetCDF support.");
2578 }
2579 
2580 #endif /* ENABLE_NETCDF */
2581 
Numeric dotprod_with_los(ConstVectorView los, const Numeric &u, const Numeric &v, const Numeric &w, const Index &atmosphere_dim)
dotprod_with_los
Definition: rte.cc:776
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
const SpeciesRecord & SpeciesData() const
The matching SpeciesRecord from species_data.
Definition: linerecord.cc:50
This header file contains all the declarations of the implemented continua and full absorption (lines...
void abs_xsec_per_speciesInit(ArrayOfMatrix &abs_xsec_per_species, const ArrayOfArrayOfSpeciesTag &tgs, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Index &abs_xsec_agenda_checked, const Verbosity &verbosity)
WORKSPACE METHOD: abs_xsec_per_speciesInit.
Definition: m_abs.cc:1468
void isotopologue_ratiosInitFromBuiltin(SpeciesAuxData &isotopologue_ratios, const Verbosity &)
WORKSPACE METHOD: isotopologue_ratiosInitFromBuiltin.
Definition: m_abs.cc:2421
This file contains basic functions to handle NetCDF data files.
void abs_xsec_per_speciesAddLines(ArrayOfMatrix &abs_xsec_per_species, const ArrayOfArrayOfSpeciesTag &tgs, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const Matrix &abs_vmrs, const ArrayOfArrayOfLineRecord &abs_lines_per_species, const ArrayOfLineshapeSpec &abs_lineshape, const SpeciesAuxData &isotopologue_ratios, const ArrayOfArrayOfLineMixingRecord &line_mixing_data, const ArrayOfArrayOfIndex &line_mixing_data_lut, const Verbosity &verbosity)
WORKSPACE METHOD: abs_xsec_per_speciesAddLines.
Definition: m_abs.cc:1531
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
Explicit construction of Arrays.
Definition: make_array.h:51
void abs_lines_per_speciesSetEmpty(ArrayOfArrayOfLineRecord &abs_lines_per_species, const ArrayOfArrayOfSpeciesTag &tgs, const Verbosity &)
WORKSPACE METHOD: abs_lines_per_speciesSetEmpty.
Definition: m_abs.cc:98
Declarations having to do with the four output streams.
bool ReadFromArtscat3Stream(istream &is, const Verbosity &verbosity)
Read one line from a stream associated with an ARTSCAT-3 file.
Definition: linerecord.cc:2064
Lineshape related specification like which lineshape to use, the normalizationfactor, and the cutoff.
Definition: absorption.h:141
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
Index Species() const
Molecular species index.
void nca_error(const int e, const String s)
Throws a runtime error for the given NetCDF error code.
Definition: nc_io.cc:671
void propmat_clearskyAddFaraday(Tensor4 &propmat_clearsky, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &rtp_vmr, const Vector &rtp_los, const Vector &rtp_mag, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyAddFaraday.
Definition: m_abs.cc:2155
const Numeric VACUUM_PERMITTIVITY
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
void abs_linesReadFromArts(ArrayOfLineRecord &abs_lines, const String &filename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_linesReadFromArts.
Definition: m_abs.cc:369
String get_tag_group_name(const ArrayOfSpeciesTag &tg)
Return the name of a tag group as a string.
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
void abs_coefCalcFromXsec(Matrix &abs_coef, ArrayOfMatrix &abs_coef_per_species, const ArrayOfMatrix &abs_xsec_per_species, const Matrix &abs_vmrs, const Vector &abs_p, const Vector &abs_t, const Verbosity &verbosity)
WORKSPACE METHOD: abs_coefCalcFromXsec.
Definition: m_abs.cc:1388
bool is_zeeman(const ArrayOfSpeciesTag &tg)
Is this a Zeeman tag group?
void abs_linesArtscat4FromArtscat3(ArrayOfLineRecord &abs_lines, const Verbosity &)
WORKSPACE METHOD: abs_linesArtscat4FromArtscat3.
Definition: m_abs.cc:114
This file contains basic functions to handle XML data files.
The Tensor7 class.
Definition: matpackVII.h:1931
void propmat_clearskyAddParticles(Tensor4 &propmat_clearsky, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &rtp_vmr, const Vector &rtp_los, const Numeric &rtp_temperature, const ArrayOfSingleScatteringData &scat_data_array, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddParticles.
Definition: m_abs.cc:2228
This file contains basic functions to handle ASCII files.
void find_xml_file(String &filename, const Verbosity &verbosity)
Find an xml file.
Definition: file.cc:464
void abs_lines_per_speciesReadFromCatalogues(ArrayOfArrayOfLineRecord &abs_lines_per_species, const ArrayOfArrayOfSpeciesTag &tgs, const ArrayOfString &filenames, const ArrayOfString &formats, const Vector &fmin, const Vector &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lines_per_speciesReadFromCatalogues.
Definition: m_abs.cc:549
const Index & Ind_lsn() const
Return the index of the normalization factor.
Definition: absorption.h:165
void propmat_clearskyZero(Tensor4 &propmat_clearsky, const Vector &f_grid, const Index &stokes_dim, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyZero.
Definition: m_abs.cc:2407
Numeric Uf() const
The upper line center frequency in Hz: If this is <0 it means no upper limit.
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:75
Numeric number_density(const Numeric &p, const Numeric &t)
number_density
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.
const Numeric ELECTRON_CHARGE
const Numeric SPEED_OF_LIGHT
This file contains the definition of Array.
const Array< IsotopologueRecord > & Isotopologue() const
Definition: absorption.h:425
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
void opt_propExtract(MatrixView ext_mat_mono_spt, VectorView abs_vec_mono_spt, const SingleScatteringData &scat_data, const Numeric za, const Numeric aa, const Numeric rtp_temperature, const Index stokes_dim, const Verbosity &verbosity)
Definition: montecarlo.cc:751
void xml_read_arts_catalogue_from_file(const String &filename, ArrayOfLineRecord &type, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
Definition: xml_io.cc:900
void abs_lineshape_per_tgDefine(ArrayOfLineshapeSpec &abs_lineshape, const ArrayOfArrayOfSpeciesTag &tgs, const ArrayOfString &shape, const ArrayOfString &normalizationfactor, const Vector &cutoff, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lineshape_per_tgDefine.
Definition: m_abs.cc:1182
Index Version() const
Return the version number.
Definition: linerecord.h:298
void abs_linesReadFromSplitArtscat(ArrayOfLineRecord &abs_lines, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_linesReadFromSplitArtscat.
Definition: m_abs.cc:410
void scat_data_array_monoCalc(ArrayOfSingleScatteringData &scat_data_array_mono, const ArrayOfSingleScatteringData &scat_data_array, const Vector &f_grid, const Index &f_index, const Verbosity &)
WORKSPACE METHOD: scat_data_array_monoCalc.
void AbsInputFromRteScalars(Vector &abs_p, Vector &abs_t, Matrix &abs_vmrs, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const Verbosity &)
WORKSPACE METHOD: AbsInputFromRteScalars.
Definition: m_abs.cc:73
void propmat_clearskyAddFromAbsCoefPerSpecies(Tensor4 &propmat_clearsky, const ArrayOfMatrix &abs_coef_per_species, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyAddFromAbsCoefPerSpecies.
Definition: m_abs.cc:2046
void abs_lines_per_speciesWriteToSplitArtscat(const String &output_file_format, const ArrayOfArrayOfLineRecord &abs_lines_per_species, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lines_per_speciesWriteToSplitArtscat.
Definition: m_abs.cc:382
void abs_cont_descriptionInit(ArrayOfString &abs_cont_names, ArrayOfString &abs_cont_options, ArrayOfVector &abs_cont_parameters, const Verbosity &verbosity)
WORKSPACE METHOD: abs_cont_descriptionInit.
Definition: m_abs.cc:2003
void abs_lines_per_speciesCreateFromLines(ArrayOfArrayOfLineRecord &abs_lines_per_species, const ArrayOfLineRecord &abs_lines, const ArrayOfArrayOfSpeciesTag &tgs, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lines_per_speciesCreateFromLines.
Definition: m_abs.cc:759
void tolower()
Convert to lower case.
Definition: mystring.h:91
Numeric Lf() const
The lower line center frequency in Hz.
A tag group can consist of the sum of several of these.
void propmat_clearskyInit(Tensor4 &propmat_clearsky, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &f_grid, const Index &stokes_dim, const Index &propmat_clearsky_agenda_checked, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyInit.
Definition: m_abs.cc:2122
Index Isotopologue() const
The index of the isotopologue species that this line belongs to.
Definition: linerecord.h:307
void AbsInputFromAtmFields(Vector &abs_p, Vector &abs_t, Matrix &abs_vmrs, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const Tensor4 &vmr_field, const Verbosity &)
WORKSPACE METHOD: AbsInputFromAtmFields.
Definition: m_abs.cc:1360
const Numeric PI
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
void abs_xsec_per_speciesAddConts(ArrayOfMatrix &abs_xsec_per_species, const ArrayOfArrayOfSpeciesTag &tgs, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const Matrix &abs_vmrs, const ArrayOfString &abs_cont_names, const ArrayOfVector &abs_cont_parameters, const ArrayOfString &abs_cont_models, const Verbosity &verbosity)
WORKSPACE METHOD: abs_xsec_per_speciesAddConts.
Definition: m_abs.cc:1806
void WriteMolTau(const Vector &f_grid, const Tensor3 &z_field, const Tensor7 &propmat_clearsky_field, const Index &atmosphere_dim, const String &filename, const Verbosity &)
WORKSPACE METHOD: WriteMolTau.
Definition: m_abs.cc:2431
void abs_speciesDefineAllInScenario(ArrayOfArrayOfSpeciesTag &tgs, Index &propmat_clearsky_agenda_checked, Index &abs_xsec_agenda_checked, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesDefineAllInScenario.
Definition: m_abs.cc:1047
void xsec_continuum_tag(MatrixView xsec, const String &name, ConstVectorView parameters, const String &model, ConstVectorView f_grid, ConstVectorView abs_p, ConstVectorView abs_t, ConstVectorView abs_n2, ConstVectorView abs_h2o, ConstVectorView abs_o2, ConstVectorView vmr, const Verbosity &verbosity)
Calculates model absorption for one continuum or full model tag.
Definition: continua.cc:16085
void abs_linesReadFromArtsObsolete(ArrayOfLineRecord &abs_lines, const String &filename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
Obsolete old ARTS catalogue reading function.
Definition: m_abs.cc:476
#define abs(x)
Definition: continua.cc:20458
void abs_o2Set(Vector &abs_o2, const ArrayOfArrayOfSpeciesTag &abs_species, const Matrix &abs_vmrs, const Verbosity &)
abs_o2Set.
Definition: m_abs.cc:1342
const Joker joker
const Numeric & Cutoff() const
Return the cutoff frequency (in Hz).
Definition: absorption.h:172
String VersionString() const
Return the version String.
Definition: linerecord.cc:34
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:32
void propmat_clearskyAddOnTheFly(Workspace &ws, Tensor4 &propmat_clearsky, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const Agenda &abs_xsec_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddOnTheFly.
Definition: m_abs.cc:2337
void abs_linesReadFromHitranPre2004(ArrayOfLineRecord &abs_lines, const String &filename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_linesReadFromHitranPre2004.
Definition: m_abs.cc:144
Implements the class MakeArray, which is a derived class of Array, allowing explicit initialization...
void abs_lineshapeDefine(ArrayOfLineshapeSpec &abs_lineshape, const String &shape, const String &normalizationfactor, const Numeric &cutoff, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lineshapeDefine.
Definition: m_abs.cc:1113
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 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
bool ReadFromJplStream(istream &is, const Verbosity &verbosity)
Read one line from a stream associated with a JPL file.
Definition: linerecord.cc:1818
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:1260
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
This can be used to make arrays out of anything.
Definition: array.h:40
void split(Array< my_basic_string< charT > > &aos, const my_basic_string< charT > &delim) const
Split string into substrings.
Definition: mystring.h:233
void open_input_file(ifstream &file, const String &name)
Open a file for reading.
Definition: file.cc:166
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:465
void abs_lines_per_speciesAddMirrorLines(ArrayOfArrayOfLineRecord &abs_lines_per_species, const Numeric &max_f, const Verbosity &)
WORKSPACE METHOD: abs_lines_per_speciesAddMirrorLines.
Definition: m_abs.cc:901
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 abs_linesReadFromMytran2(ArrayOfLineRecord &abs_lines, const String &filename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_linesReadFromMytran2.
Definition: m_abs.cc:283
Array< Array< LineRecord > > ArrayOfArrayOfLineRecord
Holds a lists of spectral line data for each tag group.
Definition: linerecord.h:1097
The class MakeVector is a special kind of Vector that can be initialized explicitly from one or more ...
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
#define ns
Definition: continua.cc:21931
#define min(a, b)
Definition: continua.cc:20460
void abs_cont_descriptionAppend(ArrayOfString &abs_cont_names, ArrayOfString &abs_cont_models, ArrayOfVector &abs_cont_parameters, const String &tagname, const String &model, const Vector &userparameters, const Verbosity &)
WORKSPACE METHOD: abs_cont_descriptionAppend.
Definition: m_abs.cc:2021
Workspace methods and template functions for supergeneric XML IO.
void abs_n2Set(Vector &abs_n2, const ArrayOfArrayOfSpeciesTag &abs_species, const Matrix &abs_vmrs, const Verbosity &)
abs_n2Set.
Definition: m_abs.cc:1315
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
void abs_linesReadFromJpl(ArrayOfLineRecord &abs_lines, const String &filename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_linesReadFromJpl.
Definition: m_abs.cc:324
#define _U_
Definition: config.h:167
const Numeric ELECTRON_MASS
void check_continuum_model(const String &name)
An auxiliary functions that checks if a given continuum model is listed in species_data.cc.
Definition: continua.cc:20236
void abs_h2oSet(Vector &abs_h2o, const ArrayOfArrayOfSpeciesTag &abs_species, const Matrix &abs_vmrs, const Verbosity &)
abs_h2oSet.
Definition: m_abs.cc:1287
void fillSpeciesAuxDataWithIsotopologueRatiosFromSpeciesData(SpeciesAuxData &sad)
Fill SpeciesAuxData with default isotopologue ratios from species data.
Definition: absorption.cc:305
#define CREATE_OUT3
Definition: messages.h:214
void mirror_los(Vector &los_mirrored, ConstVectorView los, const Index &atmosphere_dim)
mirror_los
Definition: rte.cc:2389
void abs_lines_per_speciesCompact(ArrayOfArrayOfLineRecord &abs_lines_per_species, const ArrayOfLineshapeSpec &abs_lineshape, const Vector &f_grid, const Verbosity &)
WORKSPACE METHOD: abs_lines_per_speciesCompact.
Definition: m_abs.cc:953
void WriteXML(Workspace &ws, const String &file_format, const Agenda &v, const String &f, const Index &no_clobber, const String &v_name, const String &f_name, const String &no_clobber_name, const Verbosity &verbosity)
Definition: m_xml.cc:67
void setF(Numeric new_mf)
Set the line center frequency in Hz.
Definition: linerecord.h:342
Array< LineRecord > ArrayOfLineRecord
Holds a list of spectral line data.
Definition: linerecord.h:1092
Auxiliary data for isotopologues.
Definition: absorption.h:439
This file contains header information for the dealing with command line parameters.
Index Species() const
The index of the molecular species that this line belongs to.
Definition: linerecord.h:302
#define ll
Definition: continua.cc:21705
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
void abs_xsec_agendaExecute(Workspace &ws, ArrayOfMatrix &abs_xsec_per_species, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const Matrix &abs_vmrs, const Agenda &input_agenda)
Definition: auto_md.cc:13122
#define CREATE_OUT2
Definition: messages.h:213
bool ReadFromHitran2004Stream(istream &is, const Verbosity &verbosity, const Numeric fmin=0)
Read one line from a stream associated with a HITRAN 2004 file.
Definition: linerecord.cc:558
bool ReadFromHitran2001Stream(istream &is, const Verbosity &verbosity)
Read one line from a stream associated with a HITRAN 1986-2001 file.
Definition: linerecord.cc:116
Scattering database structure and functions.
Index Isotopologue() const
Isotopologue species index.
void ext_matFromabs_vec(MatrixView ext_mat, ConstVectorView abs_vec, const Index &stokes_dim)
Derive extinction matrix from absorption vector.
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
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:50
void abs_linesReadFromHitran(ArrayOfLineRecord &abs_lines, const String &filename, const Numeric &fmin, const Numeric &fmax, const Verbosity &verbosity)
WORKSPACE METHOD: abs_linesReadFromHitran.
Definition: m_abs.cc:189
Spectral line catalog data.
Definition: linerecord.h:196
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1403
bool ReadFromMytran2Stream(istream &is, const Verbosity &verbosity)
Read one line from a stream associated with a MYTRAN2 file.
Definition: linerecord.cc:1420
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackIV.cc:402