00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00030 #include "arts.h"
00031 #include <map>
00032 #include <cfloat>
00033 #include <algorithm>
00034 #include <cmath>
00035 #include "absorption.h"
00036 #include "math_funcs.h"
00037 #include "auto_md.h"
00038 #include "messages.h"
00039
00040
00042 std::map<String, Index> SpeciesMap;
00043
00044
00045
00046
00047
00048 Numeric IsotopeRecord::CalculatePartitionFctAtTemp( Numeric
00049 temperature ) const
00050 {
00051 Numeric result = 0.;
00052 Numeric exponent = 1.;
00053
00054 ArrayOfNumeric::const_iterator it;
00055
00056 for (it=mqcoeff.begin(); it != mqcoeff.end(); it++)
00057 {
00058 result += *it * exponent;
00059 exponent *= temperature;
00060 }
00061 return result;
00062 }
00063
00064 void define_species_map()
00065 {
00066 extern const Array<SpeciesRecord> species_data;
00067 extern std::map<String, Index> SpeciesMap;
00068
00069 for ( Index i=0 ; i<species_data.nelem() ; ++i)
00070 {
00071 SpeciesMap[species_data[i].Name()] = i;
00072 }
00073 }
00074
00075
00076 ostream& operator << (ostream& os, const LineRecord& lr)
00077 {
00078
00079
00080 Index precision;
00081 switch (sizeof(Numeric)) {
00082 case sizeof(float) : precision = FLT_DIG; break;
00083 case sizeof(double) : precision = DBL_DIG; break;
00084 default: out0 << "Numeric must be double or float\n"; exit(1);
00085 }
00086
00087 os << "@"
00088 << " " << lr.Name ()
00089 << " "
00090 << setprecision(precision)
00091 << lr.F ()
00092 << " " << lr.Psf ()
00093 << " " << lr.I0 ()
00094 << " " << lr.Ti0 ()
00095 << " " << lr.Elow ()
00096 << " " << lr.Agam ()
00097 << " " << lr.Sgam ()
00098 << " " << lr.Nair ()
00099 << " " << lr.Nself ()
00100 << " " << lr.Tgam ()
00101 << " " << lr.Naux ()
00102 << " " << lr.dF ()
00103 << " " << lr.dI0 ()
00104 << " " << lr.dAgam ()
00105 << " " << lr.dSgam ()
00106 << " " << lr.dNair ()
00107 << " " << lr.dNself()
00108 << " " << lr.dPsf ();
00109
00110
00111 for ( Index i=0; i<lr.Naux(); ++i )
00112 os << " " << lr.Aux()[i];
00113
00114 return os;
00115 }
00116
00117
00126 template<class T>
00127 void extract(T& x,
00128 String& line,
00129 Index n)
00130 {
00131
00132
00133 x = T(0);
00134
00135
00136
00137
00138 istringstream item( line.substr(0,n) );
00139
00140
00141
00142
00143
00144
00145 line.erase(0,n);
00146
00147
00148
00149 item >> x;
00150 }
00151
00152
00153
00154
00155
00156
00157
00158
00159 bool LineRecord::ReadFromHitranStream(istream& is)
00160 {
00161
00162 extern Array<SpeciesRecord> species_data;
00163
00164
00165
00166
00167 const Index missing = species_data.nelem() + 100;
00168
00169
00170
00171
00172
00173
00174 static Array< Index > hspec(100);
00175
00176
00177
00178 static Array< ArrayOfIndex > hiso(100);
00179
00180
00181 static bool hinit = false;
00182
00183
00184
00185 static ArrayOfIndex warned_missing;
00186
00187 if ( !hinit )
00188 {
00189
00190
00191 hspec = missing;
00192
00193 for ( Index i=0; i<species_data.nelem(); ++i )
00194 {
00195 const SpeciesRecord& sr = species_data[i];
00196
00197
00198
00199
00200 if ( 0 < sr.Isotope()[0].HitranTag() )
00201 {
00202
00203
00204
00205
00206
00207
00208
00209 Index mo = sr.Isotope()[0].HitranTag() / 10;
00210
00211 hspec[mo] = i;
00212
00213
00214 Index n_iso = sr.Isotope().nelem();
00215 ArrayOfIndex iso_tags;
00216 iso_tags.resize(n_iso);
00217 for ( Index j=0; j<n_iso; ++j )
00218 {
00219 iso_tags[j] = sr.Isotope()[j].HitranTag();
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229 hiso[mo].resize( max(iso_tags)%10 + 1 );
00230 hiso[mo] = missing;
00231
00232
00233
00234 for ( Index j=0; j<n_iso; ++j )
00235 {
00236 if ( 0 < iso_tags[j] )
00237 {
00238
00239
00240 hiso[mo][iso_tags[j] % 10] = j;
00241 }
00242 }
00243 }
00244 }
00245
00246
00247
00248 out3 << " HITRAN index table:\n";
00249 for ( Index i=0; i<hspec.nelem(); ++i )
00250 {
00251 if ( missing != hspec[i] )
00252 {
00253
00254
00255
00256 out3 << " mo = " << i << " Species = "
00257 << setw(10) << setiosflags(ios::left)
00258 << species_data[hspec[i]].Name().c_str()
00259 << "iso = ";
00260 for ( Index j=1; j<hiso[i].nelem(); ++j )
00261 {
00262 if ( missing==hiso[i][j] )
00263 out3 << " " << "m";
00264 else
00265 out3 << " " << species_data[hspec[i]].Isotope()[hiso[i][j]].Name();
00266 }
00267 out3 << "\n";
00268 }
00269 }
00270
00271 hinit = true;
00272 }
00273
00274
00275
00276
00277
00278 String line;
00279
00280
00281 Index mo;
00282
00283
00284 bool comment = true;
00285
00286 while (comment)
00287 {
00288
00289 if (is.eof()) return true;
00290
00291
00292 if (!is) throw runtime_error ("Stream bad.");
00293
00294
00295 getline(is,line);
00296
00297
00298
00299
00300
00301 if (line.nelem() == 0 && is.eof()) return true;
00302
00303
00304
00305 if (line[line.nelem () - 1] == 13)
00306 {
00307 line.erase (line.nelem () - 1, 1);
00308 }
00309
00310
00311
00312
00313
00314 mo = 0;
00315
00316
00317 extract(mo,line,2);
00318
00319
00320
00321 if ( 0 != mo )
00322 {
00323
00324 if ( missing != hspec[mo] )
00325 {
00326 comment = false;
00327
00328
00329
00330 Numeric nChar = line.nelem() + 2;
00331 if ( nChar != 100 )
00332 {
00333 ostringstream os;
00334 os << "Invalid HITRAN 1986-2001 line data record with " << nChar <<
00335 " characters (expected: 100)." << endl << line << " n: " << line.nelem ();
00336 throw runtime_error(os.str());
00337 }
00338
00339 }
00340 else
00341 {
00342
00343
00344 if ( 0 == std::count(warned_missing.begin(),
00345 warned_missing.end(),
00346 mo) )
00347 {
00348 out0 << "Error: HITRAN mo = " << mo << " is not "
00349 << "known to ARTS.\n";
00350 warned_missing.push_back(mo);
00351 exit(1);
00352
00353
00354
00355 }
00356 }
00357 }
00358 }
00359
00360
00361
00362
00363 mspecies = hspec[mo];
00364
00365
00366 Index iso;
00367 extract(iso,line,1);
00368
00369
00370
00371
00372
00373
00374
00375 misotope = missing;
00376 if ( iso < hiso[mo].nelem() )
00377 if ( missing != hiso[mo][iso] )
00378 misotope = hiso[mo][iso];
00379
00380
00381 if (missing == misotope)
00382 {
00383 ostringstream os;
00384 os << "Species: " << species_data[mspecies].Name()
00385 << ", isotope iso = " << iso
00386 << " is unknown.";
00387 throw runtime_error(os.str());
00388 }
00389
00390
00391
00392 {
00393
00394 Numeric v;
00395
00396 extern const Numeric SPEED_OF_LIGHT;
00397
00398
00399
00400 const Numeric w2Hz = SPEED_OF_LIGHT * 100.;
00401
00402
00403 extract(v,line,12);
00404
00405
00406 mf = v * w2Hz;
00407
00408 }
00409
00410
00411 {
00412 extern const Numeric SPEED_OF_LIGHT;
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424 const Numeric hi2arts = 1e-2 * SPEED_OF_LIGHT;
00425
00426 Numeric s;
00427
00428
00429 extract(s,line,10);
00430
00431 mi0 = s * hi2arts;
00432
00433 mi0 /= species_data[mspecies].Isotope()[misotope].Abundance();
00434 }
00435
00436
00437 {
00438 Numeric r;
00439 extract(r,line,10);
00440 }
00441
00442
00443
00444 {
00445
00446
00447 Numeric gam;
00448
00449
00450 extern const Numeric ATM2PA;
00451
00452 extern const Numeric SPEED_OF_LIGHT;
00453
00454
00455 const Numeric w2Hz = SPEED_OF_LIGHT * 1e2;
00456
00457 const Numeric hi2arts = w2Hz / ATM2PA;
00458
00459
00460 extract(gam,line,5);
00461
00462
00463 magam = gam * hi2arts;
00464
00465
00466 extract(gam,line,5);
00467
00468
00469 msgam = gam * hi2arts;
00470
00471
00472 if (0==msgam)
00473 msgam = magam;
00474
00475
00476 }
00477
00478
00479
00480 {
00481
00482
00483
00484
00485 extract(melow,line,10);
00486
00487
00488 melow = wavenumber_to_joule(melow);
00489 }
00490
00491
00492
00493 {
00494
00495 extract(mnair,line,4);
00496
00497
00498 mnself = mnair;
00499
00500 }
00501
00502
00503
00504 {
00505
00506
00507 Numeric d;
00508
00509
00510 extern const Numeric ATM2PA;
00511
00512 extern const Numeric SPEED_OF_LIGHT;
00513
00514
00515 const Numeric w2Hz = SPEED_OF_LIGHT * 1e2;
00516
00517 const Numeric hi2arts = w2Hz / ATM2PA;
00518
00519
00520 extract(d,line,8);
00521
00522
00523 mpsf = d * hi2arts;
00524 }
00525
00526
00527
00528
00529 {
00530 Index eu;
00531 extract(eu,line,3);
00532 }
00533
00534
00535 {
00536 Index el;
00537 extract(el,line,3);
00538 }
00539
00540
00541 {
00542 Index eul;
00543 extract(eul,line,9);
00544 }
00545
00546
00547 {
00548 Index ell;
00549 extract(ell,line,9);
00550 }
00551
00552
00553 {
00554 Index df;
00555
00556 extract(df,line,1);
00557
00558 convHitranIERF(mdf,df);
00559 }
00560
00561
00562 {
00563 Index di0;
00564
00565 extract(di0,line,1);
00566 convHitranIERSH(mdi0,di0);
00567 }
00568
00569
00570 {
00571 Index dgam;
00572
00573 extract(dgam,line,1);
00574
00575 convHitranIERSH(mdagam,dgam);
00576
00577
00578
00579 }
00580
00581
00582
00583 mdpsf =-1;
00584
00585
00586
00587
00588
00589
00590
00591 mti0 = 296.0;
00592
00593
00594
00595 mtgam = 296.0;
00596
00597
00598 return false;
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 bool LineRecord::ReadFromHitran2004Stream(istream& is)
00610 {
00611
00612 extern Array<SpeciesRecord> species_data;
00613
00614
00615
00616
00617 const Index missing = species_data.nelem() + 100;
00618
00619
00620
00621
00622
00623
00624 static Array< Index > hspec(100);
00625
00626
00627
00628 static Array< ArrayOfIndex > hiso(100);
00629
00630
00631 static bool hinit = false;
00632
00633
00634
00635 static ArrayOfIndex warned_missing;
00636
00637 if ( !hinit )
00638 {
00639
00640
00641 hspec = missing;
00642
00643 for ( Index i=0; i<species_data.nelem(); ++i )
00644 {
00645 const SpeciesRecord& sr = species_data[i];
00646
00647
00648
00649
00650 if ( 0 < sr.Isotope()[0].HitranTag() )
00651 {
00652
00653
00654
00655
00656
00657
00658
00659 Index mo = sr.Isotope()[0].HitranTag() / 10;
00660
00661 hspec[mo] = i;
00662
00663
00664 Index n_iso = sr.Isotope().nelem();
00665 ArrayOfIndex iso_tags;
00666 iso_tags.resize(n_iso);
00667 for ( Index j=0; j<n_iso; ++j )
00668 {
00669 iso_tags[j] = sr.Isotope()[j].HitranTag();
00670 }
00671
00672
00673
00674
00675
00676
00677
00678
00679 hiso[mo].resize( max(iso_tags)%10 + 1 );
00680 hiso[mo] = missing;
00681
00682
00683
00684 for ( Index j=0; j<n_iso; ++j )
00685 {
00686 if ( 0 < iso_tags[j] )
00687 {
00688
00689
00690 hiso[mo][iso_tags[j] % 10] = j;
00691 }
00692 }
00693 }
00694 }
00695
00696
00697
00698 out3 << " HITRAN index table:\n";
00699 for ( Index i=0; i<hspec.nelem(); ++i )
00700 {
00701 if ( missing != hspec[i] )
00702 {
00703
00704
00705
00706 out3 << " mo = " << i << " Species = "
00707 << setw(10) << setiosflags(ios::left)
00708 << species_data[hspec[i]].Name().c_str()
00709 << "iso = ";
00710 for ( Index j=1; j<hiso[i].nelem(); ++j )
00711 {
00712 if ( missing==hiso[i][j] )
00713 out3 << " " << "m";
00714 else
00715 out3 << " " << species_data[hspec[i]].Isotope()[hiso[i][j]].Name();
00716 }
00717 out3 << "\n";
00718 }
00719 }
00720
00721 hinit = true;
00722 }
00723
00724
00725
00726
00727
00728 String line;
00729
00730
00731 Index mo;
00732
00733
00734 bool comment = true;
00735
00736 while (comment)
00737 {
00738
00739 if (is.eof()) return true;
00740
00741
00742 if (!is) throw runtime_error ("Stream bad.");
00743
00744
00745 getline(is,line);
00746
00747
00748
00749
00750
00751 if (line.nelem() == 0 && is.eof()) return true;
00752
00753
00754
00755 if (line[line.nelem () - 1] == 13)
00756 {
00757 line.erase (line.nelem () - 1, 1);
00758 }
00759
00760
00761
00762
00763
00764 mo = 0;
00765
00766
00767 extract(mo,line,2);
00768
00769
00770
00771 if ( 0 != mo )
00772 {
00773
00774 if ( missing != hspec[mo] )
00775 {
00776 comment = false;
00777
00778
00779
00780 Numeric nChar = line.nelem() + 2;
00781 if ( nChar != 160 )
00782 {
00783 ostringstream os;
00784 os << "Invalid HITRAN 2004 line data record with " << nChar <<
00785 " characters (expected: 160).";
00786 throw runtime_error(os.str());
00787 }
00788
00789 }
00790 else
00791 {
00792
00793
00794 if ( 0 == std::count(warned_missing.begin(),
00795 warned_missing.end(),
00796 mo) )
00797 {
00798 out0 << "Warning: HITRAN molecule number mo = " << mo << " is not "
00799 << "known to ARTS.\n";
00800 warned_missing.push_back(mo);
00801 }
00802 }
00803 }
00804 }
00805
00806
00807
00808
00809 mspecies = hspec[mo];
00810
00811
00812 Index iso;
00813 extract(iso,line,1);
00814
00815
00816
00817
00818
00819
00820
00821 misotope = missing;
00822 if ( iso < hiso[mo].nelem() )
00823 if ( missing != hiso[mo][iso] )
00824 misotope = hiso[mo][iso];
00825
00826
00827 if (missing == misotope)
00828 {
00829 ostringstream os;
00830 os << "Species: " << species_data[mspecies].Name()
00831 << ", isotope iso = " << iso
00832 << " is unknown.";
00833 throw runtime_error(os.str());
00834 }
00835
00836
00837
00838 {
00839
00840 Numeric v;
00841
00842 extern const Numeric SPEED_OF_LIGHT;
00843
00844
00845
00846 const Numeric w2Hz = SPEED_OF_LIGHT * 100.;
00847
00848
00849 extract(v,line,12);
00850
00851
00852 mf = v * w2Hz;
00853
00854 }
00855
00856
00857 {
00858 extern const Numeric SPEED_OF_LIGHT;
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870 const Numeric hi2arts = 1e-2 * SPEED_OF_LIGHT;
00871
00872 Numeric s;
00873
00874
00875 extract(s,line,10);
00876
00877 mi0 = s * hi2arts;
00878
00879 mi0 /= species_data[mspecies].Isotope()[misotope].Abundance();
00880 }
00881
00882
00883 {
00884 Numeric r;
00885 extract(r,line,10);
00886 }
00887
00888
00889
00890 {
00891
00892
00893 Numeric gam;
00894
00895
00896 extern const Numeric ATM2PA;
00897
00898 extern const Numeric SPEED_OF_LIGHT;
00899
00900
00901 const Numeric w2Hz = SPEED_OF_LIGHT * 1e2;
00902
00903 const Numeric hi2arts = w2Hz / ATM2PA;
00904
00905
00906 extract(gam,line,5);
00907
00908
00909 magam = gam * hi2arts;
00910
00911
00912 extract(gam,line,5);
00913
00914
00915 msgam = gam * hi2arts;
00916
00917
00918 if (0==msgam)
00919 msgam = magam;
00920
00921
00922 }
00923
00924
00925
00926 {
00927
00928
00929
00930
00931 extract(melow,line,10);
00932
00933
00934 melow = wavenumber_to_joule(melow);
00935 }
00936
00937
00938
00939 {
00940
00941 extract(mnair,line,4);
00942
00943
00944 mnself = mnair;
00945
00946 }
00947
00948
00949
00950 {
00951
00952
00953 Numeric d;
00954
00955
00956 extern const Numeric ATM2PA;
00957
00958 extern const Numeric SPEED_OF_LIGHT;
00959
00960
00961 const Numeric w2Hz = SPEED_OF_LIGHT * 1e2;
00962
00963 const Numeric hi2arts = w2Hz / ATM2PA;
00964
00965
00966 extract(d,line,8);
00967
00968
00969 mpsf = d * hi2arts;
00970 }
00971
00972
00973
00974
00975 {
00976 Index eu;
00977 extract(eu,line,15);
00978 }
00979
00980
00981 {
00982 Index el;
00983 extract(el,line,15);
00984 }
00985
00986
00987 {
00988 Index eul;
00989 extract(eul,line,15);
00990 }
00991
00992
00993 {
00994 Index ell;
00995 extract(ell,line,15);
00996 }
00997
00998
00999 {
01000 Index df;
01001
01002 extract(df,line,1);
01003
01004 convHitranIERF(mdf,df);
01005 }
01006
01007
01008 {
01009 Index di0;
01010
01011 extract(di0,line,1);
01012
01013 convHitranIERSH(mdi0,di0);
01014 }
01015
01016
01017 {
01018 Index dagam;
01019
01020 extract(dagam,line,1);
01021
01022 convHitranIERSH(mdagam,dagam);
01023 }
01024
01025
01026 {
01027 Index dsgam;
01028
01029 extract(dsgam,line,1);
01030
01031 convHitranIERSH(mdsgam,dsgam);
01032 }
01033
01034
01035 {
01036 Index dnair;
01037
01038 extract(dnair,line,1);
01039
01040 convHitranIERSH(mdnair,dnair);
01041 }
01042
01043
01044
01045 mdnself =-1;
01046
01047
01048 {
01049 Index dpsf;
01050
01051 extract(dpsf,line,1);
01052
01053 convHitranIERF(mdpsf,dpsf);
01054
01055 mdpsf = mdpsf / mf;
01056 }
01057
01058
01059
01060
01061
01062
01063
01064 mti0 = 296.0;
01065
01066
01067
01068 mtgam = 296.0;
01069
01070
01071 return false;
01072 }
01073
01074
01075 bool LineRecord::ReadFromMytran2Stream(istream& is)
01076 {
01077
01078 extern Array<SpeciesRecord> species_data;
01079
01080
01081
01082
01083 const Index missing = species_data.nelem() + 100;
01084
01085
01086
01087
01088
01089
01090
01091 static Array< Index > hspec(100,missing);
01092
01093
01094
01095 static Array< ArrayOfIndex > hiso(100);
01096
01097
01098 static bool hinit = false;
01099
01100
01101
01102 static ArrayOfIndex warned_missing;
01103
01104 if ( !hinit )
01105 {
01106 for ( Index i=0; i<species_data.nelem(); ++i )
01107 {
01108 const SpeciesRecord& sr = species_data[i];
01109
01110
01111
01112
01113 if ( 0 < sr.Isotope()[0].MytranTag() )
01114 {
01115
01116
01117
01118
01119
01120
01121
01122 Index mo = sr.Isotope()[0].MytranTag() / 10;
01123
01124 hspec[mo] = i;
01125
01126
01127 Index n_iso = sr.Isotope().nelem();
01128 ArrayOfIndex iso_tags;
01129 iso_tags.resize(n_iso);
01130 for ( Index j=0; j<n_iso; ++j )
01131 {
01132 iso_tags[j] = sr.Isotope()[j].MytranTag();
01133 }
01134
01135
01136
01137
01138
01139
01140
01141
01142 hiso[mo].resize( max(iso_tags)%10 + 1 );
01143 hiso[mo] = missing;
01144
01145
01146 for ( Index j=0; j<n_iso; ++j )
01147 {
01148 if ( 0 < iso_tags[j] )
01149 {
01150
01151
01152
01153 hiso[mo][iso_tags[j] % 10] = j;
01154 }
01155 }
01156 }
01157 }
01158
01159
01160
01161
01162
01163 out3 << " MYTRAN index table:\n";
01164 for ( Index i=0; i<hspec.nelem(); ++i )
01165 {
01166 if ( missing != hspec[i] )
01167 {
01168
01169
01170
01171 out3 << " mo = " << i << " Species = "
01172 << setw(10) << setiosflags(ios::left)
01173 << species_data[hspec[i]].Name().c_str()
01174 << "iso = ";
01175 for ( Index j=1; j<hiso[i].nelem(); ++j )
01176 {
01177 if ( missing==hiso[i][j] )
01178 out3 << " " << "m";
01179 else
01180 out3 << " " << species_data[hspec[i]].Isotope()[hiso[i][j]].Name();
01181 }
01182 out3 << "\n";
01183 }
01184 }
01185
01186 hinit = true;
01187 }
01188
01189
01190
01191
01192
01193 String line;
01194
01195
01196 Index mo;
01197
01198
01199 bool comment = true;
01200
01201 while (comment)
01202 {
01203
01204 if (is.eof()) return true;
01205
01206
01207 if (!is) throw runtime_error ("Stream bad.");
01208
01209
01210 getline(is,line);
01211
01212
01213
01214
01215
01216 if (line.nelem() == 0 && is.eof()) return true;
01217
01218
01219
01220
01221
01222 mo = 0;
01223
01224
01225 extract(mo,line,2);
01226
01227
01228
01229 if ( 0 != mo )
01230 {
01231
01232
01233 if ( missing != hspec[mo] ) comment = false ;
01234 else
01235 {
01236
01237
01238 if ( 0 == std::count(warned_missing.begin(),
01239 warned_missing.end(),
01240 mo) )
01241 {
01242 out0 << "Error: MYTRAN mo = " << mo << " is not "
01243 << "known to ARTS.\n";
01244 warned_missing.push_back(mo);
01245 exit(1);
01246
01247
01248
01249 }
01250 }
01251 }
01252 }
01253
01254
01255
01256
01257 mspecies = hspec[mo];
01258
01259
01260 Index iso;
01261 extract(iso,line,1);
01262
01263
01264
01265
01266
01267
01268
01269 misotope = missing;
01270 if ( iso < hiso[mo].nelem() )
01271 if ( missing != hiso[mo][iso] )
01272 misotope = hiso[mo][iso];
01273
01274
01275 if (missing == misotope)
01276 {
01277 ostringstream os;
01278 os << "Species: " << species_data[mspecies].Name()
01279 << ", isotope iso = " << iso
01280 << " is unknown.";
01281 throw runtime_error(os.str());
01282 }
01283
01284
01285
01286 {
01287
01288 Numeric v;
01289
01290
01291 extract(v,line,13);
01292
01293
01294 mf = v * 1E6;
01295
01296 }
01297
01298
01299 {
01300
01301 Numeric df;
01302 extract(df,line,8);
01303
01304 mdf = df * 1E6;
01305 }
01306
01307
01308 {
01309 extern const Numeric SPEED_OF_LIGHT;
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320 const Numeric hi2arts = 1e-2 * SPEED_OF_LIGHT;
01321
01322 Numeric s;
01323
01324
01325 extract(s,line,10);
01326
01327
01328 mi0 = s * hi2arts;
01329 }
01330
01331
01332
01333 {
01334
01335
01336 Numeric gam;
01337
01338
01339 extern const Numeric TORR2PA;
01340
01341
01342 extract(gam,line,5);
01343
01344
01345 magam = gam * 1E6 / TORR2PA;
01346
01347
01348 extract(gam,line,5);
01349
01350
01351 msgam = gam * 1E6 / TORR2PA;
01352
01353
01354 }
01355
01356
01357
01358 {
01359
01360
01361
01362
01363 extract(melow,line,10);
01364
01365
01366 melow = wavenumber_to_joule(melow);
01367 }
01368
01369
01370
01371 {
01372
01373 extract(mnair,line,4);
01374
01375
01376 extract(mnself,line,4);
01377
01378 }
01379
01380
01381
01382 {
01383
01384 extract(mtgam,line,7);
01385 }
01386
01387
01388
01389 {
01390
01391 Numeric d;
01392
01393
01394 extern const Numeric TORR2PA;
01395
01396
01397 extract(d,line,9);
01398
01399
01400 mpsf = d * 1E6 / TORR2PA;
01401 }
01402
01403
01404
01405
01406 {
01407 Index eu;
01408 extract(eu,line,3);
01409 }
01410
01411
01412 {
01413 Index el;
01414 extract(el,line,3);
01415 }
01416
01417
01418 {
01419 Index eul;
01420 extract(eul,line,9);
01421 }
01422
01423
01424 {
01425 Index ell;
01426 extract(ell,line,9);
01427 }
01428
01429 {
01430 Index di0;
01431
01432 extract(di0,line,1);
01433
01434 convMytranIER(mdi0,di0);
01435 }
01436
01437
01438 {
01439 Index dgam;
01440
01441 extract(dgam,line,1);
01442
01443 convMytranIER(mdagam,dgam);
01444 }
01445
01446
01447 {
01448 Index dnair;
01449
01450 extract(dnair,line,1);
01451
01452 convMytranIER(mdnair,dnair);
01453 }
01454
01455
01456
01457
01458
01459
01460
01461
01462 mti0 = 296.0;
01463
01464
01465
01466
01467
01468
01469
01470 return false;
01471 }
01472
01473
01474 bool LineRecord::ReadFromJplStream(istream& is)
01475 {
01476
01477 extern Array<SpeciesRecord> species_data;
01478
01479
01480
01481
01482
01483 static map<Index, SpecIsoMap> JplMap;
01484
01485
01486 static bool hinit = false;
01487
01488 if ( !hinit )
01489 {
01490
01491 out3 << " JPL index table:\n";
01492
01493 for ( Index i=0; i<species_data.nelem(); ++i )
01494 {
01495 const SpeciesRecord& sr = species_data[i];
01496
01497
01498 for ( Index j=0; j<sr.Isotope().nelem(); ++j)
01499 {
01500
01501 for (Index k=0; k<sr.Isotope()[j].JplTags().nelem(); ++k)
01502 {
01503
01504 SpecIsoMap indicies(i,j);
01505
01506 JplMap[sr.Isotope()[j].JplTags()[k]] = indicies;
01507
01508
01509
01510
01511
01512 const Index& i1 = JplMap[sr.Isotope()[j].JplTags()[k]].Speciesindex();
01513 const Index& i2 = JplMap[sr.Isotope()[j].JplTags()[k]].Isotopeindex();
01514
01515 out3 << " JPL TAG = " << sr.Isotope()[j].JplTags()[k] << " Species = "
01516 << setw(10) << setiosflags(ios::left)
01517 << species_data[i1].Name().c_str()
01518 << "iso = "
01519 << species_data[i1].Isotope()[i2].Name().c_str()
01520 << "\n";
01521 }
01522
01523 }
01524 }
01525 hinit = true;
01526 }
01527
01528
01529
01530
01531
01532 String line;
01533
01534
01535
01536 bool comment = true;
01537
01538 while (comment)
01539 {
01540
01541 if (is.eof()) return true;
01542
01543
01544 if (!is) throw runtime_error ("Stream bad.");
01545
01546
01547 getline(is,line);
01548
01549
01550
01551
01552
01553 if (line.nelem() == 0 && is.eof()) return true;
01554
01555
01556
01557
01558
01559
01560
01561
01562 Numeric v = 0.0;
01563
01564
01565 extract(v,line,13);
01566
01567
01568 if (v != 0.0)
01569 {
01570
01571 mf = v * 1E6;
01572
01573 comment = false;
01574 }
01575 }
01576
01577
01578 {
01579 Numeric df;
01580 extract(df,line,8);
01581
01582 mdf = df * 1E6;
01583 }
01584
01585
01586 {
01587
01588
01589
01590
01591
01592
01593
01594
01595 Numeric s;
01596
01597
01598 extract(s,line,8);
01599
01600
01601 s = pow(10,s);
01602
01603
01604 mi0 = s / 1E12;
01605 }
01606
01607
01608 {
01609 Index dr;
01610
01611
01612 extract(dr,line,2);
01613 }
01614
01615
01616 {
01617
01618
01619
01620
01621 extract(melow,line,10);
01622
01623
01624 melow = wavenumber_to_joule(melow);
01625 }
01626
01627
01628 {
01629 Index gup;
01630
01631
01632 extract(gup,line,3);
01633 }
01634
01635
01636 Index tag;
01637 {
01638
01639 extract(tag,line,7);
01640
01641
01642 tag = tag > 0 ? tag : -tag;
01643 }
01644
01645
01646
01647
01648 const map<Index, SpecIsoMap>::const_iterator i = JplMap.find(tag);
01649 if ( i == JplMap.end() )
01650 {
01651 ostringstream os;
01652 os << "JPL Tag: " << tag << " is unknown.";
01653 throw runtime_error(os.str());
01654 }
01655
01656 SpecIsoMap id = i->second;
01657
01658
01659
01660 mspecies = id.Speciesindex();
01661
01662
01663 misotope = id.Isotopeindex();
01664
01665
01666
01667
01668
01669
01670
01671 {
01672
01673 magam = 2.5E4;
01674
01675
01676 msgam = magam;
01677 }
01678
01679
01680
01681
01682
01683
01684
01685 {
01686 mnair = 0.75;
01687 mnself = 0.0;
01688 }
01689
01690
01691
01692
01693
01694 {
01695 mtgam = 300.0;
01696 }
01697
01698
01699
01700 {
01701 mpsf = 0.0;
01702 }
01703
01704
01705
01706
01707
01708
01709
01710
01711 mti0 = 300.0;
01712
01713
01714 return false;
01715 }
01716
01717
01718 bool LineRecord::ReadFromArtsStream(istream& is)
01719 {
01720
01721 extern Array<SpeciesRecord> species_data;
01722
01723
01724
01725
01726 static map<String, SpecIsoMap> ArtsMap;
01727
01728
01729 static bool hinit = false;
01730
01731 if ( !hinit )
01732 {
01733
01734 out3 << " ARTS index table:\n";
01735
01736 for ( Index i=0; i<species_data.nelem(); ++i )
01737 {
01738 const SpeciesRecord& sr = species_data[i];
01739
01740
01741 for ( Index j=0; j<sr.Isotope().nelem(); ++j)
01742 {
01743
01744 SpecIsoMap indicies(i,j);
01745 String buf = sr.Name()+"-"+sr.Isotope()[j].Name();
01746
01747 ArtsMap[buf] = indicies;
01748
01749
01750
01751
01752
01753 const Index& i1 = ArtsMap[buf].Speciesindex();
01754 const Index& i2 = ArtsMap[buf].Isotopeindex();
01755
01756 out3 << " Arts Identifier = " << buf << " Species = "
01757 << setw(10) << setiosflags(ios::left)
01758 << species_data[i1].Name().c_str()
01759 << "iso = "
01760 << species_data[i1].Isotope()[i2].Name().c_str()
01761 << "\n";
01762
01763 }
01764 }
01765 hinit = true;
01766 }
01767
01768
01769
01770
01771
01772 String line;
01773
01774
01775 bool comment = true;
01776
01777 while (comment)
01778 {
01779
01780 if (is.eof()) return true;
01781
01782
01783 if (!is) throw runtime_error ("Stream bad.");
01784
01785
01786 getline(is,line);
01787
01788
01789
01790
01791
01792 if (line.nelem() == 0 && is.eof()) return true;
01793
01794
01795 char c;
01796 extract(c,line,1);
01797
01798
01799 if (c == '@')
01800 {
01801 comment = false;
01802 }
01803 }
01804
01805
01806
01807 istringstream icecream(line);
01808
01809 String artsid;
01810 icecream >> artsid;
01811
01812 if (artsid.length() != 0)
01813 {
01814
01815
01816
01817 const map<String, SpecIsoMap>::const_iterator i = ArtsMap.find(artsid);
01818 if ( i == ArtsMap.end() )
01819 {
01820 ostringstream os;
01821 os << "ARTS Tag: " << artsid << " is unknown.";
01822 throw runtime_error(os.str());
01823 }
01824
01825 SpecIsoMap id = i->second;
01826
01827
01828
01829 mspecies = id.Speciesindex();
01830
01831
01832 misotope = id.Isotopeindex();
01833
01834
01835
01836 icecream >> mf;
01837
01838
01839
01840 icecream >> mpsf;
01841
01842
01843 icecream >> mi0;
01844
01845
01846
01847 icecream >> mti0;
01848
01849
01850
01851 icecream >> melow;
01852
01853
01854
01855 icecream >> magam;
01856 icecream >> msgam;
01857
01858
01859 icecream >> mnair;
01860 icecream >> mnself;
01861
01862
01863
01864 icecream >> mtgam;
01865
01866
01867 Index naux;
01868 icecream >> naux;
01869
01870
01871 maux.resize(naux);
01872
01873 for (Index i = 0; i<naux; i++)
01874 {
01875 icecream >> maux[i];
01876
01877 }
01878
01879
01880 try
01881 {
01882 icecream >> mdf;
01883 icecream >> mdi0;
01884 icecream >> mdagam;
01885 icecream >> mdsgam;
01886 icecream >> mdnair;
01887 icecream >> mdnself;
01888 icecream >> mdpsf;
01889 }
01890 catch (runtime_error x)
01891 {
01892
01893
01894
01895 mdf = -1;
01896 mdi0 = -1;
01897 mdagam = -1;
01898 mdsgam = -1;
01899 mdnair = -1;
01900 mdnself = -1;
01901 mdpsf = -1;
01902 }
01903 }
01904
01905
01906 return false;
01907 }
01908
01909
01910
01911 OneTag::OneTag(String def)
01912 {
01913
01914 extern const Array<SpeciesRecord> species_data;
01915
01916 extern std::map<String, Index> SpeciesMap;
01917
01918 String name, isoname;
01919
01920 Index n;
01921
01922
01923
01924 mlf = -1;
01925 muf = -1;
01926
01927
01928
01929
01930
01931
01932
01933 n = def.find('-');
01934 if (n != def.npos )
01935 {
01936 name = def.substr(0,n);
01937 def.erase(0,n+1);
01938 }
01939 else
01940 {
01941
01942
01943
01944 name = def;
01945 def = "";
01946 }
01947
01948
01949
01950
01951
01952 map<String, Index>::const_iterator mi = SpeciesMap.find(name);
01953 if ( mi != SpeciesMap.end() )
01954 {
01955
01956 mspecies = mi->second;
01957 }
01958 else
01959 {
01960 ostringstream os;
01961 os << "Species " << name << " is not a valid species.";
01962 throw runtime_error(os.str());
01963 }
01964
01965
01966 const SpeciesRecord& spr = species_data[mspecies];
01967
01968 if ( 0 == def.nelem() )
01969 {
01970
01971
01972
01973 misotope = spr.Isotope().nelem();
01974
01975
01976 return;
01977 }
01978
01979
01980 n = def.find('-');
01981 if (n != def.npos )
01982 {
01983 isoname = def.substr(0,n);
01984 def.erase(0,n+1);
01985 }
01986 else
01987 {
01988
01989
01990
01991 isoname = def;
01992 def = "";
01993 }
01994
01995
01996 if ( "*" == isoname )
01997 {
01998
01999 misotope = spr.Isotope().nelem();
02000 }
02001 else
02002 {
02003
02004 ArrayOfString ins;
02005 for ( Index i=0; i<spr.Isotope().nelem(); ++i )
02006 ins.push_back( spr.Isotope()[i].Name() );
02007
02008
02009
02010
02011 misotope = find( ins.begin(),
02012 ins.end(),
02013 isoname ) - ins.begin();
02014
02015
02016 if ( misotope >= ins.nelem() )
02017 {
02018 ostringstream os;
02019 os << "Isotope " << isoname << " is not a valid isotope for "
02020 << "species " << name << ".\n"
02021 << "Valid isotopes are:";
02022 for ( Index i=0; i<ins.nelem(); ++i )
02023 os << " " << ins[i];
02024 throw runtime_error(os.str());
02025 }
02026 }
02027
02028 if ( 0 == def.nelem() )
02029 {
02030
02031
02032
02033
02034 return;
02035 }
02036
02037
02038
02039
02040
02041 n = def.find('-');
02042 if (n != def.npos )
02043 {
02044
02045 String fname;
02046 fname = def.substr(0,n);
02047 def.erase(0,n+1);
02048
02049
02050 if ( "*" == fname )
02051 {
02052
02053
02054 }
02055 else
02056 {
02057
02058 istringstream is(fname);
02059 is >> mlf;
02060 }
02061 }
02062 else
02063 {
02064
02065
02066 throw runtime_error("You must either speciefy both frequency limits\n"
02067 "(at least with jokers), or none.");
02068 }
02069
02070
02071
02072
02073 if ( "*" == def )
02074 {
02075
02076
02077 }
02078 else
02079 {
02080
02081 istringstream is(def);
02082 is >> muf;
02083 }
02084 }
02085
02086
02087 String OneTag::Name() const
02088 {
02089
02090 extern const Array<SpeciesRecord> species_data;
02091
02092 const SpeciesRecord& spr = species_data[mspecies];
02093
02094 ostringstream os;
02095
02096
02097 os << spr.Name() << "-";
02098
02099
02100 if ( misotope == spr.Isotope().nelem() )
02101 {
02102
02103 os << "*-";
02104 }
02105 else
02106 {
02107 os << spr.Isotope()[misotope].Name() << "-";
02108 }
02109
02110
02111
02112
02113
02114
02115 Index precision;
02116 switch (sizeof(Numeric)) {
02117 case sizeof(float) : precision = FLT_DIG; break;
02118 case sizeof(double) : precision = DBL_DIG; break;
02119 default: out0 << "Numeric must be double or float\n"; exit(1);
02120 }
02121
02122 if ( 0 > mlf )
02123 {
02124
02125 os << "*-";
02126 }
02127 else
02128 {
02129 os << setprecision(precision);
02130 os << mlf << "-";
02131 }
02132
02133 if ( 0 > muf )
02134 {
02135
02136 os << "*";
02137 }
02138 else
02139 {
02140 os << setprecision(precision);
02141 os << muf;
02142 }
02143
02144 return os.str();
02145 }
02146
02147 ostream& operator << (ostream& os, const OneTag& ot)
02148 {
02149 return os << ot.Name();
02150 }
02151
02152
02153
02176 void get_tagindex_for_Strings(
02177 ArrayOfIndex& tags1_index,
02178 const TagGroups& tags1,
02179 const ArrayOfString& tags2_Strings )
02180 {
02181 const Index n1 = tags1.nelem();
02182 const Index n2 = tags2_Strings.nelem();
02183 TagGroups tags2;
02184 Index i1, i2, nj, j, found, ok;
02185
02186 tags1_index.resize(n2);
02187
02188 tgsDefine( tags2, tags2_Strings );
02189
02190 for ( i2=0; i2<n2; i2++ )
02191 {
02192 found = 0;
02193 for ( i1=0; (i1<n1) && !found; i1++ )
02194 {
02195 nj = tags2[i2].nelem();
02196 if ( nj == tags1[i1].nelem() )
02197 {
02198 ok = 1;
02199 for ( j=0; j<nj; j++ )
02200 {
02201 if ( tags2[i2][j].Name() != tags1[i1][j].Name() )
02202 ok = 0;
02203 }
02204 if ( ok )
02205 {
02206 found = 1;
02207 tags1_index[i2] = i1;
02208 }
02209 }
02210 }
02211 if ( !found )
02212 {
02213 ostringstream os;
02214 os << "The tag String \"" << tags2_Strings[i2] <<
02215 "\" does not match any of the given tags.\n";
02216 throw runtime_error(os.str());
02217 }
02218 }
02219 }
02220
02221
02233 void get_tag_group_index_for_tag_group( Index& tgs1_index,
02234 const TagGroups& tgs1,
02235 const Array<OneTag>& tg2 )
02236 {
02237 bool found = false;
02238
02239 for ( Index i=0;
02240 i<tgs1.nelem() && !found;
02241 ++i )
02242 {
02243
02244 if ( tg2.nelem() == tgs1[i].nelem() )
02245 {
02246 bool ok = true;
02247
02248 for ( Index j=0; j<tg2.nelem(); ++j )
02249 {
02250 if ( tg2[j].Name() != tgs1[i][j].Name() )
02251 ok = false;
02252 }
02253
02254 if ( ok )
02255 {
02256 found = true;
02257 tgs1_index = i;
02258 }
02259 }
02260 }
02261
02262 if ( !found )
02263 {
02264 ostringstream os;
02265 os << "The tag String \"" << tg2 <<
02266 "\" does not match any of the given tags.\n";
02267 throw runtime_error(os.str());
02268 }
02269 }
02270
02271
02286 String get_tag_group_name( const Array<OneTag>& tg )
02287 {
02288 String name;
02289 Index i;
02290
02291 for ( i=0; i<tg.nelem()-1; ++i )
02292 {
02293 name += tg[i].Name() + ", ";
02294 }
02295 name += tg[i].Name();
02296
02297 return name;
02298 }
02299
02300
02309 void write_lines_to_stream(ostream& os,
02310 const ArrayOfLineRecord& lines)
02311 {
02312
02313
02314 LineRecord dummy;
02315
02316
02317 os << dummy.Version() << " " << lines.nelem() << "\n";
02318
02319
02320 for ( Index i=0; i<lines.nelem(); ++i )
02321 {
02322
02323 os << lines[i] << "\n";
02324 }
02325 }
02326
02327
02329
02340 bool is_sorted( ConstVectorView x )
02341 {
02342 if( x.nelem() > 1 )
02343 {
02344 for( Index i=1; i<x.nelem(); i++ )
02345 {
02346 if( x[i] < x[i-1] )
02347 return false;
02348 }
02349 }
02350 return true;
02351 }
02352
02353
02379 void xsec_species( MatrixView xsec,
02380 ConstVectorView f_mono,
02381 ConstVectorView p_abs,
02382 ConstVectorView t_abs,
02383 ConstVectorView h2o_abs,
02384 ConstVectorView vmr,
02385 const ArrayOfLineRecord& lines,
02386 const Index ind_ls,
02387 const Index ind_lsn,
02388 const Numeric cutoff)
02389 {
02390
02391 extern const Array<LineshapeRecord> lineshape_data;
02392 extern const Array<LineshapeNormRecord> lineshape_norm_data;
02393
02394
02395 extern const Numeric SPEED_OF_LIGHT;
02396
02397
02398 extern const Numeric BOLTZMAN_CONST;
02399
02400
02401 extern const Numeric AVOGADROS_NUMB;
02402
02403
02404 extern const Numeric PLANCK_CONST;
02405
02406
02407
02408
02409
02410 static const Numeric doppler_const = sqrt( 2.0 * BOLTZMAN_CONST *
02411 AVOGADROS_NUMB) / SPEED_OF_LIGHT;
02412
02413
02414 Index nf = f_mono.nelem();
02415 Index nl = lines.nelem();
02416
02417
02418
02419
02420
02421 Vector ls(nf+1);
02422 Vector fac(nf+1);
02423
02424 bool cut = (cutoff != -1) ? true : false;
02425
02426
02427
02428 if (cut)
02429 {
02430 if ( ! is_sorted( f_mono ) )
02431 {
02432 ostringstream os;
02433 os << "If you use a lineshape function with cutoff, your\n"
02434 << "frequency grid *f_mono* must be sorted.\n"
02435 << "(Duplicate values are allowed.)";
02436 throw runtime_error(os.str());
02437 }
02438 }
02439
02440
02441 bool negative = false;
02442
02443 for (Index i = 0; !negative && i < t_abs.nelem (); i++)
02444 {
02445 if (t_abs[i] < 0.)
02446 negative = true;
02447 }
02448
02449 if (negative)
02450 {
02451 ostringstream os;
02452 os << "t_abs contains at least one negative temperature value.\n"
02453 << "This is not allowed.";
02454 throw runtime_error(os.str());
02455 }
02456
02457
02458
02459
02460 Vector f_local( nf + 1 );
02461
02462
02463
02464
02465
02466
02467
02468
02469
02470
02471 Index ii = (nf+1 < 10) ? 10 : nf+1;
02472 Vector aux(ii);
02473
02474
02475
02476
02477
02478
02479 if ( t_abs.nelem() != p_abs.nelem() )
02480 {
02481 ostringstream os;
02482 os << "Variable t_abs must have the same dimension as p_abs.\n"
02483 << "t_abs.nelem() = " << t_abs.nelem() << '\n'
02484 << "p_abs.nelem() = " << p_abs.nelem();
02485 throw runtime_error(os.str());
02486 }
02487
02488 if ( vmr.nelem() != p_abs.nelem() )
02489 {
02490 ostringstream os;
02491 os << "Variable vmr must have the same dimension as p_abs.\n"
02492 << "vmr.nelem() = " << vmr.nelem() << '\n'
02493 << "p_abs.nelem() = " << p_abs.nelem();
02494 throw runtime_error(os.str());
02495 }
02496
02497 if ( h2o_abs.nelem() != p_abs.nelem() )
02498 {
02499 ostringstream os;
02500 os << "Variable h2o_abs must have the same dimension as p_abs.\n"
02501 << "h2o_abs.nelem() = " << h2o_abs.nelem() << '\n'
02502 << "p_abs.nelem() = " << p_abs.nelem();
02503 throw runtime_error(os.str());
02504 }
02505
02506
02507
02508
02509 if ( xsec.nrows() != nf || xsec.ncols() != p_abs.nelem() )
02510 {
02511 ostringstream os;
02512 os << "Variable xsec must have dimensions [f_mono.nelem(),p_abs.nelem()].\n"
02513 << "[xsec.nrows(),xsec.ncols()] = [" << xsec.nrows()
02514 << ", " << xsec.ncols() << "]\n"
02515 << "f_mono.nelem() = " << nf << '\n'
02516 << "p_abs.nelem() = " << p_abs.nelem();
02517 throw runtime_error(os.str());
02518 }
02519
02520
02521
02522 for ( Index i=0; i<p_abs.nelem(); ++i )
02523 {
02524
02525
02526
02527 Numeric p_i = p_abs[i];
02528 Numeric t_i = t_abs[i];
02529
02530
02531
02532
02533
02534
02535 Numeric n = p_i / BOLTZMAN_CONST / t_i;
02536
02537
02538 const Numeric p_partial = p_i * vmr[i];
02539
02540
02541
02542 for ( Index l=0; l< nl; ++l )
02543 {
02544
02545
02546
02547
02548 f_local[Range(0,nf)] = f_mono;
02549
02550
02551
02552 Index nfl = nf;
02553
02554
02555
02556 Index nfls = nf;
02557
02558
02559 Numeric base=0.0;
02560
02561
02562
02563 LineRecord l_l = lines[l];
02564
02565 Numeric F0 = l_l.F();
02566
02567
02568
02569
02570
02571 Numeric intensity = l_l.I0();
02572
02573
02574 Numeric e_lower = l_l.Elow();
02575
02576
02577 Numeric e_upper = e_lower + F0 * PLANCK_CONST;
02578
02579
02580
02581
02582
02583
02584
02585 Numeric part_fct_ratio =
02586 l_l.IsotopeData().CalculatePartitionFctRatio( l_l.Ti0(),
02587 t_i );
02588
02589
02590 Numeric nom = exp(- e_lower / ( BOLTZMAN_CONST * t_i ) ) -
02591 exp(- e_upper / ( BOLTZMAN_CONST * t_i ) );
02592
02593 Numeric denom = exp(- e_lower / ( BOLTZMAN_CONST * l_l.Ti0() ) ) -
02594 exp(- e_upper / ( BOLTZMAN_CONST * l_l.Ti0() ) );
02595
02596
02597
02598
02599
02600 intensity *= part_fct_ratio * nom / denom;
02601
02602 if (lineshape_norm_data[ind_lsn].Name() == "quadratic")
02603 {
02604
02605
02606
02607
02608
02609
02610 Numeric mafac = (PLANCK_CONST * F0) / (2.000e0 * BOLTZMAN_CONST * t_i);
02611 intensity = intensity * mafac / sinh(mafac);
02612 }
02613
02614
02615
02616
02617 const Numeric theta = l_l.Tgam() / t_i;
02618 const Numeric theta_Nair = pow(theta, l_l.Nair());
02619
02620 Numeric gamma
02621 = l_l.Agam() * theta_Nair * (p_i - p_partial)
02622 + l_l.Sgam() * pow(theta, l_l.Nself()) * p_partial;
02623
02624
02625
02626 Numeric sigma = F0 * doppler_const *
02627 sqrt( t_i / l_l.IsotopeData().Mass());
02628
02629
02630
02631
02632
02633 F0 += l_l.Psf() * p_i *
02634 std::pow( theta , (Numeric).25 + (Numeric)1.5*l_l.Nair() );
02635
02636
02637
02638
02639
02640
02641
02642
02643
02644
02645
02646
02647 aux[0] = theta;
02648 aux[1] = theta_Nair;
02649 aux[2] = p_i;
02650 aux[3] = vmr[i];
02651 aux[4] = h2o_abs[i];
02652 aux[5] = l_l.Agam();
02653 aux[6] = l_l.Nair();
02654
02655 if (l_l.Naux() > 1)
02656 {
02657 aux[7] = l_l.Aux()[0];
02658 aux[8] = l_l.Aux()[1];
02659
02660 }
02661 else
02662 {
02663 aux[7] = 0.0;
02664 aux[8] = 0.0;
02665 }
02666
02667
02668
02669
02670 Index i_f_min = 0;
02671 Index i_f_max = nf-1;
02672
02673
02674 if ( cut )
02675 {
02676
02677
02678
02679
02680
02681 while ( i_f_min < nf && (F0 - cutoff) > f_mono[i_f_min] )
02682 {
02683
02684 ++i_f_min;
02685 }
02686
02687
02688
02689
02690
02691
02692
02693 while ( i_f_max >= 0 && (F0 + cutoff) < f_mono[i_f_max] )
02694 {
02695
02696 --i_f_max;
02697 }
02698
02699
02700 ++i_f_max;
02701 f_local[i_f_max] = F0 + cutoff;
02702
02703
02704 nfls = i_f_max - i_f_min + 1;
02705
02706
02707
02708
02709 nfl = nfls -1;
02710 }
02711 else
02712 {
02713
02714 }
02715
02716
02717
02718
02719
02720
02721
02722 if ( nfl > 0 )
02723 {
02724
02725 lineshape_data[ind_ls].Function()(ls,
02726 aux,F0,gamma,sigma,
02727 f_local[Range(i_f_min,nfls)],
02728 nfls);
02729
02730
02731 lineshape_norm_data[ind_lsn].Function()(fac,F0,
02732 f_local[Range(i_f_min,nfls)],
02733 t_i,nfls);
02734
02735
02736
02737 VectorView this_xsec = xsec(Range(i_f_min,nfl),i);
02738
02739
02740 VectorView this_ls = ls[Range(0,nfl)];
02741 VectorView this_fac = fac[Range(0,nfl)];
02742
02743
02744 if ( cut )
02745 {
02746
02747
02748 base = ls[nfls-1];
02749
02750
02751
02752 this_ls -= base;
02753 }
02754
02755
02756 {
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766 const Numeric factors = n * intensity * l_l.IsotopeData().Abundance();
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777 this_ls *= this_fac;
02778 this_ls *= factors;
02779 this_xsec += this_ls;
02780 }
02781 }
02782 }
02783 }
02784 }
02785
02786
02787
02788
02789
02790
02791
02792
02793
02795
02810 void refr_index_BoudourisDryAir (
02811 Vector& refr_index,
02812 ConstVectorView p_abs,
02813 ConstVectorView t_abs )
02814 {
02815 const Index n = p_abs.nelem();
02816 refr_index.resize( n );
02817
02818 assert ( n == t_abs.nelem() );
02819
02820
02821 for ( Index i=0; i<n; i++ )
02822 refr_index[i] = 1.0 + 77.593e-8 * p_abs[i] / t_abs[i];
02823 }
02824
02825
02826
02828
02842 void refr_index_Boudouris (
02843 Vector& refr_index,
02844 ConstVectorView p_abs,
02845 ConstVectorView t_abs,
02846 ConstVectorView h2o_abs )
02847 {
02848 const Index n = p_abs.nelem();
02849 refr_index.resize( n );
02850
02851 assert ( n == t_abs.nelem() );
02852 assert ( n == h2o_abs.nelem() );
02853
02854 Numeric e;
02855 Numeric p;
02856
02857 for ( Index i=0; i<n; i++ )
02858 {
02859 e = p_abs[i] * h2o_abs[i];
02860 p = p_abs[i] - e;
02861
02862 refr_index[i] = 1.0 + 77.593e-8 * p / t_abs[i] +
02863 72e-8 * e / t_abs[i] +
02864 3.754e-3 * e / (t_abs[i]*t_abs[i]);
02865 }
02866 }
02867
02879 Numeric wavenumber_to_joule(Numeric e)
02880 {
02881
02882 extern const Numeric PLANCK_CONST;
02883
02884
02885 extern const Numeric SPEED_OF_LIGHT;
02886
02887
02888 static const Numeric lower_energy_const = PLANCK_CONST * SPEED_OF_LIGHT * 1E2;
02889
02890 return e*lower_energy_const;
02891 }
02892
02893
02894
02895
02896
02897
02898
02899
02900
02901 void convHitranIERF(
02902 Numeric& mdf,
02903 const Index& df
02904 )
02905 {
02906 switch ( df )
02907 {
02908 case 0:
02909 {
02910 mdf = -1;
02911 break;
02912 }
02913 case 1:
02914 {
02915 mdf = 30*1E9;
02916 break;
02917 }
02918 case 2:
02919 {
02920 mdf = 3*1E9;
02921 break;
02922 }
02923 case 3:
02924 {
02925 mdf = 300*1E6;
02926 break;
02927 }
02928 case 4:
02929 {
02930 mdf = 30*1E6;
02931 break;
02932 }
02933 case 5:
02934 {
02935 mdf = 3*1E6;
02936 break;
02937 }
02938 case 6:
02939 {
02940 mdf = 0.3*1E6;
02941 break;
02942 }
02943 }
02944 }
02945
02946
02947
02948 void convHitranIERSH(
02949 Numeric& mdh,
02950 const Index& dh
02951 )
02952 {
02953 switch ( dh )
02954 {
02955 case 0:
02956 {
02957 mdh = -1;
02958 break;
02959 }
02960 case 1:
02961 {
02962 mdh = -1;
02963 break;
02964 }
02965 case 2:
02966 {
02967 mdh = -1;
02968 break;
02969 }
02970 case 3:
02971 {
02972 mdh = 30;
02973 break;
02974 }
02975 case 4:
02976 {
02977 mdh = 20;
02978 break;
02979 }
02980 case 5:
02981 {
02982 mdh = 10;
02983 break;
02984 }
02985 case 6:
02986 {
02987 mdh =5;
02988 break;
02989 }
02990 case 7:
02991 {
02992 mdh =2;
02993 break;
02994 }
02995 case 8:
02996 {
02997 mdh =1;
02998 break;
02999 }
03000 }
03001 mdh=mdh/100;
03002 }
03003
03004
03005
03006
03007 void convMytranIER(
03008 Numeric& mdh,
03009 const Index & dh
03010 )
03011 {
03012 switch ( dh )
03013 {
03014 case 0:
03015 {
03016 mdh = 200;
03017 break;
03018 }
03019 case 1:
03020 {
03021 mdh = 100;
03022 break;
03023 }
03024 case 2:
03025 {
03026 mdh = 50;
03027 break;
03028 }
03029 case 3:
03030 {
03031 mdh = 30;
03032 break;
03033 }
03034 case 4:
03035 {
03036 mdh = 20;
03037 break;
03038 }
03039 case 5:
03040 {
03041 mdh = 10;
03042 break;
03043 }
03044 case 6:
03045 {
03046 mdh =5;
03047 break;
03048 }
03049 case 7:
03050 {
03051 mdh = 2;
03052 break;
03053 }
03054 case 8:
03055 {
03056 mdh = 1;
03057 break;
03058 }
03059 case 9:
03060 {
03061 mdh = 0.5;
03062 break;
03063 }
03064 }
03065 mdh=mdh/100;
03066 }
03067