00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00021
00023
00044
00045
00047
00048 #include <cmath>
00049 #include "arts.h"
00050 #include "atm_funcs.h"
00051 #include "complex.h"
00052 #include "matpackI.h"
00053 #include "los.h"
00054 #include "math_funcs.h"
00055 #include "messages.h"
00056 #include "auto_wsv.h"
00057 #include "auto_md.h"
00058 extern const Numeric PI;
00059 extern const Numeric DEG2RAD;
00060 extern const Numeric RAD2DEG;
00061 extern const Numeric COSMIC_BG_TEMP;
00062 extern const Numeric SUN_TEMP;
00063 extern const Numeric PLANCK_CONST;
00064 extern const Numeric BOLTZMAN_CONST;
00065 extern const Numeric SPEED_OF_LIGHT;
00066 extern const Numeric EARTH_GRAV_CONST;
00067 extern const Numeric AVOGADROS_NUMB;
00068
00069
00071
00073
00075
00086 bool any_ground( const ArrayOfIndex& ground )
00087 {
00088 for ( Index i=0; i<ground.nelem(); i++ )
00089 {
00090 if ( ground[i] )
00091 return 1;
00092 }
00093 return 0;
00094 }
00095
00096
00097
00098
00100
00126 void los_geometric(
00127 Vector& z,
00128 Vector& psi,
00129 Numeric& l_step,
00130 const Numeric& z_plat,
00131 const Numeric& za,
00132 const Numeric& atm_limit,
00133 const Numeric& r_geoid )
00134 {
00135
00136 assert( za <= 90 );
00137
00138 Vector l;
00139 Index nz;
00140
00141
00142
00143 double a, b;
00144 double llim;
00145
00146
00147 a = r_geoid + atm_limit;
00148 b = (r_geoid+z_plat) * sin(DEG2RAD*za);
00149 llim = sqrt( a*a - b*b ) - (r_geoid+z_plat)*cos(DEG2RAD*za) ;
00150
00151
00152 if ( llim < l_step )
00153 l_step = llim*0.9999;
00154
00155
00156 linspace( l, 0, llim, l_step );
00157
00158 nz = l.nelem();
00159 z.resize( nz );
00160 psi.resize( nz );
00161
00162
00163 b = r_geoid + z_plat;
00164 a = b * b;
00165
00166 for ( Index i=0; i<nz; i++ )
00167 {
00168 z[i] = sqrt( a + l[i]*l[i] + 2.0*b*l[i]*cos(DEG2RAD*za) );
00169 psi[i] = RAD2DEG * acos( (a+z[i]*z[i]-l[i]*l[i]) / (2.0*b*z[i]) );
00170
00171
00172 if ( isnan(psi[i]) )
00173 psi[i] = 0;
00174
00175 z[i] = z[i] - r_geoid;
00176 }
00177 }
00178
00179
00180
00182
00205 void los_refraction(
00206 Vector& z,
00207 Vector& psi,
00208 Numeric& l_step,
00209 const Numeric& z_plat,
00210 const Numeric& za,
00211 const Numeric& atm_limit,
00212 const Numeric& r_geoid,
00213 const Vector& ,
00214 const Vector& z_abs,
00215 const Index& ,
00216 const Index& refr_lfac,
00217 const Vector& refr_index,
00218 const Numeric& c )
00219 {
00220
00221 assert( za <= 90 );
00222
00223
00224
00225 Index np;
00226 {
00227
00228 double a = r_geoid + atm_limit;
00229 double b = (r_geoid+z_plat) * sin(DEG2RAD*za);
00230 double llim = sqrt( a*a - b*b ) - (r_geoid+z_plat)*cos(DEG2RAD*za) ;
00231
00232
00233 if ( llim < l_step )
00234 l_step = llim*0.9999;
00235
00236 np = Index( ceil( 1.5 * ( llim/l_step + 1) ) );
00237 }
00238 Vector zv(np), pv(np);
00239
00240
00241 const double l = l_step / refr_lfac;
00242 Index i = refr_lfac;
00243 double z1;
00244 double z2 = z_plat;
00245 double rz1, rz2;
00246 double psi1;
00247 double psi2 = 0;
00248 double n1, n2;
00249 double n;
00250 double c2=c; c2 = c2 * c2;
00251 Index j;
00252 double d, e, f;
00253
00254 np = 0;
00255
00256
00257
00258
00259 const Index nz = z_abs.nelem();
00260 Index iz;
00261 for ( iz=0; (iz<nz) && (z_abs[iz]<=z2); iz++ ) {}
00262 if ( iz < nz )
00263 n2 = refr_index[iz-1] + (refr_index[iz]-refr_index[iz-1])*
00264 (z2-z_abs[iz-1])/(z_abs[iz]-z_abs[iz-1]);
00265 else
00266 n2 = 1;
00267
00268 while ( z2 <= atm_limit )
00269 {
00270
00271 z1 = z2;
00272 psi1 = psi2;
00273 n1 = n2;
00274
00275 if ( i == Index(refr_lfac) )
00276 {
00277 zv[np] = z2;
00278 pv[np] = RAD2DEG * psi2;
00279 i = 1;
00280 np++;
00281 assert( np < zv.nelem() );
00282 }
00283 else
00284 i++;
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294 for ( j=1; j<=2; j++ )
00295 {
00296
00297 if ( j == 1 )
00298 n = n1;
00299 else
00300 n = ( n1 + n2 ) / 2;
00301
00302 rz1 = z1 + r_geoid;
00303 d = rz1 * rz1;
00304 e = c2/(n*n);
00305 f = d - e;
00306
00307
00308
00309
00310 if ( f <= 0 )
00311 rz2 = sqrt( l*l + e );
00312 else
00313 rz2 = sqrt( pow( l + sqrt(f), 2 ) + e );
00314
00315 z2 = rz2 - r_geoid;
00316
00317
00318 for ( ; (iz<nz) && (z_abs[iz]<=z2); iz++ ) {}
00319 if ( iz < nz )
00320 n2 = refr_index[iz-1] + (refr_index[iz]-refr_index[iz-1])*
00321 (z2-z_abs[iz-1])/(z_abs[iz]-z_abs[iz-1]);
00322 else
00323 n2 = 1;
00324 }
00325
00326 psi2 = psi1 + acos( (d+rz2*rz2-l*l) / (2*rz1*rz2) );
00327 }
00328
00329
00330 z.resize( np );
00331 psi.resize( np );
00332 z = zv[Range(0,np)];
00333 psi = pv[Range(0,np)];
00334 }
00335
00336
00337
00339
00341
00343
00372 void los_1za(
00373 Vector& z,
00374 Vector& psi,
00375 Numeric& l_step,
00376 Index& ground,
00377 Index& start,
00378 Index& stop,
00379 Numeric& z_tan,
00380 const Numeric& z_plat,
00381 const Numeric& za,
00382 const Numeric& l_step_max,
00383 const Numeric& z_ground,
00384 const Numeric& r_geoid,
00385 const Vector& p_abs,
00386 const Vector& z_abs,
00387 const Index& refr,
00388 const Index& refr_lfac,
00389 const Vector& refr_index )
00390 {
00391 Numeric c;
00392
00393
00394 const Numeric atm_limit = last(z_abs);
00395
00396 if ( refr )
00397 {
00398 c = refr_constant( r_geoid, za, z_plat, p_abs, z_abs, atm_limit,
00399 refr_index );
00400 z_tan = ztan_refr( c, za, z_plat, z_ground, p_abs, z_abs, refr_index,
00401 r_geoid );
00402 }
00403 else
00404 z_tan = ztan_geom( za, z_plat, r_geoid );
00405
00406
00407 l_step = l_step_max;
00408
00409
00410
00411 if ( z_plat >= atm_limit )
00412 {
00413 Index nz;
00414 Numeric psi0 = 0;
00415
00416 out3 << " (z_tan = " << z_tan/1e3 << " km)";
00417
00418
00419 if ( z_tan >= atm_limit )
00420 {
00421 ground = 0;
00422 z.resize( 0 );
00423 psi.resize( 0 );
00424 nz = 1;
00425 }
00426
00427
00428 else if ( z_tan >= z_ground )
00429 {
00430 if ( !refr )
00431 {
00432 los_geometric( z, psi, l_step, z_tan, 90.0, atm_limit, r_geoid );
00433 psi0 = za - 90.0;
00434 }
00435 else
00436 {
00437 los_refraction( z, psi, l_step, z_tan, 90.0, atm_limit, r_geoid,
00438 p_abs, z_abs, refr, refr_lfac, refr_index, c );
00439
00440
00441 Numeric zmax = last( z );
00442 Numeric n = interpz( p_abs, z_abs, refr_index, zmax );
00443 Numeric theta = RAD2DEG * asin( c / ((r_geoid+zmax)*n) );
00444
00445 psi0 = theta + za - 180.0 + last(psi);
00446 }
00447
00448 ground = 0;
00449 nz = z.nelem();
00450 }
00451
00452
00453 else
00454 {
00455
00456 Numeric za_g;
00457
00458 if ( !refr )
00459 {
00460 za_g = RAD2DEG * asin( (r_geoid+z_tan) / (r_geoid+z_ground) );
00461
00462 los_geometric( z, psi, l_step, z_ground, za_g, atm_limit, r_geoid );
00463
00464 psi0 = za + za_g - 180.0;
00465 }
00466 else
00467 {
00468 za_g = RAD2DEG * asin( c / ( (r_geoid+z_ground) *
00469 n_for_z( z_ground, p_abs, z_abs, refr_index, atm_limit ) ) );
00470
00471 los_refraction( z, psi, l_step, z_ground, za_g, atm_limit, r_geoid,
00472 p_abs, z_abs, refr, refr_lfac, refr_index, c );
00473
00474
00475 Numeric zmax = last( z );
00476 Numeric n = n_for_z( zmax, p_abs, z_abs, refr_index, atm_limit );
00477 Numeric theta = RAD2DEG * asin( c / ((r_geoid+zmax)*n) );
00478
00479 psi0 = theta + za - 180.0 + last(psi);
00480 }
00481
00482 ground = 1;
00483 nz = z.nelem();
00484 }
00485
00486
00487 if ( psi0 != 0 )
00488 psi += psi0;
00489
00490 start = stop = nz - 1;
00491 }
00492
00493
00494
00495 else if ( za <= 90 )
00496 {
00497 if ( !refr )
00498 los_geometric( z, psi, l_step, z_plat, za, atm_limit, r_geoid );
00499 else
00500 los_refraction( z, psi, l_step, z_plat, za, atm_limit, r_geoid,
00501 p_abs, z_abs, refr, refr_lfac, refr_index, c );
00502 ground = 0;
00503 stop = 0;
00504 start = z.nelem() - 1;
00505 }
00506
00507
00508 else
00509 {
00510
00511
00512 double l1;
00513 double a, b;
00514
00515 out3 << " (z_tan = " << z_tan/1e3 << " km)";
00516
00517
00518 if ( z_tan >= z_ground )
00519 {
00520 if ( !refr )
00521 {
00522
00523 a = r_geoid + z_plat;
00524 b = r_geoid + z_tan;
00525 l1 = sqrt(a*a-b*b);
00526
00527
00528 if ( l1 > l_step_max/10 )
00529 {
00530
00531 stop = Index( ceil( l1 / l_step_max + 1.0 ) - 1 );
00532 l_step = l1 / Numeric(stop);
00533 }
00534
00535
00536 else
00537 {
00538 l_step = l_step_max;
00539 stop = 0;
00540 }
00541
00542 los_geometric( z, psi, l_step, z_tan, 90.0, atm_limit, r_geoid );
00543 }
00544 else
00545 {
00546
00547
00548 Numeric l = l_step / refr_lfac;
00549 los_refraction( z, psi, l, z_tan, 90.0, z_plat+l, r_geoid,
00550 p_abs, z_abs, refr, 1, refr_index, c );
00551
00552
00553
00554 {
00555 Vector ls;
00556 linspace( ls, 0, l*(z.nelem()-1) , l );
00557 l1 = interp_lin( z, ls, z_plat );
00558 }
00559
00560
00561 if ( l1 > l_step_max/10 )
00562 {
00563
00564 stop = Index( ceil( l1 / l_step_max + 1.0 ) - 1 );
00565 l_step = l1 / Numeric(stop);
00566 }
00567
00568
00569 else
00570 {
00571 l_step = l_step_max;
00572 stop = 0;
00573 }
00574
00575 los_refraction( z, psi, l_step, z_tan, 90.0, atm_limit, r_geoid,
00576 p_abs, z_abs, refr, refr_lfac, refr_index, c );
00577 }
00578
00579 ground = 0;
00580 start = z.nelem() - 1;
00581
00582
00583
00584 psi += psi[stop];
00585 }
00586
00587
00588 else
00589 {
00590 Numeric za_g;
00591
00592 if ( !refr )
00593 {
00594
00595 a = r_geoid + z_plat;
00596 b = r_geoid + z_tan;
00597 b = b * b;
00598 l1 = sqrt(a*a-b);
00599 a = r_geoid + z_ground;
00600 l1 = l1 - sqrt(a*a-b);
00601
00602
00603 if ( l1 > l_step_max/10 )
00604 {
00605
00606 stop = Index( ceil( l1 / l_step_max + 1.0 ) - 1 );
00607 l_step = l1 / Numeric(stop);
00608 }
00609
00610
00611 else
00612 {
00613 l_step = l_step_max;
00614 stop = 0;
00615 }
00616
00617 za_g = RAD2DEG * asin( (r_geoid+z_tan) / (r_geoid+z_ground) );
00618
00619 los_geometric( z, psi, l_step, z_ground, za_g, atm_limit, r_geoid );
00620 }
00621
00622 else
00623 {
00624 za_g = RAD2DEG * asin( c / ( (r_geoid+z_ground) *
00625 n_for_z( z_ground, p_abs, z_abs, refr_index, atm_limit ) ) );
00626
00627 Numeric l = l_step / refr_lfac;
00628 los_refraction( z, psi, l, z_ground, za_g, z_plat+l, r_geoid,
00629 p_abs, z_abs, refr, 1, refr_index, c );
00630
00631
00632
00633 {
00634 Vector ls;
00635 linspace( ls, 0, l*(z.nelem()-1) , l );
00636 l1 = interp_lin( z, ls, z_plat );
00637 }
00638
00639
00640 if ( l1 > l_step_max/10 )
00641 {
00642
00643 stop = Index( ceil( l1 / l_step_max + 1.0 ) - 1 );
00644 l_step = l1 / Numeric(stop);
00645 }
00646
00647
00648 else
00649 {
00650 l_step = l_step_max;
00651 stop = 0;
00652 }
00653
00654 los_refraction( z, psi, l_step, z_ground, za_g, atm_limit, r_geoid,
00655 p_abs, z_abs, refr, refr_lfac, refr_index, c );
00656 }
00657
00658 ground = 1;
00659 start = z.nelem() - 1;
00660
00661
00662
00663 psi += psi[stop];
00664 }
00665 }
00666 }
00667
00668
00669
00671
00673
00678 void y_rte (
00679 Vector& y,
00680 const Los& los,
00681 ConstVectorView f_mono,
00682 ConstVectorView y_space,
00683 const ArrayOfMatrix& source,
00684 const ArrayOfMatrix& trans,
00685 ConstVectorView e_ground,
00686 const Numeric& t_ground,
00687 const Index& z_start,
00688 const Index& z_end )
00689 {
00690
00691 const Index n=los.start.nelem();
00692 const Index nf=f_mono.nelem();
00693 Vector y_tmp(nf);
00694 Index iy0=0;
00695
00696 out2 << " Integrating the radiative transfer eq. with emission.\n";
00697
00698 assert (z_start >= 0);
00699 assert (z_end < n);
00700
00701
00702 y.resize( nf* (z_end - z_start + 1) );
00703
00704
00705
00706 Vector y_ground(f_mono.nelem());
00707 if ( any_ground(los.ground) )
00708 {
00709 if ( t_ground <= 0 )
00710 throw runtime_error(
00711 "There are intersections with the ground, but the ground\n"
00712 "temperature is set to be <=0 (are dummy values used?).");
00713 if ( e_ground.nelem() != nf )
00714 throw runtime_error(
00715 "There are intersections with the ground, but the frequency and\n"
00716 "ground emission vectors have different lengths (are dummy values\n"
00717 "used?).");
00718 out2 << " There are intersections with the ground.\n";
00719 planck( y_ground, f_mono, t_ground );
00720 }
00721
00722
00723 out3 << " Zenith angle nr: ";
00724 for ( Index i = z_start; i <= z_end; i++ )
00725 {
00726 if ( (i%20)==0 )
00727 out3 << "\n ";
00728 out3 << " " << i; cout.flush();
00729
00730
00731 rte( y_tmp, los.start[i], los.stop[i], trans[i],
00732 source[i], y_space, los.ground[i], e_ground, y_ground);
00733
00734
00735 y[Range(iy0,nf)] = y_tmp;
00736
00737
00738 iy0 += nf;
00739 }
00740 out3 << "\n";
00741 }
00742
00743
00748 void y_tau (
00749 Vector& y,
00750 const Los& los,
00751 const ArrayOfMatrix& trans,
00752 ConstVectorView e_ground,
00753 const Index& z_start,
00754 const Index& z_end )
00755 {
00756
00757 const Index n=los.start.nelem();
00758 const Index nf=trans[0].nrows();
00759 Index iy, iy0=0;
00760 Vector y_tmp;
00761
00762 assert (z_start >= 0);
00763 assert (z_end < n);
00764
00765 out2 << " Calculating optical thicknesses.\n";
00766
00767
00768 y.resize( nf*n );
00769 y = 1.0;
00770
00771
00772 if ( any_ground(los.ground) )
00773 {
00774 if ( e_ground.nelem() != nf )
00775 throw runtime_error(
00776 "There are intersections with the ground, but the frequency and\n"
00777 "ground emission vectors have different lengths (are dummy values\n"
00778 "used?).");
00779 out2 << " There are intersections with the ground.\n";
00780 }
00781
00782
00783 out3 << " Zenith angle nr: ";
00784 for ( Index i = z_start; i <= z_end; i++ )
00785 {
00786 if ( (i%20)==0 )
00787 out3 << "\n ";
00788 out3 << " " << i; cout.flush();
00789
00790
00791 bl( y_tmp, los.start[i], los.stop[i], trans[i], los.ground[i], e_ground );
00792
00793
00794 for ( iy=0; iy<nf; iy++ )
00795 y[iy0+iy] = -log( y_tmp[iy] );
00796
00797
00798 iy0 += nf;
00799 }
00800 out3 << "\n";
00801 }
00802
00803
00804
00805
00807
00809
00816 void r_geoidStd( Numeric& r_geoid )
00817 {
00818 extern const Numeric EARTH_RADIUS;
00819 r_geoid = EARTH_RADIUS;
00820 }
00821
00822
00829 void r_geoidWGS84(
00830 Numeric& r_geoid,
00831 const Numeric& latitude,
00832 const Numeric& obsdirection )
00833 {
00834 check_if_in_range( -90, 90, latitude, "latitude" );
00835 check_if_in_range( -360, 360, obsdirection, "obsdirection" );
00836
00837 const Numeric rq = 6378.138e3, rp = 6356.752e3;
00838 Numeric a, b, rns, rew;
00839
00840
00841 a = cos(latitude*DEG2RAD);
00842 b = sin(latitude*DEG2RAD);
00843 rns = rq*rq*rp*rp/pow(rq*rq*a*a+rp*rp*b*b,(Numeric)1.5);
00844 rew = rq*rq/sqrt(rq*rq*a*a+rp*rp*b*b);
00845
00846
00847 a = cos(obsdirection*DEG2RAD);
00848 b = sin(obsdirection*DEG2RAD);
00849 r_geoid = 1/(a*a/rns+b*b/rew);
00850 }
00851
00852
00853
00860 void groundOff(
00861 Numeric& z_ground,
00862 Numeric& t_ground,
00863 Vector& e_ground,
00864 const Vector& z_abs )
00865 {
00866 z_ground = z_abs[0];
00867 t_ground = 0;
00868 e_ground.resize( 0 );
00869 }
00870
00871
00872
00879 void groundSet(
00880 Numeric& z_ground,
00881 Numeric& t_ground,
00882 Vector& e_ground,
00883 const Vector& p_abs,
00884 const Vector& t_abs,
00885 const Vector& z_abs,
00886 const Vector& f_mono,
00887 const Numeric& z,
00888 const Numeric& e )
00889 {
00890 check_if_in_range( 0, 1, e, "e" );
00891 z_ground = z;
00892 t_ground = interpz( p_abs, z_abs, t_abs, z );
00893 e_ground.resize( f_mono.nelem() );
00894 e_ground = e;
00895 }
00896
00897
00898
00905 void groundAtBottom(
00906 Numeric& z_ground,
00907 Numeric& t_ground,
00908 Vector& e_ground,
00909 const Vector& t_abs,
00910 const Vector& z_abs,
00911 const Vector& f_mono,
00912 const Numeric& e )
00913 {
00914 check_if_in_range( 0, 1, e, "e" );
00915 z_ground = z_abs[0];
00916 t_ground = t_abs[0];
00917 e_ground.resize( f_mono.nelem() );
00918 e_ground = e;
00919 }
00920
00921
00922
00929 void groundFlatSea(
00930 Numeric& z_ground,
00931 Numeric& t_ground,
00932 Vector& e_ground,
00933 const Vector& p_abs,
00934 const Vector& t_abs,
00935 const Vector& z_abs,
00936 const Vector& f_mono,
00937 const Vector& za_pencil,
00938 const Numeric& z_plat,
00939 const Numeric& r_geoid,
00940 const Index& refr,
00941 const Vector& refr_index,
00942 const String& pol,
00943 const Numeric& t_skin )
00944 {
00945 if( z_abs[0] > 0 || z_abs[z_abs.nelem()-1] < 0 )
00946 throw runtime_error( "The WSV *z_abs* must span 0 m." );
00947 if( min(f_mono) < 10e9 || max(f_mono) > 1000e9 )
00948 throw runtime_error(
00949 "Frequencies below 10 GHz or above 1000 GHz are not allowed." );
00950
00951 z_ground = z_abs[0];
00952
00953 if( t_skin > 0 )
00954 { t_ground = t_skin; }
00955 else
00956 { t_ground = interpz( p_abs, z_abs, t_abs, 0 ); }
00957
00958
00959 Numeric n_plat = 1, n_ground = 1;
00960 if( refr )
00961 {
00962 n_plat = n_for_z( z_plat, p_abs, z_abs, refr_index, last(z_abs) );
00963 n_ground = n_for_z( z_ground, p_abs, z_abs, refr_index, last(z_abs) );
00964 }
00965 Numeric sintheta = n_plat * (r_geoid+z_plat) * sin(DEG2RAD*max(za_pencil))
00966 / ( n_ground * (r_geoid+z_ground) );
00967
00968
00969 if( sintheta > 1 )
00970 { sintheta = 1; }
00971
00972 const Numeric costheta = sqrt( 1 - sintheta*sintheta );
00973
00974
00975
00976
00977
00978
00979
00980
00981 const Numeric theta = 1 - 300 / t_ground;
00982 const Numeric e0 = 77.66 - 103.3 * theta;
00983 const Numeric e1 = 0.0671 * e0;
00984 const Numeric f1 = 20.2 + 146.4 * theta + 316 * theta * theta;
00985 const Numeric e2 = 3.52;
00986 const Numeric f2 = 39.8 * f1;
00987
00988 const Numeric n1 = 1;
00989
00990
00991
00992
00993 e_ground.resize( f_mono.nelem() );
00994
00995 for( Index i=0; i<f_mono.nelem(); i++ )
00996 {
00997 const Complex ifGHz( 0.0, f_mono[i]/1e9 );
00998
00999 const Complex eps = e2 + (e1-e2) / (Numeric(1.0)-ifGHz/f2) +
01000 (e0-e1) / (Numeric(1.0)-ifGHz/f1);
01001 const Complex n2 = sqrt( eps );
01002 Complex a,b;
01003
01004 if( pol == "v" )
01005 {
01006 a = n2 * costheta;
01007 b = n1 * cos( asin( n1 * sintheta / n2.real() ) );
01008 }
01009 else if( pol == "h" )
01010 {
01011 a = n1 * costheta;
01012 b = n2 * cos( asin( n1 * sintheta / n2.real() ) );
01013 }
01014 else
01015 throw runtime_error(
01016 "The keyword argument *pol* must be \"v\" or \"h\"." );
01017
01018
01019 const Numeric r = pow( abs( ( a - b ) / ( a + b ) ), Numeric(2.0) );
01020
01021 e_ground[i] = 1 - r;
01022 }
01023 }
01024
01025
01026
01033 void emissionOn( Index& emission )
01034 {
01035 emission = 1;
01036 }
01037
01038
01045 void emissionOff( Index& emission )
01046 {
01047 emission = 0;
01048 }
01049
01050
01051
01058 void losCalc( Los& los,
01059 Vector& z_tan,
01060 const Numeric& z_plat,
01061 const Vector& za,
01062 const Numeric& l_step,
01063 const Vector& p_abs,
01064 const Vector& z_abs,
01065 const Index& refr,
01066 const Index& refr_lfac,
01067 const Vector& refr_index,
01068 const Numeric& z_ground,
01069 const Numeric& r_geoid )
01070 {
01071 Index n = za.nelem();
01072
01073
01074 check_if_bool( refr, "refr" );
01075 if ( z_ground < z_abs[0] )
01076 throw runtime_error(
01077 "There is a gap between the ground and the lowest absorption altitude.");
01078 if ( z_plat < z_ground )
01079 throw runtime_error("Your platform altitude is below the ground.");
01080 if ( z_plat < z_abs[0] )
01081 throw runtime_error(
01082 "The platform cannot be below the lowest absorption altitude.");
01083 if ( refr && ( p_abs.nelem() != refr_index.nelem() ) )
01084 throw runtime_error(
01085 "Refraction is turned on, but the length of refr_index does not match\n"
01086 "the length of p_abs. Are dummy vales used?");
01087 if ( refr && ( refr_lfac < 1 ) )
01088 throw runtime_error(
01089 "Refraction is turned on, but the refraction length factor is < 1. \n"
01090 "Are dummy vales used?");
01091
01092
01093 los.p.resize( n );
01094 los.psi.resize( n );
01095 los.z.resize( n );
01096 los.l_step.resize( n );
01097 los.ground.resize( n );
01098 los.start.resize( n );
01099 los.stop.resize( n );
01100 z_tan.resize( n );
01101
01102
01103 if ( refr == 0 )
01104 out2 << " Calculating line of sights WITHOUT refraction.\n";
01105 else if ( refr == 1 )
01106 out2 << " Calculating line of sights WITH refraction.\n";
01107
01108 out3 << " z_plat: " << z_plat/1e3 << " km\n";
01109
01110
01111 for ( Index i=0; i<n; i++ )
01112 {
01113 out3 << " za: " << za[i] << " degs.";
01114
01115 los_1za( los.z[i], los.psi[i], los.l_step[i], los.ground[i], los.start[i],
01116 los.stop[i], z_tan[i], z_plat, za[i], l_step, z_ground, r_geoid,
01117 p_abs, z_abs, refr, refr_lfac, refr_index );
01118 out3 << "\n";
01119
01120
01121 los.p[i].resize( los.z[i].nelem() );
01122 z2p( los.p[i], z_abs, p_abs, los.z[i] );
01123 }
01124 }
01125
01126
01127
01134 void zaFromZtan(
01135
01136 Vector& za,
01137 const String& ,
01138
01139 const Vector& z_tan,
01140 const Numeric& z_plat,
01141 const Vector& p_abs,
01142 const Vector& z_abs,
01143 const Index& refr,
01144 const Vector& refr_index,
01145 const Numeric& r_geoid,
01146 const Numeric& z_ground )
01147 {
01148 check_lengths( p_abs, "p_abs", z_abs, "z_abs" );
01149 if ( refr )
01150 check_lengths( p_abs, "p_abs", refr_index, "refr_index" );
01151
01152 const Numeric atm_limit = last(z_abs);
01153 const Index nz = z_tan.nelem();
01154
01155 za.resize(nz);
01156
01157 for (Index i=0; i<nz; i++)
01158 {
01159 if (z_tan[i]>z_plat)
01160 throw runtime_error(
01161 "One tangent altitude is larger than the platform altitude");
01162
01163
01164 if (!refr)
01165 za[i] = 90 + RAD2DEG*acos ( (r_geoid + z_tan[i]) / (r_geoid + z_plat) );
01166
01167
01168 else
01169 {
01170 Numeric nz_plat = n_for_z(z_plat,p_abs,z_abs,refr_index,atm_limit);
01171 if (z_tan[i]>=0)
01172 {
01173
01174 Numeric nza = n_for_z(z_tan[i],p_abs,z_abs,refr_index,atm_limit);
01175 Numeric c = (r_geoid + z_tan[i]) * nza;
01176 za[i] = 180 - RAD2DEG * asin( c / nz_plat / (r_geoid + z_plat));
01177 }
01178 else
01179 {
01180
01181 Numeric ze = RAD2DEG * acos((r_geoid + z_tan[i]) / r_geoid);
01182
01183 Numeric nze = n_for_z(z_ground,p_abs,z_abs,refr_index,atm_limit);
01184 Numeric c = r_geoid * sin(DEG2RAD * (90-ze)) * nze;
01185 za[i] = 180 - RAD2DEG * asin( c / nz_plat / (r_geoid + z_plat));
01186 }
01187 }
01188 }
01189 }
01190
01191
01192
01199 void zaFromDeltat(
01200
01201 Vector& za,
01202
01203 const String& ,
01204
01205 const Numeric& z_plat,
01206 const Vector& p_abs,
01207 const Vector& z_abs,
01208 const Numeric& l_step,
01209 const Index& refr,
01210 const Index& refr_lfac,
01211 const Vector& refr_index,
01212 const Numeric& r_geoid,
01213 const Numeric& z_ground,
01214
01215 const Numeric& delta_t,
01216 const Vector& z_tan_lim )
01217
01218 {
01219
01220 check_lengths( p_abs, "p_abs", z_abs, "z_abs" );
01221 if ( refr ) {
01222 check_lengths( p_abs, "p_abs", refr_index, "refr_index" );
01223 }
01224 if ( z_tan_lim[0] > z_tan_lim[1] )
01225 throw runtime_error(
01226 "The lower tangent latitude is larger than the higher (in z_tan_lim).");
01227
01228
01229 if (!refr)
01230 {
01231
01232 Vector phi(2);
01233 Vector za_lim(2);
01234 String zastr = "za_lim";
01235
01236 zaFromZtan(za_lim, zastr, z_tan_lim, z_plat, p_abs, z_abs, refr,
01237 refr_index, r_geoid, z_ground);
01238
01239 phi[0] = za_lim[0] - 90;
01240 phi[1] = za_lim[1] - 90;
01241
01242 const Numeric ang_step = RAD2DEG * delta_t *
01243 sqrt (EARTH_GRAV_CONST / pow(r_geoid + z_plat,3));
01244
01245 if (((phi[0]-phi[1])/ang_step) <=0)
01246 throw runtime_error("The time given is too large to have any cross-link in the given altitude range");
01247
01248 const Index n=Index(floor((phi[0]-phi[1])/ang_step));
01249
01250 za.resize(n);
01251 for ( Index j=0;j<n;j++ )
01252 za[j] = 90 + phi[0] - (j * ang_step);
01253 }
01254
01255
01256 else
01257 {
01258 const Index ztanstep = 100;
01259 const Index n=Index(floor((z_tan_lim[1]-z_tan_lim[0])/ztanstep));
01260
01261 Vector z_tan_1(n);
01262 Vector z_tan_2(n);
01263 za.resize(n);
01264
01265
01266 for ( Index j=0;j<n;j++ )
01267 z_tan_1[j]=z_tan_lim[0]+j*ztanstep;
01268
01269
01270 String za_str = "za";
01271 zaFromZtan(za, za_str, z_tan_1, z_plat, p_abs, z_abs, refr, refr_index,
01272 r_geoid, z_ground);
01273
01274
01275 Los los;
01276 losCalc(los,z_tan_2,z_plat,za,l_step,p_abs,z_abs,refr,refr_lfac,
01277 refr_index,z_ground,r_geoid);
01278
01279 Vector psizb(n);
01280 for ( Index j=0;j<n;j++ )
01281 psizb[j]=los.psi[j][0];
01282
01283 const Numeric psitop = psizb[n-1];
01284 const Numeric psibot = psizb[0];
01285
01286
01287 const Numeric ang_step = RAD2DEG * delta_t *
01288 sqrt (EARTH_GRAV_CONST / pow(r_geoid + z_plat,3));
01289
01290
01291 if (((psibot-psitop)/ang_step)<=0)
01292 throw runtime_error("The time given is too large to have any cross-link in the given altitude range");
01293 const Index np=Index(floor((psibot-psitop)/ang_step));
01294
01295 Vector psit(np);
01296 for ( Index j=0;j<np;j++ )
01297 psit[j]=psibot-j*ang_step;
01298
01299
01300
01301 z_tan_1.resize(np);
01302 for ( Index j=0;j<np;j++ )
01303 {
01304 z_tan_1[j]=interp_lin(psizb,z_tan_2,psit[j]);
01305 }
01306
01307
01308 zaFromZtan(za, za_str, z_tan_1, z_plat, p_abs, z_abs, refr, refr_index,
01309 r_geoid, z_ground);
01310 }
01311 }
01312
01313
01314
01321 void sourceCalc(
01322 ArrayOfMatrix& source,
01323 const Index& emission,
01324 const Los& los,
01325 const Vector& p_abs,
01326 const Vector& t_abs,
01327 const Vector& f_mono )
01328 {
01329 check_if_bool( emission, "emission" );
01330 check_lengths( p_abs, "p_abs", t_abs, "t_abs" );
01331
01332 if ( emission == 0 )
01333 {
01334 out2 << " Setting the source array to be empty.\n";
01335 source.resize( 0 );
01336 }
01337
01338 else
01339 {
01340 Vector tlos;
01341 const Index nza=los.start.nelem();
01342 const Index nf=f_mono.nelem();
01343 Index nlos;
01344 Matrix b;
01345 Index iv, ilos;
01346
01347 out2 << " Calculating the source function for LTE and no scattering.\n";
01348
01349
01350 source.resize(nza);
01351
01352
01353
01354
01355
01356 out3 << " Zenith angle nr: ";
01357 for (Index i=0; i<nza; i++ )
01358 {
01359 if ( (i%20)==0 )
01360 out3 << "\n ";
01361 out3 << " " << i; cout.flush();
01362
01363 if ( los.p[i].nelem() > 0 )
01364 {
01365 nlos = los.p[i].nelem();
01366 tlos.resize( nlos );
01367 interpp( tlos, p_abs, t_abs, los.p[i] );
01368 b.resize( nf, nlos );
01369 planck( b, f_mono, tlos );
01370 source[i].resize(nf,nlos-1);
01371 for ( ilos=0; ilos<(nlos-1); ilos++ )
01372 {
01373 for ( iv=0; iv<nf; iv++ )
01374 source[i](iv,ilos) = ( b(iv,ilos) + b(iv,ilos+1) ) / 2.0;
01375 }
01376 }
01377 }
01378 out3 << "\n";
01379 }
01380 }
01381
01382
01383
01390 void transCalc(
01391 ArrayOfMatrix& trans,
01392 const Los& los,
01393 const Vector& p_abs,
01394 const Matrix& abs )
01395 {
01396 check_length_ncol( p_abs, "p_abs", abs, "abs" );
01397
01398
01399 const Index n = los.start.nelem();
01400 const Index nf = abs.nrows();
01401 Index np;
01402 Index row, col;
01403 Matrix abs2 ;
01404 Numeric w;
01405
01406 out2 << " Calculating transmissions WITHOUT scattering.\n";
01407
01408
01409 trans.resize(n);
01410
01411
01412
01413
01414 out3 << " Zenith angle nr: ";
01415 for (Index i=0; i<n; i++ )
01416 {
01417 if ( (i%20)==0 )
01418 out3 << "\n ";
01419 out3 << " " << i; cout.flush();
01420
01421 np = los.p[i].nelem();
01422 if ( np > 0 )
01423 {
01424 abs2.resize( nf, np );
01425 interpp( abs2, p_abs, abs, los.p[i] );
01426 trans[i].resize( nf, np-1 );
01427 w = -0.5*los.l_step[i];
01428 for ( row=0; row<nf; row++ )
01429 {
01430 for ( col=0; col<(np-1); col++ )
01431 trans[i](row,col) = exp( w * ( abs2(row,col)+abs2(row,col+1) ) );
01432 }
01433 }
01434 }
01435 out3 << "\n";
01436 }
01437
01438
01439
01446 void y_spaceStd(
01447 Vector& y_space,
01448 const Vector& f,
01449 const String& choice )
01450 {
01451 y_space.resize( f.nelem() );
01452
01453 if ( choice == "zero" )
01454 {
01455 y_space = 0.0;
01456 out2 << " Setting y_space to zero.\n";
01457 }
01458 else if ( choice == "cbgr" )
01459 {
01460 planck( y_space, f, COSMIC_BG_TEMP );
01461 out2 << " Setting y_space to cosmic background radiation.\n";
01462 }
01463 else if ( choice == "sun" )
01464 {
01465 planck( y_space, f, SUN_TEMP );
01466 out2 << " Setting y_space to blackbody radiation corresponding to "
01467 << "the Sun temperature\n";
01468 }
01469 else
01470 throw runtime_error(
01471 "Possible choices for y_space are \"zero\", \"cbgr\" and \"sun\".");
01472
01473 }
01474
01475
01476
01483 void yCalc (
01484 Vector& y,
01485 const Index& emission,
01486 const Los& los,
01487 const Vector& f_mono,
01488 const Vector& y_space,
01489 const ArrayOfMatrix& source,
01490 const ArrayOfMatrix& trans,
01491 const Vector& e_ground,
01492 const Numeric& t_ground )
01493 {
01494
01495
01496
01497 check_if_bool( emission, "emission" );
01498 if ( los.p.nelem() != trans.nelem() )
01499 throw runtime_error(
01500 "The number of zenith angles of *los* and *trans* are different.");
01501 if ( emission )
01502 {
01503 check_lengths( f_mono, "f_mono", y_space, "y_space" );
01504 if ( los.p.nelem() != source.nelem() )
01505 throw runtime_error(
01506 "The number of zenith angles of *los* and *source* are different.");
01507 }
01508
01509 if ( emission == 0 )
01510 y_tau( y, los, trans, e_ground, 0, los.start.nelem () - 1 );
01511 else
01512 y_rte( y, los, f_mono, y_space, source, trans, e_ground, t_ground,
01513 0, los.start.nelem () - 1 );
01514 }
01515
01522 void sourcetransyCalcSaveMemory(
01523 Vector& y,
01524 const Index& emission,
01525 const Los& los,
01526 const Vector& p_abs,
01527 const Vector& t_abs,
01528 const Vector& f_mono,
01529 const Matrix& abs,
01530 const Vector& y_space,
01531 const Vector& e_ground,
01532 const Numeric& t_ground,
01533 const Index& f_chunksize)
01534 {
01535
01536 const Index n=los.start.nelem();
01537 const Index nf=f_mono.nelem();
01538
01539 assert (f_chunksize > 0);
01540
01541 Index chunksize;
01542 if (f_chunksize > nf) chunksize=nf;
01543 else chunksize=f_chunksize;
01544
01545
01546 y.resize( nf*n );
01547
01548 for (Index i = 0; i < nf / chunksize + 1; i++)
01549 {
01550 Index nf_local;
01551
01552 if ((i+1) * chunksize <= nf)
01553 nf_local = chunksize;
01554 else
01555 nf_local = nf % (i * chunksize);
01556
01557 if (! nf_local) break;
01558
01559 Range r (i * chunksize, nf_local);
01560
01561 ConstVectorView f_monolocal (f_mono [r]);
01562 ConstVectorView e_groundlocal (e_ground [r]);
01563 ConstVectorView y_spacelocal (y_space [r]);
01564 ConstMatrixView abslocal (abs (r, joker));
01565
01566
01567 ArrayOfMatrix sourcelocal;
01568 sourceCalc(
01569 sourcelocal,
01570
01571 emission,
01572 los,
01573 p_abs,
01574 t_abs,
01575 f_monolocal);
01576
01577
01578 ArrayOfMatrix translocal;
01579 transCalc(
01580 translocal,
01581
01582 los,
01583 p_abs,
01584 abslocal);
01585
01586
01587
01588
01589 check_if_bool( emission, "emission" );
01590 if ( los.p.nelem() != translocal.nelem() )
01591 throw runtime_error( "The number of zenith angles of *los* and "
01592 "*translocal* are different.");
01593 if ( emission )
01594 {
01595 if ( los.p.nelem() != sourcelocal.nelem() )
01596 throw runtime_error( "The number of zenith angles of *los* "
01597 "and *sourcelocal* are different.");
01598 }
01599
01600 for (Index j = 0; j < n; j++)
01601 {
01602
01603 Vector ylocal;
01604
01605 if ( emission == 0 )
01606 y_tau( ylocal, los, translocal, e_groundlocal, j, j );
01607 else
01608 y_rte( ylocal, los, f_monolocal, y_spacelocal, sourcelocal,
01609 translocal, e_groundlocal, t_ground, j, j );
01610
01611
01612 y[Range(i * chunksize + j * nf, nf_local)] = ylocal;
01613 }
01614 }
01615
01616 }
01617
01624 void yTB (
01625 Vector& y,
01626 const Vector& f_mono,
01627 const Vector& za_pencil )
01628 {
01629 out2 << " Converts the values of *y* to Planck temperatures.\n";
01630 invplanck( y, f_mono, za_pencil );
01631 }
01632
01633
01634
01641 void yTRJ (
01642 Vector& y,
01643 const Vector& f_mono,
01644 const Vector& za_pencil )
01645 {
01646 out2 << " Converts the values of *y* to Rayleigh-Jean temperatures.\n";
01647 invrayjean( y, f_mono, za_pencil );
01648 }
01649
01650
01651
01658 void MatrixTRJ (
01659 Matrix& kout,
01660
01661 const String& kout_name,
01662
01663 const Vector& f_mono,
01664 const Vector& za_pencil,
01665
01666 const Matrix& kin,
01667
01668 const String& kin_name)
01669 {
01670 out2 << " Converts the values of *" << kin_name << "* to Rayleigh-Jean "
01671 << "temperatures,\n and stores the result in *" << kout_name
01672 << "*.\n";
01673
01674
01675 if ( kout.nrows()!=kin.nrows() || kout.ncols()!=kin.ncols() )
01676 kout.resize(kin.nrows(),kin.ncols());
01677
01678 for ( Index i=0; i<kin.ncols(); i++ )
01679 {
01680 kout(Range(joker),i) = kin(Range(joker),i);
01681 invrayjean( kout(Range(joker),i), f_mono, za_pencil );
01682 }
01683 }
01684
01685
01686
01693 void MatrixTB (
01694 Matrix& kout,
01695
01696 const String& kout_name,
01697
01698 const Vector& f_mono,
01699 const Vector& za_pencil,
01700
01701 const Matrix& kin,
01702
01703 const String& kin_name)
01704 {
01705 out2 << " Converts the values of *" << kin_name << "* to Planck "
01706 << "temperatures,\n and stores the result in *" << kout_name
01707 << "*.\n";
01708
01709
01710 if (kout.nrows()!=kin.nrows() ||
01711 kout.ncols()!=kin.ncols())
01712 kout.resize(kin.nrows(),kin.ncols());
01713
01714 Vector y(kin.nrows());
01715 for ( Index i=0; i<kin.ncols(); i++ )
01716 {
01717 y = kin(Range(joker),i);
01718 invplanck( y, f_mono, za_pencil );
01719 kout(Range(joker),i) = y;
01720 }
01721 }
01722
01723
01724
01725
01732 void CoolingRates(
01733 Matrix& coolrate,
01734 const Numeric& lstep0,
01735 const Vector& p_abs,
01736 const Vector& z_abs,
01737 const Vector& t_abs,
01738 const Vector& f_mono,
01739 const Matrix& absorption,
01740 const Vector& za_pencil,
01741 const Index& refr,
01742 const Index& refr_lfac,
01743 const Vector& refr_index,
01744 const Numeric& r_geoid,
01745 const Numeric& z_ground,
01746 const Vector& e_ground,
01747 const Numeric& t_ground,
01748 const Vector& p_coolrate,
01749 const Numeric& lstep_limit )
01750 {
01751
01752 const Index nf = f_mono.nelem();
01753 const Index nza = za_pencil.nelem();
01754 const Index np = p_coolrate.nelem();
01755
01756
01757 if ( za_pencil[0] != 0 || za_pencil[nza-1] != 180 )
01758 throw runtime_error(
01759 "The given zenith angle grid must start with 0 and end with 180" );
01760
01761
01762 coolrate.resize( nf, np );
01763 coolrate = 0;
01764
01765
01766 Vector z_coolrate( np ), t_coolrate( np );
01767 Matrix abs_coolrate( nf, np );
01768 interpp( z_coolrate, p_abs, z_abs, p_coolrate );
01769 interpp( t_coolrate, p_abs, t_abs, p_coolrate );
01770 interpp( abs_coolrate, p_abs, absorption, p_coolrate );
01771
01772
01773
01774 Matrix intgrmatrix( nf, nza, 0 );
01775
01776
01777 const Index emission = 1;
01778 Vector y_space( nf );
01779 y_spaceStd( y_space, f_mono, "cbgr" );
01780
01781
01782 Los los;
01783 Vector z_tan, y;
01784 ArrayOfMatrix source, trans;
01785
01786
01787
01788 for( Index ip=0; ip<np; ip++ )
01789 {
01790
01791
01792
01793
01794
01795 for( Index iza=1; iza<nza-1; iza++ )
01796 {
01797 Numeric lstep = lstep0 / abs( cos( DEG2RAD*za_pencil[iza] ) );
01798 if( lstep > lstep_limit || abs( za_pencil[iza] - 90 ) < 0.01 )
01799 lstep = lstep_limit;
01800
01801 losCalc( los, z_tan, z_coolrate[ip], Vector(1,za_pencil[iza]),
01802 lstep, p_abs, z_abs, refr, refr_lfac, refr_index,
01803 z_ground, r_geoid );
01804 sourceCalc( source, emission, los, p_abs, t_abs, f_mono );
01805 transCalc( trans, los, p_abs, absorption );
01806 yCalc ( y, emission, los, f_mono, y_space, source, trans,
01807 e_ground, t_ground );
01808
01809 const Numeric sinv = sin( DEG2RAD*za_pencil[iza] );
01810
01811 for( Index iv=0; iv<nf; iv++ )
01812 {
01813 Numeric bb_eff;
01814 if( los.stop[0] == 0)
01815 { bb_eff = source[0](iv,los.stop[0]); }
01816 else
01817 { bb_eff = source[0](iv,los.stop[0]-1); }
01818
01819 intgrmatrix(iv,iza) = (bb_eff-y[iv])*abs_coolrate(iv,ip)*sinv;
01820 }
01821 }
01822
01823
01824
01825 const Numeric cp = 1.00e3;
01826 const Numeric rho = 28.8 * p_coolrate[ip] /
01827 ( t_coolrate[ip] * BOLTZMAN_CONST * AVOGADROS_NUMB );
01828 const Numeric factor1 = 2 * PI / (cp*rho);
01829
01830 for( Index iza=1; iza<nza; iza++ )
01831 {
01832 const Numeric factor2 = DEG2RAD*(za_pencil[iza]-za_pencil[iza-1]);
01833
01834 for( Index iv=0; iv<nf; iv++ )
01835 {
01836 coolrate(iv,ip) += factor1 * factor2 *
01837 ( intgrmatrix(iv,iza) + intgrmatrix(iv,iza-1) ) / 2;
01838 }
01839 }
01840 }
01841 }