00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00023
00025
00037
00038
00040
00041 #include <math.h>
00042 #include <stdexcept>
00043 #include "arts.h"
00044 #include "matpackI.h"
00045 #include "messages.h"
00046 #include "math_funcs.h"
00047 #include "make_vector.h"
00048
00049 extern const Numeric DEG2RAD;
00050 extern const Numeric RAD2DEG;
00051 extern const Numeric PLANCK_CONST;
00052 extern const Numeric SPEED_OF_LIGHT;
00053 extern const Numeric BOLTZMAN_CONST;
00054
00055
00056
00058
00060
00062
00075 void planck (
00076 MatrixView B,
00077 ConstVectorView f,
00078 ConstVectorView t )
00079 {
00080
00081 static const double a = 2.0*PLANCK_CONST/(SPEED_OF_LIGHT*SPEED_OF_LIGHT);
00082 static const double b = PLANCK_CONST/BOLTZMAN_CONST;
00083
00084 const Index n_f = f.nelem();
00085 const Index n_t = t.nelem();
00086 Index i_f, i_t;
00087 Numeric c, d;
00088
00089 assert( n_f==B.nrows() );
00090 assert( n_t==B.ncols() );
00091
00092 for ( i_f=0; i_f<n_f; i_f++ )
00093 {
00094 c = a * f[i_f]*f[i_f]*f[i_f];
00095 d = b * f[i_f];
00096 for ( i_t=0; i_t<n_t; i_t++ )
00097 B(i_f,i_t) = c / (exp(d/t[i_t]) - 1.0);
00098 }
00099 }
00100
00101
00102
00104
00114 void planck (
00115 VectorView B,
00116 ConstVectorView f,
00117 Numeric t )
00118 {
00119
00120 static const double a = 2.0*PLANCK_CONST/(SPEED_OF_LIGHT*SPEED_OF_LIGHT);
00121 static const double b = PLANCK_CONST/BOLTZMAN_CONST;
00122 const double c = b/t;
00123
00124 assert( B.nelem()==f.nelem() );
00125
00126 for ( Index i=0; i<f.nelem(); i++ )
00127 {
00128 B[i] = a * f[i]*f[i]*f[i] / ( exp( f[i]*c ) - 1.0 );
00129 }
00130 }
00131
00132
00133
00135
00145 void invplanck (
00146 VectorView y,
00147 ConstVectorView f,
00148 ConstVectorView za )
00149 {
00150 const Index nf = f.nelem();
00151 const Index nza = za.nelem();
00152 const Index ny = y.nelem();
00153 Index i0;
00154
00155
00156 const double a = PLANCK_CONST/BOLTZMAN_CONST;
00157 const double b = 2*PLANCK_CONST/(SPEED_OF_LIGHT*SPEED_OF_LIGHT);
00158 double c,d;
00159
00160
00161 if ( max(y) > 1e-4 )
00162 throw runtime_error("The spectrum cannot be in expected intensity unit "
00163 "(impossible value detected).");
00164
00165 if ( nf*nza != ny )
00166 {
00167 ostringstream os;
00168 os << "The length of *y* does not match *f_mono* and *za_pencil*.\n"
00169 << "y.nelem(): " << y.nelem() << "\n"
00170 << "Should be f_mono.nelem()*za_pencil.nelem(): "
00171 << f.nelem() * za.nelem() << "\n"
00172 << "f_mono.nelem(): " << f.nelem() << "\n"
00173 << "za_pencil.nelem(): " << za.nelem();
00174 throw runtime_error(os.str());
00175 }
00176
00177 for ( Index i=0; i<nf; i++ )
00178 {
00179 c = a*f[i];
00180 d = b*f[i]*f[i]*f[i];
00181 for ( Index j=0; j<nza; j++ )
00182 {
00183 i0 = j*nf + i;
00184 y[i0] = c / ( log(d/y[i0]+1) );
00185 }
00186 }
00187 }
00188
00189
00190
00192
00202 void invrayjean (
00203 VectorView y,
00204 ConstVectorView f,
00205 ConstVectorView za )
00206 {
00207 const Index nf = f.nelem();
00208 const Index nza = za.nelem();
00209 const Index ny = y.nelem();
00210 Index i0;
00211
00212
00213
00214
00215
00216 double b;
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230 if ( nf*nza != ny )
00231 {
00232 ostringstream os;
00233 os << "The length of *y* does not match *f_mono* and *za_pencil*.\n"
00234 << "y.nelem(): " << y.nelem() << "\n"
00235 << "Should be f_mono.nelem()*za_pencil.nelem(): "
00236 << f.nelem() * za.nelem() << "\n"
00237 << "f_mono.nelem(): " << f.nelem() << "\n"
00238 << "za_pencil.nelem(): " << za.nelem();
00239 throw runtime_error(os.str());
00240 }
00241
00242 for ( Index i=0; i<nf; i++ )
00243 {
00244
00245
00246
00247 b = SPEED_OF_LIGHT / ((double)f[i] * (double)f[i])
00248 * SPEED_OF_LIGHT / (2 * BOLTZMAN_CONST);
00249 for ( Index j=0; j<nza; j++ )
00250 {
00251 i0 = j*nf + i;
00252 y[i0] = b * (double)y[i0];
00253 }
00254 }
00255 }
00256
00257
00258
00260
00270 Numeric number_density (
00271 Numeric p,
00272 Numeric t )
00273 {
00274 assert( 0!=t );
00275 return p/t/BOLTZMAN_CONST;
00276 }
00277
00278
00279
00281
00291 Vector number_density (
00292 ConstVectorView p,
00293 ConstVectorView t )
00294 {
00295 assert( p.nelem()==t.nelem() );
00296
00297
00298
00299 Vector dummy(p);
00300
00301
00302 dummy /= t;
00303 dummy /= BOLTZMAN_CONST;
00304
00305 return dummy;
00306 }
00307
00308
00309
00311
00322 Numeric g_of_z (
00323 Numeric r_geoid,
00324 Numeric g0,
00325 Numeric z )
00326 {
00327 return g0 * pow( r_geoid/(r_geoid+z), 2 );
00328 }
00329
00331
00340 Numeric g_of_lat (Numeric latitude)
00341 {
00342 extern const Numeric EARTH_EQUATORIAL_RADIUS;
00343 extern const Numeric EARTH_POLAR_RADIUS;
00344 extern const Numeric DEG2RAD;
00345
00346
00347 check_if_in_range( -90, 90, latitude, "latitude" );
00348
00349
00350 Numeric latG = atan( (EARTH_EQUATORIAL_RADIUS*EARTH_EQUATORIAL_RADIUS) /
00351 (EARTH_POLAR_RADIUS*EARTH_POLAR_RADIUS) * tan(latitude*DEG2RAD));
00352
00353
00354 return 9.7803267714 * ( ( 1.0 + 0.0019318514 * sin(latG) * sin(latG) ) /
00355 sqrt( 1.0 - 0.00669438 * sin(latG) * sin(latG) ) );
00356 }
00357
00358
00359
00360
00361
00363
00365
00367
00386 void rte_iterate (
00387 VectorView y,
00388 const Index start_index,
00389 const Index stop_index,
00390 ConstMatrixView tr,
00391 ConstMatrixView s,
00392 const Index n_f )
00393 {
00394 Index i_f;
00395 Index i_z;
00396 Index i_step;
00397
00398 if ( start_index >= stop_index )
00399 i_step = -1;
00400
00401 else
00402 i_step = 1;
00403
00404 for ( i_z=start_index; i_z!=(stop_index+i_step); i_z+=i_step )
00405 {
00406 for ( i_f=0; i_f<n_f; i_f++ )
00407 y[i_f] = y[i_f]*tr(i_f,i_z) + s(i_f,i_z) * ( 1.0-tr(i_f,i_z) );
00408 }
00409 }
00410
00411
00412
00414
00433 void rte (
00434 VectorView y,
00435 const Index start_index,
00436 const Index stop_index,
00437 ConstMatrixView tr,
00438 ConstMatrixView s,
00439 ConstVectorView y_space,
00440 const Index ground,
00441 ConstVectorView e_ground,
00442 ConstVectorView y_ground )
00443 {
00444 const Index n_f = tr.nrows();
00445 Index i_f;
00446 Index i_break;
00447 Index i_start;
00448
00449
00450 y = y_space;
00451
00452
00453
00454
00455 if ( start_index > 0 )
00456 {
00457
00458 if ( ground > 0 )
00459 i_break = ground-1;
00460 else
00461 i_break = 0;
00462
00463
00464 rte_iterate( y, start_index-1, i_break, tr, s, n_f );
00465
00466
00467
00468
00469 if ( !(stop_index==0 && ground==0) )
00470 {
00471
00472 i_start = 0;
00473 i_break = stop_index - 1;
00474
00475
00476
00477 if ( ground > 0 )
00478 {
00479 for ( i_f=0; i_f<n_f; i_f++ )
00480 y[i_f] = y[i_f]*(1.0-e_ground[i_f]) + y_ground[i_f]*e_ground[i_f];
00481
00482 if ( ground > 1 )
00483 {
00484 i_start = ground - 2;
00485 i_break = 0;
00486 }
00487 }
00488
00489
00490 rte_iterate( y, i_start, i_break, tr, s, n_f );
00491
00492 }
00493 }
00494 }
00495
00496
00497
00499
00518 void bl_iterate (
00519 VectorView y,
00520 const Index start_index,
00521 const Index stop_index,
00522 ConstMatrixView tr,
00523 const Index n_f )
00524 {
00525 Index i_f;
00526 Index i_z;
00527 Index i_step;
00528
00529 if ( start_index >= stop_index )
00530 i_step = -1;
00531 else
00532 i_step = 1;
00533
00534 for ( i_z=start_index; i_z!=(stop_index+i_step); i_z+=i_step )
00535 {
00536 for ( i_f=0; i_f<n_f; i_f++ )
00537 y[i_f] *= tr(i_f,i_z);
00538 }
00539 }
00540
00541
00542
00544
00560 void bl (
00561 Vector& y,
00562 const Index start_index,
00563 const Index stop_index,
00564 ConstMatrixView tr,
00565 const Index ground,
00566 ConstVectorView e_ground )
00567 {
00568 if ( start_index < stop_index )
00569 throw runtime_error("The start index cannot be "
00570 "smaller than the stop index." );
00571
00572 const Index nf = tr.nrows();
00573 Index iy;
00574
00575
00576 y.resize( nf );
00577 y = 1.0;
00578
00579
00580 if ( stop_index > 0 )
00581 {
00582 bl_iterate( y, 0, stop_index-1, tr, nf );
00583 y *= y;
00584 }
00585
00586
00587 if ( start_index != stop_index )
00588 bl_iterate( y, stop_index, start_index-1, tr, nf );
00589
00590
00591 if ( ground > 0 )
00592 {
00593 for ( iy=0; iy<nf; iy++ )
00594 y[iy] *= ( 1.0 - e_ground[iy] );
00595 }
00596 }
00597
00598
00599
00601
00603
00605
00621 void z2p(
00622 VectorView p,
00623 ConstVectorView z0,
00624 ConstVectorView p0,
00625 ConstVectorView z )
00626 {
00627 assert( p.nelem()==z.nelem() );
00628 if ( z.nelem() > 0 )
00629 {
00630
00631 Vector logp0(p0.nelem());
00632 transform( logp0, log, p0 );
00633
00634 interp_lin_vector( p, z0, logp0, z );
00635 transform( p, exp, p );
00636 }
00637 }
00638
00639
00640
00642
00658 void interpp(
00659 VectorView x,
00660 ConstVectorView p0,
00661 ConstVectorView x0,
00662 ConstVectorView p )
00663 {
00664 assert( x.nelem()==p.nelem() );
00665
00666
00667 Vector logp0(p0.nelem());
00668 Vector logp(p.nelem());
00669 transform( logp0, log, p0 );
00670 transform( logp, log, p );
00671
00672 interp_lin_vector( x, logp0, x0, logp );
00673 }
00674
00675 void interpp_cloud(
00676 VectorView x,
00677 ConstVectorView p0,
00678 ConstVectorView x0,
00679 ConstVectorView p )
00680 {
00681 assert( x.nelem()==p.nelem() );
00682
00683 interp_lin_vector( x, p0, x0, p );
00684 }
00685
00686
00688
00705 void interpp(
00706 MatrixView A,
00707 ConstVectorView p0,
00708 ConstMatrixView A0,
00709 ConstVectorView p )
00710 {
00711 assert( A.nrows() == A0.nrows() );
00712 assert( A.ncols() == p.nelem() );
00713
00714
00715 Vector logp0(p0.nelem());
00716 Vector logp(p.nelem());
00717 transform( logp0, log, p0 );
00718 transform( logp, log, p );
00719
00720 interp_lin_matrix( A, logp0, A0, logp );
00721 }
00722
00723
00724
00726
00739 Numeric interpp(
00740 ConstVectorView p0,
00741 ConstVectorView x0,
00742 const Numeric p )
00743 {
00744
00745 Vector logp0(p0.nelem());
00746 transform( logp0, log, p0 );
00747
00748 return interp_lin( logp0, x0, log(p) );
00749 }
00750
00751
00752
00754
00774 void interpz(
00775 VectorView x,
00776 ConstVectorView p0,
00777 ConstVectorView z0,
00778 ConstVectorView x0,
00779 ConstVectorView z )
00780 {
00781 assert( x.nelem()==z.nelem() );
00782 Vector p(z.nelem());
00783 z2p( p, z0, p0, z );
00784 interpp( x, p0, x0, p );
00785 }
00786
00787
00788
00790
00810 Numeric interpz(
00811 ConstVectorView p0,
00812 ConstVectorView z0,
00813 ConstVectorView x0,
00814 const Numeric z )
00815 {
00816 Vector x(1);
00817 MakeVector Z(z);
00818 interpz( x, p0, z0, x0, Z );
00819 return x[0];
00820 }
00821
00822
00823
00825
00827
00829
00839 Numeric ztan_geom(
00840 const Numeric za,
00841 const Numeric z_plat,
00842 const Numeric r_geoid )
00843 {
00844 Numeric z_tan;
00845 if ( za >= 90 )
00846 z_tan = (r_geoid+z_plat)*sin(DEG2RAD*za) - r_geoid;
00847 else
00848 z_tan = 999e3;
00849 return z_tan;
00850 }
00851
00852
00853
00855
00871 Numeric n_for_z(
00872 const Numeric z,
00873 ConstVectorView p_abs,
00874 ConstVectorView z_abs,
00875 ConstVectorView refr_index,
00876 const Numeric atm_limit )
00877
00878 {
00879 if ( z > atm_limit )
00880 return 1.0;
00881 else
00882 return interpz( p_abs, z_abs, refr_index, z );
00883 }
00884
00885
00887
00908 Numeric refr_constant(
00909 const Numeric r_geoid,
00910 const Numeric za,
00911 const Numeric z_plat,
00912 ConstVectorView p_abs,
00913 ConstVectorView z_abs,
00914 const Numeric atm_limit,
00915 ConstVectorView refr_index )
00916 {
00917 Numeric n_plat = n_for_z( z_plat, p_abs, z_abs, refr_index, atm_limit );
00918
00919 return (r_geoid + z_plat) * sin(DEG2RAD*za) * n_plat;
00920 }
00921
00922
00923
00925
00940 Numeric ztan_refr(
00941 const Numeric c,
00942 const Numeric za,
00943 const Numeric z_plat,
00944 const Numeric z_ground,
00945 ConstVectorView p_abs,
00946 ConstVectorView z_abs,
00947 ConstVectorView refr_index,
00948 const Numeric r_geoid )
00949 {
00950 const Numeric atm_limit = last(z_abs);
00951 if ( za < 90 )
00952 return ztan_geom( za, z_plat, r_geoid );
00953 else
00954 {
00955 const Index n = z_abs.nelem();
00956 Index i;
00957
00958 for ( i=(n-1); (i>=0) && (r_geoid+z_abs[i])*refr_index[i]>c; i-- )
00959 {
00960 if ( z_abs[i] <= z_ground )
00961 {
00962 Numeric n_ground = n_for_z(z_ground,p_abs,z_abs,refr_index,atm_limit);
00963 Numeric theta = RAD2DEG*asin(c/n_ground/(r_geoid+z_ground));
00964 return ztan_geom( 180-theta, z_ground, r_geoid );
00965 }
00966 }
00967 if ( i == (n-1) )
00968 return ztan_geom( za, z_plat, r_geoid );
00969 else
00970 {
00971 Vector zs(2), cs(2);
00972 zs[0] = z_abs[i];
00973 zs[1] = z_abs[i+1];
00974 cs[0] = (r_geoid+z_abs[i])*refr_index[i];
00975 cs[1] = (r_geoid+z_abs[i+1])*refr_index[i+1];
00976 return interp_lin( cs, zs, c );
00977 }
00978 }
00979 }
00980
00982
00998 void e_eq_water( VectorView e_eq,
00999 ConstVectorView t )
01000 {
01001 Index cur_elem;
01002 double a, b, c, d, e, t_cur;
01003
01004
01005 assert( e_eq.nelem() == t.nelem() );
01006
01007
01008 a = -6096.9385;
01009 b = 21.2409642;
01010 c = -2.711193e-2;
01011 d = 1.673952e-5;
01012 e = 2.433502;
01013
01014 for ( cur_elem = 0; cur_elem < t.nelem(); cur_elem++ )
01015 {
01016 t_cur = t[cur_elem];
01017 e_eq[cur_elem] = exp( a / t_cur + b + c * t_cur + d * t_cur * t_cur + e * log( t_cur ) );
01018 }
01019 }
01020
01022
01038 void e_eq_ice( VectorView e_eq,
01039 ConstVectorView t )
01040 {
01041 Index cur_elem;
01042 double a, b, c, d, e, t_cur;
01043
01044
01045 assert( e_eq.nelem() == t.nelem() );
01046
01047
01048 a = -6024.5282;
01049 b = 29.32707;
01050 c = 1.0613868e-2;
01051 d = -1.3198825e-5;
01052 e = -0.49382577;
01053
01054 for ( cur_elem = 0; cur_elem < t.nelem(); cur_elem++ )
01055 {
01056 t_cur = t[cur_elem];
01057 e_eq[cur_elem] = exp( a / t_cur + b + c * t_cur + d * t_cur * t_cur + e * log( t_cur ) );
01058 }
01059 }