icedb  version 0.5.1
Snow particle scattering database API
refractImpl.cpp
Go to the documentation of this file.
1 #define _SCL_SECURE_NO_WARNINGS // Issue with linterp
2 #pragma warning( disable : 4244 ) // ICEDB_LOG_WARNINGs C4244 and C4267: size_t to int and int <-> _int64
3 #pragma warning( disable : 4267 )
4 #include <cmath>
5 #include <complex>
6 #include <valarray>
7 #include <mutex>
8 #include "../icedb/refract/refract.hpp"
9 #include "../icedb/refract/refractBase.hpp"
10 #include "../icedb/zeros.hpp"
11 #include "../icedb/units/units.hpp"
12 #include "../private/linterp.h"
13 #include "../icedb/error.hpp"
14 #include "../icedb/logging.hpp"
15 
16 namespace {
17  std::mutex m_setup;
18 
19  enum class hanelAmedium { WATER_RE, WATER_IM, ICE_RE, ICE_IM, NACL_RE, NACL_IM, SEASALT_RE, SEASALT_IM };
20  std::shared_ptr<InterpMultilinear<1, double > > setupHanelA(hanelAmedium m)
21  {
22  std::lock_guard<std::mutex> lock(m_setup);
23  static std::shared_ptr<InterpMultilinear<1, double > > res_water_re, res_water_im,
24  res_ice_re, res_ice_im, res_nacl_re, res_nacl_im, res_seasalt_re, res_seasalt_im;
25  if (!res_water_re) {
26 
27  const size_t nWvlengths = 90;
28  const double Dwavelengths[nWvlengths] = { 0.2, 0.25, 0.3, 0.337, 0.4, 0.488, 0.515,
29  0.55, 0.633, 0.694, 0.86, 1.06, 1.3, 1.536, 1.8, 2, 2.25, 2.5,
30  2.7, 3, 3.2, 3.392, 3.5, 3.75, 4, 4.5, 5, 5.5, 6, 6.2, 6.5,
31  7.2, 7.9, 8.2, 8.5, 8.7, 9, 9.2, 9.5, 9.8, 10, 10.591, 11,
32  11.5, 12.5, 13, 14, 14.8, 15, 16.4, 17.2, 18, 18.5, 20, 21.3,
33  22.5, 25, 27.9, 30, 35, 40, 45, 50, 55, 60, 65, 70, 80, 90,
34  100, 110, 120, 135, 150, 165, 180, 200, 250, 300, 400, 500,
35  750, 1000, 1500, 2000, 3000, 5000, 10000, 20000, 30000 };
36  static std::vector<double > wavelengths(nWvlengths);
37  for (size_t i = 0; i < nWvlengths; ++i) wavelengths[i] = Dwavelengths[i];
38  const double water_re[nWvlengths] = { 1.396, 1.362, 1.349, 1.345, 1.339, 1.335,
39  1.334, 1.333, 1.332, 1.331, 1.329, 1.326, 1.323, 1.318, 1.312, 1.306,
40  1.292, 1.261, 1.188, 1.371, 1.478, 1.422, 1.4, 1.369, 1.351, 1.332,
41  1.325, 1.298, 1.265, 1.363, 1.339, 1.312, 1.294, 1.286, 1.278, 1.272,
42  1.262, 1.255, 1.243, 1.229, 1.218, 1.179, 1.153, 1.126, 1.123, 1.146,
43  1.21, 1.258, 1.27, 1.346, 1.386, 1.423, 1.443, 1.48, 1.491, 1.506, 1.531,
44  1.549, 1.551, 1.532, 1.519, 1.536, 1.587, 1.645, 1.703, 1.762, 1.821,
45  1.92, 1.979, 2.037, 2.06, 2.082, 2.094, 2.106, 2.109, 2.113, 2.117,
46  2.12, 2.121, 2.142, 2.177, 2.291, 2.437, 2.562, 2.705, 3.013, 3.627,
47  4.954, 6.728, 7.682 };
48  const double water_im[nWvlengths] = { 1.10E-07, 3.35E-08, 1.60E-08, 8.45E-09,
49  1.86E-09, 9.69E-10, 1.18E-09, 1.96E-09, 1.46E-08, 3.05E-08, 3.29E-07,
50  4.18E-06, 3.69E-05, 9.97E-05, 1.15E-04, 1.10E-03, 3.90E-04, 0.00174,
51  0.019, 0.272, 0.0924, 0.0204, 0.0094, 0.0035, 0.0046, 0.0134, 0.0124,
52  0.0116, 0.107, 0.088, 0.0392, 0.0321, 0.0339, 0.0351, 0.0367, 0.0379,
53  0.0399, 0.0415, 0.0444, 0.0479, 0.0508, 0.0674, 0.0968, 0.142, 0.259,
54  0.305, 0.37, 0.396, 0.402, 0.427, 0.429, 0.426, 0.421, 0.393, 0.379,
55  0.37, 0.356, 0.339, 0.328, 0.336, 0.385, 0.449, 0.514, 0.551, 0.587,
56  0.582, 0.576, 0.52, 0.49, 0.46, 0.44, 0.42, 0.41, 0.4, 0.41, 0.42,
57  0.43, 0.46, 0.49, 0.53, 0.58, 0.65, 0.73, 0.96, 1.19, 1.59, 2.14, 2.79,
58  2.87, 2.51 };
59  //static std::vector<double > water_re(nWvlengths);
60  //for (size_t i = 0; i < nWvlengths; ++i) water_re[i] = water_re[i], water_im[i]);
61  const double ice_re[nWvlengths] = { 1.394, 1.351, 1.334, 1.326, 1.32, 1.313, 1.312,
62  1.311, 1.308, 1.306, 1.303, 1.3, 1.295, 1.29, 1.282, 1.273, 1.256, 1.225,
63  1.163, 1.045, 1.652, 1.51, 1.453, 1.391, 1.361, 1.34, 1.327, 1.299, 1.296,
64  1.313, 1.32, 1.318, 1.313, 1.306, 1.291, 1.282, 1.269, 1.261, 1.245, 1.219,
65  1.197, 1.098, 1.093, 1.176, 1.387, 1.472, 1.569, 1.579, 1.572, 1.531, 1.534,
66  1.522, 1.51, 1.504, 1.481, 1.455, 1.414, 1.358, 1.325, 1.226, 1.202, 1.299,
67  1.629, 1.767, 1.585, 1.748, 1.869, 1.903, 1.856, 1.832, 1.821, 1.819, 1.819,
68  1.823, 1.829, 1.832, 1.827, 1.807, 1.795, 1.786, 1.783, 1.782, 1.781, 1.782,
69  1.782, 1.782, 1.783, 1.784, 1.784, 1.784 };
70  const double ice_im[nWvlengths] = { 1.50E-08, 8.60E-09, 5.50E-09, 4.50E-09, 2.71E-09,
71  1.75E-09, 2.19E-09, 3.11E-09, 1.09E-08, 2.62E-08, 2.15E-07, 1.96E-06,
72  1.32E-05, 6.10E-04, 1.13E-04, 1.61E-03, 2.13E-04, 7.95E-04, 0.00293, 0.429,
73  0.283, 0.0401, 0.0161, 0.007, 0.01, 0.0287, 0.012, 0.0217, 0.0647, 0.0683,
74  0.0559, 0.0544, 0.0479, 0.039, 0.0391, 0.04, 0.0429, 0.0446, 0.0459, 0.047,
75  0.051, 0.131, 0.239, 0.36, 0.422, 0.389, 0.283, 0.191, 0.177, 0.125, 0.107,
76  0.0839, 0.076, 0.067, 0.0385, 0.0291, 0.0299, 0.049, 0.065, 0.155, 0.344,
77  0.601, 0.543, 0.42, 0.39, 0.49, 0.399, 0.235, 0.165, 0.139, 0.126, 0.12,
78  0.108, 0.0962, 0.0846, 0.065, 0.0452, 0.0222, 0.0153, 0.0125, 0.0106, 0.008,
79  0.0065, 0.0049, 0.004, 0.003, 0.0021, 0.0013, 7.90E-04, 5.90E-04 };
80  //static std::vector<std::complex<double> > ice(nWvlengths);
81  //for (size_t i = 0; i < nWvlengths; ++i) ice[i] = std::complex<double>(ice_re[i], ice_im[i]);
82  const double nacl_re[nWvlengths] = { 1.79, 1.655, 1.607, 1.587, 1.567, 1.553, 1.55,
83  1.547, 1.542, 1.539, 1.534, 1.531, 1.529, 1.528, 1.527, 1.527, 1.526, 1.525,
84  1.525, 1.524, 1.524, 1.523, 1.523, 1.522, 1.522, 1.52, 1.519, 1.517, 1.515,
85  1.515, 1.513, 1.51, 1.507, 1.505, 1.504, 1.503, 1.501, 1.5, 1.498, 1.496,
86  1.495, 1.491, 1.488, 1.484, 1.476, 1.471, 1.462, 1.454, 1.451, 1.435, 1.425,
87  1.414, 1.406, 1.382, 1.36, 1.33, 1.27, 1.17, 1.08, 0.78, 0.58, 0.27, 0.14,
88  0.31, 4.52, 5.28, 3.92, 3.17, 2.87, 2.74, 2.64, 2.59, 2.54, 2.5, 2.48, 2.47,
89  2.45, 2.44, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43 };
90  const double nacl_im[nWvlengths] = { 3.10E-09, 2.30E-09, 1.50E-09, 8.70E-10, 3.80E-10,
91  1.10E-10, 4.90E-11, 6.80E-11, 1.10E-10, 1.50E-10, 2.40E-10, 3.50E-10,
92  4.80E-10, 6.10E-10, 7.50E-10, 8.60E-10, 1.00E-09, 1.10E-09, 1.20E-09,
93  1.40E-09, 1.50E-09, 1.60E-09, 1.65E-09, 1.80E-09, 1.80E-09, 1.80E-09,
94  1.70E-09, 2.60E-09, 4.90E-09, 5.80E-09, 7.20E-09, 1.00E-08, 1.40E-08,
95  1.50E-08, 1.60E-08, 1.70E-08, 1.90E-08, 2.00E-08, 3.00E-08, 4.40E-08,
96  5.30E-08, 8.00E-08, 1.30E-07, 3.30E-07, 1.40E-06, 2.80E-06, 8.80E-06,
97  2.30E-05, 2.70E-05, 7.60E-05, 1.30E-04, 2.00E-04, 2.90E-04, 6.20E-04,
98  9.90E-04, 0.0014, 0.0035, 0.01, 0.026, 0.14, 0.66, 1.08, 1.99, 3.46, 6.94,
99  0.761, 0.271, 0.123, 0.0968, 0.087, 0.079, 0.077, 0.072, 0.064, 0.056, 0.052,
100  0.047, 0.041, 0.03, 0.027, 0.024, 0.012, 0.008, 0.0061, 0.0047, 0.0029,
101  0.0024, 5.60E-04, 4.10E-04, 2.60E-04 };
102  //static std::vector<std::complex<double> > nacl(nWvlengths);
103  //for (size_t i = 0; i < nWvlengths; ++i) nacl[i] = std::complex<double>(nacl_re[i], nacl_im[i]);
104  const double seasalt_re[nWvlengths] = { 1.51, 1.51, 1.51, 1.51, 1.5, 1.5, 1.5,
105  1.5, 1.49, 1.49, 1.48, 1.47, 1.47, 1.46, 1.45, 1.45, 1.44, 1.43, 1.4,
106  1.61, 1.49, 1.48, 1.48, 1.47, 1.48, 1.49, 1.47, 1.42, 1.41, 1.6, 1.46,
107  1.42, 1.4, 1.42, 1.48, 1.6, 1.65, 1.61, 1.58, 1.56, 1.54, 1.5, 1.48,
108  1.48, 1.42, 1.41, 1.41, 1.43, 1.45, 1.56, 1.74, 1.78, 1.77, 1.76, 1.76,
109  1.76, 1.76, 1.77, 1.77, 1.76, 1.74, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
110  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
111  const double seasalt_im[nWvlengths] = { 1.00E-04, 5.00E-06, 2.00E-06, 4.00E-07,
112  3.00E-08, 2.00E-08, 1.00E-08, 1.00E-08, 2.00E-08, 1.00E-07, 3.00E-06,
113  2.00E-04, 4.00E-04, 6.00E-04, 8.00E-04, 0.001, 0.002, 0.004, 0.007,
114  0.01, 0.003, 0.002, 0.0016, 0.0014, 0.0014, 0.0014, 0.0025, 0.0036,
115  0.011, 0.022, 0.005, 0.007, 0.013, 0.02, 0.026, 0.03, 0.028, 0.026,
116  0.018, 0.016, 0.015, 0.014, 0.014, 0.014, 0.016, 0.018, 0.023, 0.03,
117  0.035, 0.09, 0.12, 0.13, 0.135, 0.152, 0.165, 0.18, 0.205, 0.275, 0.3,
118  0.5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
119  0, 0, 0, 0, 0, 0, 0, 0 };
120  //static std::vector<std::complex<double> > seasalt(nWvlengths);
121  //for (size_t i = 0; i < nWvlengths; ++i) seasalt[i] = std::complex<double>(seasalt_re[i], seasalt_im[i]);
122 
123  static std::vector< std::vector<double>::iterator > grid_iter_list;
124  grid_iter_list.push_back(wavelengths.begin());
125 
126  // the size of the grid in each dimension
127  static std::array<size_t, 1> grid_sizes;
128  grid_sizes[0] = wavelengths.size();
129  // total number of elements
130  static size_t num_elements = grid_sizes[0];// *grid_sizes[1];
131 
132 
133  res_water_re = std::shared_ptr<InterpMultilinear<1, double > >
134  (new InterpMultilinear<1, double >
135  (grid_iter_list.begin(), grid_sizes.begin(), water_re, water_re + num_elements));
136  res_water_im = std::shared_ptr<InterpMultilinear<1, double > >
137  (new InterpMultilinear<1, double >
138  (grid_iter_list.begin(), grid_sizes.begin(), water_im, water_im + num_elements));
139  res_ice_re = std::shared_ptr<InterpMultilinear<1, double > >
140  (new InterpMultilinear<1, double >
141  (grid_iter_list.begin(), grid_sizes.begin(), ice_re, ice_re + num_elements));
142  res_ice_im = std::shared_ptr<InterpMultilinear<1, double > >
143  (new InterpMultilinear<1, double >
144  (grid_iter_list.begin(), grid_sizes.begin(), ice_im, ice_im + num_elements));
145  res_nacl_re = std::shared_ptr<InterpMultilinear<1, double > >
146  (new InterpMultilinear<1, double >
147  (grid_iter_list.begin(), grid_sizes.begin(), nacl_re, nacl_re + num_elements));
148  res_nacl_im = std::shared_ptr<InterpMultilinear<1, double > >
149  (new InterpMultilinear<1, double >
150  (grid_iter_list.begin(), grid_sizes.begin(), nacl_im, nacl_im + num_elements));
151  res_seasalt_re = std::shared_ptr<InterpMultilinear<1, double > >
152  (new InterpMultilinear<1, double >
153  (grid_iter_list.begin(), grid_sizes.begin(), seasalt_re, seasalt_re + num_elements));
154  res_seasalt_im = std::shared_ptr<InterpMultilinear<1, double > >
155  (new InterpMultilinear<1, double >
156  (grid_iter_list.begin(), grid_sizes.begin(), seasalt_im, seasalt_im + num_elements));
157  }
158  if (m == hanelAmedium::WATER_RE) return res_water_re;
159  else if (m == hanelAmedium::WATER_IM) return res_water_im;
160  else if (m == hanelAmedium::ICE_RE) return res_ice_re;
161  else if (m == hanelAmedium::ICE_IM) return res_ice_im;
162  else if (m == hanelAmedium::NACL_RE) return res_nacl_re;
163  else if (m == hanelAmedium::NACL_IM) return res_nacl_im;
164  else if (m == hanelAmedium::SEASALT_RE) return res_seasalt_re;
165  else if (m == hanelAmedium::SEASALT_IM) return res_seasalt_im;
166  else ICEDB_throw(icedb::error::error_types::xUnimplementedFunction)
167  .add<std::string>("Reason", "Bad enum input for hanelAmedium. "
168  "This is not a "
169  "case that the Hanel dielectric code can handle.");
170  return nullptr;
171  }
172  enum class hanelBmedium { SAND_O_RE, SAND_O_IM, SAND_E_RE, SAND_E_IM, DUST_LIKE_RE, DUST_LIKE_IM };
173  std::shared_ptr<InterpMultilinear<1, double > > setupHanelB(hanelBmedium m)
174  {
175  std::lock_guard<std::mutex> lock(m_setup);
176  static std::shared_ptr<InterpMultilinear<1, double> > res_sand_o_re, res_sand_o_im,
177  res_sand_e_re, res_sand_e_im, res_dust_re, res_dust_im;
178  if (!res_sand_o_re) {
179 
180  const size_t nWvlengths = 68, nWvlengthsDust = 61;
181  const double D_wavelengths[nWvlengths] = { 0.2, 0.25, 0.3, 0.337, 0.4, 0.488,
182  0.515, 0.55, 0.633, 0.694, 0.86, 1.06, 1.3, 1.536, 1.8, 2, 2.25, 2.5,
183  2.7, 3, 3.2, 3.392, 3.5, 3.75, 4, 4.5, 5, 5.5, 6, 6.2, 6.5, 7.2, 7.9,
184  8.2, 8.5, 8.7, 9, 9.2, 9.5, 9.8, 10, 10.591, 11, 11.5, 12.5, 13, 14,
185  14.8, 15, 16.4, 17.2, 18, 18.5, 20, 21.3, 22.5, 25, 27.9, 30, 35, 40,
186  50, 60, 80, 100, 150, 200, 300 };
187  static std::vector<double > wavelengths(nWvlengths);
188  for (size_t i = 0; i < nWvlengths; ++i) wavelengths[i] = D_wavelengths[i];
189  const double sand_o_re[nWvlengths] = { 1.665, 1.686, 1.677, 1.665, 1.656,
190  1.667, 1.666, 1.664, 1.655, 1.65, 1.635, 1.627, 1.623, 1.62, 1.615,
191  1.611, 1.607, 1.602, 1.597, 1.591, 1.585, 1.579, 1.576, 1.566, 1.555,
192  1.535, 1.511, 1.466, 1.42, 1.389, 1.342, 1.161, 0.662, 0.175, 0.356,
193  0.854, 0.221, 0.612, 4.274, 2.933, 2.615, 2.18, 2.015, 1.844, 1.564,
194  2.119, 1.715, 1.727, 1.669, 1.163, 2.011, 1.584, 1.346, 0.145, 0.346,
195  5.37, 0.946, 2.162, 3.446, 2.655, 2.491, 2.446, 2.388, 2.347, 2.326,
196  2.31, 2.299, 2.296 };
197  const double sand_o_im[nWvlengths] = { 0.147, 0.117, 0.09, 0.0778, 0.0335,
198  0.0109, 0.00809, 0.00474, 5.22E-04, 5.61E-05, 2.46E-04, 1.92E-06,
199  6.51E-07, 3.63E-07, 2.67E-07, 3.00E-07, 3.31E-07, 8.58E-07, 1.19E-06,
200  2.60E-05, 7.13E-06, 6.78E-06, 8.66E-06, 5.75E-05, 7.10E-05, 4.86E-04,
201  0.00539, 0.00533, 0.00656, 0.00764, 0.00571, 0.0124, 0.0922, 0.632,
202  1.82, 2.19, 2.06, 4.07, 0.356, 0.0789, 0.0474, 0.0214, 0.0172, 0.019,
203  2.02, 0.0515, 0.0394, 0.0381, 0.0348, 0.233, 0.655, 0.0912, 0.0603,
204  0.886, 2.59, 1.35, 1.84, 0.382, 0.353, 0.0243, 0.0107, 0.00613, 0.00428,
205  0.00494, 0.00187, 0.00143, 0.00121, 0.00106 };
206  //static std::vector<std::complex<double> > sand_o(nWvlengths);
207  //for (size_t i = 0; i < nWvlengths; ++i) sand_o[i] = std::complex<double>(sand_o_re[i], sand_o_im[i]);
208  const double sand_e_re[nWvlengths] = { 1.665, 1.686, 1.677, 1.665, 1.656,
209  1.667, 1.666, 1.664, 1.655, 1.65, 1.635, 1.627, 1.623, 1.62, 1.615, 1.611,
210  1.607, 1.602, 1.597, 1.591, 1.585, 1.579, 1.576, 1.566, 1.555, 1.535,
211  1.511, 1.466, 1.42, 1.389, 1.342, 1.16, 0.639, 0.155, 0.239, 1.557, 0.251,
212  1.59, 3.733, 2.794, 2.533, 2.151, 2.003, 1.859, 1.305, 2.524, 1.762, 1.58,
213  1.533, 0.955, 1.404, 0.489, 0.229, 1.389, 2.853, 2.698, 1.559, 4.4, 4.196,
214  2.836, 2.606, 2.495, 2.436, 2.389, 2.37, 2.352, 2.347, 2.343 };
215  const double sand_e_im[nWvlengths] = { 0.147, 0.117, 0.09, 0.0778, 0.0335,
216  0.0109, 0.00809, 0.00474, 5.22E-04, 5.61E-05, 2.46E-04, 1.92E-06, 6.51E-07,
217  3.63E-07, 2.67E-07, 3.00E-07, 3.31E-07, 1.09E-06, 1.59E-06, 8.69E-06,
218  4.53E-06, 6.12E-06, 8.66E-06, 3.58E-05, 7.10E-05, 3.45E-04, 0.00427, 0.00503,
219  0.00468, 0.00764, 0.00717, 0.0146, 0.0938, 0.693, 1.74, 0.709, 2.56, 5.85,
220  0.207, 0.0629, 0.0405, 0.0198, 0.016, 0.0159, 0.123, 0.291, 0.0229, 0.0293,
221  0.0333, 0.481, 0.153, 0.249, 1.1, 4.57, 0.448, 0.295, 0.154, 0.915, 1.41,
222  0.0369, 0.0151, 0.00509, 0.00319, 0.00192, 0.0014, 8.24E-04, 5.88E-04, 3.43E-04 };
223  //static std::vector<std::complex<double> > sand_e(nWvlengths);
224  //for (size_t i = 0; i < nWvlengths; ++i) sand_e[i] = std::complex<double>(sand_e_re[i], sand_e_im[i]);
225  const double dust_re[nWvlengthsDust] = { 1.53, 1.53, 1.53, 1.53, 1.53, 1.53,
226  1.53, 1.53, 1.53, 1.53, 1.52, 1.52, 1.46, 1.4, 1.33, 1.26, 1.22, 1.18, 1.18,
227  1.16, 1.22, 1.26, 1.28, 1.27, 1.26, 1.26, 1.25, 1.22, 1.15, 1.14, 1.13, 1.4,
228  1.15, 1.13, 1.3, 1.4, 1.7, 1.72, 1.73, 1.74, 1.75, 1.62, 1.62, 1.59, 1.51,
229  1.47, 1.52, 1.57, 1.57, 1.6, 1.63, 1.64, 1.64, 1.68, 1.77, 1.9, 1.97, 1.89,
230  1.8, 1.9, 2.1 };
231  const double dust_im[nWvlengthsDust] = { 0.07, 0.03, 0.008, 0.008, 0.008, 0.008,
232  0.008, 0.008, 0.008, 0.008, 0.008, 0.008, 0.008, 0.008, 0.008, 0.008, 0.009,
233  0.009, 0.013, 0.012, 0.01, 0.013, 0.011, 0.011, 0.012, 0.014, 0.016, 0.021,
234  0.037, 0.039, 0.042, 0.055, 0.04, 0.074, 0.09, 0.1, 0.14, 0.15, 0.162, 0.162,
235  0.162, 0.12, 0.105, 0.1, 0.09, 0.1, 0.085, 0.1, 0.1, 0.1, 0.1, 0.115, 0.12,
236  0.22, 0.28, 0.28, 0.24, 0.32, 0.42, 0.5, 0.6 };
237  //static std::vector<std::complex<double> > dust(nWvlengthsDust);
238  //for (size_t i = 0; i < nWvlengthsDust; ++i) dust[i] = std::complex<double>(dust_re[i], dust_im[i]);
239 
240  static std::vector< std::vector<double>::iterator > grid_iter_list;
241  grid_iter_list.push_back(wavelengths.begin());
242 
243  // the size of the grid in each dimension
244  static std::array<size_t, 1> grid_sizes;
245  grid_sizes[0] = wavelengths.size();
246  // total number of elements
247  static size_t num_elements = grid_sizes[0];// *grid_sizes[1];
248 
249  static std::array<size_t, 1> grid_sizes_small;
250  grid_sizes_small[0] = nWvlengthsDust;
251  // total number of elements
252  static size_t num_elements_small = grid_sizes_small[0];// *grid_sizes[1];
253 
254  res_sand_o_re = std::shared_ptr<InterpMultilinear<1, double > >
255  (new InterpMultilinear<1, double >
256  (grid_iter_list.begin(), grid_sizes.begin(), sand_o_re, sand_o_re + num_elements));
257  res_sand_e_re = std::shared_ptr<InterpMultilinear<1, double > >
258  (new InterpMultilinear<1, double >
259  (grid_iter_list.begin(), grid_sizes.begin(), sand_e_re, sand_e_re + num_elements));
260  res_dust_re = std::shared_ptr<InterpMultilinear<1, double > >
261  (new InterpMultilinear<1, double >
262  (grid_iter_list.begin(), grid_sizes_small.begin(), dust_re, dust_re + num_elements_small));
263  res_sand_o_im = std::shared_ptr<InterpMultilinear<1, double > >
264  (new InterpMultilinear<1, double >
265  (grid_iter_list.begin(), grid_sizes.begin(), sand_o_im, sand_o_im + num_elements));
266  res_sand_e_im = std::shared_ptr<InterpMultilinear<1, double > >
267  (new InterpMultilinear<1, double >
268  (grid_iter_list.begin(), grid_sizes.begin(), sand_e_im, sand_e_im + num_elements));
269  res_dust_im = std::shared_ptr<InterpMultilinear<1, double > >
270  (new InterpMultilinear<1, double >
271  (grid_iter_list.begin(), grid_sizes_small.begin(), dust_im, dust_im + num_elements_small));
272  }
273  if (m == hanelBmedium::SAND_O_RE) return res_sand_o_re;
274  else if (m == hanelBmedium::SAND_O_IM) return res_sand_o_im;
275  else if (m == hanelBmedium::SAND_E_RE) return res_sand_e_re;
276  else if (m == hanelBmedium::SAND_E_IM) return res_sand_e_im;
277  else if (m == hanelBmedium::DUST_LIKE_RE) return res_dust_re;
278  else if (m == hanelBmedium::DUST_LIKE_IM) return res_dust_im;
279  else ICEDB_throw(icedb::error::error_types::xUnimplementedFunction)
280  .add<std::string>("Reason", "Bad enum input for hanelBmedium. "
281  "This is not a "
282  "case that the Hanel dielectric code can handle.");
283  return nullptr;
284  }
285 }
286 
287 // Water complex refractive index
288 // from Liu's mcx.f
289 // LIEBE, HUFFORD AND MANABE, INT. J. IR & MM WAVES V.12, pp.659-675
290 // (1991); Liebe et al, AGARD Conf. Proc. 542, May 1993.
291 // Valid from 0 to 1000 GHz. freq in GHz, temp in K
292 void icedb::refract::implementations::mWaterLiebe(double f, double t, std::complex<double> &m)
293 {
294  if (f < 0 || f > 1000 || t > 273.15)
295  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
296  .add<double>("Frequency (GHz)", f)
297  .add<double>("Temperature (K)", t)
298  .add<std::string>("Reason", "Allowed freq. range (GHz) is (0,1000), and allowed temp. range is <= 273.15 K.");
299 
300  double theta1 = 1.0 - (300.0 / t);
301  double eps0 = 77.66 - (103.3*theta1);
302  double eps1 = .0671*eps0;
303  double eps2 = 3.52;
304  double fp = (316.*theta1 + 146.4)*theta1 + 20.20;
305  double fs = 39.8*fp;
306  using namespace std;
307  complex<double> eps;
308  eps = std::complex<double>(eps0 - eps1, 0) / std::complex<double>(1.0, f / fp)
309  + std::complex<double>(eps1 - eps2, 0) / std::complex<double>(1.0, f / fs)
310  + std::complex<double>(eps2, 0);
311  m = sqrt(eps);
312  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2,
313  "mWaterLiebe result for freq: " << f << " temp " << t << " is " << m);
314 }
315 
316 void icedb::refract::implementations::mWaterFreshMeissnerWentz(double f, double tK, std::complex<double> &m)
317 {
318  if (f < 0 || f > 500 || tK < 273.15)
319  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
320  .add<double>("Frequency (GHz)", f)
321  .add<double>("Temperature (K)", tK)
322  .add<std::string>("Reason", "Allowed freq. range (GHz) is (0,500), and temp. must be >= 273.15 K.");
323 
324  const double as[11] = {
325  5.7230, 0.022379, -0.00071237, 5.0478,
326  -0.070315, 0.00060059, 3.6143, 0.028841,
327  0.13652, 0.0014825, 0.00024166
328  };
329 
330  double tC = tK - 273.15;
331  if (tC < -20 || tC > 40)
332  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
333  .add<double>("Frequency (GHz)", f)
334  .add<double>("Temperature (K)", tK)
335  .add<std::string>("Reason", "Allowed temp. range (C) is (-20,40)");
336 
337 
338  // static dielectric constant for pure water (Stogryn)
339  double es = (37088.6 - (82.168*tC)) / (tC + 421.854);
340 
341  double e1 = as[0] + (as[1] * tC) + (as[2] * tC*tC);
342  double nu1 = (45 + tC) / (as[3] + (tC*as[4]) + (tC*tC*as[5]));
343  double einf = as[6] + (tC*as[7]);
344  double nu2 = (45 + tC) / (as[8] + (tC*as[9]) + (tC*tC*as[10]));
345  double sigma = 0;
346  // vacuum electric permittivity
347  const double oneover2pie0 = 17.97510; // GHz m / S
348 
349  using namespace std;
350  complex<double> eps;
351  eps = std::complex<double>(es - e1, 0) / (complex<double>(1, f / nu1));
352  eps += std::complex<double>(e1 - einf, 0) / (complex<double>(1, f / nu2));
353  eps += std::complex<double>(einf, -sigma * oneover2pie0 / f); // -sigma / (f*2.*pi*e0));
354  m = sqrt(eps);
355  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2,
356  "mWaterFreshMeissnerWentz result for freq: " << f << " GHz temp " << tK << " K is " << m);
357 }
358 
359 
360 void icedb::refract::implementations::mIceMatzler(double f, double t, std::complex<double> &m)
361 {
362  if (f < 0 || f > 1000 || t > 273.15)
363  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
364  .add<double>("Frequency (GHz)", f)
365  .add<double>("Temperature (K)", t)
366  .add<std::string>("Reason", "Allowed freq. range (GHz) is (0,1000), and allowed temp. range is <= 273.15 K.");
367 
368  double er = 0;
369  if (t>243.0)
370  er = 3.1884 + 9.1e-4*(t - 273.0);
371  else
372  er = 3.1611 + 4.3e-4*(t - 243.0);
373  // Imaginary part
374  double theta = 300.0 / (t)-1.0;
375  double alpha = (0.00504 + 0.0062*theta)*exp(-22.1*theta);
376  double dbeta = exp(-9.963 + 0.0372*(t - 273.16));
377  const double B1 = 0.0207;
378  const double B2 = 1.16e-11;
379  const double b = 335;
380  double betam = B1 / t*exp(b / t) / pow((exp(b / t) - 1.0), 2) + B2*pow(f, 2.0);
381  double beta = betam + dbeta;
382  double ei = alpha / f + beta*f;
383  std::complex<double> e(er, -ei);
384  m = sqrt(e);
385  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2,
386  "mIceMatzler result for freq: " << f << " GHz and temp " << t << " K is " << m);
387 }
388 
389 void icedb::refract::implementations::mIceWarren(double f, double t, std::complex<double> &m)
390 {
391  if (f < 0.167 || f > 8600 || t > -1 || t < -60)
392  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
393  .add<double>("Frequency (GHz)", f)
394  .add<double>("Temperature (degC)", t)
395  .add<std::string>("Reason", "Allowed freq. range (GHz) is (0.167,8600), and allowed temp. range is -60 - -1 deg C.");
396 
397  // Warren table 2 is used for interpolation
398  static bool setup = false;
399  static std::vector<double> tempCs, wavelengths;
400  static std::vector<double > vals_re, vals_im;
401  static std::shared_ptr<InterpMultilinear<2, double > > interp_ML_warren_re, interp_ML_warren_im;
402  auto setupWarren = [&]()
403  {
404  std::lock_guard<std::mutex> lock(m_setup);
405  if (setup) return;
406 
407  tempCs.push_back(-1); tempCs.push_back(-5); tempCs.push_back(-20); tempCs.push_back(-60);
408  wavelengths.reserve(62);
409  vals_re.reserve(600); vals_im.reserve(600);
410 
411  // First wavelength is at tbl[0], last is at tbl[549]
412  static const double tbl[] = {
413  .1670, 1.8296, 8.30e-2, 1.8296, 8.30e-2, 1.8296, 8.30e-2, 1.8296, 8.30e-2,
414  .1778, 1.8236, 6.90e-2, 1.8236, 6.90e-2, 1.8236, 6.90e-2, 1.8236, 6.90e-2,
415  .1884, 1.8315, 5.70e-2, 1.8315, 5.70e-2, 1.8315, 5.70e-2, 1.8315, 5.70e-2,
416  .1995, 1.8275, 4.56e-2, 1.8275, 4.56e-2, 1.8275, 4.56e-2, 1.8275, 4.45e-2,
417  .2113, 1.8222, 3.79e-2, 1.8222, 3.79e-2, 1.8222, 3.79e-2, 1.8222, 3.55e-2,
418  .2239, 1.8172, 3.14e-2, 1.8172, 3.14e-2, 1.8172, 3.14e-2, 1.8172, 2.91e-2,
419  .2371, 1.8120, 2.62e-2, 1.8120, 2.62e-2, 1.8120, 2.62e-2, 1.8120, 2.44e-2,
420  .2512, 1.8070, 2.24e-2, 1.8070, 2.24e-2, 1.8070, 2.19e-2, 1.8070, 1.97e-2,
421  .2661, 1.8025, 1.96e-2, 1.8025, 1.96e-2, 1.8025, 1.88e-2, 1.8025, 1.67e-2,
422  .2818, 1.7983, 1.76e-2, 1.7983, 1.76e-2, 1.7983, 1.66e-2, 1.7983, 1.40e-2,
423  .2985, 1.7948, 1.67e-2, 1.7948, 1.67e-2, 1.7948, 1.54e-2, 1.7948, 1.26e-2,
424  .3162, 1.7921, 1.62e-2, 1.7921, 1.60e-2, 1.7921, 1.47e-2, 1.7921, 1.08e-2,
425  .3548, 1.7884, 1.55e-2, 1.7884, 1.50e-2, 1.7884, 1.35e-2, 1.7884, 8.90e-3,
426  .3981, 1.7860, 1.47e-2, 1.7860, 1.40e-2, 1.7860, 1.25e-2, 1.7860, 7.34e-3,
427  .4467, 1.7843, 1.39e-2, 1.7843, 1.31e-2, 1.7843, 1.15e-2, 1.7843, 6.40e-3,
428  .5012, 1.7832, 1.32e-2, 1.7832, 1.23e-2, 1.7832, 1.06e-2, 1.7832, 5.60e-3,
429  .5623, 1.7825, 1.25e-2, 1.7825, 1.15e-2, 1.7825, 9.77e-3, 1.7825, 5.00e-3,
430  .6310, 1.7820, 1.18e-2, 1.7820, 1.08e-2, 1.7820, 9.01e-3, 1.7820, 4.52e-3,
431  .7943, 1.7817, 1.06e-2, 1.7817, 9.46e-3, 1.7816, 7.66e-3, 1.7815, 3.68e-3,
432  1.000, 1.7816, 9.54e-3, 1.7816, 8.29e-3, 1.7814, 6.52e-3, 1.7807, 2.99e-3,
433  1.259, 1.7819, 8.56e-3, 1.7819, 7.27e-3, 1.7816, 5.54e-3, 1.7801, 2.49e-3,
434  2.500, 1.7830, 6.21e-3, 1.7830, 4.91e-3, 1.7822, 3.42e-3, 1.7789, 1.55e-3,
435  5.000, 1.7843, 4.49e-3, 1.7843, 3.30e-3, 1.7831, 2.10e-3, 1.7779, 9.61e-4,
436  10.00, 1.7852, 3.24e-3, 1.7852, 2.22e-3, 1.7838, 1.29e-3, 1.7773, 5.95e-4,
437  20.00, 1.7862, 2.34e-3, 1.7861, 1.49e-3, 1.7839, 7.93e-4, 1.7772, 3.69e-4,
438  32.00, 1.7866, 1.88e-3, 1.7863, 1.14e-3, 1.7840, 5.70e-4, 1.7772, 2.67e-4,
439  35.00, 1.7868, 1.74e-3, 1.7864, 1.06e-3, 1.7840, 5.35e-4, 1.7772, 2.51e-4,
440  40.00, 1.7869, 1.50e-3, 1.7865, 9.48e-4, 1.7840, 4.82e-4, 1.7772, 2.29e-4,
441  45.00, 1.7870, 1.32e-3, 1.7865, 8.50e-4, 1.7840, 4.38e-4, 1.7772, 2.11e-4,
442  50.00, 1.7870, 1.16e-3, 1.7865, 7.66e-4, 1.7840, 4.08e-4, 1.7772, 1.96e-4,
443  60.00, 1.7871, 8.80e-4, 1.7865, 6.30e-4, 1.7839, 3.50e-4, 1.7772, 1.73e-4,
444  70.00, 1.7871, 6.95e-4, 1.7865, 5.20e-4, 1.7838, 3.20e-4, 1.7772, 1.55e-4,
445  90.00, 1.7872, 4.64e-4, 1.7865, 3.84e-4, 1.7837, 2.55e-4, 1.7772, 1.31e-4,
446  111.0, 1.7872, 3.40e-4, 1.7865, 2.96e-4, 1.7837, 2.12e-4, 1.7772, 1.13e-4,
447  120.0, 1.7872, 3.11e-4, 1.7865, 2.70e-4, 1.7837, 2.00e-4, 1.7772, 1.06e-4,
448  130.0, 1.7872, 2.94e-4, 1.7865, 2.52e-4, 1.7837, 1.86e-4, 1.7772, 9.90e-5,
449  140.0, 1.7872, 2.79e-4, 1.7865, 2.44e-4, 1.7837, 1.75e-4, 1.7772, 9.30e-5,
450  150.0, 1.7872, 2.70e-4, 1.7865, 2.36e-4, 1.7837, 1.66e-4, 1.7772, 8.73e-5,
451  160.0, 1.7872, 2.64e-4, 1.7865, 2.30e-4, 1.7837, 1.56e-4, 1.7772, 8.30e-5,
452  170.0, 1.7872, 2.58e-4, 1.7865, 2.28e-4, 1.7837, 1.49e-4, 1.7772, 7.87e-5,
453  180.0, 1.7872, 2.52e-4, 1.7865, 2.25e-4, 1.7837, 1.44e-4, 1.7772, 7.50e-5,
454  200.0, 1.7872, 2.49e-4, 1.7865, 2.20e-4, 1.7837, 1.35e-4, 1.7772, 6.83e-5,
455  250.0, 1.7872, 2.54e-4, 1.7865, 2.16e-4, 1.7837, 1.21e-4, 1.7772, 5.60e-5,
456  290.0, 1.7872, 2.64e-4, 1.7865, 2.17e-4, 1.7837, 1.16e-4, 1.7772, 4.96e-5,
457  320.0, 1.7872, 2.74e-4, 1.7865, 2.20e-4, 1.7837, 1.16e-4, 1.7772, 4.55e-5,
458  350.0, 1.7872, 2.89e-4, 1.7665, 2.25e-4, 1.7837, 1.17e-4, 1.7772, 4.21e-5,
459  380.0, 1.7872, 3.05e-4, 1.7865, 2.32e-4, 1.7837, 1.20e-4, 1.7772, 3.91e-5,
460  400.0, 1.7872, 3.15e-4, 1.7865, 2.39e-4, 1.7837, 1.23e-4, 1.7772, 3.76e-5,
461  450.0, 1.7872, 3.46e-4, 1.7865, 2.60e-4, 1.7837, 1.32e-4, 1.7772, 3.40e-5,
462  500.0, 1.7872, 3.82e-4, 1.7865, 2.86e-4, 1.7837, 1.44e-4, 1.7772, 3.10e-5,
463  600.0, 1.7872, 4.62e-4, 1.7865, 3.56e-4, 1.7837, 1.68e-4, 1.7772, 2.64e-5,
464  640.0, 1.7872, 5.00e-4, 1.7865, 3.83e-4, 1.7837, 1.80e-4, 1.7772, 2.51e-5,
465  680.0, 1.7872, 5.50e-4, 1.7865, 4.15e-4, 1.7837, 1.90e-4, 1.7772, 2.43e-5,
466  720.0, 1.7872, 5.95e-4, 1.7865, 4.45e-4, 1.7837, 2.09e-4, 1.7772, 2.39e-5,
467  760.0, 1.7872, 6.47e-4, 1.7865, 4.76e-4, 1.7837, 2.16e-4, 1.7772, 2.37e-5,
468  800.0, 1.7872, 6.92e-4, 1.7865, 5.08e-4, 1.7837, 2.29e-4, 1.7772, 2.38e-5,
469  840.0, 1.7872, 7.42e-4, 1.7865, 5.40e-4, 1.7837, 2.40e-4, 1.7772, 2.40e-5,
470  900.0, 1.7872, 8.20e-4, 1.7865, 5.86e-4, 1.7837, 2.60e-4, 1.7772, 2.46e-5,
471  1000., 1.7872, 9.70e-4, 1.7865, 6.78e-4, 1.7837, 2.92e-4, 1.7772, 2.66e-5,
472  2000., 1.7872, 1.95e-3, 1.7865, 1.28e-3, 1.7837, 6.10e-4, 1.7772, 4.45e-5,
473  5000., 1.7872, 5.78e-3, 1.7865, 3.55e-3, 1.7840, 1.02e-3, 1.7772, 8.70e-5,
474  8600., 1.7880, 9.70e-3, 1.7872, 5.60e-3, 1.7845, 1.81e-3, 1.7780, 1.32e-4
475  };
476 
477  for (size_t i = 0; i < 558; ++i)
478  {
479  if (i % 9 == 0) wavelengths.push_back(tbl[i]);
480  else if (i % 9 == 1 || i % 9 == 3 || i % 9 == 5 || i % 9 == 7) vals_re.push_back(tbl[i]);
481  else vals_im.push_back(tbl[i]);
482  }
483 
484  static std::vector< std::vector<double>::iterator > grid_iter_list;
485  grid_iter_list.push_back(tempCs.begin());
486  grid_iter_list.push_back(wavelengths.begin());
487 
488  // the size of the grid in each dimension
489  static std::array<size_t, 2> grid_sizes;
490  grid_sizes[0] = tempCs.size();
491  grid_sizes[1] = wavelengths.size();
492  // total number of elements
493  static size_t num_elements = grid_sizes[0] * grid_sizes[1];
494 
495  // construct the interpolator. the last two arguments are pointers to the underlying data
496  interp_ML_warren_re = std::shared_ptr<InterpMultilinear<2, double > >(new InterpMultilinear<2, double >
497  (grid_iter_list.begin(), grid_sizes.begin(), vals_re.data(), vals_re.data() + num_elements));
498  interp_ML_warren_im = std::shared_ptr<InterpMultilinear<2, double > >(new InterpMultilinear<2, double >
499  (grid_iter_list.begin(), grid_sizes.begin(), vals_im.data(), vals_im.data() + num_elements));
500 
501  setup = true;
502  };
503 
504 
505  setupWarren();
506 
507  // interpolate one value
508  double wvlen = icedb::units::conv_spec("GHz", "mm").convert(f);
509  std::array<double, 2> args = { wvlen, t - 273.15 };
510  m = std::complex<double>(interp_ML_warren_re->interp(args.begin()),
511  -1.0 * interp_ML_warren_im->interp(args.begin()));
512  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mIceWarren result for freq: " << f << " GHz, temp " << t << " K is " << m);
513 }
514 
515 
516 void icedb::refract::implementations::mWaterHanel(double lambda, std::complex<double> &m)
517 {
518  if (lambda < 0.2 || lambda > 30000)
519  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
520  .add<double>("Wavelength (um)", lambda)
521  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,30000).");
522 
523  array<double, 1> args = { lambda };
524  m = std::complex<double>(setupHanelA(hanelAmedium::WATER_RE)->interp(args),
525  -1.0 * setupHanelA(hanelAmedium::WATER_IM)->interp(args));
526  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mWaterHanel result for lambda: " << lambda << " is " << m);
527 }
528 
529 void icedb::refract::implementations::mIceHanel(double lambda, std::complex<double> &m)
530 {
531  if (lambda < 0.2 || lambda > 30000)
532  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
533  .add<double>("Wavelength (um)", lambda)
534  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,30000).");
535 
536  array<double, 1> args = { lambda };
537  m = std::complex<double>(setupHanelA(hanelAmedium::ICE_RE)->interp(args),
538  -1.0 * setupHanelA(hanelAmedium::ICE_IM)->interp(args));
539  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mIceHanel result for lambda: " << lambda << " is " << m);
540 }
541 
542 void icedb::refract::implementations::mNaClHanel(double lambda, std::complex<double> &m)
543 {
544  if (lambda < 0.2 || lambda > 30000)
545  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
546  .add<double>("Wavelength (um)", lambda)
547  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,30000).");
548 
549  //double wvlen = icedb::units::conv_spec("GHz", "mm").convert(f);
550  array<double, 1> args = { lambda };
551  m = std::complex<double>(setupHanelA(hanelAmedium::NACL_RE)->interp(args),
552  -1.0 * setupHanelA(hanelAmedium::NACL_IM)->interp(args));
553  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mNaClHanel result for lambda: " << lambda << " is " << m);
554 }
555 
556 void icedb::refract::implementations::mSeaSaltHanel(double lambda, std::complex<double> &m)
557 {
558  if (lambda < 0.2 || lambda > 30000)
559  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
560  .add<double>("Wavelength (um)", lambda)
561  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,30000).");
562  //double wvlen = icedb::units::conv_spec("GHz", "mm").convert(f);
563  array<double, 1> args = { lambda };
564  m = std::complex<double>(setupHanelA(hanelAmedium::SEASALT_RE)->interp(args),
565  -1.0 * setupHanelA(hanelAmedium::SEASALT_IM)->interp(args));
566  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mSeaSaltHanel result for lambda: " << lambda << " is " << m);
567 }
568 
569 void icedb::refract::implementations::mDustHanel(double lambda, std::complex<double> &m)
570 {
571  if (lambda < 0.2 || lambda > 300)
572  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
573  .add<double>("Wavelength (um)", lambda)
574  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,300).");
575  //double wvlen = icedb::units::conv_spec("GHz", "mm").convert(f);
576  array<double, 1> args = { lambda };
577  m = std::complex<double>(setupHanelB(hanelBmedium::DUST_LIKE_RE)->interp(args),
578  -1.0 * setupHanelB(hanelBmedium::DUST_LIKE_IM)->interp(args));
579  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mDustHanel result for lambda: " << lambda << " is " << m);
580 }
581 
582 void icedb::refract::implementations::mSandOHanel(double lambda, std::complex<double> &m)
583 {
584  if (lambda < 0.2 || lambda > 300)
585  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
586  .add<double>("Wavelength (um)", lambda)
587  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,300).");
588  //double wvlen = icedb::units::conv_spec("GHz", "mm").convert(f);
589  array<double, 1> args = { lambda };
590  m = std::complex<double>(setupHanelB(hanelBmedium::SAND_O_RE)->interp(args),
591  -1.0 * setupHanelB(hanelBmedium::SAND_O_IM)->interp(args));
592  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mSandOHanel result for lambda: " << lambda << " is " << m);
593 }
594 
595 void icedb::refract::implementations::mSandEHanel(double lambda, std::complex<double> &m)
596 {
597  if (lambda < 0.2 || lambda > 300)
598  ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
599  .add<double>("Wavelength (um)", lambda)
600  .add<std::string>("Reason", "Allowed wavelength range (um) is (0.2,300).");
601  //double wvlen = icedb::units::conv_spec("GHz", "mm").convert(f);
602  array<double, 1> args = { lambda };
603  m = std::complex<double>(setupHanelB(hanelBmedium::SAND_E_RE)->interp(args),
604  -1.0 * setupHanelB(hanelBmedium::SAND_E_IM)->interp(args));
605  ICEDB_log("refract", icedb::logging::ICEDB_LOG_DEBUG_2, "mSandEHanel result for lambda: " << lambda << " is " << m);
606 }
607 
anonymous_namespace{refractImpl.cpp}::hanelAmedium::WATER_IM
anonymous_namespace{refractImpl.cpp}::setupHanelA
std::shared_ptr< InterpMultilinear< 1, double > > setupHanelA(hanelAmedium m)
Definition: refractImpl.cpp:20
icedb::refract::implementations::mIceWarren
void mIceWarren(double f, double t, std::complex< double > &m)
Ice complex refractive index for microwave/uv.
Definition: refractImpl.cpp:389
icedb::refract::implementations::mWaterHanel
void mWaterHanel(double lambda, std::complex< double > &m)
Water complex refractive index for ir/vis.
Definition: refractImpl.cpp:516
icedb::refract::implementations::mIceMatzler
void mIceMatzler(double f, double t, std::complex< double > &m)
Definition: refractImpl.cpp:360
anonymous_namespace{refractImpl.cpp}::hanelBmedium::SAND_E_RE
anonymous_namespace{refractImpl.cpp}::hanelBmedium::SAND_O_RE
anonymous_namespace{refractImpl.cpp}::m_setup
std::mutex m_setup
Definition: refractImpl.cpp:17
anonymous_namespace{refractImpl.cpp}::hanelAmedium
hanelAmedium
Definition: refractImpl.cpp:19
anonymous_namespace{refractImpl.cpp}::hanelBmedium::DUST_LIKE_RE
std
STL namespace.
anonymous_namespace{refractImpl.cpp}::hanelBmedium::SAND_O_IM
icedb::refract::implementations::mNaClHanel
void mNaClHanel(double lambda, std::complex< double > &m)
Sodium chloride refractive index for ir/vis.
Definition: refractImpl.cpp:542
icedb::refract::implementations::mSandEHanel
void mSandEHanel(double lambda, std::complex< double > &m)
Sand E-ray refractive index for ir/vis (birefringent)
Definition: refractImpl.cpp:595
icedb::refract::implementations::mSandOHanel
void mSandOHanel(double lambda, std::complex< double > &m)
Sand O-ray refractvie index for ir/vis (birefringent)
Definition: refractImpl.cpp:582
anonymous_namespace{refractImpl.cpp}::hanelAmedium::NACL_IM
anonymous_namespace{refractImpl.cpp}::hanelAmedium::SEASALT_IM
ICEDB_throw
#define ICEDB_throw(x)
Definition: error.hpp:88
anonymous_namespace{refractImpl.cpp}::hanelAmedium::ICE_RE
anonymous_namespace{refractImpl.cpp}::hanelAmedium::SEASALT_RE
anonymous_namespace{refractImpl.cpp}::hanelAmedium::WATER_RE
icedb::units::conv_spec
Perform interconversions between frequency, wavelength and wavenumber (GHz, Hz, m, cm, um, cm^-1, m^-1)
Definition: units.hpp:56
icedb::refract::implementations::mDustHanel
void mDustHanel(double lambda, std::complex< double > &m)
Dust-like particle refractive index for ir/vis.
Definition: refractImpl.cpp:569
icedb::refract::implementations::mWaterLiebe
void mWaterLiebe(double f, double t, std::complex< double > &m)
Definition: refractImpl.cpp:292
icedb::units::converter::convert
virtual double convert(double inVal) const
Definition: units.cpp:52
icedb::refract::implementations::mWaterFreshMeissnerWentz
void mWaterFreshMeissnerWentz(double f, double t, std::complex< double > &m)
Definition: refractImpl.cpp:316
anonymous_namespace{refractImpl.cpp}::hanelAmedium::ICE_IM
anonymous_namespace{refractImpl.cpp}::hanelBmedium::DUST_LIKE_IM
anonymous_namespace{refractImpl.cpp}::setupHanelB
std::shared_ptr< InterpMultilinear< 1, double > > setupHanelB(hanelBmedium m)
Definition: refractImpl.cpp:173
ICEDB_log
#define ICEDB_log(c, p, x)
Definition: logging.hpp:45
anonymous_namespace{refractImpl.cpp}::hanelAmedium::NACL_RE
anonymous_namespace{refractImpl.cpp}::hanelBmedium
hanelBmedium
Definition: refractImpl.cpp:172
anonymous_namespace{refractImpl.cpp}::hanelBmedium::SAND_E_IM
icedb::refract::implementations::mIceHanel
void mIceHanel(double lambda, std::complex< double > &m)
Ice complex refractive index for ir/vis.
Definition: refractImpl.cpp:529
icedb::refract::implementations::mSeaSaltHanel
void mSeaSaltHanel(double lambda, std::complex< double > &m)
Sea salt refractive index for ir/vis.
Definition: refractImpl.cpp:556
Generated by   doxygen 1.8.13