ARTS 2.5.11 (git: 725533f0)
refraction.cc
Go to the documentation of this file.
1/*===========================================================================
2 === File description
3 ===========================================================================*/
4
13/*===========================================================================
14 === External declarations
15 ===========================================================================*/
16
17#include "refraction.h"
18#include <cmath>
19#include "arts_constants.h"
20#include "auto_md.h"
21#include "matpack_complex.h"
22#include "geodetic.h"
23#include "interpolation.h"
24#include "special_interp.h"
25#include "check_input.h"
26#include "arts_conversions.h"
27
28inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
29inline constexpr Numeric RAD2DEG=Conversion::rad2deg(1);
30inline constexpr Numeric TEMP_0_C=Constant::temperature_at_0c;
31
32/*===========================================================================
33 === The functions (in alphabetical order)
34 ===========================================================================*/
35
37
57void complex_n_water_liebe93(Matrix& complex_n,
58 const Vector& f_grid,
59 const Numeric& t) {
60 chk_if_in_range("t", t, TEMP_0_C - 40, TEMP_0_C + 100);
61 chk_if_in_range("min of f_grid", min(f_grid), 10e9, 1000e9);
62 chk_if_in_range("max of f_grid", max(f_grid), 10e9, 1000e9);
63
64 const Index nf = f_grid.nelem();
65
66 complex_n.resize(nf, 2);
67
68 // Implementation following epswater93.m (by C. Mätzler), part of Atmlab,
69 // but numeric values strictly following the paper version (146, not 146.4)
70 const Numeric theta = 1 - 300 / t;
71 const Numeric e0 = 77.66 - 103.3 * theta;
72 const Numeric e1 = 0.0671 * e0;
73 const Numeric f1 = 20.2 + 146 * theta + 316 * theta * theta;
74 const Numeric e2 = 3.52;
75 const Numeric f2 = 39.8 * f1;
76
77 for (Index iv = 0; iv < nf; iv++) {
78 const Complex ifGHz(0.0, f_grid[iv] / 1e9);
79
80 Complex n = sqrt(e2 + (e1 - e2) / (Numeric(1.0) - ifGHz / f2) +
81 (e0 - e1) / (Numeric(1.0) - ifGHz / f1));
82
83 complex_n(iv, 0) = n.real();
84 complex_n(iv, 1) = n.imag();
85 }
86}
87
89
107void complex_n_ice_matzler06(Matrix& complex_n,
108 const Vector& f_grid,
109 const Numeric& t) {
110 chk_if_in_range("t", t, 20., 280.);
111 chk_if_in_range("min of f_grid", min(f_grid), 10e6, 3000e9);
112 chk_if_in_range("max of f_grid", max(f_grid), 10e6, 3000e9);
113
114 const Index nf = f_grid.nelem();
115
116 complex_n.resize(nf, 2);
117
118 // some parametrization constants
119 const Numeric B1 = 0.0207;
120 const Numeric B2 = 1.16e-11;
121 const Numeric b = 335.;
122
123 const Numeric deltabeta = exp(-9.963 + 0.0372 * (t - 273));
124 const Numeric ebdt = exp(b / t);
125 const Numeric betam = (B1 / t) * ebdt / ((ebdt - 1.) * (ebdt - 1.));
126
127 const Numeric theta = 300. / t - 1;
128 const Numeric alfa = (0.00504 + 0.0062 * theta) * exp(-22.1 * theta);
129 const Numeric reps = 3.1884 + 9.1e-4 * (t - 273);
130
131 for (Index iv = 0; iv < nf; iv++) {
132 Numeric f = f_grid[iv] / 1e9;
133 Numeric beta = betam + B2 * f * f + deltabeta;
134 Numeric ieps = alfa / f + beta * f;
135
136 Complex eps(reps, ieps);
137 Complex n = sqrt(eps);
138 complex_n(iv, 0) = n.real();
139 complex_n(iv, 1) = n.imag();
140 }
141}
142
144
171 Numeric& refr_index_air,
172 Numeric& refr_index_air_group,
173 const Agenda& refr_index_air_agenda,
174 ConstVectorView p_grid,
175 ConstVectorView refellipsoid,
176 ConstTensor3View z_field,
177 ConstTensor3View t_field,
178 ConstTensor4View vmr_field,
179 ConstVectorView f_grid,
180 const Numeric& r) {
181 Numeric rtp_pressure, rtp_temperature;
182 Vector rtp_vmr;
183
184 // Pressure grid position
185 ArrayOfGridPos gp(1);
186 gridpos(gp, z_field(joker, 0, 0), Vector(1, r - refellipsoid[0]));
187
188 // Altitude interpolation weights
189 Matrix itw(1, 2);
190 interpweights(itw, gp);
191
192 // Pressure
193 Vector dummy(1);
194 itw2p(dummy, p_grid, gp, itw);
195 rtp_pressure = dummy[0];
196
197 // Temperature
198 interp(dummy, itw, t_field(joker, 0, 0), gp);
199 rtp_temperature = dummy[0];
200
201 // VMR
202 const Index ns = vmr_field.nbooks();
203 //
204 rtp_vmr.resize(ns);
205 //
206 for (Index is = 0; is < ns; is++) {
207 interp(dummy, itw, vmr_field(is, joker, 0, 0), gp);
208 rtp_vmr[is] = dummy[0];
209 }
210
212 refr_index_air,
213 refr_index_air_group,
214 rtp_pressure,
215 rtp_temperature,
216 rtp_vmr,
217 Vector{f_grid},
218 refr_index_air_agenda);
219}
220
222
251 Numeric& refr_index_air,
252 Numeric& refr_index_air_group,
253 const Agenda& refr_index_air_agenda,
254 ConstVectorView p_grid,
255 ConstVectorView lat_grid,
256 ConstVectorView refellipsoid,
257 ConstTensor3View z_field,
258 ConstTensor3View t_field,
259 ConstTensor4View vmr_field,
260 ConstVectorView f_grid,
261 const Numeric& r,
262 const Numeric& lat) {
263 Numeric rtp_pressure, rtp_temperature;
264 Vector rtp_vmr;
265
266 // Determine the geometric altitudes at *lat*
267 const Index np = p_grid.nelem();
268 Vector z_grid(np);
269 ArrayOfGridPos gp_lat(1);
270 //
271 gridpos(gp_lat, lat_grid, ExhaustiveConstVectorView{lat});
272 z_at_lat_2d(z_grid, p_grid, lat_grid, z_field(joker, joker, 0), gp_lat[0]);
273
274 // Determine the ellipsoid radius at *lat*
275 const Numeric rellips = refell2d(refellipsoid, lat_grid, gp_lat[0]);
276
277 // Altitude (equal to pressure) grid position
278 ArrayOfGridPos gp_p(1);
279 gridpos(gp_p, z_grid, Vector(1, r - rellips));
280
281 // Altitude interpolation weights
282 Matrix itw(1, 2);
283 Vector dummy(1);
284 interpweights(itw, gp_p);
285
286 // Pressure
287 itw2p(dummy, p_grid, gp_p, itw);
288 rtp_pressure = dummy[0];
289
290 // Temperature
291 itw.resize(1, 4);
292 interpweights(itw, gp_p, gp_lat);
293 interp(dummy, itw, t_field(joker, joker, 0), gp_p, gp_lat);
294 rtp_temperature = dummy[0];
295
296 // VMR
297 const Index ns = vmr_field.nbooks();
298 //
299 rtp_vmr.resize(ns);
300 //
301 for (Index is = 0; is < ns; is++) {
302 interp(dummy, itw, vmr_field(is, joker, joker, 0), gp_p, gp_lat);
303 rtp_vmr[is] = dummy[0];
304 }
305
307 refr_index_air,
308 refr_index_air_group,
309 rtp_pressure,
310 rtp_temperature,
311 rtp_vmr,
312 Vector{f_grid},
313 refr_index_air_agenda);
314}
315
344 Numeric& refr_index_air,
345 Numeric& refr_index_air_group,
346 const Agenda& refr_index_air_agenda,
347 ConstVectorView p_grid,
348 ConstVectorView lat_grid,
349 ConstVectorView lon_grid,
350 ConstVectorView refellipsoid,
351 ConstTensor3View z_field,
352 ConstTensor3View t_field,
353 ConstTensor4View vmr_field,
354 ConstVectorView f_grid,
355 const Numeric& r,
356 const Numeric& lat,
357 const Numeric& lon) {
358 Numeric rtp_pressure, rtp_temperature;
359 Vector rtp_vmr;
360
361 // Determine the geometric altitudes at *lat* and *lon*
362 const Index np = p_grid.nelem();
363 Vector z_grid(np);
364 ArrayOfGridPos gp_lat(1), gp_lon(1);
365 //
366 gridpos(gp_lat, lat_grid, ExhaustiveConstVectorView{lat});
367 gridpos(gp_lon, lon_grid, ExhaustiveConstVectorView{lon});
369 z_grid, p_grid, lat_grid, lon_grid, z_field, gp_lat[0], gp_lon[0]);
370
371 // Determine the elipsoid radius at *lat*
372 const Numeric rellips = refell2d(refellipsoid, lat_grid, gp_lat[0]);
373
374 // Altitude (equal to pressure) grid position
375 ArrayOfGridPos gp_p(1);
376 gridpos(gp_p, z_grid, Vector(1, r - rellips));
377
378 // Altitude interpolation weights
379 Matrix itw(1, 2);
380 Vector dummy(1);
381 interpweights(itw, gp_p);
382
383 // Pressure
384 itw2p(dummy, p_grid, gp_p, itw);
385 rtp_pressure = dummy[0];
386
387 // Temperature
388 itw.resize(1, 8);
389 interpweights(itw, gp_p, gp_lat, gp_lon);
390 interp(dummy, itw, t_field, gp_p, gp_lat, gp_lon);
391 rtp_temperature = dummy[0];
392
393 // VMR
394 const Index ns = vmr_field.nbooks();
395 //
396 rtp_vmr.resize(ns);
397 //
398 for (Index is = 0; is < ns; is++) {
399 interp(
400 dummy, itw, vmr_field(is, joker, joker, joker), gp_p, gp_lat, gp_lon);
401 rtp_vmr[is] = dummy[0];
402 }
403
405 refr_index_air,
406 refr_index_air_group,
407 rtp_pressure,
408 rtp_temperature,
409 rtp_vmr,
410 Vector{f_grid},
411 refr_index_air_agenda);
412}
413
415
441 Numeric& refr_index_air,
442 Numeric& refr_index_air_group,
443 Numeric& dndr,
444 const Agenda& refr_index_air_agenda,
445 ConstVectorView p_grid,
446 ConstVectorView refellipsoid,
447 ConstTensor3View z_field,
448 ConstTensor3View t_field,
449 ConstTensor4View vmr_field,
450 ConstVectorView f_grid,
451 const Numeric& r) {
453 refr_index_air,
454 refr_index_air_group,
455 refr_index_air_agenda,
456 p_grid,
457 refellipsoid,
458 z_field,
459 t_field,
460 vmr_field,
461 f_grid,
462 r);
463
464 const Numeric n0 = refr_index_air;
465 Numeric dummy;
466
468 refr_index_air,
469 dummy,
470 refr_index_air_agenda,
471 p_grid,
472 refellipsoid,
473 z_field,
474 t_field,
475 vmr_field,
476 f_grid,
477 r + 1);
478
479 dndr = refr_index_air - n0;
480
481 refr_index_air = n0;
482}
483
485
518 Numeric& refr_index_air,
519 Numeric& refr_index_air_group,
520 Numeric& dndr,
521 Numeric& dndlat,
522 const Agenda& refr_index_air_agenda,
523 ConstVectorView p_grid,
524 ConstVectorView lat_grid,
525 ConstVectorView refellipsoid,
526 ConstTensor3View z_field,
527 ConstTensor3View t_field,
528 ConstTensor4View vmr_field,
529 ConstVectorView f_grid,
530 const Numeric& r,
531 const Numeric& lat) {
533 refr_index_air,
534 refr_index_air_group,
535 refr_index_air_agenda,
536 p_grid,
537 lat_grid,
538 refellipsoid,
539 z_field,
540 t_field,
541 vmr_field,
542 f_grid,
543 r,
544 lat);
545
546 const Numeric n0 = refr_index_air;
547 Numeric dummy;
548
550 refr_index_air,
551 dummy,
552 refr_index_air_agenda,
553 p_grid,
554 lat_grid,
555 refellipsoid,
556 z_field,
557 t_field,
558 vmr_field,
559 f_grid,
560 r + 1,
561 lat);
562
563 dndr = refr_index_air - n0;
564
565 const Numeric dlat = 1e-4;
566
568 refr_index_air,
569 dummy,
570 refr_index_air_agenda,
571 p_grid,
572 lat_grid,
573 refellipsoid,
574 z_field,
575 t_field,
576 vmr_field,
577 f_grid,
578 r,
579 lat + dlat);
580
581 dndlat = (refr_index_air - n0) / (DEG2RAD * dlat * r);
582
583 refr_index_air = n0;
584}
585
587
624 Numeric& refr_index_air,
625 Numeric& refr_index_air_group,
626 Numeric& dndr,
627 Numeric& dndlat,
628 Numeric& dndlon,
629 const Agenda& refr_index_air_agenda,
630 ConstVectorView p_grid,
631 ConstVectorView lat_grid,
632 ConstVectorView lon_grid,
633 ConstVectorView refellipsoid,
634 ConstTensor3View z_field,
635 ConstTensor3View t_field,
636 ConstTensor4View vmr_field,
637 ConstVectorView f_grid,
638 const Numeric& r,
639 const Numeric& lat,
640 const Numeric& lon) {
642 refr_index_air,
643 refr_index_air_group,
644 refr_index_air_agenda,
645 p_grid,
646 lat_grid,
647 lon_grid,
648 refellipsoid,
649 z_field,
650 t_field,
651 vmr_field,
652 f_grid,
653 r,
654 lat,
655 lon);
656
657 const Numeric n0 = refr_index_air;
658 Numeric dummy;
659
661 refr_index_air,
662 dummy,
663 refr_index_air_agenda,
664 p_grid,
665 lat_grid,
666 lon_grid,
667 refellipsoid,
668 z_field,
669 t_field,
670 vmr_field,
671 f_grid,
672 r + 1,
673 lat,
674 lon);
675
676 dndr = refr_index_air - n0;
677
678 const Numeric dlat = 1e-4;
679
681 refr_index_air,
682 dummy,
683 refr_index_air_agenda,
684 p_grid,
685 lat_grid,
686 lon_grid,
687 refellipsoid,
688 z_field,
689 t_field,
690 vmr_field,
691 f_grid,
692 r,
693 lat + dlat,
694 lon);
695
696 dndlat = (refr_index_air - n0) / (DEG2RAD * dlat * r);
697
698 const Numeric dlon = 1e-4;
699
701 refr_index_air,
702 dummy,
703 refr_index_air_agenda,
704 p_grid,
705 lat_grid,
706 lon_grid,
707 refellipsoid,
708 z_field,
709 t_field,
710 vmr_field,
711 f_grid,
712 r,
713 lat,
714 lon + dlon);
715
716 dndlon = (refr_index_air - n0) / (DEG2RAD * dlon * r * cos(DEG2RAD * lat));
717
718 refr_index_air = n0;
719}
720
722 const Index& only_valid_range,
723 const Numeric& frequency,
724 const Numeric& temperature,
725 const Numeric& density) {
726 //convert frequency to wavelength
727 const Numeric wavelength = Conversion::freq2wavelen(frequency) * 1e6; // [µm]
728
729 //Reference values
730 const Numeric T_star = 273.15; // [K]
731 const Numeric rho_star = 1000; // [kg/m^3]
732 const Numeric lambda_star = 0.589; // [µm]
733
734 //check input
735 bool T_ok = (temperature > T_star - 12) && (temperature < T_star + 500);
736 bool rho_ok = (density > 0) && (density < 1060);
737 bool wvl_ok = (wavelength > 0.2) && (wavelength < 1.9);
738
739 if (only_valid_range) {
740 ARTS_USER_ERROR_IF(!T_ok,
741 "Refractive index is calculated outside range of "
742 "validity \n",
743 "Temperature must be between ",
744 T_star - 12,
745 "K and ",
746 T_star + 500,
747 "K\n"
748 "Desired temperature: ",
749 temperature,
750 "K \n")
751
752 ARTS_USER_ERROR_IF(!rho_ok,
753 "Refractive index is calculated outside range of "
754 "validity \n",
755 "Density must be between ",
756 0,
757 "kg m^-3 and ",
758 1060,
759 "kg m^-3\n"
760 "Desired density: ",
761 density,
762 "kg m^-3 \n")
763
764 Numeric frq_upper_limit = Conversion::wavelen2freq(0.2 * 1e-6);
765 Numeric frq_lower_limit = Conversion::wavelen2freq(1.9 * 1e-6);
766
767 ARTS_USER_ERROR_IF(!wvl_ok,
768 "Refractive index is calculated outside range of "
769 "validity \n",
770 "Frequency must be between ",
771 frq_lower_limit,
772 "Hz and ",
773 frq_upper_limit,
774 "Hz\n"
775 "Desired density: ",
776 frequency,
777 "Hz \n")
778 }
779
780 //coefficients
781 const Numeric a0 = 0.244257733;
782 const Numeric a1 = 9.74634476e-3;
783 const Numeric a2 = -3.73234996e-3;
784 const Numeric a3 = 2.68678472e-4;
785 const Numeric a4 = 1.58920570e-3;
786 const Numeric a5 = 2.45934259e-3;
787 const Numeric a6 = 0.900704920;
788 const Numeric a7 = -1.66626219e-2;
789
790 const Numeric lambda_uv = 0.2292020;
791 const Numeric lambda_ir = 5.432937;
792
793 //normalize input variables
794 Numeric T_bar = temperature / T_star;
795 Numeric rho_bar = density / rho_star;
796 Numeric lambda_bar = wavelength / lambda_star;
797
798 //set up right hands side of eq A1
799 Numeric rhs = a0;
800 rhs += a1 * rho_bar;
801 rhs += a2 * T_bar;
802 rhs += a3 * lambda_bar * lambda_bar * T_bar;
803 rhs += a4 / (lambda_bar * lambda_bar);
804 rhs += a5 / (lambda_bar * lambda_bar - lambda_uv * lambda_uv);
805 rhs += a6 / (lambda_bar * lambda_bar - lambda_ir * lambda_ir);
806 rhs += a7 * rho_bar * rho_bar;
807
808 Complex a;
809 a = rho_bar * rhs;
810
811 //solve Eq A1 (see paper in documentation after n
812 Complex n_cmplx = sqrt(-2 * a - 1);
813 n_cmplx /= sqrt(a - 1);
814
815 //refractive index
816 n = - real(n_cmplx);
817}
base max(const Array< base > &x)
Max function.
Definition array.h:128
base min(const Array< base > &x)
Min function.
Definition array.h:144
Constants of physical expressions as constexpr.
Common ARTS conversions.
void refr_index_air_agendaExecute(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Vector &f_grid, const Agenda &input_agenda)
Definition auto_md.cc:27633
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
The Agenda class.
Workspace class.
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition geodetic.cc:1287
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
constexpr Numeric temperature_at_0c
Global constant, Temperature in Celsius of 0 Kelvin.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr auto freq2wavelen(auto x) noexcept
Conversion from Hz to wavelength.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
constexpr auto wavelen2freq(auto x) noexcept
Conversion from wavelength to Hz.
void refr_gradients_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, Numeric &dndlon, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat, const Numeric &lon)
refr_gradients_3d
void get_refr_index_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat, const Numeric &lon)
void complex_n_ice_matzler06(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_ice_matzler06
void get_refr_index_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r)
get_refr_index_1d
void refr_gradients_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat)
refr_gradients_2d
constexpr Numeric DEG2RAD
Definition refraction.cc:28
void get_refr_index_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat)
get_refr_index_2d
constexpr Numeric RAD2DEG
Definition refraction.cc:29
void refr_gradients_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r)
refr_gradients_1d
void complex_n_water_liebe93(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_water_liebe93
Definition refraction.cc:57
void refractive_index_water_and_steam_VisNIR(Numeric &n, const Index &only_valid_range, const Numeric &frequency, const Numeric &temperature, const Numeric &density)
Refractive index of water and steam for the optical and near infrared .
constexpr Numeric TEMP_0_C
Definition refraction.cc:30
Refraction functions.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void z_at_latlon(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, const GridPos &gp_lat, const GridPos &gp_lon)
Returns the geomtrical altitudes of p_grid for one latitude and one longitude.
void z_at_lat_2d(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstMatrixView z_field, const GridPos &gp_lat)
Returns the geomtrical altitudes of p_grid for one latitude.
Header file for special_interp.cc.
#define a
#define b